Merge branch 'development' into ppm_full_well_balancing

Exec/hydro_tests/gresho_vortex/problem_initialize_state_data.H

@@ -73,7 +73,7 @@ void problem_initialize_state_data (int i, int j, int k,
                 }
 
                 u_tot += u_phi;
-                reint += p/(gamma_const - 1.0_rt);
+                reint += p/(eos_rp::eos_gamma - 1.0_rt);
             }
         }
     }

Exec/radiation_tests/Rad2Tshock/problem_initialize.H

@@ -14,8 +14,8 @@ void problem_initialize ()
 
     eos_state.rho = problem::rho0;
     eos_state.T = problem::T0;
-    for (int n = 0; n < NumSpec; n++) {
-        eos_state.xn[n] = 0.0_rt;
+    for (auto & X : eos_state.xn) {
+        X = 0.0_rt;
     }
     eos_state.xn[0]  = 1.0_rt;
 

Exec/radiation_tests/Rad2Tshock/problem_initialize_rad_data.H

@@ -12,6 +12,9 @@ void problem_initialize_rad_data (int i, int j, int k,
                                   const GeometryData& geomdata)
 {
 
+    amrex::ignore_unused(nugroup);
+    amrex::ignore_unused(dnugroup);
+
     const Real* dx = geomdata.CellSize();
     const Real* problo = geomdata.ProbLo();
 

Exec/radiation_tests/Rad2Tshock/problem_initialize_state_data.H

@@ -16,7 +16,7 @@ void problem_initialize_state_data (int i, int j, int k,
     // Provides the simulation to be run in the x,y,or z direction
     // where length direction is the length side in a square prism
 
-    Real length_cell;
+    Real length_cell{};
     if (problem::idir == 1) {
         length_cell = problo[0] + dx[0] * (static_cast<Real>(i) + 0.5_rt);
     } else if (problem::idir == 2) {

Exec/science/flame_wave/ci-benchmarks/job_info_params.txt

@@ -95,7 +95,7 @@
     castro.sdc_quadrature = 0
     castro.sdc_extra = 0
     castro.sdc_solver = 1
-    castro.use_axisymmetric_geom_source = 1
+    castro.use_geom_source = 1
     castro.add_sdc_react_source_to_advection = 1
     castro.hydro_memory_footprint_ratio = -1
     castro.fixed_dt = -1

Exec/science/subch_planar/Problem_Derive.cpp

@@ -23,9 +23,8 @@ void ca_dergradpoverp(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*nco
 
 #if AMREX_SPACEDIM == 3
     amrex::Error("3D not supported");
-#endif
-#if AMREX_SPACEDIM == 1
-    amrex::Error("1D not supported");
+#elif AMREX_SPACEDIM == 1
+    return; // Skip for 1D
 #endif
 
     Real dxinv = 1.0_rt / dx[0];
@@ -226,9 +225,8 @@ void ca_dergradpoverp1(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*nc
 
 #if AMREX_SPACEDIM == 3
     amrex::Error("3D not supported");
-#endif
-#if AMREX_SPACEDIM == 1
-    amrex::Error("1D not supported");
+#elif AMREX_SPACEDIM == 1
+    return; // Skip for 1D
 #endif
 
     Real dxinv = 1.0_rt / dx[0];
@@ -434,9 +432,8 @@ void ca_dergradpx(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/
 
 #if AMREX_SPACEDIM == 3
     amrex::Error("3D not supported");
-#endif
-#if AMREX_SPACEDIM == 1
-    amrex::Error("1D not supported");
+#elif AMREX_SPACEDIM == 1
+    return; // Skip for 1D
 #endif
 
     Real dxinv = 1.0_rt / dx[0];
@@ -545,9 +542,8 @@ void ca_dergradpy(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/
 
 #if AMREX_SPACEDIM == 3
     amrex::Error("3D not supported");
-#endif
-#if AMREX_SPACEDIM == 1
-    amrex::Error("1D not supported");
+#elif AMREX_SPACEDIM == 1
+    return; // Skip for 1D
 #endif
 
     Real dxinv = 1.0_rt / dx[0];
@@ -658,9 +654,8 @@ void ca_dergradrhooverrho(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /
 
 #if AMREX_SPACEDIM == 3
     amrex::Error("3D not supported");
-#endif
-#if AMREX_SPACEDIM == 1
-    amrex::Error("1D not supported");
+#elif AMREX_SPACEDIM == 1
+    return; // Skip for 1D
 #endif
 
     Real dxinv = 1.0_rt / dx[0];
@@ -719,4 +714,4 @@ void ca_dergradrhooverrho(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /
 
     });
 
-}
+}

Exec/science/subch_planar/problem_initialize_state_data.H

@@ -25,7 +25,9 @@ void problem_initialize_state_data (int i, int j, int k,
     Real z = problo[2] + dx[2] * (static_cast<Real>(k) + 0.5_rt);
 #endif
 
-#if AMREX_SPACEDIM == 2
+#if AMREX_SPACEDIM == 1
+    Real height = x;
+#elif AMREX_SPACEDIM == 2
     Real height = y;
 #else
     Real height = z;

Source/driver/Castro.cpp

@@ -737,7 +737,7 @@ Castro::Castro (Amr&            papa,
       }
       radiation->regrid(level, grids, dmap);
 
-      rad_solver.reset(new RadSolve(parent, level, grids, dmap));
+      rad_solver = std::make_unique<RadSolve>(parent, level, grids, dmap);
     }
 #endif
 
@@ -751,7 +751,7 @@ Castro::Castro (Amr&            papa,
 Castro::~Castro ()  // NOLINT(modernize-use-equals-default)
 {
 #ifdef RADIATION
-    if (radiation != 0) {
+    if (radiation != nullptr) {
       //radiation->cleanup(level);
       radiation->close(level);
     }

Source/driver/Castro_util.H

@@ -157,10 +157,12 @@ AMREX_GPU_HOST_DEVICE AMREX_INLINE
 Real volume(const int& i, const int& j, const int& k,
             const GeometryData& geomdata)
 {
+    //
     // Given 3D cell-centered indices (i,j,k), return the volume of the zone.
-    // Note that Castro has no support for angular coordinates, so
-    // this function only provides Cartesian in 1D/2D/3D, Cylindrical (R-Z)
-    // in 2D, and Spherical in 1D.
+    // this function only provides Cartesian in 1D/2D/3D,
+    // Cylindrical (R-Z) in 2D,
+    // and Spherical in 1D and 2D (R-THETA).
+    //
 
     amrex::ignore_unused(i);
     amrex::ignore_unused(j);
@@ -210,8 +212,7 @@ Real volume(const int& i, const int& j, const int& k,
 
         vol = dx[0] * dx[1];
 
-    }
-    else {
+    } else if (coord == 1) {
 
         // Cylindrical
 
@@ -220,6 +221,20 @@ Real volume(const int& i, const int& j, const int& k,
 
         vol = M_PI * (r_l + r_r) * dx[0] * dx[1];
 
+    } else {
+
+        // Spherical
+
+        constexpr Real twoThirdsPi = 2.0_rt * M_PI / 3.0_rt;
+
+        Real r_l = geomdata.ProbLo()[0] + static_cast<Real>(i) * dx[0];
+        Real r_r = geomdata.ProbLo()[0] + static_cast<Real>(i+1) * dx[0];
+        Real theta_l = geomdata.ProbLo()[1] + static_cast<Real>(j) * dx[1];
+        Real theta_r = geomdata.ProbLo()[1] + static_cast<Real>(j+1) * dx[1];
+
+        vol = twoThirdsPi * (std::cos(theta_l) - std::cos(theta_r)) * dx[0] *
+            (r_r * r_r + r_r * r_l + r_l * r_l);
+
     }
 
 #else
@@ -290,8 +305,7 @@ Real area(const int& i, const int& j, const int& k,
             a = dx[0];
         }
 
-    }
-    else  {
+    } else  if (coord == 1) {
 
         // Cylindrical
 
@@ -304,6 +318,24 @@ Real area(const int& i, const int& j, const int& k,
             a = 2.0_rt * M_PI * r * dx[0];
         }
 
+    } else {
+
+        // Spherical
+
+        if (idir == 0) {
+            Real r = geomdata.ProbLo()[0] + static_cast<Real>(i) * dx[0];
+            Real theta_l = geomdata.ProbLo()[1] + static_cast<Real>(j) * dx[1];
+            Real theta_r = geomdata.ProbLo()[1] + static_cast<Real>(j+1) * dx[1];
+
+            a = 2.0_rt * M_PI * r * r * (std::cos(theta_l) - std::cos(theta_r));
+        }
+        else {
+            Real r = geomdata.ProbLo()[0] + (static_cast<Real>(i) + 0.5_rt) * dx[0];
+            Real theta = geomdata.ProbLo()[1] + static_cast<Real>(j) * dx[1];
+
+            a = 2.0_rt * M_PI * std::sin(theta) * r * dx[0];
+        }
+
     }
 
 #else

Source/driver/_cpp_parameters

@@ -272,8 +272,9 @@ sdc_extra                    int           0
 # which SDC nonlinear solver to use?  1 = Newton, 2 = VODE, 3 = VODE for first iter
 sdc_solver                   int           1
 
-# for 2-d axisymmetry, do we include the geometry source terms from Bernand-Champmartin?
-use_axisymmetric_geom_source bool           1
+# Do we include geometry source terms due to local unit vectors in non-Cartesian Coord?
+# We currently support R-Z cylinderical 2D (Bernand-Champmartin) and R-THETA spherical 2D
+use_geom_source              bool          1
 
 # for simplified-SDC, do we add the reactive source prediction to the interface states
 # used in the advective source construction?

Source/driver/math.H

@@ -3,6 +3,11 @@
 
 #include <AMReX_REAL.H>
 #include <AMReX_Array.H>
+#include <AMReX_BLassert.H>
+#include <cmath>
+
+using namespace amrex::literals;
+
 
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
 void
@@ -16,4 +21,11 @@ cross_product(amrex::GpuArray<amrex::Real, 3> const& a,
 
 }
 
+
+AMREX_GPU_HOST_DEVICE AMREX_INLINE
+amrex::Real cot(amrex::Real x) {
+
+    AMREX_ASSERT(x != 0.0_rt || x != M_PI);
+    return std::cos(x) / std::sin(x);
+}
 #endif

Source/gravity/Gravity.H

@@ -567,12 +567,6 @@ public:
                                      const amrex::Vector<amrex::MultiFab*>& res,
                                      amrex::Real time);
 
-  static inline amrex::Real
-  get_const_grav() {
-    return gravity::const_grav;
-  }
-
-
 };
 
 ///

Source/gravity/Gravity.cpp

@@ -881,11 +881,17 @@ Gravity::get_old_grav_vector(int level, MultiFab& grav_vector, Real time)
     MultiFab grav(grids[level], dmap[level], AMREX_SPACEDIM, ng);
     grav.setVal(0.0,ng);
 
-    if (gravity::gravity_type == "ConstantGrav") {
+    const Geometry& geom = parent->Geom(level);
 
-       // Set to constant value in the AMREX_SPACEDIM direction and zero in all others.
+    if (gravity::gravity_type == "ConstantGrav") {
 
-       grav.setVal(gravity::const_grav,AMREX_SPACEDIM-1,1,ng);
+        if (AMREX_SPACEDIM == 2 && geom.Coord() == 2) {
+            // 2D spherical r-theta, we want g in the radial direction
+            grav.setVal(gravity::const_grav, 0, 1, ng);
+        } else {
+            // Set to constant value in the AMREX_SPACEDIM direction and zero in all others.
+            grav.setVal(gravity::const_grav, AMREX_SPACEDIM-1, 1, ng);
+        }
 
     } else if (gravity::gravity_type == "MonopoleGrav") {
 
@@ -895,7 +901,6 @@ Gravity::get_old_grav_vector(int level, MultiFab& grav_vector, Real time)
 
     } else if (gravity::gravity_type == "PoissonGrav") {
 
-       const Geometry& geom = parent->Geom(level);
        amrex::average_face_to_cellcenter(grav, amrex::GetVecOfConstPtrs(grad_phi_prev[level]), geom);
        grav.mult(-1.0, ng); // g = - grad(phi)
 
@@ -952,11 +957,17 @@ Gravity::get_new_grav_vector(int level, MultiFab& grav_vector, Real time)
 
     MultiFab grav(grids[level],dmap[level],AMREX_SPACEDIM,ng);
     grav.setVal(0.0,ng);
+    const Geometry& geom = parent->Geom(level);
 
     if (gravity::gravity_type == "ConstantGrav") {
 
-       // Set to constant value in the AMREX_SPACEDIM direction
-       grav.setVal(gravity::const_grav,AMREX_SPACEDIM-1,1,ng);
+        if (AMREX_SPACEDIM == 2 && geom.Coord() == 2) {
+            // 2D spherical r-theta, we want g in the radial direction
+            grav.setVal(gravity::const_grav, 0, 1, ng);
+        } else {
+            // Set to constant value in the AMREX_SPACEDIM direction
+            grav.setVal(gravity::const_grav, AMREX_SPACEDIM-1, 1, ng);
+        }
 
     } else if (gravity::gravity_type == "MonopoleGrav") {
 
@@ -967,7 +978,6 @@ Gravity::get_new_grav_vector(int level, MultiFab& grav_vector, Real time)
 
     } else if (gravity::gravity_type == "PoissonGrav") {
 
-        const Geometry& geom = parent->Geom(level);
         amrex::average_face_to_cellcenter(grav, amrex::GetVecOfConstPtrs(grad_phi_curr[level]), geom);
         grav.mult(-1.0, ng); // g = - grad(phi)