Skip to content

Commit

Permalink
Merge branch 'development' of ssh://github.com/AMReX-Astro/Castro int…
Browse files Browse the repository at this point in the history 
…o development
  • Loading branch information
@zingale
zingale committed Jan 10, 2024
2 parents a739019 + b14f494 commit e1eecac
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 e1eecac

Please sign in to comment.