fix scaling of the source correction for the shock detection (#2711)

There was an extra dx scaling in here.
This also updates the source correction to follow the same logic as use_pslope
and makes it a runtime option.

Source/driver/_cpp_parameters

@@ -414,6 +414,9 @@ disable_shock_burning        int           0
 # shock detection threshold for grad{P} / P
 shock_detection_threshold    Real          0.6666666666666666666666_rt
 
+# do we subtract off the hydrostatic pressure when evaluating a shock?
+shock_detection_include_sources   int      1
+
 # initial guess for the temperature when inverting the EoS (e.g. when
 # calling eos_input_re)
 T_guess                     Real           1.e8

Source/hydro/Castro_ctu_hydro.cpp

@@ -283,7 +283,7 @@ Castro::construct_ctu_hydro_source(Real time, Real dt)
       });
 
       if (hybrid_riemann == 1 || compute_shock) {
-        shock(obx, q_arr, src_q_arr, shk_arr);
+        shock(obx, q_arr, old_src_arr, shk_arr);
       }
       else {
         amrex::ParallelFor(obx,

Source/hydro/Castro_mol_hydro.cpp

@@ -143,7 +143,7 @@ Castro::construct_mol_hydro_source(Real time, Real dt, MultiFab& A_update)
         src_q.resize(qbx, NQSRC);
         Array4<Real> const src_q_arr = src_q.array();
 
-        if (hybrid_riemann == 1 || compute_shock || (sdc_order == 2)) {
+        if (sdc_order == 2) {
 
             amrex::ParallelFor(qbx,
             [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
@@ -162,7 +162,7 @@ Castro::construct_mol_hydro_source(Real time, Real dt, MultiFab& A_update)
 
 
         if (hybrid_riemann == 1 || compute_shock) {
-            shock(obx, q_arr, src_q_arr, shk_arr);
+            shock(obx, q_arr, source_in_arr, shk_arr);
         }
         else {
             amrex::ParallelFor(obx,

Source/hydro/advection_util.cpp

@@ -61,7 +61,7 @@ Castro::ctoprim(const Box& bx,
 void
 Castro::shock(const Box& bx,
               Array4<Real const> const& q_arr,
-              Array4<Real const> const& q_src_arr,
+              Array4<Real const> const& U_src_arr,
               Array4<Real> const& shk) {
 
   // This is a basic multi-dimensional shock detection algorithm.
@@ -136,14 +136,23 @@ Castro::shock(const Box& bx,
     // We subtract off the hydrostatic force, since the pressure that
     // balances that is not available to make a shock.
     // We'll use a centered diff for the pressure gradient.
-    Real dP_x = 0.5_rt * (q_arr(i+1,j,k,QPRES) - q_arr(i-1,j,k,QPRES)) - dx[0] * q_src_arr(i,j,k,QU);
+    Real dP_x = 0.5_rt * (q_arr(i+1,j,k,QPRES) - q_arr(i-1,j,k,QPRES));
+    if (shock_detection_include_sources == 1) {
+        dP_x += -0.25_rt * dx[0] * (U_src_arr(i+1,j,k,UMX) + 2.0_rt * U_src_arr(i,j,k,UMX) + U_src_arr(i-1,j,k,UMX));
+    }
     Real dP_y = 0.0_rt;
     Real dP_z = 0.0_rt;
 #if AMREX_SPACEDIM >= 2
-    dP_y = 0.5_rt * (q_arr(i,j+1,k,QPRES) - q_arr(i,j-1,k,QPRES)) - dx[1] * q_src_arr(i,j,k,QV);
+    dP_y = 0.5_rt * (q_arr(i,j+1,k,QPRES) - q_arr(i,j-1,k,QPRES));
+    if (shock_detection_include_sources == 1) {
+        dP_y += -0.25_rt * dx[1] * (U_src_arr(i,j+1,k,UMY) + 2.0_rt * U_src_arr(i,j,k,UMY) + U_src_arr(i,j-1,k,UMY));
+    }
 #endif
 #if AMREX_SPACEDIM == 3
-    dP_z = 0.5_rt * (q_arr(i,j,k+1,QPRES) - q_arr(i,j,k-1,QPRES)) - dx[2] * q_src_arr(i,j,k,QW);
+    dP_z = 0.5_rt * (q_arr(i,j,k+1,QPRES) - q_arr(i,j,k-1,QPRES));
+    if (shock_detection_include_sources == 1) {
+        dP_z += -0.25_rt * dx[2] * (U_src_arr(i,j,k+1,UMZ) + 2.0_rt * U_src_arr(i,j,k,UMZ) + U_src_arr(i,j,k-1,UMZ));
+    }
 #endif
 
     //Real gradPdx_over_P = std::sqrt(dP_x * dP_x + dP_y * dP_y + dP_z * dP_z) / q_arr(i,j,k,QPRES);
@@ -157,7 +166,7 @@ Castro::shock(const Box& bx,
                                   dP_y * q_arr(i,j,k,QV) +
                                   dP_z * q_arr(i,j,k,QW)) / vel;
     }
-    gradPdx_over_P /= (q_arr(i,j,k,QPRES) / std::max(dx[0], std::max(dx[1], dx[2])));
+    gradPdx_over_P /= q_arr(i,j,k,QPRES);
 
     if (gradPdx_over_P > castro::shock_detection_threshold && div_u < 0.0_rt) {
       shk(i,j,k) = 1.0_rt;