From 1816f1a9d2d8937ca2f0ac0bd579b4969832c3ad Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Thu, 18 Jan 2024 09:00:33 -0500 Subject: [PATCH] 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 | 3 +++ Source/hydro/Castro_ctu_hydro.cpp | 2 +- Source/hydro/Castro_mol_hydro.cpp | 4 ++-- Source/hydro/advection_util.cpp | 19 ++++++++++++++----- 4 files changed, 20 insertions(+), 8 deletions(-) diff --git a/Source/driver/_cpp_parameters b/Source/driver/_cpp_parameters index add5643322..8325fc86b6 100644 --- 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 index 96d64c4f55..c1a29ae0c2 100644 --- 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 index aadc6a41d8..3fddd8d072 100644 --- 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 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 index 3fe363592a..b40f8ee8af 100644 --- 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 const& q_arr, - Array4 const& q_src_arr, + Array4 const& U_src_arr, Array4 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;