add r-theta geometric source terms

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
@@ -116,7 +116,7 @@
     castro.use_post_step_regrid = 0
     castro.max_subcycles = 10
     castro.sdc_iters = 2
-    castro.stopping_criterion_field = 
+    castro.stopping_criterion_field =
     castro.stopping_criterion_value = 1e+200
     castro.dtnuc_e = 1e+200
     castro.dtnuc_X = 1e+200
@@ -192,12 +192,12 @@
 [*] problem.H_star = 2000
 [*] problem.atm_delta = 100
     problem.fuel1_name = helium-4
-    problem.fuel2_name = 
-    problem.fuel3_name = 
-    problem.fuel4_name = 
+    problem.fuel2_name =
+    problem.fuel3_name =
+    problem.fuel4_name =
 [*] problem.ash1_name = nickel-56
-    problem.ash2_name = 
-    problem.ash3_name = 
+    problem.ash2_name =
+    problem.ash3_name =
     problem.fuel1_frac = 1
     problem.fuel2_frac = 0
     problem.fuel3_frac = 0

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,8 @@
 
 #include <AMReX_REAL.H>
 #include <AMReX_Array.H>
+#include <AMReX_BLassert.H>
+#include <cmath>
 
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
 void
@@ -16,4 +18,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,38 @@
-#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 geomtric 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) {
+  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,6 +54,7 @@ Castro::construct_old_geom_source(MultiFab& source, MultiFab& state_in, Real tim
 #endif
     }
 
+#endif
 }
 
 
@@ -55,14 +63,15 @@ void
 Castro::construct_new_geom_source(MultiFab& source, MultiFab& state_old, MultiFab& state_new, Real time, Real dt)
 {
 
-  if (geom.Coord() != 1) {
-      return;
-  }
+    if (geom.Coord() == 0) {
+        return;
+    }
 
-  if (use_axisymmetric_geom_source == 0) {
-      return;
-  }
+    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 +81,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(time, dt, state_in, geom_src);
+  } else {
+      fill_RTheta_geom_source(time, dt, state_in, geom_src);
+  }
 
   Real mult_factor = -0.5;
 
@@ -84,7 +97,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_in, geom_src);
+  } else {
+      fill_RTheta_geom_source(time, dt, state_in, geom_src);
+  }
 
   MultiFab::Saxpy(source, mult_factor, geom_src, 0, 0, source.nComp(), 0);   // NOLINT(readability-suspicious-call-argument)
 
@@ -104,17 +121,17 @@ 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 ([[maybe_unused]] Real time, [[maybe_unused]] 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
 
@@ -148,3 +165,51 @@ Castro::fill_geom_source ([[maybe_unused]] Real time, [[maybe_unused]] Real dt,
     });
   }
 }
+
+
+
+void
+Castro::fill_RTheta_geom_source ([[maybe_unused]] Real time, [[maybe_unused]] 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.
+
+  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);
+
+    });
+  }
+}