From e65b3b9c2b50c38d0c2ed381516a24e6867df9cd Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sun, 17 Dec 2023 14:21:24 -0500 Subject: [PATCH] add a no-limiting option to PPM (#2670) This is useful for convergence testing. --- Source/driver/_cpp_parameters | 4 +- Source/hydro/ppm.H | 117 +++++++++++++++++++--------------- 2 files changed, 67 insertions(+), 54 deletions(-) diff --git a/Source/driver/_cpp_parameters b/Source/driver/_cpp_parameters index 774d6916ac..4d1dde2581 100644 --- a/Source/driver/_cpp_parameters +++ b/Source/driver/_cpp_parameters @@ -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 diff --git a/Source/hydro/ppm.H b/Source/hydro/ppm.H index 9070a48bc1..1fe14199fb 100644 --- a/Source/hydro/ppm.H +++ b/Source/hydro/ppm.H @@ -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]); + + } }