make the shock threshold runtime settable

AMReX-Astro · Jan 7, 2024 · 4123355 · 4123355
1 parent 69d7001
commit 4123355
Show file tree

Hide file tree

Showing 2 changed files with 5 additions and 3 deletions.
diff --git a/Source/driver/_cpp_parameters b/Source/driver/_cpp_parameters
@@ -408,8 +408,12 @@ react_rho_min                Real          0.0
 react_rho_max                Real          1.e200
 
 # disable burning inside hydrodynamic shock regions
+# note: requires compiling with `USE_SHOCK_VAR=TRUE`
 disable_shock_burning        int           0
 
+# shock detection threshold for grad{P} / P
+shock_detection_threshold    Real          0.6666666666666666666666_rt
+
 # initial guess for the temperature when inverting the EoS (e.g. when
 # calling eos_input_re)
 T_guess                     Real           1.e8

diff --git a/Source/hydro/advection_util.cpp b/Source/hydro/advection_util.cpp
@@ -68,8 +68,6 @@ Castro::shock(const Box& bx,
   // we look for |grad P . dx| / P > 2/3 and div u < 0
   // This is basically the method in Gronow et al. 2020
 
-  constexpr Real eps = 2.0_rt / 3.0_rt;
-
   const auto dx = geom.CellSizeArray();
   const int coord_type = geom.Coord();
 
@@ -150,7 +148,7 @@ Castro::shock(const Box& bx,
 
     Real gradPdx_over_P = std::sqrt(dP_x * dP_x + dP_y * dP_y + dP_z * dP_z) / q_arr(i,j,k,QPRES);
 
-    if (gradPdx_over_P > eps && div_u < 0.0_rt) {
+    if (gradPdx_over_P > castro::shock_detection_threshold && div_u < 0.0_rt) {
       shk(i,j,k) = 1.0_rt;
     } else {
       shk(i,j,k) = 0.0_rt;