add geometric source terms for 2D spherical geometry

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

@@ -94,7 +94,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

Source/driver/_cpp_parameters

@@ -268,8 +268,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/sources/Castro_geom.cpp

@@ -1,31 +1,43 @@
-#include "Castro.H"
+#include <Castro.H>
+#include <math.H>
 
 using namespace amrex;
 
 ///
-/// this adds the geometric source term for axisymmetric
-/// coordinates as described in Bernand-Champmartin.  This only
-/// applies to axisymmetric geometry.
+/// this adds the geometric source terms for non-Cartesian Coordinates
+/// This includes 2D Cylindrical (R-Z) coordinate as described in Bernand-Champmartin
+/// and 2D Spherical (R-THETA) coordinate.
 ///
+
 void
 Castro::construct_old_geom_source(MultiFab& source, MultiFab& state_in, Real time, Real dt)
 {
 
-  if (geom.Coord() != 1) {
+  amrex::ignore_unused(source);
+  amrex::ignore_unused(state_in);
+  amrex::ignore_unused(time);
+  amrex::ignore_unused(dt);
+
+  if (geom.Coord() == 0) {
       return;
   }
 
-  if (use_axisymmetric_geom_source == 0) {
+  if (use_geom_source == 0) {
       return;
   }
 
+#if AMREX_SPACEDIM == 2
   const Real strt_time = ParallelDescriptor::second();
 
   MultiFab geom_src(grids, dmap, source.nComp(), 0);
 
   geom_src.setVal(0.0);
 
-  fill_geom_source(time, dt, state_in, geom_src);
+  if (geom.Coord() == 1) {
+      fill_RZ_geom_source(time, dt, state_in, geom_src);
+  } else {
+      fill_RTheta_geom_source(time, dt, state_in, geom_src);
+  }
 
   Real mult_factor = 1.0;
 
@@ -47,22 +59,31 @@ Castro::construct_old_geom_source(MultiFab& source, MultiFab& state_in, Real tim
 #endif
     }
 
+#endif
 }
 
 
 
 void
-Castro::construct_new_geom_source(MultiFab& source, MultiFab& state_old, MultiFab& state_new, Real time, Real dt)
+Castro::construct_new_geom_source(MultiFab& source, MultiFab& state_old,
+                                  MultiFab& state_new, Real time, Real dt)
 {
 
-  if (geom.Coord() != 1) {
+  amrex::ignore_unused(source);
+  amrex::ignore_unused(state_old);
+  amrex::ignore_unused(state_new);
+  amrex::ignore_unused(time);
+  amrex::ignore_unused(dt);
+
+  if (geom.Coord() == 0) {
       return;
   }
 
-  if (use_axisymmetric_geom_source == 0) {
+  if (use_geom_source == 0) {
       return;
   }
 
+#if AMREX_SPACEDIM == 2
   const Real strt_time = ParallelDescriptor::second();
 
   MultiFab geom_src(grids, dmap, source.nComp(), 0);
@@ -72,7 +93,11 @@ Castro::construct_new_geom_source(MultiFab& source, MultiFab& state_old, MultiFa
   // Subtract off the old-time value first
   Real old_time = time - dt;
 
-  fill_geom_source(old_time, dt, state_old, geom_src);
+  if (geom.Coord() == 1) {
+      fill_RZ_geom_source(old_time, dt, state_old, geom_src);
+  } else {
+      fill_RTheta_geom_source(old_time, dt, state_old, geom_src);
+  }
 
   Real mult_factor = -0.5;
 
@@ -84,7 +109,11 @@ Castro::construct_new_geom_source(MultiFab& source, MultiFab& state_old, MultiFa
 
   mult_factor = 0.5;
 
-  fill_geom_source(time, dt, state_new, geom_src);
+  if (geom.Coord() == 1) {
+      fill_RZ_geom_source(time, dt, state_new, geom_src);
+  } else {
+      fill_RTheta_geom_source(time, dt, state_new, geom_src);
+  }
 
   MultiFab::Saxpy(source, mult_factor, geom_src, 0, 0, source.nComp(), 0);   // NOLINT(readability-suspicious-call-argument)
 
@@ -104,20 +133,22 @@ Castro::construct_new_geom_source(MultiFab& source, MultiFab& state_old, MultiFa
 #endif
     }
 
-
+#endif
 }
 
 
 
 void
-Castro::fill_geom_source ([[maybe_unused]] Real time, [[maybe_unused]] Real dt,
-                          MultiFab& cons_state, MultiFab& geom_src)
+Castro::fill_RZ_geom_source (Real time, Real dt, MultiFab& cons_state, MultiFab& geom_src)
 {
 
-  // Compute the geometric source for axisymmetric coordinates
+  // Compute the geometric source for axisymmetric coordinates (R-Z)
   // resulting from taking the divergence of (rho U U) in cylindrical
   // coordinates.  See the paper by Bernard-Champmartin
 
+  amrex::ignore_unused(time);
+  amrex::ignore_unused(dt);
+
   auto dx = geom.CellSizeArray();
   auto prob_lo = geom.ProbLoArray();
 
@@ -148,3 +179,53 @@ Castro::fill_geom_source ([[maybe_unused]] Real time, [[maybe_unused]] Real dt,
     });
   }
 }
+
+
+
+void
+Castro::fill_RTheta_geom_source (Real time, Real dt, MultiFab& cons_state, MultiFab& geom_src)
+{
+
+  // Compute the geometric source resulting from taking the divergence of (rho U U)
+  // in Spherical 2D (R-Theta) coordinate.
+
+  amrex::ignore_unused(time);
+  amrex::ignore_unused(dt);
+
+  auto dx = geom.CellSizeArray();
+  auto prob_lo = geom.ProbLoArray();
+
+#ifdef _OPENMP
+#pragma omp parallel
+#endif
+  for (MFIter mfi(geom_src, TilingIfNotGPU()); mfi.isValid(); ++mfi) {
+
+    const Box& bx = mfi.tilebox();
+
+    Array4<Real const> const U_arr = cons_state.array(mfi);
+
+    Array4<Real> const src = geom_src.array(mfi);
+
+    amrex::ParallelFor(bx,
+    [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+    {
+
+      // Cell-centered Spherical Radius and Theta
+      Real r = prob_lo[0] + (static_cast<Real>(i) + 0.5_rt)*dx[0];
+      Real theta = prob_lo[1] + (static_cast<Real>(j) + 0.5_rt)*dx[1];
+
+      // radial momentum: F = rho (v_theta**2 + v_phi**2) / r
+      src(i,j,k,UMX) = (U_arr(i,j,k,UMY) * U_arr(i,j,k,UMY) +
+                        U_arr(i,j,k,UMZ) * U_arr(i,j,k,UMZ)) / (U_arr(i,j,k,URHO) * r);
+
+      // Theta momentum F = rho v_phi**2 cot(theta) / r - rho v_r v_theta / r
+      src(i,j,k,UMY) = (U_arr(i,j,k,UMZ) * U_arr(i,j,k,UMZ) * cot(theta) -
+                        U_arr(i,j,k,UMX) * U_arr(i,j,k,UMY)) / (U_arr(i,j,k,URHO) * r);
+
+      // Phi momentum: F = - rho v_r v_phi / r - rho v_theta v_phi cot(theta) / r
+      src(i,j,k,UMZ) = (- U_arr(i,j,k,UMY) * U_arr(i,j,k,UMZ) * cot(theta) -
+                        U_arr(i,j,k,UMX) * U_arr(i,j,k,UMZ)) / (U_arr(i,j,k,URHO) * r);
+
+    });
+  }
+}

Source/sources/Castro_sources.H

@@ -330,15 +330,27 @@
 
 
 ///
-/// Fill ``ext_src`` with axisymmetric geometry sources
+/// Fill ``ext_src`` with 2D cylindrical R-Z geometry sources
 ///
 /// @param time     current time
 /// @param dt       timestep
 /// @param S        state
 /// @param ext_src  MultiFab to fill with sources
 ///
-    void fill_geom_source(amrex::Real time, amrex::Real dt,
-                          amrex::MultiFab& cons_state, amrex::MultiFab& geom_src);
+    void fill_RZ_geom_source(amrex::Real time, amrex::Real dt,
+                             amrex::MultiFab& cons_state, amrex::MultiFab& geom_src);
+
+
+///
+/// Fill ``ext_src`` with 2D spherical R-Theta geometry sources
+///
+/// @param time     current time
+/// @param dt       timestep
+/// @param S        state
+/// @param ext_src  MultiFab to fill with sources
+///
+    void fill_RTheta_geom_source(amrex::Real time, amrex::Real dt,
+                                 amrex::MultiFab& cons_state, amrex::MultiFab& geom_src);
 
 
 ///