Merge branch 'development' into add_sdc_normalization

AMReX-Astro · Oct 12, 2023 · db30a84 · db30a84
2 parents d860c8f + f8446cd
commit db30a84
Show file tree

Hide file tree

Showing 20 changed files with 311 additions and 901 deletions.
diff --git a/CHANGES.md b/CHANGES.md
@@ -1,5 +1,8 @@
 # 23.11
 
+  * The ca_rad_source hook in Fortran has been removed. The existing
+    problem_rad_source() hook in C++ can be used instead. (#2626)
+
   * The compile option USE_AUX_UPDATE has been removed. If you want to
     manually update the auxiliary parameters, you can use an external
     source term or you can use the problem post-timestep hook. (#2614)

diff --git a/Docs/source/reactions.rst b/Docs/source/reactions.rst
@@ -57,8 +57,32 @@ in the inputs file. Reactions are enabled by setting::
 the efficiency).
 
 It is possible to set the maximum and minimum temperature and density for allowing
-reactions to occur in a zone using the parameters ``castro.react_T_min``,
-``castro.react_T_max``, ``castro.react_rho_min`` and ``castro.react_rho_max``.
+reactions to occur in a zone using the parameters:
+
+* ``castro.react_T_min`` and ``castro.react_T_max`` for temperature
+
+* ``castro.react_rho_min`` and ``castro.react_rho_max`` for density
+
+.. index:: castro.disable_shock_burning, USE_SHOCK_VAR
+
+Burning can also be disabled inside shocks.  This requires that the code be
+compiled with::
+
+  USE_SHOCK_VAR = TRUE
+
+in the ``GNUmakefile``.  This will allocate storage for a shock flag in the conserved
+state array.  This flag is computed via a multidimensional shock detection algorithm
+that looks for compression (:math:`\nabla \cdot \ub < 0`) along with a pressure jump
+in the direction of compression.  The runtime parameter::
+
+   castro.disable_shock_burning = 1
+
+will skip reactions in a zone where we've detected a shock.
+
+.. note::
+
+   Both the compilation with ``USE_SHOCK_VAR = TRUE`` and the runtime parameter
+   ``castro.disable_shock_burning = 1`` are needed to turn off burning in shocks.
 
 Reactions Flowchart
 ===================

diff --git a/Source/driver/Castro_setup.cpp b/Source/driver/Castro_setup.cpp
@@ -685,15 +685,15 @@ Castro::variableSetUp ()
     std::cout << "Radiation::nGroups = " << Radiation::nGroups << std::endl;
   }
 
-  char rad_name[10];
+  std::string rad_name;
   if (!Radiation::do_multigroup) {
     desc_lst
       .setComponent(Rad_Type, Rad, "rad", bc,
                     genericBndryFunc);
   }
   else {
       for (int i = 0; i < Radiation::nGroups; i++) {
-        sprintf(rad_name, "rad%d", i);
+        rad_name = "rad" + std::to_string(i);
         desc_lst
           .setComponent(Rad_Type, i, rad_name, bc,
                         genericBndryFunc);

diff --git a/Source/radiation/CastroRad_1d.f90 b/Source/radiation/CastroRad_1d.f90
diff --git a/Source/radiation/CastroRad_2d.f90 b/Source/radiation/CastroRad_2d.f90
diff --git a/Source/radiation/CastroRad_3d.f90 b/Source/radiation/CastroRad_3d.f90
diff --git a/Source/radiation/HABEC.H b/Source/radiation/HABEC.H
@@ -367,6 +367,222 @@ namespace HABEC
             }
         });
     }
+
+    AMREX_INLINE
+    void hdterm (Array4<Real> const dterm,
+                 Array4<Real> const er,
+                 const Box& reg,
+                 int cdir, int bct, Real bcl,
+                 Array4<Real const> const bcval,
+                 Array4<int const> const mask,
+                 Array4<Real const> const d,
+                 const Real* dx)
+    {
+        Real h;
+
+        // Index shift for whichever edge we're checking.
+        // Negative means we're looking at the lo edge,
+        // positive means we're looking at the hi edge.
+        // The first group is for cell centers, the second
+        // group is for cell edges.
+
+        int icp = 0;
+        int jcp = 0;
+        int kcp = 0;
+
+        int iep = 0;
+        int jep = 0;
+        int kep = 0;
+
+#if AMREX_SPACEDIM == 1
+        if (cdir == 0) {
+            h = dx[0];
+            icp = -1;
+        }
+        else if (cdir == 1) {
+            h = dx[0];
+            icp = 1;
+            iep = 1;
+        }
+#elif AMREX_SPACEDIM == 2
+        if (cdir == 0) {
+            h = dx[0];
+            icp = -1;
+        }
+        else if (cdir == 2) {
+            h = dx[0];
+            icp = 1;
+            iep = 1;
+        }
+        else if (cdir == 1) {
+            h = dx[1];
+            jcp = -1;
+        }
+        else if (cdir == 3) {
+            h = dx[1];
+            jcp = 1;
+            jep = 1;
+        }
+#else
+        if (cdir == 0) {
+            h = dx[0];
+            icp = -1;
+        }
+        else if (cdir == 3) {
+            h = dx[0];
+            icp = 1;
+            iep = 1;
+        }
+        else if (cdir == 1) {
+            h = dx[1];
+            jcp = -1;
+        }
+        else if (cdir == 4) {
+            h = dx[1];
+            jcp = 1;
+            jep = 1;
+        }
+        else if (cdir == 2) {
+            h = dx[2];
+            kcp = -1;
+        }
+        else if (cdir == 5) {
+            h = dx[2];
+            kcp = 1;
+            kep = 1;
+        }
+#endif
+
+        amrex::LoopOnCpu(reg, [=] (int i, int j, int k) noexcept
+        {
+            if (mask.contains(i+icp,j+jcp,k+kcp)) {
+                if (mask(i+icp,j+jcp,k+kcp) > 0) {
+                    if (bct == LO_DIRICHLET) {
+                        Real d_sign = 1.0_rt;
+                        if (iep != 0 || jep != 0 || kep != 0) {
+                            // right edge
+                            d_sign = -1.0_rt;
+                        }
+                        dterm(i+iep,j+jep,k+kep) = d(i+iep,j+jep,k+kep) * d_sign * (er(i,j,k) - bcval(i+icp,j+jcp,k+kcp)) / (0.5_rt * h + bcl);
+                    }
+                }
+            }
+        });
+    }
+
+    AMREX_INLINE
+    void hdterm3 (Array4<Real> const dterm,
+                  Array4<Real const> const er,
+                  const Box& reg,
+                  int cdir, int bctype,
+                  Array4<int const> const tf,
+                  Real bcl,
+                  Array4<Real const> const bcval,
+                  Array4<int const> const mask,
+                  Array4<Real const> const d,
+                  const Real* const dx)
+    {
+        Real h;
+
+        // Index shift for whichever edge we're checking.
+        // Negative means we're looking at the lo edge,
+        // positive means we're looking at the hi edge.
+        // The first group is for cell centers, the second
+        // group is for cell edges.
+
+        int icp = 0;
+        int jcp = 0;
+        int kcp = 0;
+
+        int iep = 0;
+        int jep = 0;
+        int kep = 0;
+
+#if AMREX_SPACEDIM == 1
+        if (cdir == 0) {
+            h = dx[0];
+            icp = -1;
+        }
+        else if (cdir == 1) {
+            h = dx[0];
+            icp = 1;
+            iep = 1;
+        }
+#elif AMREX_SPACEDIM == 2
+        if (cdir == 0) {
+            h = dx[0];
+            icp = -1;
+        }
+        else if (cdir == 2) {
+            h = dx[0];
+            icp = 1;
+            iep = 1;
+        }
+        else if (cdir == 1) {
+            h = dx[1];
+            jcp = -1;
+        }
+        else if (cdir == 3) {
+            h = dx[1];
+            jcp = 1;
+            jep = 1;
+        }
+#else
+        if (cdir == 0) {
+            h = dx[0];
+            icp = -1;
+        }
+        else if (cdir == 3) {
+            h = dx[0];
+            icp = 1;
+            iep = 1;
+        }
+        else if (cdir == 1) {
+            h = dx[1];
+            jcp = -1;
+        }
+        else if (cdir == 4) {
+            h = dx[1];
+            jcp = 1;
+            jep = 1;
+        }
+        else if (cdir == 2) {
+            h = dx[2];
+            kcp = -1;
+        }
+        else if (cdir == 5) {
+            h = dx[2];
+            kcp = 1;
+            kep = 1;
+        }
+#endif
+
+        amrex::LoopOnCpu(reg, [=] (int i, int j, int k) noexcept
+        {
+            if (mask.contains(i+icp,j+jcp,k+kcp)) {
+                if (mask(i+icp,j+jcp,k+kcp) > 0) {
+                    int bct;
+                    if (bctype == -1) {
+                        bct = tf(i+icp,j+jcp,k+kcp);
+                    }
+                    else {
+                        bct = bctype;
+                    }
+                    if (bct == LO_DIRICHLET) {
+                        Real d_sign = 1.0_rt;
+                        if (iep != 0 || jep != 0 || kep != 0) {
+                            // right edge
+                            d_sign = -1.0_rt;
+                        }
+                        dterm(i+iep,j+jep,k+kep) = d(i+iep,j+jep,k+kep) * d_sign * (er(i,j,k) - bcval(i+icp,j+jcp,k+kcp)) / (0.5_rt * h + bcl);
+                    }
+                    else if (bct == LO_NEUMANN && bcval(i+icp,j+jcp,k+kcp) == 0.0_rt) {
+                        dterm(i+iep,j+jep,k+kep) = 0.0_rt;
+                    }
+                }
+            }
+        });
+    }
 } // namespace HABEC
 
 #endif