From b14f49459e68e0157a1b3532fe38b66102741000 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Wed, 10 Jan 2024 11:54:35 -0500 Subject: [PATCH] a simpler / better shock detection algorithm (#2704) The original algorithm was based on one in Flash but it had some normalizations in it that did not make sense. This now looks at | (grad P - rho g) . v / |v| | / (P / dx) > eps div{U} < 0 This projects grad P in the direction of the velocity (v / |v|) and also removes the hydrostatic pressure, since only the pressure in excess of that is available to make a shock. Here eps defaults to 2/3 and can be set at runtime via castro.shock_detection_threshold --- Source/driver/_cpp_parameters | 4 ++ Source/hydro/Castro_ctu_hydro.cpp | 26 ++++---- Source/hydro/Castro_hydro.H | 1 + Source/hydro/Castro_mol_hydro.cpp | 71 +++++++++++---------- Source/hydro/advection_util.cpp | 101 +++++++----------------------- 5 files changed, 79 insertions(+), 124 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/Castro_ctu_hydro.cpp b/Source/hydro/Castro_ctu_hydro.cpp index 985527f437..96d64c4f55 100644 --- a/Source/hydro/Castro_ctu_hydro.cpp +++ b/Source/hydro/Castro_ctu_hydro.cpp @@ -267,17 +267,6 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) bool compute_shock = false; #endif - if (hybrid_riemann == 1 || compute_shock) { - shock(obx, q_arr, shk_arr); - } - else { - amrex::ParallelFor(obx, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - shk_arr(i,j,k) = 0.0; - }); - } - // get the primitive variable hydro sources src_q.resize(qbx3, NQSRC); @@ -293,6 +282,17 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) hydro::src_to_prim(i, j, k, dt, U_old_arr, q_arr, old_src_arr, src_corr_arr, src_q_arr); }); + if (hybrid_riemann == 1 || compute_shock) { + shock(obx, q_arr, src_q_arr, shk_arr); + } + else { + amrex::ParallelFor(obx, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + shk_arr(i,j,k) = 0.0; + }); + } + // work on the interface states qxm.resize(obx, NQ); @@ -1143,8 +1143,8 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) flux_arr(i,j,k,USHK) = 0.e0; #endif #ifdef NSE_NET - flux_arr(i,j,k,UMUP) = 0.e0; - flux_arr(i,j,k,UMUN) = 0.e0; + flux_arr(i,j,k,UMUP) = 0.e0; + flux_arr(i,j,k,UMUN) = 0.e0; #endif }); diff --git a/Source/hydro/Castro_hydro.H b/Source/hydro/Castro_hydro.H index bad17cb707..647abc10ac 100644 --- a/Source/hydro/Castro_hydro.H +++ b/Source/hydro/Castro_hydro.H @@ -141,6 +141,7 @@ /// void shock(const amrex::Box& bx, amrex::Array4 const& q_arr, + amrex::Array4 const& q_src_arr, amrex::Array4 const& shk); /// diff --git a/Source/hydro/Castro_mol_hydro.cpp b/Source/hydro/Castro_mol_hydro.cpp index 05c3687eb9..aadc6a41d8 100644 --- a/Source/hydro/Castro_mol_hydro.cpp +++ b/Source/hydro/Castro_mol_hydro.cpp @@ -129,6 +129,29 @@ Castro::construct_mol_hydro_source(Real time, Real dt, MultiFab& A_update) }); } + +#ifdef SHOCK_VAR + bool compute_shock = true; +#else + bool compute_shock = false; +#endif + + // for well-balancing and shock detection, we need to + // primitive variable source terms + + const Box& qbx = amrex::grow(bx, NUM_GROW_SRC); + src_q.resize(qbx, NQSRC); + Array4 const src_q_arr = src_q.array(); + + if (hybrid_riemann == 1 || compute_shock || (sdc_order == 2)) { + + amrex::ParallelFor(qbx, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + hydro::src_to_prim(i, j, k, dt, uin_arr, q_arr, source_in_arr, src_q_arr); + }); + } + // get the interface states and shock variable shk.resize(obx, 1); @@ -137,21 +160,16 @@ Castro::construct_mol_hydro_source(Real time, Real dt, MultiFab& A_update) // Multidimensional shock detection // Used for the hybrid Riemann solver -#ifdef SHOCK_VAR - bool compute_shock = true; -#else - bool compute_shock = false; -#endif if (hybrid_riemann == 1 || compute_shock) { - shock(obx, q_arr, shk_arr); + shock(obx, q_arr, src_q_arr, shk_arr); } else { - amrex::ParallelFor(obx, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - shk_arr(i,j,k) = 0.0; - }); + amrex::ParallelFor(obx, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + shk_arr(i,j,k) = 0.0; + }); } const Box& xbx = amrex::surroundingNodes(bx, 0); @@ -251,11 +269,11 @@ Castro::construct_mol_hydro_source(Real time, Real dt, MultiFab& A_update) idir, true); if (do_hydro == 0) { - amrex::ParallelFor(nbx, NUM_STATE, - [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept - { - f_avg_arr(i,j,k,n) = 0.0; - }); + amrex::ParallelFor(nbx, NUM_STATE, + [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept + { + f_avg_arr(i,j,k,n) = 0.0; + }); } @@ -389,8 +407,8 @@ Castro::construct_mol_hydro_source(Real time, Real dt, MultiFab& A_update) flux_arr(i,j,k,n) = 0.0; #endif #ifdef NSE_NET - } else if (n == UMUP || n == UMUN) { - flux_arr(i,j,k,n) = 0.0; + } else if (n == UMUP || n == UMUN) { + flux_arr(i,j,k,n) = 0.0; #endif } else { @@ -454,18 +472,6 @@ Castro::construct_mol_hydro_source(Real time, Real dt, MultiFab& A_update) dq.resize(obx, NQ); auto dq_arr = dq.array(); - // for well-balancing, we need to primitive variable - // source terms - const Box& qbx = amrex::grow(bx, NUM_GROW_SRC); - src_q.resize(qbx, NQSRC); - Array4 const src_q_arr = src_q.array(); - - amrex::ParallelFor(qbx, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - hydro::src_to_prim(i, j, k, dt, uin_arr, q_arr, source_in_arr, src_q_arr); - }); - mol_plm_reconstruct(obx, idir, q_arr, flatn_arr, src_q_arr, dq_arr, @@ -506,12 +512,11 @@ Castro::construct_mol_hydro_source(Real time, Real dt, MultiFab& A_update) flux_arr(i,j,k,USHK) = 0.e0; #endif #ifdef NSE_NET - flux_arr(i,j,k,UMUP) = 0.e0; - flux_arr(i,j,k,UMUN) = 0.e0; + flux_arr(i,j,k,UMUP) = 0.e0; + flux_arr(i,j,k,UMUN) = 0.e0; #endif }); - // apply artificial viscosity apply_av(nbx, idir, div_arr, uin_arr, flux_arr); diff --git a/Source/hydro/advection_util.cpp b/Source/hydro/advection_util.cpp index 40752a6e0d..c75bf51268 100644 --- a/Source/hydro/advection_util.cpp +++ b/Source/hydro/advection_util.cpp @@ -61,18 +61,12 @@ Castro::ctoprim(const Box& bx, void Castro::shock(const Box& bx, Array4 const& q_arr, + Array4 const& q_src_arr, Array4 const& shk) { // This is a basic multi-dimensional shock detection algorithm. - // This implementation follows Flash, which in turn follows - // AMRA and a Woodward (1995) (supposedly -- couldn't locate that). - // - // The spirit of this follows the shock detection in Colella & - // Woodward (1984) - // - - constexpr Real small = 1.e-10_rt; - constexpr Real eps = 0.33e0_rt; + // we look for |grad P . dx| / P > 2/3 and div u < 0 + // This is basically the method in Gronow et al. 2020 const auto dx = geom.CellSizeArray(); const int coord_type = geom.Coord(); @@ -137,81 +131,32 @@ Castro::shock(const Box& bx, #endif } - // find the pre- and post-shock pressures in each direction - Real px_pre; - Real px_post; - Real e_x; - - if (q_arr(i+1,j,k,QPRES) - q_arr(i-1,j,k,QPRES) < 0.0_rt) { - px_pre = q_arr(i+1,j,k,QPRES); - px_post = q_arr(i-1,j,k,QPRES); - } else { - px_pre = q_arr(i-1,j,k,QPRES); - px_post = q_arr(i+1,j,k,QPRES); - } - - // use compression to create unit vectors for the shock direction - e_x = std::pow(q_arr(i+1,j,k,QU) - q_arr(i-1,j,k,QU), 2); - - Real py_pre; - Real py_post; - Real e_y; - -#if (AMREX_SPACEDIM >= 2) - if (q_arr(i,j+1,k,QPRES) - q_arr(i,j-1,k,QPRES) < 0.0_rt) { - py_pre = q_arr(i,j+1,k,QPRES); - py_post = q_arr(i,j-1,k,QPRES); - } else { - py_pre = q_arr(i,j-1,k,QPRES); - py_post = q_arr(i,j+1,k,QPRES); - } - - e_y = std::pow(q_arr(i,j+1,k,QV) - q_arr(i,j-1,k,QV), 2); - -#else - py_pre = 0.0_rt; - py_post = 0.0_rt; - e_y = 0.0_rt; + // now compute (grad P - rho g) . dx + // 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_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); #endif - - Real pz_pre; - Real pz_post; - Real e_z; - -#if (AMREX_SPACEDIM == 3) - if (q_arr(i,j,k+1,QPRES) - q_arr(i,j,k-1,QPRES) < 0.0_rt) { - pz_pre = q_arr(i,j,k+1,QPRES); - pz_post = q_arr(i,j,k-1,QPRES); - } else { - pz_pre = q_arr(i,j,k-1,QPRES); - pz_post = q_arr(i,j,k+1,QPRES); - } - - e_z = std::pow(q_arr(i,j,k+1,QW) - q_arr(i,j,k-1,QW), 2); - -#else - pz_pre = 0.0_rt; - pz_post = 0.0_rt; - - e_z = 0.0_rt; +#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); #endif - Real denom = 1.0_rt / (e_x + e_y + e_z + small); - - e_x = e_x * denom; - e_y = e_y * denom; - e_z = e_z * denom; - - // project the pressures onto the shock direction - Real p_pre = e_x * px_pre + e_y * py_pre + e_z * pz_pre; - Real p_post = e_x * px_post + e_y * py_post + e_z * pz_post; + //Real gradPdx_over_P = std::sqrt(dP_x * dP_x + dP_y * dP_y + dP_z * dP_z) / q_arr(i,j,k,QPRES); + Real vel = std::sqrt(q_arr(i,j,k,QU) * q_arr(i,j,k,QU) + + q_arr(i,j,k,QV) * q_arr(i,j,k,QV) + + q_arr(i,j,k,QW) * q_arr(i,j,k,QW)); - // test for compression + pressure jump to flag a shock - // this avoid U = 0, so e_x, ... = 0 - Real pjump = p_pre == 0 ? 0.0_rt : eps - (p_post - p_pre) / p_pre; + Real gradPdx_over_P = std::abs(dP_x * q_arr(i,j,k,QU) + + 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]))); - if (pjump < 0.0 && 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;