Second-order immersed-boundary wall gradients: two probes with a quadratic wall-normal profile

Summary

MFC's immersed-boundary wall gradients are first order, and that first order eventually contaminates every derived quantity — wall heat flux, Nusselt number, wall shear, and the species fluxes that drive surface reactions. This issue records the diagnosis, the measurements behind it, and a concrete implementation that reaches second order, with the degradation path it needs in order to be safe in dense geometries.

The problem, measured

Steady annulus between two isothermal cylinders (exact Kirchhoff solution), aligned grid, T1 = 600 K, at 10/20/40/80 cells per inner radius. Convergence orders for successive grid pairs, current code:

quantity	10 -> 20	20 -> 40	40 -> 80
T field L1	1.91	1.45	1.12
T field Linf	0.90	0.78	0.88
wall heat flux Q	1.59	1.53	0.93
local sector flux Linf	1.20	0.29	3.03

Every order decays toward 1 as the grid refines. A coarse-pair estimate looks second order; the asymptotic behaviour is first.

Why: two independent O(h) terms, and fixing one is not enough

The wall gradient is formed as a one-sided difference between the wall value and an image-point value,
g = (u_I - u_w)/d_I, with d_I = O(h). Its error has two parts:

Truncation of the linear wall-normal profile: (d_I/2) * u_nn = O(h).
Interpolation error divided by the probe distance: e_I/d_I = O(h^2)/O(h) = O(h).

So a second-order wall gradient needs both a wall-normal polynomial of degree >= 2 and interpolation at the probes with error O(h^3). Raising either alone leaves the other at first order. This is stated explicitly in the literature:

Jost & Glockner (2021): a Neumann condition with p = 2 Lagrange interpolation converges at first order; p = 3 is needed for second.
Yousefzadeh & Battiato (2019), §3: "the bilinear interpolation approximates the value of the derivative at the mirror point with first order".
Ezra, Peles & Kozak (2025), Int. J. Heat Mass Transf. 252:127314: standard ghost-cell extrapolation gives a second-order solution but only a first-order Nusselt number.

Our own measurements match this decomposition. Replacing inverse-distance weighting with a TSC (quadratic B-spline) kernel — which fixes term 2 only — improves exactly the quantities that interpolation limits and leaves the rest untouched:

quantity	IDW + clamp	TSC	`ib_ip_min_dist` (#1958)
T field L1	1.893e-04, p=1.91	9.123e-05, p=2.21	1.410e-04, p=2.13
T field Linf	1.933e-03, p=0.90	4.162e-04, p=2.09	5.608e-04, p=1.75
wall heat flux Q	5.782e-04, p=1.59	1.885e-04, p=2.47	5.895e-04, p=2.10
local sector flux Linf	1.026e-01, p=1.20	1.025e-01, p=1.19	1.031e-01, p=1.20
interface overshoot (low)	4.921e-03, p=0.16	4.890e-03, p=0.16	4.993e-03, p=0.11

(values at 20 cells/R, orders from 10 -> 20, single estimates)

Two things to read from that table. TSC clearly helps where interpolation is the limit. And the local sector flux is identical to within 0.6% across three quite different interpolation schemes, all at p ~ 1.2 — so that quantity is limited by something none of them touch, consistent with the truncation term plus how the surface cuts cells.

Proposed implementation

Two probes per wall point with a quadratic wall-normal profile, evaluated with the TSC kernel, used only for wall gradients. The ghost fill and therefore the WENO stencils are left alone, so positivity risk stays confined to diagnostics and boundary-condition inputs rather than entering the solver.

1. Probe placement

For each ghost point, take the wall foot at |levelset| along the outward normal, then place probes at fixed multiples of the local spacing, s1 = h and s2 = 2h, where h = (dx dy dz)^(1/num_dims). Fixed distances rather than |levelset|-derived ones are what keep the fit conditioned when a ghost point lies almost on the surface.

2. Interpolation at the probes

TSC (quadratic B-spline) weights on the 3x3x3 stencil centred on the cell nearest each probe. 1D weights, per direction, for fractional offset f in [-0.5, 0.5]:

w(1) = 0.5*(1.5 - |f + 1|)^2
w(2) = 0.75 - f^2
w(3) = 0.5*(1.5 - |f - 1|)^2

Keep every weight: the kernel is a partition of unity and reproduces a linear field exactly, and zeroing solid cells then renormalising destroys both properties — which is how inverse-distance weighting loses linear exactness and moves the effective sample off the image point.

3. The fit

dphi/dn|_w = [s2^2 (phi_1 - phi_w) - s1^2 (phi_2 - phi_w)] / [s1 s2 (s2 - s1)]

For a Newton surface balance the wall value is the unknown, so return the two coefficients of the form that is linear in it:

dphi/dn|_w = a - b*phi_w,    a = (s2^2 phi_1 - s1^2 phi_2)/(s1 s2 (s2 - s1)),    b = (s1 + s2)/(s1 s2)

Setting a = phi_IP/d, b = 1/d reproduces the existing one-sided difference exactly, which makes the change verifiable as a no-op before it is switched on.

4. Graceful degradation (required, not optional)

Each probe reports the fraction of its stencil weight carried by genuine fluid cells. Tiers:

tier	condition	gradient	order
1	both probes have sufficient fluid weight	quadratic, two probes	2
2	outer probe is mostly solid-bound	linear, one probe	1
3	inner probe also unusable	leave existing behaviour	-

This is load-bearing rather than defensive. A s2 = 2h probe with a 27-point stencil reaches about 3h into the fluid, so it fails in gaps under roughly 3 cells, at rims, and between close-packed bodies. With a soft-sphere contact model (collision_model = 1), force exists only while bodies overlap, so during every contact the gap has closed, there are no fluid cells between the surfaces, and the contact patch is necessarily tier 3. Record the tier per ghost point and report the tally, so a run cannot claim an order it did not achieve over that surface.

Verification

Annulus against the exact Kirchhoff solution at 10/20/40/80 cells/R: report orders for the field (L1, Linf), the integrated wall flux, the local sector flux, and the interface overshoots separately. A single metric hides the failure — the integrated flux can improve while the local flux does not move.
Confirm the tier mix is ~100% tier 1 on an isolated convex body, so that any order shortfall is attributable to the scheme rather than to degradation.
Re-run with a = phi_IP/d, b = 1/d and confirm results match the current code to round-off, separating plumbing errors from the effect of the fit.
Report the tier mix on a dense packing, where tier 2 and 3 are expected; that fraction bounds the claim.

Scope notes

The species flux rho D dY_k/dn in the surface balance uses the same first-order one-sided difference and has the same problem. It is num_species probe pairs per ghost point, so the cost should be measured on the thermal path first.
The non-convergent interface overshoot (p ~ 0.11-0.16 across all three schemes) is not addressed by this change and may need a sharper interface treatment.
ib_ip_min_dist (2nd order IBM #1958) is complementary, not an alternative: it improves stencil conditioning but keeps O(h^2) interpolation, so it cannot by itself produce a second-order wall gradient.

Reference

Ezra, Peles & Kozak (2025), Int. J. Heat Mass Transf. 252:127314, is the only work found that demonstrates a second-order Nusselt number in compressible flow with a ghost-cell IBM, using two image points plus TSC interpolation; with shocks their method stays better than first order where the standard treatment falls below it.

Filed with Claude Code.

An additional argument for the probe-based gradient, in the regime where it matters most.

s_blend_ghost_state (src/simulation/m_ibm.fpp:2592) limits the thermal extrapolation so the ghost temperature stays inside the window the thermodynamic fits are valid over (T_surface_min = 200 K, T_surface_max = 5000 K in src/common/m_constants.fpp):

theta_T = theta_T_max
if (T_s < T_IP) theta_T = min(theta_T, 0.9*(T_s - T_surface_min)/(T_IP - T_s))
if (T_s > T_IP) theta_T = min(theta_T, 0.9*(T_surface_max - T_s)/(T_s - T_IP))
T_g = T_s + theta_T*(T_s - T_IP)

That clamp is necessary — the code comment gives the reason, that a cold wall in hot gas otherwise evaluates the NASA polynomials far below their T_low. But the cap it imposes depends strongly on the temperature ratio:

case	T_s	T_IP	cap on theta_T
modest ratio (annulus, T1 = 600 K, either wall)	300 / 600	310 / 590	inert
T1 = 1500 K, cold wall	300	400	0.900
cold wall in hot gas	300	1500	0.075
hot reacting surface	1800	1700	inert

(theta_T_max is itself computed per ghost point by f_face_theta, so 1.0 above is only a nominal reference; the 0.075 is an absolute cap that binds whatever that function returns.)

Two consequences.

It is not the cause of the order decay reported in the issue. In the annulus study quoted above (T1 = 600 K) the clamp is inert at both walls, so the decay of every order toward 1 is the O(h) truncation term rather than any clipping. Worth stating explicitly, since the clamp is an obvious suspect.
The two-probe fit is not subject to it. The gradient proposed here is formed from values interpolated in the fluid at s1 = h and s2 = 2h; it never reads the clamped ghost temperature. So in a configuration where the extrapolation is capped around 0.075 — a cold body suddenly exposed to hot gas, i.e. an ignition transient — the probe-based gradient is the only one of the options that is not suppressed by roughly a factor of 13. The ghost fill remains clamped, as it must for thermodynamic validity, but the quantity that drives the surface energy balance no longer inherits that limit.

This also suggests a sharper test than the annulus for whichever implementation lands: a cold body stepped into hot gas, comparing the surface heat flux from the clamped ghost extrapolation against the probe-based gradient. The annulus has a modest temperature ratio and will not exercise this at all.

Written with Claude Code.

Three-grid update: the thesis holds, and one claim in the issue was wrong.

The orders in the issue body rest on a single 10 -> 20 pair. Adding 40 cells/R gives a second estimate, and it matters:

line	quantity	IDW + clamp	TSC	`ib_ip_min_dist`
field	L1	1.91 / 1.46	2.21 / 1.85	2.13 / 1.79
field	Linf	0.90 / 0.78	2.09 / 1.66	1.75 / 1.83
wall gradient	local sector flux Linf	1.19 / 0.29	1.19 / 0.29	1.20 / 0.29
integrated	wall heat flux Q	1.59 / 1.53	2.47 / 2.09	2.10 / 1.78
interface	overshoot (low)	-	0.16 / 2.33	0.11 / 2.37

This supports the central argument. Every scheme's order decays under refinement: TSC 2.21 -> 1.85 (L1), 2.09 -> 1.66 (Linf), 2.47 -> 2.09 (Q); ib_ip_min_dist 2.13 -> 1.79; IDW 1.91 -> 1.46. None of these is asymptotically second order. That is exactly what the two-term decomposition predicts — all three improve the interpolation error O(h^2)/d = O(h) and none removes the truncation of the linear wall-normal profile, (d/2) u_nn = O(h). Only a quadratic normal profile removes that term, which is the proposal here. A single coarse pair flatters all three schemes.

The ranking is unchanged: TSC leads at both grid pairs on every quantity interpolation governs (Q at 20 -> 40 is 2.09 against 1.78 and 1.53).

Correction to the issue body. It states that the interface overshoot sits at p ~ 0.16 and is not converging, and lists that as something this change does not address. That is wrong — over 20 -> 40 it converges at 2.33 (TSC) and 2.37 (ib_ip_min_dist). The stall was an artifact of the single coarse pair. The "Scope notes" bullet about a non-convergent interface overshoot should be disregarded.

One clarification on the local wall flux. The three schemes give the same order sequence to three figures (1.19 / 0.29, 1.19 / 0.29, 1.20 / 0.29). That quantity is computed by the analysis script from conductive face fluxes in the near-wall gas band and, in its own words, "needs no ghost values" — so it measures near-wall fluid resolution rather than the ghost-derived gradient, and no interpolation change can move it. It is also a Linf over 16 angular sectors, i.e. a max over discrete bins, so its pairwise order is not a smooth function of h and the 0.29 should not be read as an asymptotic rate. The ghost-derived integrated flux Q is the quantity that tracks the wall-gradient treatment, and that is where the schemes separate.

Written with Claude Code.

Correction to the verification plan in this issue, and a scope warning worth attaching before anyone implements it.

The verification plan tests the wrong quantity. Step 1 proposes the annulus against the exact Kirchhoff solution and reports the integrated wall flux Q. But Q there is read from D/ib_heat.dat, i.e. the face fluxes the solver actually accumulates, and those are set by the ghost fill. This proposal deliberately does not change the ghost fill — it changes how wall gradients are computed for q_w, wall shear, species fluxes and the surface balance. On a non-reacting conduction case such as the annulus, the change does not fire at all.

So the annulus is the right case for comparing interpolation schemes (which is what the tables above do) and the wrong case for verifying this change. Verification should instead be:

A reacting-surface case, comparing the wall heat flux formed from the two probes against the existing one-sided difference, at several resolutions.
A cold body stepped into hot gas, which is also the only configuration that exercises the thermal-extrapolation clamp (see the comment above): the annulus temperature ratio leaves that clamp inert.
The tier mix on an isolated convex body (expect ~100% tier 1) and on a dense packing (where tiers 2 and 3 are expected), so the achieved order is never claimed beyond the surface fraction that earned it.

The scope warning. If the near-wall fluxes themselves must improve — drag, Nu, the conducted flux — then this proposal is not sufficient, and the next step is a bounded higher-order ghost fill: minmod(quadratic, linear), Barth-Jespersen clipping to the stencil and wall values, or WENO-type extrapolation with nonlinear weights that fall back near shocks (Tan & Shu 2010). That carries positivity risk which this proposal does not: Seshadri & De (2019/2021), §4.6, report quadratic interpolation near shocks producing negative pressure or density.

And a limit that order alone may not reach. Başkaya et al. (2024) assessed ghost-cell IBMs on cold isothermal and ablating walls with gas-surface chemistry. Two independent ghost-cell implementations gave wall heat fluxes that were clearly wrong — one reports 0.167 MW/m² against 7.19 MW/m² from a cut-cell method (Table 7), with the same disagreement visible in their Fig. 18 — while a cut-cell method at the same reconstruction order matched body-fitted results. They attribute the discrepancy to conservation errors at steep wall gradients rather than to reconstruction order (§4.4, §5).

If that holds for a hot gasifying wall, then raising the wall-gradient order is necessary but not sufficient, and the remaining error is structural to the ghost-cell approach. That does not argue against this change — second-order gradients are still needed for surface kinetics, and this is the cheap and low-risk part — but it does argue for checking a reacting hot-wall case against an independent reference before treating the resulting wall fluxes as trustworthy.

Written with Claude Code.

Correction to my three-grid comment above: it claimed more than the data shows, and one of the metrics does not measure what I implied.

1. "None of these is asymptotically second order" was an extrapolation, not a result. Fitting e = a h^2 + b h (h = 1/N, N cells per inner radius) to the annulus errors:

scheme	grids	O(h) component	O(h) overtakes O(h^2) at	single power law p (L1 / Q)
IDW + clamp	4	clear	about 17-27 cells/R	1.49 / 1.37
TSC	3	none detectable	-	2.03 / 2.28
`ib_ip_min_dist`	3	small	about 260-360 cells/R	1.96 / 1.94

Only the current code is shown to drift toward first order. For TSC and ib_ip_min_dist the pairwise orders do fall (TSC L1 2.21 -> 1.85), but three grids cannot distinguish converging to 2 from drifting lower. These fits are crude (5-10% rms residual).

2. Q does not measure the local wall gradient. Q_inner_err and Q_outer_err agree to about five digits (TSC, 20 cells/R: -1.88538e-4 against -1.88544e-4). Q is fixed by steady-state conservation, so it is a global integral. The sector flux needs no ghost values, as noted before. So none of the annulus metrics measures the ghost-derived local wall gradient, and this study can neither confirm nor refute a first-order wall-gradient error. The two-O(h)-terms decomposition above remains the literature's analysis, untested here.

3. A measurement that does speak to boundary-condition enforcement. The gas-temperature error against the exact solution, extrapolated to the wall, is the effective wall temperature the fluid sees. Orders for 10 -> 20 / 20 -> 40:

	IDW + clamp	TSC	`ib_ip_min_dist`
wall error, inner	0.91 / 1.81	2.75 / 2.63	2.03 / 1.76
wall error, outer	1.79 / 1.46	2.31 / 1.92	2.15 / 1.87
wall error, inner, value at 20 cells/R	2.7e-4	4.1e-5	3.5e-4

With TSC or ib_ip_min_dist the wall temperature is enforced at about second order over 10-40 cells/R; the current code is nearer 1.5 and several times larger. Since both alternatives keep a linear wall-normal profile and still remove the first-order component, the likelier culprit is the clamp / stencil-centroid mismatch rather than truncation of the linear profile. That is a hypothesis, not isolated.

4. Part of the decay is not wall-related. The interior (more than 0.25 R from both walls) also loses order for every scheme (IDW 2.08 / 1.35, TSC 2.02 / 1.79, ib_ip_min_dist 2.03 / 1.66), so the wall-gradient argument cannot be the whole explanation. For TSC the error is about 88% a uniform signed offset at every grid; ib_ip_min_dist's offset has the opposite sign, so it is not a shared reference mismatch.

Net effect on this proposal. The case for a quadratic wall-normal profile rests on the cited literature, not on these measurements. For temperature with TSC there is currently no evidence in MFC that profile truncation limits the order. Not yet tested: velocity / no-slip enforcement (an exact-solution Taylor-Couette annulus would be the analogue), and the ghost temperature itself, which is never written out (ghost_points%T_cond is in memory only).

The verification plan in the previous comment still stands.

Written with Claude Code.

IBM: opt-in two-probe quadratic wall treatment (ib_wall_order = 2)

Navigation Menu

Repository navigation

Second-order immersed-boundary wall gradients: two probes with a quadratic wall-normal profile #2084

Summary

The problem, measured

Why: two independent O(h) terms, and fixing one is not enough

Proposed implementation

1. Probe placement

2. Interpolation at the probes

3. The fit

4. Graceful degradation (required, not optional)

Verification

Scope notes

Reference

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Second-order immersed-boundary wall gradients: two probes with a quadratic wall-normal profile #2084

Description

Summary

The problem, measured

Why: two independent O(h) terms, and fixing one is not enough

Proposed implementation

1. Probe placement

2. Interpolation at the probes

3. The fit

4. Graceful degradation (required, not optional)

Verification

Scope notes

Reference

Activity

sbryngelson commented on Oct 15, 2026

sbryngelson commented on Oct 15, 2026

sbryngelson commented on Oct 15, 2026

sbryngelson commented on Oct 15, 2026

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions