Repository navigation
IBM image-point interpolation: bilinear in 2D, trilinear in 3D, and accuracy verification #2026
Description
@danieljvickers — tracking the bilinear interpolation work you are already doing, including the corresponding trilinear treatment in 3D and verification alongside #1958.
Problem
The second-order wall extrapolation in #1958 fixes the placement of the velocity boundary condition, but s_compute_interpolation_coeffs in src/simulation/m_ibm.fpp still uses normalized inverse-distance-squared image-point weights. Those weights preserve constants but generally do not reproduce linear fields. Moving the image point farther into the fluid improves stencil availability without changing this property.
For an entirely fluid 2D stencil with corners (0,0), (h,0), (0,h), (h,h), an image point at (0.2h,0.3h), and u=y, the current arithmetic gives u_IP = 0.23422065216744875h instead of 0.3h. The error is 0.06577934783255125h. An arithmetic reproduction at five refinements (h=1/16 through 1/256) halves the error each time. This is a local interpolation check, not a full MFC convergence study and not a claim that every solver observable must converge at first order.
Scope
- Use bilinear image-point interpolation in 2D and its trilinear counterpart in 3D for fully fluid stencils. Trilinear is the 3D analogue, not a higher-order replacement for bilinear in 2D.
- Decide whether the existing inverse-distance method should remain selectable for compatibility/comparison. If so, document the interpolation choice independently of
ib_second_order_vel. - Define reconstruction and fallback behavior when the nominal stencil intersects solid cells or another body. Simply masking solid corners and renormalizing bilinear/trilinear weights generally loses linear reproduction. A larger fluid stencil or boundary-constrained reconstruction needs sufficient independent points and a conditioning check.
- Preserve the existing protection for coincident image points, empty stencils, and overlapping bodies; verify GPU/MPI compatibility.
- Coordinate with s_ibm_correct_state writes q_prim inside the loop that reads q_prim at image-point stencils: read-after-write race across ghost points #1901: interpolation should read a consistent input state when ghost states are updated in parallel.
Verification
- Check constant and affine-field reproduction to roundoff on fully fluid 2D and 3D stencils, including nonuniform grids and points near corners/faces.
- Measure quadratic-field interpolation error under refinement; expect second-order error for a regular, fully fluid multilinear stencil.
- Exercise oblique/curved walls, different subcell offsets, partially masked stencils, and the documented fallback. Do not infer accuracy of the fallback from the fully fluid test.
- Combine interpolation with 2nd order IBM #1958's wall extrapolation in a smooth wall-bounded convergence case, measuring near-wall velocity and surface quantities as well as bulk error.
- Include a prescribed moving-wall case that crosses cells, with separate spatial and temporal refinement and multiple grid offsets. This verifies reconstruction under motion; geometric conservation and the body RK update remain separate work (Moving immersed boundary: add geometric-conservation-law (GCL) treatment to eliminate spurious pressure / CFL blow-up #1636, Moving IBM: RK state weights and centroid/angle stage velocities are incorrect #2025).
Mittal et al. (2008) use bilinear/trilinear image-point interpolation as part of their sharp-interface construction: author PDF. Achieving the intended overall order also depends on boundary reconstruction, stencil quality, geometry, and time integration.
Reproducer
This standard-library Python snippet reproduces the current interpolation formula on the all-fluid stencil:
from math import log2
previous_error = None
for n in (16, 32, 64, 128, 256):
h = 1 / n
x, y = 0.2 * h, 0.3 * h
points = [(0, 0), (h, 0), (0, h), (h, h)]
raw = [1 / ((a - x)**2 + (b - y)**2) for a, b in points]
weights = [w / sum(raw) for w in raw]
value = sum(w * b for w, (a, b) in zip(weights, points))
error = abs(value - y)
order = None if previous_error is None else log2(previous_error / error)
print(n, error, order)
previous_error = errorSource reviewed: #1958 head 82bd69ffabe040a8efde9e9985c527b0605a445c, base a6add01e35239cf9e87de38972e968ad1799c04c. This issue and the arithmetic check were prepared with Codex assistance; no new full-solver convergence run is claimed here.
Activity
sbryngelson commentedon Oct 13, 2026
Two things for whoever picks this up: a documentation trap, and a measurement that changes the rationale.
1. The docstring already claims this is implemented.
m_ibm.fpp:1028on master readsbut the coefficients it uses are inverse-distance-weighted, with non-fluid contributions dropped:
So reading the source suggests this issue is already done. Worth fixing the comment independently of the change itself.
2. The cap, not IDW, is what costs the order -- and that changes why this issue matters. I had assumed from Luo et al. 2017 Part II (IDW is first order because it drops non-fluid contributions) that IDW was the cause here. A replica of the ghost-value algebra says otherwise. Steady conduction in an annulus, inner circle an isothermal ghost-cell IB, the ghost treatment replaying
m_ibm.fppexactly, with interpolation and the mirror-weight capd = max(ip_dist, |levelset|)as independent switches:With the cap the wall gradient does not converge at all over a factor of 8 in resolution, and its Linf grows (8.16e-02 -> 1.40e-01). The reason is that the cap's error is scale-invariant:
ghost_geometryon the same geometry has the cap binding at 17.9 % of face ghost points with|levelset| - ip_distup to 0.050 dx againstip_distof 0.503..0.905 dx, so the relative error(|ls| - ip_dist)/ip_distis ~10 % at every resolution. The cap pinsthetato 1 as though the image point sat at|levelset|, while the IDW centroid actually sits nearer the wall.IDW-without-cap and bilinear come out equivalent in order. So the value of this issue is not that IDW is too low order -- it is that a bilinear sample makes the cap unnecessary: in the replica, bilinear's cap-on and cap-off results are bit-identical, i.e. the cap never binds. That matters because deleting the cap on its own is not a fix; tried directly, it moved 32 of 43 IBM tests by up to 7.5e+04 relative, since
theta > 1amplifies the ghost value. The cap is load-bearing for stability, and this change removes the need for it.Caveats. This is a Laplace replica (steady, pure conduction, 2D circle, no WENO or flow); it isolates the ghost-value algebra and nothing else, and it is not a measurement in MFC. The replica's bilinear moves body-side nodes to their wall projections carrying
T_wall; a different treatment would change its constants, though not the cap conclusion, which does not depend on the bilinear column. Also worth setting expectations: the wall gradient stays first order even with bilinear (Linf orders 0.81, 0.96, 0.75), so this change should be expected to buy a converging gradient and a second-order field, not a second-order gradient.The concrete check once it is implemented: run a curved isothermal wall at 40 -> 80 cells per radius and confirm the wall-gradient error actually falls, since that is the quantity that is currently frozen.
Prepared with Claude Code.