IBM: extrapolate hot/reacting-wall ghost state from where the image point is actually sampled #2062

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.

This PR changes CFD results at hot or reacting IB walls; the verification is below.

Summary

Fixes the ghost state at isothermal (thermal_bc = 1) and reacting (surface_reaction = 1) IB walls in chemistry runs (#1821). Before this fix, near-surface ghosts were too far from Twall, and the surface-reaction solve used the wrong gradient length.

Mechanism

s_blend_ghost_state builds the ghost value as the linear mirror phi_g = phi_s + theta*(phi_s - phi_IP) with theta = 1. That is only correct if phi_IP is sampled |levelset| in front of the wall. s_compute_interpolation_coeffs uses inverse-distance weights that drop solid stencil cells. For a ghost point near the surface, the mirror point lies inside the ghost's own (solid) cell, so phi_IP actually comes from fluid cells about 1 dx farther out. The ghost then gets roughly 2 T_s - T(~1 dx) where it should get about T_s. At a hot wall in cold gas that is several hundred K too hot. s_ibm_correct_state also passes |levelset| as the gas-side gradient length d to the surface solve. At these points that overestimates the gradient, so the surface fluxes are too large as well.

Fix

• ghost_point%ip_dist stores the wall distance of the weighted stencil centroid sum c_i x_i, which is where the image-point value is actually sampled. It is computed in s_compute_interpolation_coeffs by scalar accumulation. An array-constructor version gave garbage inside the offload loop on MI300A.
• s_ibm_correct_state sets d = max(ip_dist, |levelset|).
  • d is the surface-solve gradient length.
  • |levelset|/d is passed to s_blend_ghost_state as theta_max, the linear-extrapolation weight. It is at most 1 and equals 1 for an ideal mirror. The existing temperature-window and species-positivity limits still apply on top.
• Adiabatic / inert IB walls are unaffected: with T_s = T_IP and Ys_s = Ys_IP, theta has no effect.

How it was found

A 3D Mach-2 shock hitting a curtain of 413 moving reacting spheres (12 cells/D, 1800 K walls) went NaN at about step 200. Isolation runs ruled out motion, collisions, the wall-temperature model, surface chemistry, the transport time step, and CFL; halving CFL gave NaN at the same physical time and cell. The failure needed a hot isothermal-type wall in 3D. The NaN cells were in narrow gaps (3–5 cells) between neighbouring spheres. In those gaps, first fluid cells sitting on the surface in staircase notches had over-hot ghosts on 2–3 faces and were heated to 3000–4080 K before blowing up. Ghosts with levelset ≈ 0 held about 2850 K against an 1800 K wall, and the mean ghost T did not fall toward Twall as levelset → 0. With the fix, the same case runs to t_stop, and a 5x longer variant also completes. At the old failure time, the hottest near-wall fluid cell is about 2190 K.

Verification

The verification runs (1 node x 4 MI300A, Tuolumne) were done on a development branch that has this fix on top of master plus unmerged IB thermal features (thermal layer, lumped wall temperature). This PR ports only the ghost-distance change. The resolved Q1/Q2 rod cases use only thermal_bc = 1 and surface_reaction, which master has. Case A uses --thermal-layer 0.1, and the 2D curtain uses thermal_bc = 3; both are from the unmerged features. OLD = before the fix, NEW = with the fix.

Implied wall temperature. For each fluid/ghost face pair, T is interpolated linearly in levelset to the zero crossing. This is the wall temperature the discrete profile actually imposes. Results are front-half means.

Resolved rod (Mach 0.2, Re_D = 471, isothermal 1800 K wall, T∞ = 496 K; Q1 non-reacting, Q2 with surface reactions):

The OLD error decays as about h^0.8. NEW's is about 15x smaller and decays as about h^2.6.

Case A (Mach 1.5, Re_D = 5.4e4, reacting, 1800 K; boundary layer unresolved at 2–4 cells):

OLD imposes a wall 85–210 K too hot at every resolution tested. NEW converges onto Twall.

Heat flux, resolved non-reacting rod (Q1). Q_cons is the heat the IB injects, from a global energy balance; two box sizes agree to 0.5%. Q_fit is the gas-side gradient flux.

• NEW's heat input varies by 0.2% over 40–80 cells/D.
• NEW's gas-side flux matches the energy it injects (Q_fit/Q_cons = 1.00–1.01). OLD's does not (0.92–0.95).
• The OLD−NEW gap shrinks with dx, as expected for an error in the ghost value alone.

Carbon flux, resolved reacting rod (Q2). Regression rate is the mean over t* = 2.5–4.

NEW is grid-converged to 0.3%. OLD drifts by +4.9% and moves away from NEW as the grid is refined. This is consistent with the bug: the gradient-length error is an O(1) factor at the affected ghost points, not O(dx). An offline replica of the weights on the case-A rod shows that 23–27% of ghost points have d/|levelset| > 2 at every resolution from 40 to 200 cells/D.

Case A burning rate and heat flux.

The burning rate drops by 19–24%. The wall heat flux through the staircase changes much less (0–10%).

2D shock-curtain smoke run (38 particles, 12 cells/D, 1800 K): total carbon loss goes from 4.952e-3 to 4.344e-3 kg/s/m (−12.3%). Mean near-wall gas is 80–220 K cooler.

What this does not address

• Case A is still under-resolved at 200 cells/D under both builds. Its burning rate rises by about 4% from 150 to 200 cells/D, and extrapolating to zero cell size is unreliable. Its absolute rates are therefore limited by the boundary-layer resolution, independent of this fix.
• The fix does not remove the near-wall gas maxima of 2700–3100 K in the 2D curtain smoke run. Ghost maxima of about 3100 K remain too. These are deep ghosts (theta ≈ 1) next to cold gas, where the full mirror is correct. The overshoot comes from the hot-wall / cold-gas start plus shock compression at 12 cells/D, not from the misplaced image point.
• Under shock compression, fluid cells right at a 3D staircase wall can still sit above Twall (about 2200–2600 K in the curtain case). This comes from the first-order staircase boundary.

Tests

CPU, ./mfc.sh test --no-mpi --no-gpu -j 24:

• --only IBM: 61/61 pass, unchanged.
• --only Chemistry: 23/23 pass after the golden changes below.

Goldens:

• F52F0D4C (2D -> Example -> ibm_reacting_surface, regenerated). This is a 1200 K reacting carbon cylinder; the suite runs it on a 26x26 grid (about 4 cells across the rod) for 1.27 ms. It goes through exactly the changed path, so its solution moves throughout the wake. The largest absolute change is 9.29e4 in energy (rel 1.78, a ghost cell). The largest relative change is 2.1e2, at a near-zero value. This is the only existing golden in the IBM and Chemistry suites that moved. ip_dist is computed for every IB ghost point but is only read in the chemistry branch of s_ibm_correct_state, so other goldens should not move. The full suite was not run locally.
• 204838D0 (2D -> Chemistry -> Isothermal Wall -> IBM Hot Cylinder, new). This is a 1500 K non-reacting IB cylinder in 300 K h2o2 gas, 50x50, 20 steps, tol 1e-6, with no new mechanism. It pins the thermal ghost extrapolation, which had no registered case. Run against the pre-fix source, it fails with energy differing by up to 5.1e4 (12% of the field range) and density by up to 2.2e-2 (7%) in near-wall and ghost cells. A golden can detect a change in this path but cannot by itself assert that the wall temperature is correct. That check is the Timp convergence study above.

./mfc.sh format and ./mfc.sh precheck pass.

A fix for this is already written. The pull request was closed only to reduce the number of simultaneously open PRs from this author, not because the work was rejected or found incorrect.

• Closed PR: IBM: extrapolate hot/reacting-wall ghost state from where the image point is actually sampled #1969
• Branch: fix-ib-ghost-distance on the sbryngelson fork, at 76cde939da
• Reopening needs no force-push and no rebase: gh pr reopen 1969 -R MFlowCode/MFC

Tracked with the pr-to-do label so the whole set is discoverable together.

Filed with Claude Code.