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

[sbryngelson]

Contribution Policy

We do not accept pull requests generated primarily by AI without genuine understanding or real-world usage context.

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.

---

Acknowledgement

• I confirm this PR meets the above expectations and reflects my own understanding and real-world context.

---

PR template credit: junegunn

---

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):

run	cells/D	Timp − Tw OLD, mean (rms) [K]	Timp − Tw NEW, mean (rms) [K]
Q1	40	+38.8 (49.1)	+2.5 (3.3)
Q1	60	+26.4 (35.9)	+0.9 (1.4)
Q1	80	+22.6 (29.7)	+0.4 (0.8)
Q2 (reacting)	40 / 60 / 80	+37.0 / +25.1 / +21.5	+2.6 / +1.0 / +0.4

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):

cells/D	Timp − Tw OLD, mean (rms) [K]	Timp − Tw NEW, mean (rms) [K]
40	+180 (365)	−115 (234)
60	+212 (355)	−53 (185)
100	+145 (211)	−11 (94)
150	+109 (149)	+4 (36)
200	+85 (113)	+4 (17)

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.

cells/D	40	60	80
Q_cons OLD [W/m]	2363	2346	2340
Q_cons NEW [W/m]	2319	2316	2315
NEW/OLD − 1	−1.9%	−1.3%	−1.1%
Q_fit OLD / NEW [W/m]	2237 / 2324	2183 / 2339	2149 / 2316

• 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.

cells/D	40	60	80
OLD [µm/s]	12.88	13.10	13.51
NEW [µm/s]	11.26	11.27	11.24
NEW/OLD − 1	−12.6%	−14.0%	−16.8%

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.

cells/D	regression OLD / NEW [µm/s]	NEW/OLD − 1	Q_face OLD / NEW [kW/m]
40	92.3 / 70.3	−23.9%	17.0 / 15.3
60	109.4 / 83.5	−23.7%	20.2 / 18.7
100	121.7 / 95.1	−21.9%	19.0 / 18.6
150	130.4 / 105.5	−19.1%	19.2 / 18.9
200	136.0 / 109.9	−19.2%	18.6 / 18.5

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.

[github-actions]

Lines of Code

File	Lines	Diff
src/simulation/m_ibm.fpp	1861	+20
src/common/m_derived_types.fpp	483	+1

Directory	Lines	Diff
common	10449	+1
simulation	28465	+20
total	47447	+21

[sbryngelson]

Validation from a particle-resolved IB campaign on Tuolumne (MI300A); build = master a6add01 + #1953, #1955, #1956, #1957, #1962 and this PR. These runs check the fixed thermal_bc = 1 wall against heat-transfer references; they were not repeated without this PR.

Nusselt number, isothermal wall (Tw/T∞ = 1.1, 330 K in 300 K moist air, M = 0.2, gas reactions off)

case	cells/D	Nu (MFC)	reference	error
2D cylinder, Re 20	20 / 40	2.456 / 2.428	2.409 (Lange et al. 1998)	+1.9 / +0.8 %
2D cylinder, Re 40	20 / 40 / 80	3.303 / 3.243 / 3.204	3.281	+0.7 / −1.1 / −2.3 %
2D cylinder, Re 100	20 / 40	5.298 / 5.151	5.128	+3.3 / +0.5 %
sphere, Re 100	20	7.15	7.29 (Ranz–Marshall), 6.55 (Whitaker)	−2.0 / +9.1 %

Churchill–Bernstein gives the same picture (within 3 % at 40 cells/D).

• Nu = hD/k_f, h = (q̄_w − q̄_w,0)/(Tw − T∞). q_w = −k(Tw) ∂T/∂n from a quadratic least-squares fit to fluid cells within 3.5 dx of the wall; q̄_w,0 is the same run with Tw = T∞, which removes viscous heating (adiabatic-wall T ≈ 301.5 K, 5–6 % of the 30 K difference at M = 0.2). k_f = k(315 K) from Cantera with the run's mechanism; correlations at Pr_f = 0.69.
• Domain [−15D, 25D] × [−16D, 16D] (3.1 % blockage, uncorrected); sphere [−5D, 15D] × [−6D, 6D]².

Hot wall, reacting (2D cylinder, Tw = 1400–3000 K in 1000 K air, 16 cells/D, thermal_layer 0.25 D). The gas within 3 dx of the wall never exceeds Tw (maximum 0.13 % below, e.g. 2596.9 K at 2600 K). The carbon flux from the surface solve, which uses the corrected gradient length, is within 2 % of a kinetic/film-theory series estimate across the kinetic → diffusion transition (details on #1956).

Caveats:

• ΔT = 30 K is mild, so the Nu runs mainly show the fix does not bias an ordinary isothermal wall; the large pre-fix errors in the PR body were at 1800 K.
• Nu at Re 40 drifts down with refinement (observed order 0.59, Richardson −4.7 % vs Lange). A likely contributor is the first-order no-slip wall (IB no-slip ghost velocity is set to the wall value, not reflected (first-order wall) #2014), whose coarse-grid slip raises heat transfer; not confirmed.
• A linear fit within 1.5 dx changes the uncorrected Nu by 1.4–5.4 % (wall-gradient fit uncertainty).

Made with Claude Code.