Compute the shock flag once -- at the start of the timestep

Source/hydro/Castro_ctu.cpp

@@ -9,9 +9,6 @@ using namespace amrex;
 
 void
 Castro::consup_hydro(const Box& bx,
-#ifdef SHOCK_VAR
-                     Array4<Real const> const& shk,
-#endif
                      Array4<Real> const& U_new,
                      Array4<Real> const& flux0,
                      Array4<Real const> const& qx,
@@ -68,11 +65,6 @@ Castro::consup_hydro(const Box& bx,
 
       U_new(i,j,k,n) = U_new(i,j,k,n) - dt * pdu;
 
-#ifdef SHOCK_VAR
-    } else if (n == USHK) {
-      U_new(i,j,k,USHK) = shk(i,j,k);
-#endif
-
     } else if (n == UMX) {
       // Add gradp term to momentum equation -- only for axisymmetric
       // coords (and only for the radial flux).

Source/hydro/Castro_ctu_hydro.cpp

@@ -1193,9 +1193,6 @@ Castro::construct_ctu_hydro_source(Real time, Real dt)  // NOLINT(readability-co
 #endif
 
       consup_hydro(bx,
-#ifdef SHOCK_VAR
-                   shk_arr,
-#endif
                    update_arr,
                    flx_arr, qx_arr,
 #if AMREX_SPACEDIM >= 2

Source/hydro/Castro_hydro.H

@@ -711,9 +711,6 @@
 /// @param dt      timestep
 ///
     void consup_hydro(const amrex::Box& bx,
-#ifdef SHOCK_VAR
-                      amrex::Array4<amrex::Real const> const& shk,
-#endif
                       amrex::Array4<amrex::Real> const& U_new,
                       amrex::Array4<amrex::Real> const& flux0,
                       amrex::Array4<amrex::Real const> const& qx,

Source/sources/Castro_sources.H

@@ -57,7 +57,7 @@
 /// @param dt       the timestep to advance (e.g., go from time to
 ///                    time + dt)
 ///
-    advance_status do_old_sources (amrex::Real time, amrex::Real dt);
+    advance_status do_old_sources (amrex::Real time, amrex::Real dt, bool apply_to_state = true);
 
 
 ///

Source/sources/Castro_sources.cpp

@@ -151,7 +151,7 @@ Castro::do_old_sources(
     // Optionally print out diagnostic information about how much
     // these source terms changed the state.
 
-    if (print_update_diagnostics) {
+    if (apply_to_state && print_update_diagnostics) {
       bool is_new = false;
       print_all_source_changes(dt, is_new);
     }
@@ -174,7 +174,7 @@ Castro::do_old_sources(
 }
 
 advance_status
-Castro::do_old_sources (Real time, Real dt)
+Castro::do_old_sources (Real time, Real dt, bool apply_to_state)
 {
     advance_status status {};
 
@@ -192,7 +192,8 @@ Castro::do_old_sources (Real time, Real dt)
 #ifdef MHD
                    Bx_old, By_old, Bz_old,
 #endif
-                   old_source, Sborder, S_new, time, dt);
+                   old_source, Sborder, S_new, time, dt,
+                   apply_to_state);
 
     return status;
 }
@@ -537,6 +538,71 @@ Castro::pre_advance_operators (Real time, Real dt)  // NOLINT(readability-conver
     construct_old_gravity(time);
 #endif
 
+#ifdef SHOCK_VAR
+    // we want to compute the shock flag that will be used
+    // (optionally) in disabling reactions in shocks.  We compute this
+    // only once per timestep using the time-level n data.
+
+    // We need to compute the old sources -- note that we will
+    // recompute the old sources after the burn, so this is done here
+    // only for evaluating the shock flag.
+
+    const bool apply_to_state{false};
+    do_old_sources(time, dt, apply_to_state);
+
+    MultiFab& S_old = get_old_data(State_Type);
+    MultiFab& old_source = get_old_data(Source_Type);
+
+    FArrayBox shk(The_Async_Arena());
+    FArrayBox q(The_Async_Arena()), qaux(The_Async_Arena());
+
+    for (MFIter mfi(Sborder, TilingIfNotGPU()); mfi.isValid(); ++mfi) {
+        const Box& bx = mfi.tilebox();
+        const Box& obx = mfi.growntilebox(1);
+
+        shk.resize(bx, 1);
+#ifdef RADIATION
+        q.resize(obx, NQ);
+#else
+        q.resize(obx, NQTHERM);
+#endif
+        qaux.resize(obx, NQAUX);
+
+        Array4<Real> const shk_arr = shk.array();
+        Array4<Real> const q_arr = q.array();
+        Array4<Real> const qaux_arr = qaux.array();
+
+        Array4<Real const> const Sborder_old_arr = Sborder.array(mfi);
+        Array4<Real> const S_old_arr = S_old.array(mfi);
+        Array4<Real> const old_src_arr = old_source.array(mfi);
+
+        ctoprim(obx, time, Sborder_old_arr, q_arr, qaux_arr);
+
+        shock(bx, q_arr, old_src_arr, shk_arr);
+
+        // now store it in S_old -- we'll fillpatch into Sborder in a bit
+
+        // Note: we still compute the shock flag in the hydro for the
+        // hybrid-Riemann (for now) since that version will have seen the
+        // effect of the burning in the first dt), but that version is
+        // never stored in State_Type
+
+        amrex::ParallelFor(bx,
+        [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+        {
+            S_old_arr(i,j,k,USHK) = shk_arr(i,j,k);
+        });
+
+    }
+
+    // we only computed the shock flag on the interior, but the first
+    // burn needs ghost cells, so FillPatch just the shock flag
+
+    if (Sborder.nGrow() > 0) {
+      AmrLevel::FillPatch(*this, Sborder, Sborder.nGrow(), time, State_Type, USHK, 1, USHK);
+    }
+
+#endif
 
     // If we are Strang splitting the reactions, do the old-time contribution now.