Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

[WIP] an alternate shock detection algorithm -- this uses the divu projection #2884

Closed
wants to merge 3 commits into from
Closed

[WIP] an alternate shock detection algorithm -- this uses the divu projection #2884

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
8fc2515
an alternate shock detection algorithm -- this uses the divu projection
zingale Jun 25, 2024
8db7ccb
Merge branch 'development' into divu_shock_flag
zingale Jul 2, 2024
17c0344
Merge branch 'development' into divu_shock_flag
zingale Jul 9, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
124 changes: 66 additions & 58 deletions Source/hydro/advection_util.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -82,97 +82,105 @@ Castro::shock(const Box& bx,
amrex::ParallelFor(bx,
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
Real div_u = 0.0_rt;
Real div_u = 0.0_rt;
Real du_x{};
Real dv_y{};
Real dw_z{};

// construct div{U}
if (coord_type == 0) {
// construct div{U}
if (coord_type == 0) {

// Cartesian
div_u += 0.5_rt * (q_arr(i+1,j,k,QU) - q_arr(i-1,j,k,QU)) * dxinv;
// Cartesian
du_x = 0.5_rt * (q_arr(i+1,j,k,QU) - q_arr(i-1,j,k,QU)) * dxinv;
#if (AMREX_SPACEDIM >= 2)
div_u += 0.5_rt * (q_arr(i,j+1,k,QV) - q_arr(i,j-1,k,QV)) * dyinv;
dv_y = 0.5_rt * (q_arr(i,j+1,k,QV) - q_arr(i,j-1,k,QV)) * dyinv;
#endif
#if (AMREX_SPACEDIM == 3)
div_u += 0.5_rt * (q_arr(i,j,k+1,QW) - q_arr(i,j,k-1,QW)) * dzinv;
dw_z = 0.5_rt * (q_arr(i,j,k+1,QW) - q_arr(i,j,k-1,QW)) * dzinv;
#endif

#if AMREX_SPACEDIM <= 2
} else if (coord_type == 1) {
} else if (coord_type == 1) {

// r-z
Real rc = (i + 0.5_rt) * dx[0];
Real rm = (i - 1 + 0.5_rt) * dx[0];
Real rp = (i + 1 + 0.5_rt) * dx[0];
// r-z
Real rc = (i + 0.5_rt) * dx[0];
Real rm = (i - 1 + 0.5_rt) * dx[0];
Real rp = (i + 1 + 0.5_rt) * dx[0];

#if (AMREX_SPACEDIM == 1)
div_u += 0.5_rt * (rp * q_arr(i+1,j,k,QU) - rm * q_arr(i-1,j,k,QU)) / (rc * dx[0]);
#endif
du_x = 0.5_rt * (rp * q_arr(i+1,j,k,QU) -
rm * q_arr(i-1,j,k,QU)) / (rc * dx[0]);
#if (AMREX_SPACEDIM == 2)
div_u += 0.5_rt * (rp * q_arr(i+1,j,k,QU) - rm * q_arr(i-1,j,k,QU)) / (rc * dx[0]) +
0.5_rt * (q_arr(i,j+1,k,QV) - q_arr(i,j-1,k,QV)) * dyinv;
dv_y = 0.5_rt * (q_arr(i,j+1,k,QV) - q_arr(i,j-1,k,QV)) * dyinv;
#endif
#endif

#if AMREX_SPACEDIM == 1
} else if (coord_type == 2) {
} else if (coord_type == 2) {

// 1-d spherical
Real rc = (i + 0.5_rt) * dx[0];
Real rm = (i - 1 + 0.5_rt) * dx[0];
Real rp = (i + 1 + 0.5_rt) * dx[0];
// 1-d spherical
Real rc = (i + 0.5_rt) * dx[0];
Real rm = (i - 1 + 0.5_rt) * dx[0];
Real rp = (i + 1 + 0.5_rt) * dx[0];

div_u += 0.5_rt * (rp * rp * q_arr(i+1,j,k,QU) - rm * rm * q_arr(i-1,j,k,QU)) / (rc * rc * dx[0]);
du_x = 0.5_rt * (rp * rp * q_arr(i+1,j,k,QU) -
rm * rm * q_arr(i-1,j,k,QU)) / (rc * rc * dx[0]);
#endif

#ifndef AMREX_USE_GPU

} else {
amrex::Error("ERROR: invalid coord_type in shock");
} else {
amrex::Error("ERROR: invalid coord_type in shock");
#endif
}
}

div_u = du_x + dv_y + dw_z;

// 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));
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));
// 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));
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;
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));
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));
}
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));
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));
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);
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));
Real cdu_x = std::min(du_x, 0.0);
Real cdv_y = std::min(dv_y, 0.0);
Real cdw_z = std::min(dw_z, 0.0);

Real gradPdx_over_P{0.0_rt};
if (vel != 0.0) {
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);
Real divu_mag = std::sqrt(cdu_x * cdu_x +
cdv_y * cdv_y +
cdw_z * cdw_z);

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;
}
Real gradPdx_over_P{0.0_rt};
if (divu_mag != 0.0) {
gradPdx_over_P = std::abs(std::abs(dP_x) * cdu_x +
std::abs(dP_y) * cdv_y +
std::abs(dP_z) * cdw_z) / divu_mag;
}
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;
} else {
shk(i,j,k) = 0.0_rt;
}
});

}
Expand Down