Skip to content

Commit 30273f0

Browse files
BenWibkingBenWibking
committed
fix PLM convergence
1 parent fa6db2b commit 30273f0

3 files changed

Lines changed: 35 additions & 15 deletions

File tree

‎src/QuokkaSimulation.hpp‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -2652,10 +2652,10 @@ void QuokkaSimulation<problem_t>::hydroFluxFunction(amrex::MultiFab &primVar_mf,
26522652
ng_reconstruct, 2);
26532653
}
26542654
} else if (reconstructionOrder_ == 2) {
2655-
HyperbolicSystem<problem_t>::template ReconstructStatesPLM<DIR, SlopeLimiter::minmod>(primVar_mf, leftState, rightState, ng_reconstruct, nvars);
2655+
HyperbolicSystem<problem_t>::template ReconstructStatesPLM<DIR, SlopeLimiter::sweby>(primVar_mf, leftState, rightState, ng_reconstruct, nvars);
26562656
if constexpr (Physics_Traits<problem_t>::is_mhd_enabled) {
2657-
HyperbolicSystem<problem_t>::template ReconstructStatesPLM<DIR, SlopeLimiter::minmod>(cc_bfield_perp_comps_mf, leftState_bfield,
2658-
rightState_bfield, ng_reconstruct, 2);
2657+
HyperbolicSystem<problem_t>::template ReconstructStatesPLM<DIR, SlopeLimiter::sweby>(cc_bfield_perp_comps_mf, leftState_bfield,
2658+
rightState_bfield, ng_reconstruct, 2);
26592659
}
26602660
} else if (reconstructionOrder_ == 1) {
26612661
HyperbolicSystem<problem_t>::template ReconstructStatesConstant<DIR>(primVar_mf, leftState, rightState, ng_reconstruct, nvars);
@@ -3151,7 +3151,7 @@ void QuokkaSimulation<problem_t>::fluxFunction(amrex::Array4<const amrex::Real>
31513151
x1ReconstructRange, nvars);
31523152
} else if (radiationReconstructionOrder_ == 2) {
31533153
// PLM and donor cell are interface-centered kernels
3154-
HyperbolicSystem<problem_t>::template ReconstructStatesPLM<DIR, SlopeLimiter::MC>(primVar.array(), x1LeftState.array(), x1RightState.array(),
3154+
HyperbolicSystem<problem_t>::template ReconstructStatesPLM<DIR, SlopeLimiter::sweby>(primVar.array(), x1LeftState.array(), x1RightState.array(),
31553155
x1ReconstructRange, nvars);
31563156
} else if (radiationReconstructionOrder_ == 1) {
31573157
HyperbolicSystem<problem_t>::template ReconstructStatesConstant<DIR>(primVar.array(), x1LeftState.array(), x1RightState.array(),

‎src/hydro/mhd_system.hpp‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -911,13 +911,13 @@ void MHDSystem<problem_t>::ReconstructTo(FluxDir dir, arrayconst_t &cState, arra
911911
} else if (reconstructionOrder == 2) {
912912
switch (dir) {
913913
case FluxDir::X1:
914-
MHDSystem<problem_t>::template ReconstructStatesPLM<FluxDir::X1, SlopeLimiter::minmod>(cState, lState, rState, box_r, 1);
914+
MHDSystem<problem_t>::template ReconstructStatesPLM<FluxDir::X1, SlopeLimiter::sweby>(cState, lState, rState, box_r, 1);
915915
break;
916916
case FluxDir::X2:
917-
MHDSystem<problem_t>::template ReconstructStatesPLM<FluxDir::X2, SlopeLimiter::minmod>(cState, lState, rState, box_r, 1);
917+
MHDSystem<problem_t>::template ReconstructStatesPLM<FluxDir::X2, SlopeLimiter::sweby>(cState, lState, rState, box_r, 1);
918918
break;
919919
case FluxDir::X3:
920-
MHDSystem<problem_t>::template ReconstructStatesPLM<FluxDir::X3, SlopeLimiter::minmod>(cState, lState, rState, box_r, 1);
920+
MHDSystem<problem_t>::template ReconstructStatesPLM<FluxDir::X3, SlopeLimiter::sweby>(cState, lState, rState, box_r, 1);
921921
break;
922922
}
923923
} else if (reconstructionOrder == 1) {

‎src/hyperbolic_system.hpp‎

Lines changed: 28 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -36,7 +36,7 @@ enum redoFlag { none = 0, redo = 1 };
3636
} // namespace quokka
3737

3838
// Define enum for slope limiter type
39-
enum SlopeLimiter { minmod = 0, MC };
39+
enum SlopeLimiter { minmod = 0, sweby, MC };
4040

4141
using array_t = amrex::Array4<amrex::Real> const;
4242
using arrayconst_t = amrex::Array4<const amrex::Real> const;
@@ -47,25 +47,45 @@ template <typename problem_t> class HyperbolicSystem
4747
public:
4848
template <SlopeLimiter limiter> AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE static auto SlopeFunc(amrex::Real x, amrex::Real y) -> amrex::Real
4949
{
50-
static_assert(limiter == SlopeLimiter::minmod || limiter == SlopeLimiter::MC, "Invalid slope limiter specified.");
50+
static_assert(limiter == SlopeLimiter::minmod || limiter == SlopeLimiter::sweby || limiter == SlopeLimiter::MC,
51+
"Invalid slope limiter specified.");
5152
if constexpr (limiter == SlopeLimiter::minmod) {
52-
return minmod(x, y);
53+
return Sweby(x, y, 1.0);
54+
}
55+
if constexpr (limiter == SlopeLimiter::sweby) {
56+
return Sweby(x, y, 1.5);
5357
}
5458
if constexpr (limiter == SlopeLimiter::MC) {
55-
return MC(x, y);
59+
return Sweby(x, y, 2.0);
5660
}
5761
}
5862

5963
[[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE static auto MC(double a, double b) -> double
6064
{
61-
return 0.5 * (sgn(a) + sgn(b)) * std::min({0.5 * std::abs(a + b), 2.0 * std::abs(a), 2.0 * std::abs(b)});
65+
return Sweby(a, b, 2.0);
6266
}
6367

6468
[[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE static auto minmod(double a, double b) -> double
6569
{
6670
return 0.5 * (sgn(a) + sgn(b)) * std::min(std::abs(a), std::abs(b));
6771
}
6872

73+
[[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE static auto minmod3(double a, double b, double c) -> double
74+
{
75+
if ((sgn(a) == sgn(b)) && (sgn(a) == sgn(c))) {
76+
return sgn(a) * std::min({std::abs(a), std::abs(b), std::abs(c)});
77+
}
78+
return 0.0;
79+
}
80+
81+
[[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE static auto Sweby(double a, double b, double sigma) -> double
82+
{
83+
const double sigma_a = sigma * a;
84+
const double sigma_b = sigma * b;
85+
const double centered = 0.5 * (a + b);
86+
return minmod3(sigma_a, centered, sigma_b);
87+
}
88+
6989
[[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE static auto median(double a, double b, double c) -> double
7090
{
7191
return std::max(std::min(a, b), std::min(std::max(a, b), c));
@@ -255,7 +275,7 @@ HyperbolicSystem<problem_t>::ReconstructStatesPLM(quokka::Array4View<amrex::Real
255275
// permute array indices according to dir
256276
auto [i, j, k] = quokka::reorderMultiIndex<DIR>(i_in, j_in, k_in);
257277

258-
// Unlike PPM, PLM with MC or minmod limiters is TVD.
278+
// Unlike PPM, PLM with Sweby limiters (sigma=1 midmod, sigma=1.5 Sweby, sigma=2 MC) is TVD.
259279
// (There are no spurious oscillations, *except* in the slow-moving shock problem,
260280
// which can produce unphysical oscillations even when using upwind Godunov fluxes.)
261281
// However, most tests fail when using PLM reconstruction because
@@ -273,8 +293,8 @@ HyperbolicSystem<problem_t>::ReconstructStatesPLM(quokka::Array4View<amrex::Real
273293
// (This converges at second order in spatial resolution.)
274294
const auto lslope = HyperbolicSystem<problem_t>::template SlopeFunc<limiter>(q(i, j, k, n) - q(i - 1, j, k, n), q(i - 1, j, k, n) - q(i - 2, j, k, n));
275295
const auto rslope = HyperbolicSystem<problem_t>::template SlopeFunc<limiter>(q(i + 1, j, k, n) - q(i, j, k, n), q(i, j, k, n) - q(i - 1, j, k, n));
276-
leftState(i, j, k, n) = q(i - 1, j, k, n) + 0.25 * lslope; // NOLINT
277-
rightState(i, j, k, n) = q(i, j, k, n) - 0.25 * rslope; // NOLINT
296+
leftState(i, j, k, n) = q(i - 1, j, k, n) + 0.5 * lslope; // NOLINT
297+
rightState(i, j, k, n) = q(i, j, k, n) - 0.5 * rslope; // NOLINT
278298
}
279299

280300
template <typename problem_t>

0 commit comments

Comments
 (0)