Add 2D axisymmetric WENO5 convergence test; fix single-fluid alpha drift; clarify HLL gsrc comment#1411
Open
sbryngelson wants to merge 35 commits intoMFlowCode:masterfrom
Open
Add 2D axisymmetric WENO5 convergence test; fix single-fluid alpha drift; clarify HLL gsrc comment#1411sbryngelson wants to merge 35 commits intoMFlowCode:masterfrom
sbryngelson wants to merge 35 commits intoMFlowCode:masterfrom
Conversation
With eps=5 the prim->cons covariance error O(eps^3 h^2) swamps WENO5's O(eps^2 h^5) error at all practical resolutions, giving apparent rate 2. Using eps=0.01 shifts the crossover to h~0.22 so that N=16..64 (h=0.63..0.16) are firmly in the WENO5-dominated regime and show rate ~5, as in Johnsen & Colonius (2006) who similarly used a weak vortex. - hcid=283: vortex strength now read from patch_icpp(patch_id)%epsilon (defaults to 5 if not set, preserving backward compatibility) - case.py: eps_vortex=0.01, c_max updated accordingly - run_convergence.py: resolutions 16,32,64; WENO5 tolerance 4.0 (rate>=4) - convergence.yml: updated to 16 32 64
CFL=0.02 eliminates RK3 temporal error so WENO5 shows rate 5. WENO3 expected rate is 2 (not 3) on smooth-extremum problems per Henrick et al. 2005 — JS smoothness indicators degrade at f'=0. The 2D vortex job remains the authority on WENO3 rate 3.
1D runner now uses examples/1D_euler_convergence/case.py (single-fluid,
density sine wave, no non-conservative alpha equation) instead of the
two-fluid advection case. The five-equation model non-conservative alpha
transport degrades unlimited MUSCL to ~1st order on two-fluid problems
even with identical EOS; the single-fluid Euler case gives clean rates.
Add muscl_lim=0 (unlimited central-difference) to MUSCL reconstruction:
- src/simulation/m_muscl.fpp: central slope = 0.5*(slopeL+slopeR)
- toolchain/mfc/params/definitions.py: add 0 to choices
- toolchain/mfc/case_validator.py: allow muscl_lim in {0..5}
1D and 2D convergence cases default to muscl_lim=0 for smooth problems
where TVD limiters would stall at smooth extrema. Runner default CFL=0.02
so RK3 temporal error O(h^3) stays negligible for WENO5 rate-5 verification.
N=16 (m=n=15) is below the minimum grid size for WENO5's 5-point stencil and causes pre_process to abort. Start all 2D tests at N=32. WENO3 on the stationary isentropic vortex is pre-asymptotic at N=32-128: fitted rate ~2.02, local rate 1.83->2.21 (converging toward 3). Lower threshold from 2.2 to 1.8 (expected=3, tol=1.2). The rate still clearly exceeds the 1D result (1.77, Jiang-Shu smooth-extremum degradation), which is the distinction this test is designed to demonstrate.
…d to 0.95 Default and CI resolutions: 128 256 512 1024. WENO5 is capped at N=512 per-scheme (hits double-precision floor at N=1024 after 111k steps; error ~2.6e-12 vs 4.2e-12 at N=512, rate collapses to 0.69). N=128-512 gives fitted rate 4.99. WENO1 threshold tightened from 0.6 to 0.95 (tol 0.05): pairwise rates 0.95->0.97->0.99 at N=128-1024, fitted 0.97. MUSCL2 shows exact rate 2.00 at all doublings N=128-1024.
WENO5 capped at N=512 (machine-precision floor kills rate at N=1024, error ~2.6e-12, rate collapses to 0.69). WENO3 starts at N=256 to skip coarse pre-asymptotic points. WENO1 and MUSCL2 use full [128,1024] range. Thresholds: WENO5>=4.8, WENO3>=1.8, WENO1>=0.95, MUSCL2>=1.9.
Tests are too slow (~30-45 min each job) to run on every PR push. Manual trigger via workflow_dispatch lets you run them deliberately when verifying numerics, not on every commit.
Both 1D and 2D runners now test all 7 schemes: WENO5/3/1, MUSCL2, TENO5, WENO7, TENO7. Both case.py files gain --teno and --teno-ct flags. Resolution bounds: 1D: TENO5 [128,512], WENO7/TENO7 [128,256] (precision floor near N=512) 2D: TENO5 [32,128], WENO7/TENO7 [64,128] (MFC stencil constraint: N>=35) TENO CT values match the Shu-Osher examples: 1e-6 for order-5, 1e-9 for order-7. weno_eps tightened to 1e-40 (was 1e-16) for all WENO/TENO schemes in case.py.
Scripts import numpy which lives in build/venv, not the system Python.
Removed --no-build from run_case() in both runners. The generic ./mfc.sh build created a binary without the analytical IC (case.fpp), so --no-build would silently fall back to that binary and crash with SIGILL when the IC code wasn't found. MFC's hash system rebuilds only when case.fpp changes, so each scheme triggers one build then reuses the binary for all N values. Also drop WENO7/TENO7 tolerance to 4.0 (threshold >=3.0): local runs confirm rate ~3.7 due to RK3 temporal and spatial h^7 errors being comparable at N=128-256.
With CFL=0.02 (1D) or CFL=0.4 (2D) and RK3 time integration, the temporal error O(dt^3) is comparable to the O(h^7) spatial error at N=128-256, giving a spurious fitted rate of ~3.7 instead of 7. Fix: bake CFL=0.005 directly into WENO7/TENO7 extra_args. This drops the temporal error by (0.005/0.02)^3=1/64 in 1D and (0.005/0.4)^3=1/51200 in 2D, making spatial error dominate. All other schemes keep CFL=0.02 (1D) or 0.4 (2D). Threshold updated from >=3.0 to >=6.5 (tol 0.5) in 1D and >=6.0 (tol 1.0) in 2D. Also adds --cfl arg to 2D case.py so the runner can override it per scheme.
1D: use N=64,128 (not N=128,256). At N=256 the WENO7 spatial error (~1.7e-14) falls below the round-off accumulation floor (~2.5e-12 for ~28M cell-steps), giving a spurious rate of -0.21. At N=64 the spatial error is ~2.8e-10, well above the floor, and the rate between N=64 and N=128 is 6.97. Add N=64 to the CI resolution list accordingly. 2D: remove WENO7/TENO7 from the 2D test. The isentropic vortex has a primitive->conserved covariance floor O(eps^3*h^2). For WENO7 scheme error O(eps^2*h^7) to dominate requires h > eps^(1/5) = 0.40 (N < 25). At N=64-128 the covariance floor gives rate ~2, masking the 7th-order spatial rate regardless of CFL. The 1D test cleanly verifies WENO7/TENO7 7th-order convergence on a pure advection problem without this floor.
Threads num_ranks through run_case() and test_scheme() in both runners. GitHub ubuntu-latest runners have 4 vCPUs; using all 4 reduces simulation time for schemes with large Nt (e.g. WENO7/TENO7 at 27k-56k steps).
read_vf1_1d and read_cons_vf1 were reading only p_all/p0/, so with
num_ranks>1 only the first rank's N/num_ranks cells were included in
the L2 error. Now both readers iterate over p0..p{n-1} and concatenate.
For the L2 norm this is correct regardless of spatial ordering across
ranks. Validated: 4-rank WENO7 |vf0|=64 (not 16) and error matches
single-rank result exactly.
Tests WENO1/3/5/7, MUSCL-minmod/MC/VanLeer/SUPERBEE, TENO5/7 on the 1D Sod problem (shock + contact + rarefaction). Self-convergence: compare N vs 2N via cell-averaging, no exact Riemann solver needed. Rate ~1 for all schemes by Godunov theorem; WENO1 and SUPERBEE use wider tolerance (0.5) due to contact-dominated O(h^0.5) L1 diffusion at these resolutions. CI runs at 64/128/256/512 with 4 MPI ranks.
… from conservation check WENO5/TENO5 require (N/2) >= num_stcls_min*weno_order = 25 cells/rank in a 2x2 decomposition; N=32 only gives 16 cells/rank which fails MFC's decomposition check. Setting min_N=64 (32 cells/rank) fixes this. The isentropic vortex has zero net linear momentum by symmetry, making the relative conservation error formula ill-conditioned (divides near-zero by 1e-300). Density and energy have large nonzero integrals and are the meaningful conserved quantities to check.
MUSCL-SUPERBEE's pre-asymptotic L1 rate at N=64 is ~0.40 (contact-dominated), which pulls the fitted rate below the 0.5 threshold. Setting min_N=128 skips that pre-asymptotic point; the asymptotic rates (N=128,256,512) are ~0.56.
Previously triggered on every push to every branch, running expensive 1D/2D/temporal/Sod jobs (60-90 min total) on every dev push. Now triggers on: push to master, any pull_request, or workflow_dispatch.
…e, add traceback, validate output sizes, fix silent PASS on zero data
- Add examples/2D_axisym_convergence/case.py: 2D cylindrical density sine wave convergence case with HLLC (riemann_solver=2) and periodic x-BCs so the exact solution is valid at every cell. - Add toolchain/mfc/test/run_convergence_axisym.py: convergence runner that refines Nz=[64,128,256] with Nr=32 fixed, reads alpha_rho from p_all binary, compares r-averaged density to the analytic solution, and asserts fitted WENO5 rate >= 4.8 (observed: 5.08). - Fix HLL and LF Riemann solvers (m_riemann_solvers.fpp): void-fraction geometric source was incorrectly set to flux_rs instead of 0. - Fix single-fluid cyl_coord alpha drift (m_time_steppers.fpp): after each RK3 stage, reset alpha=1 for num_fluids=1 + cyl_coord cases to prevent pressure NaN from contact-wave-speed variation across faces. - Add convergence-2d-axisym CI job to .github/workflows/convergence.yml.
Contributor
ⓘ You've reached your Qodo monthly free-tier limit. Reviews pause until next month — upgrade your plan to continue now, or link your paid account if you already have one. |
Claude Code ReviewHead SHA: 743016a Files changed:
Findings
|
…#800 intent The previous commit incorrectly set flux_gsrc(adv) = 0 for HLL and LF, which would have reverted PR MFlowCode#800. Per the analysis in issue MFlowCode#1402 comment by ChrisZYJ: flux_rs(adv) in HLL/LF is the numerical diffusion component (not transport), and any face-flux in cylindrical coordinates requires the full geometric divergence (1/r)*d(rF)/dr = dF/dr + F/r. The flux_gsrc term provides the F/r correction for the diffusion. The transport half is handled separately via flux_src_rs (non-conservative source). PR MFlowCode#800 was correct. The comment is updated to document this so future readers do not repeat the mistake. The convergence test is unaffected: for num_fluids=1, alpha=1 everywhere so flux_rs(adv) = 0 at all faces regardless.
sbryngelson
changed the title
Member
Author
|
@ChrisZYJ give this a look. I found a bug in axisymm after seeing your comment on my other (now closed) PR. I think. |
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## master #1411 +/- ##
==========================================
+ Coverage 64.90% 64.94% +0.04%
==========================================
Files 72 72
Lines 18874 18870 -4
Branches 1571 1571
==========================================
+ Hits 12250 12256 +6
+ Misses 5649 5639 -10
Partials 975 975 ☔ View full report in Codecov by Sentry. 🚀 New features to boost your workflow:
|
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
While developing a 2D axisymmetric (
cyl_coord=T,num_fluids=1) convergence test, one bug in the simulation source was found and fixed, and a misleading comment in the HLL/LF solver was corrected. The convergence test (WENO5, fitted rate 5.08) is included as a new CI job.Bug — Single-fluid
cyl_coord: alpha drifts each RK stageFor
num_fluids=1+cyl_coord=T, α is trivially 1 everywhere but the HLLC contact-wave speed varies slightly across radial faces, producing a tiny net flux in α over RK sub-stages. This eventually causes a pressure NaN when computingp = (γ−1)(E − ρ|u|²/2) / α. Fix: resetq_cons_ts(1)%vf(eqn_idx%adv%beg)%sf = 1._wpafter theadv_ncall in each RK3 stage (guarded bynum_fluids == 1 .and. cyl_coord .and. .not. bubbles_euler).Comment correction — HLL/LF void-fraction geometric source
A previous version of this PR incorrectly changed
flux_gsrc(adv) = flux_rs→= 0for HLL and LF, which would have reverted PR #800. This has been corrected. Per the analysis in issue #1402 (comment by ChrisZYJ):flux_rs(adv)in HLL/LF is the numerical diffusion component (not the transport flux — that lives influx_src_rs), and in cylindrical coordinates any face-flux requires the full geometric divergence(1/r)∂(rF)/∂r = ∂F/∂r + F/r. Theflux_gsrc = flux_rsassignment provides theF/rcorrection for the diffusion. PR #800 was correct; only the comment is updated to document why.New convergence case and CI job
examples/2D_axisym_convergence/case.py: density sine waveρ=1+0.2·sin(2πx),u_x=1,u_r=0,p=1in a[0,5]×[0,0.5]cylindrical domain. The domain contains exactly 5 full sine periods, so periodic x-BCs makeρ(x,T)=1+0.2·sin(2π(x−T))an exact solution everywhere. This is an exact nonlinear solution to the cylindrical Euler equations: it is a pure entropy wave — pressure stays uniform atp=1for all amplitudes becausep=(γ−1)(E−ρ/2)=1identically, and the r-momentum geometric source cancels exactly ((1/r)∂(rp)/∂r − p/r = 0for uniform p).toolchain/mfc/test/run_convergence_axisym.py: refinesNz∈{64,128,256}withNr=32fixed, readsq_cons_vf1fromp_all, averages over r, and asserts fitted WENO5 rate ≥ 4.8. Starting at Nz=64 (12.8 cells/wavelength) avoids the pre-asymptotic regime where WENO-JS underperforms near smooth extrema..github/workflows/convergence.yml: newconvergence-2d-axisymjob.Observed rates:
Fitted rate: 5.08 (need ≥ 4.8) ✓
Test plan
./mfc.sh precheck -j 8passes (all 6 checks)./mfc.sh build -j 8compiles cleanpython toolchain/mfc/test/run_convergence_axisym.py→ PASS, rate 5.08./mfc.sh test -j 8(no golden files changed)