Skip to content

Commit

Permalink
add a no-limiting option to PPM (#2670)
Browse files Browse the repository at this point in the history 
This is useful for convergence testing.
  • Loading branch information
@zingale
zingale authored Dec 17, 2023
1 parent 1613d3d commit e65b3b9
Show file tree
Hide file tree
Showing 2 changed files with 67 additions and 54 deletions.
4 changes: 3 additions & 1 deletion Source/driver/_cpp_parameters
View file
Original file line number Diff line number Diff line change
Expand Up @@ -84,9 +84,11 @@ hybrid_hydro int 0
# reconstruction type:
# 0: piecewise linear;
# 1: classic Colella \& Woodward ppm;
# 2: extrema-preserving ppm
ppm_type int 1

# do we limit the ppm parabola?
ppm_do_limiting int 1

# For MHD + PLM, do we limit on characteristic or primitive variables
mhd_limit_characteristic int 1

Expand Down
117 changes: 64 additions & 53 deletions Source/hydro/ppm.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -57,84 +57,95 @@ ppm_reconstruct(const Real* s,
const Real flatn, Real& sm, Real& sp) {


// Compute van Leer slopes
if (ppm_do_limiting) {

Real dsl = 2.0_rt * (s[im1] - s[im2]);
Real dsr = 2.0_rt * (s[i0] - s[im1]);
// Compute van Leer slopes

Real dsvl_l = 0.0_rt;
if (dsl*dsr > 0.0_rt) {
Real dsc = 0.5_rt * (s[i0] - s[im2]);
dsvl_l = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc),amrex::min(std::abs(dsl),std::abs(dsr)));
}
Real dsl = 2.0_rt * (s[im1] - s[im2]);
Real dsr = 2.0_rt * (s[i0] - s[im1]);

dsl = 2.0_rt * (s[i0] - s[im1]);
dsr = 2.0_rt * (s[ip1] - s[i0]);
Real dsvl_l = 0.0_rt;
if (dsl*dsr > 0.0_rt) {
Real dsc = 0.5_rt * (s[i0] - s[im2]);
dsvl_l = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc),amrex::min(std::abs(dsl),std::abs(dsr)));
}

Real dsvl_r = 0.0_rt;
if (dsl*dsr > 0.0_rt) {
Real dsc = 0.5_rt * (s[ip1] - s[im1]);
dsvl_r = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc),amrex::min(std::abs(dsl),std::abs(dsr)));
}
dsl = 2.0_rt * (s[i0] - s[im1]);
dsr = 2.0_rt * (s[ip1] - s[i0]);

// Interpolate s to edges
Real dsvl_r = 0.0_rt;
if (dsl*dsr > 0.0_rt) {
Real dsc = 0.5_rt * (s[ip1] - s[im1]);
dsvl_r = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc),amrex::min(std::abs(dsl),std::abs(dsr)));
}

sm = 0.5_rt * (s[i0] + s[im1]) - (1.0_rt/6.0_rt) * (dsvl_r - dsvl_l);
// Interpolate s to edges

// Make sure sedge lies in between adjacent cell-centered values
sm = 0.5_rt * (s[i0] + s[im1]) - (1.0_rt/6.0_rt) * (dsvl_r - dsvl_l);

sm = amrex::max(sm, amrex::min(s[i0], s[im1]));
sm = amrex::min(sm, amrex::max(s[i0], s[im1]));
// Make sure sedge lies in between adjacent cell-centered values

sm = amrex::max(sm, amrex::min(s[i0], s[im1]));
sm = amrex::min(sm, amrex::max(s[i0], s[im1]));

// Compute van Leer slopes

dsl = 2.0_rt * (s[i0] - s[im1]);
dsr = 2.0_rt * (s[ip1] - s[i0]);
// Compute van Leer slopes

dsvl_l = 0.0_rt;
if (dsl*dsr > 0.0_rt) {
Real dsc = 0.5_rt * (s[ip1] - s[im1]);
dsvl_l = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc), amrex::min(std::abs(dsl),std::abs(dsr)));
}
dsl = 2.0_rt * (s[i0] - s[im1]);
dsr = 2.0_rt * (s[ip1] - s[i0]);

dsl = 2.0_rt * (s[ip1] - s[i0]);
dsr = 2.0_rt * (s[ip2] - s[ip1]);
dsvl_l = 0.0_rt;
if (dsl*dsr > 0.0_rt) {
Real dsc = 0.5_rt * (s[ip1] - s[im1]);
dsvl_l = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc), amrex::min(std::abs(dsl),std::abs(dsr)));
}

dsvl_r = 0.0_rt;
if (dsl*dsr > 0.0_rt) {
Real dsc = 0.5_rt * (s[ip2] - s[i0]);
dsvl_r = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc), amrex::min(std::abs(dsl),std::abs(dsr)));
}
dsl = 2.0_rt * (s[ip1] - s[i0]);
dsr = 2.0_rt * (s[ip2] - s[ip1]);

// Interpolate s to edges
dsvl_r = 0.0_rt;
if (dsl*dsr > 0.0_rt) {
Real dsc = 0.5_rt * (s[ip2] - s[i0]);
dsvl_r = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc), amrex::min(std::abs(dsl),std::abs(dsr)));
}

sp = 0.5_rt * (s[ip1] + s[i0]) - (1.0_rt/6.0_rt) * (dsvl_r - dsvl_l);
// Interpolate s to edges

// Make sure sedge lies in between adjacent cell-centered values
sp = 0.5_rt * (s[ip1] + s[i0]) - (1.0_rt/6.0_rt) * (dsvl_r - dsvl_l);

sp = amrex::max(sp, amrex::min(s[ip1], s[i0]));
sp = amrex::min(sp, amrex::max(s[ip1], s[i0]));
// Make sure sedge lies in between adjacent cell-centered values

sp = amrex::max(sp, amrex::min(s[ip1], s[i0]));
sp = amrex::min(sp, amrex::max(s[ip1], s[i0]));

// Flatten the parabola

sm = flatn * sm + (1.0_rt - flatn) * s[i0];
sp = flatn * sp + (1.0_rt - flatn) * s[i0];
// Flatten the parabola

// Modify using quadratic limiters -- note this version of the limiting comes
// from Colella and Sekora (2008), not the original PPM paper.
sm = flatn * sm + (1.0_rt - flatn) * s[i0];
sp = flatn * sp + (1.0_rt - flatn) * s[i0];

if ((sp - s[i0]) * (s[i0] - sm) <= 0.0_rt) {
sp = s[i0];
sm = s[i0];
// Modify using quadratic limiters -- note this version of the limiting comes
// from Colella and Sekora (2008), not the original PPM paper.

} else if (std::abs(sp - s[i0]) >= 2.0_rt * std::abs(sm - s[i0])) {
sp = 3.0_rt * s[i0] - 2.0_rt * sm;
if ((sp - s[i0]) * (s[i0] - sm) <= 0.0_rt) {
sp = s[i0];
sm = s[i0];

} else if (std::abs(sm - s[i0]) >= 2.0_rt * std::abs(sp - s[i0])) {
sm = 3.0_rt * s[i0] - 2.0_rt * sp;
}
} else if (std::abs(sp - s[i0]) >= 2.0_rt * std::abs(sm - s[i0])) {
sp = 3.0_rt * s[i0] - 2.0_rt * sm;

} else if (std::abs(sm - s[i0]) >= 2.0_rt * std::abs(sp - s[i0])) {
sm = 3.0_rt * s[i0] - 2.0_rt * sp;
}

} else {

// unlimited PPM reconstruction (Eq. 1.9 in the PPM paper)

sm = (7.0_rt/12.0_rt) * (s[i0] + s[im1]) - (1.0_rt/12.0_rt) * (s[im2] + s[ip1]);
sp = (7.0_rt/12.0_rt) * (s[ip1] + s[i0]) - (1.0_rt/12.0_rt) * (s[im1] + s[ip2]);

}

}

Expand Down

0 comments on commit e65b3b9

Please sign in to comment.