From 41233551319b48c1bdfcff0132390dc1d2e21889 Mon Sep 17 00:00:00 2001
From: Michael Zingale <michael.zingale@stonybrook.edu>
Date: Sun, 7 Jan 2024 10:43:03 -0500
Subject: [PATCH] make the shock threshold runtime settable

---
 Source/driver/_cpp_parameters   | 4 ++++
 Source/hydro/advection_util.cpp | 4 +---
 2 files changed, 5 insertions(+), 3 deletions(-)

diff --git a/Source/driver/_cpp_parameters b/Source/driver/_cpp_parameters
index 4d1dde2581..add5643322 100644
--- 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
index 2816db4d70..d43a14fbed 100644
--- 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;