Bound chemistry face states with a scaling limiter (NaN at finite-cylinder rims)

[sbryngelson]

Contribution Policy

All contributions are expected to demonstrate:

• A clear understanding of the codebase
• Alignment with product direction
• Thoughtful reasoning behind changes
• Evidence of real-world usage or hands-on experience with the problem

If these expectations are not met, we would prefer to implement the changes ourselves rather than spend time reviewing low-effort submissions.

---

Fixes #1964

Root cause

This is not a levelset or image-point bug. I checked every ghost point of the repro cylinder offline (same marker, levelset, image-point and stencil logic as master). Every image point lands in fluid, and every interpolation stencil has at least 4 fluid points (cm3c_cases/debug_E/cyl_rim/tools/geom.py).

The failure happens like this:

1. At the cap/side rim, the ghost column under the first side-fluid cell (axial depth 0.5 dx) mirrors different faces. Layer 1 (radial depth 0.16 dx) is nearest the side and mirrors the expansion: p = 1.27e5. Layers 2–4 are nearest the cap and mirror the stagnation region in front of it: p = 2.3–2.8e5.
2. With a hot isothermal wall (thermal_bc = 1, Twall = 1800 K), the ghost densities in that column are 0.14, 0.24, 0.27 and 0.28 kg/m^3, next to fluid at 1.05. The column is non-monotone.
3. WENO picks the smooth all-ghost sub-stencil and extrapolates the inward-rising ghost density out to the fluid/ghost face. The face density undershoots to 0.059, so face T = p·W/(ρR) ≈ 7600 K.
4. Above the NASA-7 fit (3500 K for O2/CO2 in this mechanism), cp_O2/R evaluates to −10. HLLC then gets a NaN sound speed and returns a NaN flux on that face. Only the fluid cell goes NaN, because IB-cell RHS is zeroed.

Instrumented trace from a restart at step 21 of the repro: at stage 3 of step 22, face (12, 17.5, 24) has ρ_R = 0.059 and p_R = 1.31e5, and the y-flux there is NaN. Non-reacting finite rods (no chemistry, adiabatic wall) run cleanly. With chemistry, an adiabatic IB (thermal_bc = 0) or an ambient-temperature wall (Twall = 350 K) also runs cleanly. The hot isothermal wall is the ingredient that drives this to NaN.

A convex corner has no single ghost value that suits both faces. A real fix at the IB level needs a new corner treatment (e.g., per-direction ghost states). This PR bounds the reconstructed face states where the NaN is produced.

#1969 alone does not fix this. #1969 extrapolates the hot-wall ghost from where the image point is actually sampled. On the repro it cools the shallowest rim ghost layer (step 21: 3176 K → 2063 K, ρ 0.137 → 0.193) but leaves layers 2–4 at ~3120 K, so the column is still non-monotone. The NaN is delayed, not removed: the 3D repro NaNs at step 26 instead of 22 (cell (12, 18, 21), the same rim position), and the 2D rectangle case at step 70 instead of 24.

Change

• m_rhs.fpp: new s_bound_face_states(id), run after reconstruction when chemistry is on. It is a Zhang–Shu scaling limiter. Each cell's two face states (all primitive variables) are scaled toward the cell state by one θ ∈ [0, 1]. θ is the largest value that keeps face ρ, p and R_mix above 1e-12 of the cell value and face T = p/(ρ R_mix) ≤ T_surface_max. T_surface_max is 5000 K, the thermo window m_ibm already uses. ρ, p and R_mix are all linear in θ. So the positivity bound is closed form, and the T bound is the first root of a quadratic. That makes θ continuous in the data, with no iteration. θ = 1 (no change) wherever the faces are already in bounds. There is no lower T bound: an earlier 200 K bound moved 1D -> Chemistry -> MultiComponent Diffusion. The range and the GPU loop style follow the species-bounds limiter (Keep species mass fractions in [0, 1] by construction instead of clipping #1963), so it runs over the range the reconstruction filled, with θ shared by both faces of a cell.
• This replaces the previous version of this PR, which reverted an offending face to first order. That was a hard switch on one face; this is bounded, continuous and keeps the reconstruction wherever it is admissible.
• New test 2D -> Chemistry -> IBM -> Hot Wall Rectangle Corners (4646053A). It uses h2o2.yaml, 72x48 cells, 30 steps, ~3 s on CPU. It is the 2D analog: Mach 1.5 air onto an 1800 K isothermal IB rectangle. On master it aborts at step 30 (ICFL is NaN).

Evidence (CPU, CCE, serial --no-mpi)

• 3D repro (x-aligned finite cylinder, D = 2e-4, L = 3D, 12 cells/D, 72x48x48, Mach 1.5 axial air, viscous, no-slip, chemistry with reactions off, surface reactions off, Twall = 1800 K, fixed dt = 5.31e-9; run in two restart segments of 150 steps):
  • RED on master: NaN at step 22 at cell (12, 17, 24), the first fluid cell outside the side, 0.5 dx behind the upstream cap.
  • GREEN with this PR: all 300 steps, exit 0, no NaN. The limiter acts on 154 distinct cells (144 IB ghost cells, 10 fluid cells), all at the rims. Minimum θ is 0.11; unlimited face T reached 7e6 K in ghost cells and 1.5e4 K in fluid cells. Near the rim, the hottest fluid cell at step 300 is 2553 K. The previous first-order revert gave 2852 K there.
• 2D rectangle, 300 steps (same setup as the test): master NaNs at step 24; this PR runs all 300 steps. The limiter acts on 20 cells, all at the corners (19 IB ghost cells, 1 fluid cell); minimum θ is 0.27. In the 30-step test it acts on 14 ghost cells only.
• ./mfc.sh test --no-mpi --no-gpu:
  • -o IBM: 62 passed. A69D2D28 failed once with a transient No such file or directory: '0' and passed on rerun.
  • -o Cylinder: 8 passed.
  • -o Chemistry: 23 passed.
• Goldens regenerated:
  • 4646053A: the new test. Its golden was made with the revert and differs from this limiter's result.
  • B317404C (2D -> IBM -> Vieille Burn Rate) and E5B66084 (2D -> Example -> ibm_burning_grain): the same case. On master, a blowing ghost cell inside the grain reconstructs a negative face density (ρ_face = −0.059 next to a cell ρ = 0.020, so face T < 0). HLLC is fed that state on master. The revert only touched that face, and the change stayed under the 1e-3 test tolerance. The limiter scales both faces of that cell (θ = 0.20). Fluid momentum moves by up to 3.7% relative where |value| > 1e-6, and the IB force (small in magnitude) changes a lot. With a per-face θ instead, I checked that B317404C stays within its tolerance. I kept the shared per-cell θ for consistency with Keep species mass fractions in [0, 1] by construction instead of clipping #1963. Reviewers may prefer per-face θ.
  • No other golden moved.
• ./mfc.sh precheck passes.

Not verified

• GPU builds (OpenACC/OpenMP) and MPI runs of the new kernel.
• Compilers other than CCE. I set the new test's tolerance to 1e-5 instead of 1e-10: its corner faces amplify roundoff, and CI lanes for Keep species mass fractions in [0, 1] by construction instead of clipping #1963's similar test differed from the CCE golden by up to 6e-7 relative.
• The original reacting-surfaces rod case (surface reactions on). That is on a fork branch, not master.
• Flows that are legitimately above 5000 K get their face states scaled toward the cell state. Cell states are not limited.
• Rim fluid still heats to ~2500 K, above Twall. That is the corner ghost treatment, which this PR does not change.

[sbryngelson]

I checked this PR independently against the final commit 7b2c337 (CPU, CCE, --no-mpi):

• 3D repro, GREEN. The x-aligned finite cylinder with the 1800 K isothermal wall runs all 300 steps with exit 0 (run time 3587 s). Master NaNs at step 22.
  • At step 300 no cell near the body is NaN or non-physical.
  • Near-wall fluid peaks at 2852 K against the 1800 K wall. That is coarse-grid mirror-ghost heating and is not addressed here.
• One correction to the body. The 2D 200-step run and the count of "28 face states at the corners" came from an earlier version of the check that also had the 200 K lower bound. The final version was checked through the new 30-step test only:
  • master: NaN;
  • this PR: passes -o Chemistry (23/23).

Developed with Claude Code.

[sbryngelson]

Note on interaction with #1969: in the GREEN run here, near-wall gas still reaches ~2850 K next to the 1800 K isothermal wall. That over-heating comes from the existing hot-wall ghost extrapolation (image-point sample distance), which #1969 fixes; with #1969 the rim ghost states are less extreme, so this face-state guard should trigger less often. The two are complementary: #1969 corrects the ghost values, this PR guards the reconstruction against out-of-range face states that remain at sharp corners. (Developed with Claude Code.)

[github-actions]

Lines of Code

File	Lines	Diff
src/simulation/m_rhs.fpp	2033	+65

Directory	Lines	Diff
simulation	28568	+65
total	47579	+65

[sbryngelson]

Closing this for now to bring the number of open PRs from me down to four, so the ones that stay open get real review attention rather than being lost in a queue of 35.

No work is lost and nothing needs re-doing. The branch fix-cylinder-rim remains on the fork at b7b8ddbbff, this PR's diff and discussion stay here, and the head commit is additionally archived in a local verified bundle. Reopening takes no force-push and no rebase:

gh pr reopen 1975 -R MFlowCode/MFC

I intend to reopen these in small batches as the open slots free up. Apologies for the notification noise.

Opened and closed with Claude Code.

[sbryngelson]

Root-cause context for the NaN this PR limits, from a separate investigation.

This PR bounds chemistry face states with a scaling limiter because of NaNs at finite-cylinder rims. Chasing an apparently different near-IB NaN — a sphere in Mach 0.9 flow at Re_D = 100 — produced a quantitative criterion that may explain why a limiter was needed, and suggests the limiter is treating a symptom of case setup rather than a solver defect.

Measured failure map. Same case, same code, varying only Mach and cells-per-diameter (1 rank, minimal 0.44M-cell domain, reproduces identically on CPU and GPU):

Mach	20 cells/D	24	30	40	48
0.6	—	clean	clean	NaN step 9	NaN step 6
0.9	clean	NaN step 14	NaN step 7	—	—
1.2	—	NaN step 6	—	—	—

Failure gets worse with refinement, which is the signature of something other than ordinary under-resolution.

What the state does. Instrumenting the ghost publish to report extremes over all ghost points, every RK stage: the minimum ghost pressure falls to 6.9 kPa and the minimum ghost density to 0.072, against ambient 101325 Pa and 1.168 kg/m3 — roughly 1/15 and 1/16 of ambient. The adjacent fluid follows it down smoothly, and pressure and velocity go NaN while density is still positive and finite (0.94), i.e. the temperature solve fails at near-vacuum rather than any field being corrupted.

The criterion. Re_D fixes the body size, so demanding both a Mach number and Re_D = 100 forces a very small body: D = Re*mu/(rho*U) = 5.06e-06 m here, giving dx = 2.11e-07 m at 24 cells/D. The near-wall viscous stress then is

tau ~ mu * 2U/dx,   dx = D/N,   D = Re*mu/(rho*U)   =>   tau/p ~ 2*gamma*M^2*N/Re

(mu cancels). Over the map above that parameter runs 0.30 to 0.97 — the wall viscous stress is 30 to 100 percent of ambient pressure in every one of these runs. With the local pressure already down about 20% from the startup expansion, a stress that size drives internal energy through zero in a single stage.

Consistent with this, an ablation at the failing point: --inviscid runs clean to step 45 where the baseline dies at 14; disabling reactions changes nothing; disabling weno_Re_flux changes nothing. Only removing the viscous terms avoids it.

Why this may matter here. If the cylinder-rim NaN is also in a regime where 2*gamma*M^2*N/Re approaches 1, then the chemistry face states are going out of range because the near-wall state is being driven to vacuum, and a limiter bounds the consequence rather than the cause. Worth checking that parameter for the motivating case before relying on the limiter — and worth knowing that in our case the cure was choosing a Mach/Re/resolution combination that is not in that band, which also turned out to be the physically correct choice for comparison against reference data.

The parameter does not separate our cases perfectly (M0.6/40 fails at 0.40 while M0.9/20 is clean at 0.45), so it is an indicator, not a sharp threshold.

Written with Claude Code.