Skip to content

Commit

Permalink
a simpler / better shock detection algorithm (#2704)
Browse files Browse the repository at this point in the history 
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
  • Loading branch information
@zingale
zingale authored Jan 10, 2024
1 parent 9213839 commit b14f494
Show file tree
Hide file tree
Showing 5 changed files with 79 additions and 124 deletions.
4 changes: 4 additions & 0 deletions Source/driver/_cpp_parameters
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
26 changes: 13 additions & 13 deletions Source/hydro/Castro_ctu_hydro.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand All @@ -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);
Expand Down Expand Up @@ -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
});

Expand Down
1 change: 1 addition & 0 deletions Source/hydro/Castro_hydro.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -141,6 +141,7 @@
///
void shock(const amrex::Box& bx,
amrex::Array4<amrex::Real const> const& q_arr,
amrex::Array4<amrex::Real const> const& q_src_arr,
amrex::Array4<amrex::Real> const& shk);

///
Expand Down
71 changes: 38 additions & 33 deletions Source/hydro/Castro_mol_hydro.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<Real> 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);
Expand All @@ -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);
Expand Down Expand Up @@ -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;
});

}

Expand Down Expand Up @@ -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 {

Expand Down Expand Up @@ -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<Real> 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,
Expand Down Expand Up @@ -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);

Expand Down
101 changes: 23 additions & 78 deletions Source/hydro/advection_util.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -61,18 +61,12 @@ 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& 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();
Expand Down Expand Up @@ -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;
Expand Down

0 comments on commit b14f494

Please sign in to comment.