better balancing

AMReX-Astro · Jan 12, 2024 · 6f7cae4 · 6f7cae4
1 parent f951227
commit 6f7cae4
Show file tree

Hide file tree

Showing 4 changed files with 19 additions and 10 deletions.
diff --git a/Source/driver/_cpp_parameters b/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

diff --git a/Source/hydro/Castro_ctu_hydro.cpp b/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,

diff --git a/Source/hydro/Castro_mol_hydro.cpp b/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,

diff --git a/Source/hydro/advection_util.cpp b/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,17 +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_arr(i,j,k,QRHO) * 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_arr(i,j,k,QRHO) * 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_arr(i,j,k,QRHO) * 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);