Skip to content

Commit e61113d

Browse files
BenWibkingpre-commit-ci[bot]
BenWibking
and
pre-commit-ci[bot]
authored
Fix PLM convergence (#1552)
### Description There was a factor-of-two bug in the PLM slopes that preventing PLM convergence at order 2. I have also switched to use the Sweby $\beta=1.5$ limiter, which is sharper than minmod. (Minmod only converges at order ~1.5 due to severe extrema clipping. MC converges at almost exactly order 2.) It now converges at nearly order 2, as it should: ``` Running Richardson convergence test for HydroWave: Resolution Error Norm ---------- ---------- 128 5.680154e-09 256 1.544953e-09 512 4.469956e-10 1024 1.255872e-10 2048 3.500900e-11 Reached maximum resolution (nx = 2048) without achieving the target error 2.000e-11 Convergence Rate Analysis: Resolution Pair Observed Rate Expected Rate --------------- ------------- ------------- 128 -> 256 1.88 2.0 256 -> 512 1.79 2.0 512 -> 1024 1.83 2.0 1024 -> 2048 1.84 2.0 Overall convergence rate: 1.84 Expected rate: 2.0 ``` <img width="1200" height="800" alt="hydro_wave_convergence_limiters" src="https://github.com/user-attachments/assets/361acc88-db5d-4682-a9f7-09dd64948b77" /> ### Related issues N/A ### Checklist _Before this pull request can be reviewed, all of these tasks should be completed. Denote completed tasks with an `x` inside the square brackets `[ ]` in the Markdown source below:_ - [x] I have added a description (see above). - [x] I have added a link to any related issues (if applicable; see above). - [x] I have read the [Contributing Guide](https://github.com/quokka-astro/quokka/blob/development/CONTRIBUTING.md). - [ ] I have added tests for any new physics that this PR adds to the code. - [ ] *(For quokka-astro org members)* I have manually triggered the GPU tests with the magic comment `/azp run`. --------- Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
1 parent fa6db2b commit e61113d

4 files changed

Lines changed: 96 additions & 63 deletions

File tree

‎scripts/python/hydro_wave_convergence_analysis.py‎

Lines changed: 60 additions & 44 deletions
Original file line numberDiff line numberDiff line change
@@ -32,6 +32,13 @@ class FitResult:
3232
sse: float
3333

3434

35+
@dataclass
36+
class ConvergenceRun:
37+
label: str
38+
data: Sequence[ConvergenceDatum]
39+
fit: FitResult
40+
41+
3542
def read_convergence_csv(path: Path) -> List[ConvergenceDatum]:
3643
if not path.exists():
3744
raise FileNotFoundError(f"Could not find convergence file: {path}")
@@ -161,39 +168,32 @@ def geometric_space(x_min: float, x_max: float, num: int) -> List[float]:
161168
return [math.exp(log_min + (log_max - log_min) * i / (num - 1)) for i in range(num)]
162169

163170

164-
def make_plot(
165-
data: Sequence[ConvergenceDatum],
166-
fit: FitResult,
167-
output_path: Path,
168-
) -> None:
171+
def make_plot(runs: Sequence[ConvergenceRun], output_path: Path) -> None:
169172
import matplotlib
170173

171174
matplotlib.use("Agg")
172175

173176
import matplotlib.pyplot as plt
174177

175-
nx_values = [item.nx for item in data]
176-
error_values = [item.error for item in data]
177-
178-
nx_min = min(nx_values)
179-
nx_max = max(nx_values)
180-
nx_plot = geometric_space(nx_min, nx_max, 200)
181-
182-
model_values = []
183-
asymptotic_values = []
184-
for nx in nx_plot:
185-
asymptotic = fit.amplitude * (nx ** (-fit.order))
186-
model_values.append(math.sqrt(asymptotic * asymptotic + fit.floor * fit.floor))
187-
asymptotic_values.append(asymptotic)
188-
189178
fig, ax = plt.subplots(figsize=(6.0, 4.0))
190-
ax.loglog(nx_values, error_values, "o", label="Simulation data")
191-
ax.loglog(nx_plot, model_values, "-", label="Power law + floor fit")
192-
ax.loglog(nx_plot, asymptotic_values, "--", label="Asymptotic power law")
193-
194-
x_lo, x_hi = ax.get_xlim()
195-
ax.axhline(fit.floor, color="gray", linestyle="--", label="Roundoff floor")
196-
ax.set_xlim(x_lo, x_hi)
179+
color_cycle = plt.rcParams["axes.prop_cycle"].by_key()["color"]
180+
for run_index, run in enumerate(runs):
181+
color = color_cycle[run_index % len(color_cycle)]
182+
nx_values = [item.nx for item in run.data]
183+
error_values = [item.error for item in run.data]
184+
185+
nx_plot = geometric_space(min(nx_values), max(nx_values), 200)
186+
model_values = []
187+
asymptotic_values = []
188+
for nx in nx_plot:
189+
asymptotic = run.fit.amplitude * (nx ** (-run.fit.order))
190+
model_values.append(math.sqrt(asymptotic * asymptotic + run.fit.floor * run.fit.floor))
191+
asymptotic_values.append(asymptotic)
192+
193+
ax.loglog(nx_values, error_values, "o", color=color, label="_nolegend_")
194+
ax.loglog(nx_plot, model_values, "-", color=color, label=f"{run.label} fit")
195+
ax.loglog(nx_plot, asymptotic_values, "--", color=color, label="_nolegend_")
196+
ax.axhline(run.fit.floor, color=color, linestyle=":", alpha=0.5, label="_nolegend_")
197197

198198
ax.set_xlabel("Cells per wavelength (N_x)")
199199
ax.set_ylabel("L1 error norm")
@@ -211,8 +211,15 @@ def main() -> None:
211211
parser.add_argument(
212212
"--csv",
213213
type=Path,
214-
default=Path("tests") / "hydro_wave_convergence.csv",
215-
help="Path to the hydro wave convergence CSV file (default: tests/hydro_wave_convergence.csv).",
214+
action="append",
215+
default=[],
216+
help="Path to a hydro wave convergence CSV file (repeatable).",
217+
)
218+
parser.add_argument(
219+
"--label",
220+
action="append",
221+
default=[],
222+
help="Optional label for a CSV file (repeatable, matches --csv order).",
216223
)
217224
parser.add_argument(
218225
"--output",
@@ -227,29 +234,38 @@ def main() -> None:
227234
)
228235
args = parser.parse_args()
229236

230-
data = read_convergence_csv(args.csv)
231-
fit = fit_power_law_with_floor(data)
232-
orders = compute_observed_orders(data)
237+
csv_paths = args.csv or [Path("tests") / "hydro_wave_convergence.csv"]
238+
if args.label and len(args.label) != len(csv_paths):
239+
raise ValueError("Provide the same number of --label entries as --csv entries.")
240+
241+
labels = args.label or [path.stem for path in csv_paths]
242+
runs: List[ConvergenceRun] = []
243+
244+
for label, csv_path in zip(labels, csv_paths):
245+
data = read_convergence_csv(csv_path)
246+
fit = fit_power_law_with_floor(data)
247+
orders = compute_observed_orders(data)
248+
runs.append(ConvergenceRun(label=label, data=data, fit=fit))
233249

234-
print("Hydro wave convergence statistics:\n")
235-
print(f"{'nx':>8} {'dx':>14} {'error':>16}")
236-
for row in data:
237-
print(f"{row.nx:8d} {row.dx:14.6e} {row.error:16.6e}")
250+
print(f"Hydro wave convergence statistics for {label} ({csv_path}):\n")
251+
print(f"{'nx':>8} {'dx':>14} {'error':>16}")
252+
for row in data:
253+
print(f"{row.nx:8d} {row.dx:14.6e} {row.error:16.6e}")
238254

239-
print("\nObserved refinement ratios and orders (successive levels):")
240-
print(f"{'nx':>8} {'ratio':>10} {'order':>10}")
241-
for nx, ratio, order in orders:
242-
print(f"{nx:8d} {ratio:10.3f} {order:10.3f}")
255+
print("\nObserved refinement ratios and orders (successive levels):")
256+
print(f"{'nx':>8} {'ratio':>10} {'order':>10}")
257+
for nx, ratio, order in orders:
258+
print(f"{nx:8d} {ratio:10.3f} {order:10.3f}")
243259

244-
print("\nPower-law + floor fit parameters:")
245-
print(f" amplitude (A): {fit.amplitude:.6e}")
246-
print(f" order (p): {fit.order:.3f}")
247-
print(f" floor (f): {fit.floor:.6e}")
260+
print("\nPower-law + floor fit parameters:")
261+
print(f" amplitude (A): {fit.amplitude:.6e}")
262+
print(f" order (p): {fit.order:.3f}")
263+
print(f" floor (f): {fit.floor:.6e}\n")
248264

249265
if args.no_plot:
250266
print("\nSkipping plot generation (per --no-plot).")
251267
else:
252-
make_plot(data, fit, args.output)
268+
make_plot(runs, args.output)
253269
print(f"\nSaved log-log plot to {args.output}")
254270

255271

‎src/QuokkaSimulation.hpp‎

Lines changed: 5 additions & 5 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,8 +3151,8 @@ 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(),
3155-
x1ReconstructRange, nvars);
3154+
HyperbolicSystem<problem_t>::template ReconstructStatesPLM<DIR, SlopeLimiter::sweby>(primVar.array(), x1LeftState.array(), x1RightState.array(),
3155+
x1ReconstructRange, nvars);
31563156
} else if (radiationReconstructionOrder_ == 1) {
31573157
HyperbolicSystem<problem_t>::template ReconstructStatesConstant<DIR>(primVar.array(), x1LeftState.array(), x1RightState.array(),
31583158
x1ReconstructRange, nvars);

‎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 & 11 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 class 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,42 @@ 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

59-
[[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE static auto MC(double a, double b) -> double
60-
{
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)});
62-
}
63+
[[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE static auto MC(double a, double b) -> double { return Sweby(a, b, 2.0); }
6364

6465
[[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE static auto minmod(double a, double b) -> double
6566
{
6667
return 0.5 * (sgn(a) + sgn(b)) * std::min(std::abs(a), std::abs(b));
6768
}
6869

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

258-
// Unlike PPM, PLM with MC or minmod limiters is TVD.
275+
// Unlike PPM, PLM with Sweby limiters (sigma=1 minmod, sigma=1.5 Sweby, sigma=2 MC) is TVD.
259276
// (There are no spurious oscillations, *except* in the slow-moving shock problem,
260277
// which can produce unphysical oscillations even when using upwind Godunov fluxes.)
261278
// However, most tests fail when using PLM reconstruction because
@@ -273,8 +290,8 @@ HyperbolicSystem<problem_t>::ReconstructStatesPLM(quokka::Array4View<amrex::Real
273290
// (This converges at second order in spatial resolution.)
274291
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));
275292
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
293+
leftState(i, j, k, n) = q(i - 1, j, k, n) + 0.5 * lslope; // NOLINT
294+
rightState(i, j, k, n) = q(i, j, k, n) - 0.5 * rslope; // NOLINT
278295
}
279296

280297
template <typename problem_t>

0 commit comments

Comments
 (0)