Skip to content

IBM: curvature and body-acceleration terms in the ghost pressure of slip and moving walls - #2209

Open
sbryngelson wants to merge 7 commits into
MFlowCode:masterfrom
sbryngelson:fix/ib-slip-wall-ghost-state
Open

sbryngelson wants to merge 7 commits into
MFlowCode:masterfrom
sbryngelson:fix/ib-slip-wall-ghost-state

Conversation

@sbryngelson

@sbryngelson sbryngelson commented Nov 12, 2026 •

Copy link
Copy Markdown
Member

Four commits on origin/master 9649d50 (after #2199): the curvature term, the acceleration term, a test and docs commit, and a golden commit (twelve regenerated goldens, one new).

Summary

The ghost pressure of an immersed boundary comes from the wall-normal momentum balance, dp/dn = -rho (a_w . n) + rho sum_i kappa_i u_i^2, n out of the body (#2169 added the centripetal term of a spinning wall). Two terms were missing:

  1. Curvature at slip walls. Fluid sliding along a curved slip wall needs the pressure to rise away from a convex body, dp/dn = rho u_t^2 / R for a circle or sphere. A slip ghost point mirrored the image-point pressure, so a steady vortex around a slip cylinder, an exact Euler solution, was not held in balance. Now s_compute_ghost_point_pressure takes the image-point velocity, forms its tangential part relative to the wall and adds rho kappa u_t^2. kappa is known for circles and spheres (geometry 2, 8: 1/R) and for the curved side of a cylinder (geometry 10: 1/R, none along its axis); every other shape has kappa = 0 and keeps today's condition. No-slip walls are unchanged (their fluid is held at the wall velocity).
  2. Acceleration of a prescribed body. patch_ib%accel and %angular_accel hold (vel - step_vel)/dt and the angular counterpart at the first Runge-Kutta stage. The ghost pressure adds the translational acceleration (two-way bodies keep force/mass) and alpha x r_w to the existing omega x (omega x r_w). The carve-out that kept the image-point pressure for prescribed bodies of any kind (the centroid acceleration of a hinged body was missing) is removed, so every moving body takes the balance. A body of constant velocity has zero acceleration and an unchanged ghost state (bitwise).

The wall velocity expression is now one function, f_wall_velocity, shared by the ghost velocity (no-slip and slip) and the pressure. The chemistry branch that keeps the image-point composition took its temperature from the image-point pressure and the ghost density; it now uses the ghost pressure, so T_ghost = T_IP (identical wherever the pressure is not modified, i.e. every fixed no-slip wall).

Verification

NVIDIA A100, nvfortran 25.11, --gpu acc, 1 rank; "master" is 9649d50, "this PR" is the branch. Inviscid gamma = 1.4, WENO5, HLLC, SSP-RK3, adaptive CFL 0.5, extrapolation boundaries.

Exact steady vortex around a fixed slip circle

u_theta = M c0 R / r with isentropic Bernoulli pressure and density from the stagnation state (p0 = 1e5 Pa, rho0 = 1.2 kg/m3), tangent to the wall, so an exact steady Euler solution; the domain is +-3 R, the run is 0.4 (half width)/c0 from the exact field. Error is the rms over the cells in the first wall-distance band of (p - p_exact) / (0.5 rho0 (M c0)^2), (u_theta - u_theta,exact) / (M c0) and u_r / (M c0) (zero in the exact solution). 10 / 20 / 40 cells/R; observed orders between pairs. "TSC" is ib_ip_interp = 2.

M Scheme Pressure Order Tangential velocity Order Radial velocity Order
0.3 master 7.79e-02 / 3.39e-02 / 1.83e-02 1.2, 0.9 3.03e-02 / 1.76e-02 / 9.07e-03 0.8, 1.0 2.98e-02 / 1.46e-02 / 8.07e-03 1.0, 0.9
0.3 this PR 3.28e-02 / 1.58e-02 / 6.38e-03 1.1, 1.3 3.19e-02 / 2.00e-02 / 1.07e-02 0.7, 0.9 1.31e-02 / 6.73e-03 / 2.61e-03 1.0, 1.4
0.3 master, TSC 7.42e-02 / 3.62e-02 / 1.89e-02 1.0, 0.9 3.00e-02 / 1.70e-02 / 8.17e-03 0.8, 1.1 2.90e-02 / 1.56e-02 / 8.46e-03 0.9, 0.9
0.3 this PR, TSC 1.23e-02 / 4.53e-03 / 1.97e-03 1.4, 1.2 3.26e-02 / 2.25e-02 / 1.39e-02 0.5, 0.7 2.92e-03 / 1.29e-03 / 5.37e-04 1.2, 1.3
0.6 master 7.23e-02 / 2.98e-02 / 1.53e-02 1.3, 1.0 3.77e-02 / 1.89e-02 / 8.35e-03 1.0, 1.2 4.81e-02 / 2.40e-02 / 1.33e-02 1.0, 0.9
0.6 this PR 2.58e-02 / 1.43e-02 / 5.11e-03 0.8, 1.5 3.37e-02 / 1.97e-02 / 9.95e-03 0.8, 1.0 1.92e-02 / 9.45e-03 / 3.74e-03 1.0, 1.3
0.6 master, TSC 6.66e-02 / 3.17e-02 / 1.54e-02 1.1, 1.0 3.67e-02 / 1.94e-02 / 8.07e-03 0.9, 1.3 4.78e-02 / 2.61e-02 / 1.40e-02 0.9, 0.9
0.6 this PR, TSC 6.49e-03 / 2.74e-03 / 1.19e-03 1.2, 1.2 3.59e-02 / 2.50e-02 / 1.54e-02 0.5, 0.7 2.21e-03 / 1.11e-03 / 3.62e-04 1.0, 1.6

At 40 cells/R the first-cell pressure error falls 2.9x (M 0.3) and 3.0x (M 0.6) with the default image point and 9.6x and 12.9x with TSC; the spurious flow through the wall falls 3.1x and 3.6x, and 15.8x and 38.7x with TSC. The pressure order stays near 1 to 1.5. The tangential velocity error does not improve; at 40 cells/R it is 1.2x (default) and 1.7x (TSC) larger at M 0.3, and 1.2x and 1.9x larger at M 0.6. A slip ghost copies the image point's tangential velocity, which is first order where it varies with distance (here du_t/dn = -u_t/R), and the old pressure error had partly masked that. It is not addressed here (see the follow-up below).

The same vortex around a z-axis cylinder (geometry = 10, turned 90 degrees about y, 25 periodic cells in z) reproduces the 2D circle: first-cell pressure error 3.28e-2 (circle) and 3.28e-2 (cylinder) at 10 cells/R and 1.58e-2 and 1.62e-2 at 20, against 7.68e-2 and 3.40e-2 on master. The sphere branch (same kappa as the circle) has no exact-solution test.

Accelerating prescribed cylinder (added mass)

A D = 1 circle with moving_ibm = 1, vel(1) = a t (a = 1e-3 c^2/D), in inviscid gas at rest (rho = 1, c = 1), slip wall, domain +-6 D. Incompressible potential flow gives the surface pressure p = p0 + rho a R^2 cos(theta)/r (net force -rho pi R^2 a). Error is the rms of (p - p_exact) / (rho a R) in wall-distance bands at t = 4 D/c, with the centroid taken from the saved IB state.

Band (cells) master, 16 / 32 / 64 cells/D this PR, 16 / 32 / 64 cells/D
0-1 1.37e-01 / 4.23e-02 / 5.43e-02 1.18e-01 / 3.82e-02 / 5.14e-02
1-2 8.93e-02 / 3.30e-02 / 3.93e-02 7.70e-02 / 2.98e-02 / 3.78e-02
2-4 5.85e-02 / 2.83e-02 / 3.23e-02 4.82e-02 / 2.51e-02 / 3.20e-02

The acceleration term lowers the near-wall pressure error by 14%, 10% and 5% in the first band (the curvature term alone changes it by less than 2%: u is small here). The gain is small and shrinks with resolution, and the error does not converge: it stalls at about 3 to 5% of the signal between 32 and 64 cells/D, a floor that comes from something else in this test (acoustic tail of the impulsive start in a +-6 D box, the compressibility correction, and the discrete body position), which I did not separate. A viscous no-slip variant (Re 1000) is dominated by the boundary layer (error 0.24 to 0.31 in both builds, no change with the term) and says nothing about it.

Not conclusive

  • Flow past a fixed slip circle started from the incompressible potential field (M = 0.1, +-5 D box, run to t = 1 D/U, errors against Cp = 1 - 4 sin^2 theta on the front half): the first-band Cp error is 1.5e-1 / 9.5e-2 / 5.9e-2 on master and 1.3e-1 / 8.0e-2 / 5.3e-2 with this PR at 10 / 20 / 40 cells/R, and the second-band error stalls at 3.6e-2 in both. The error is dominated by something other than the wall pressure, so this test does not discriminate.
  • Mach 1.5 shock over a free slip cylinder (examples/2D_mibm_shock_cylinder with the grid set to 250 / 500 / 1000 cells across the 4 x 4 domain, t = 1.5): the cylinder's x centroid at the end is 0.5564 / 0.5570 / 0.5613 on master and 0.5592 / 0.5581 / 0.5618 with this PR. The difference between builds is below the grid error. All fields finite. No external reference values were used.

Tests and goldens

  • New test IBM -> Slip Circle -> Cross Flow (a slip circle in 100 m/s gas at 1e5 Pa). The acceleration term is exercised by the existing IBM -> Prescribed Kinematics -> Pitch Ramp tests (analytic prescribed velocities cannot be given to the test harness, which builds one generic binary).
  • Twelve goldens change, all intended: ten slip circle and sphere tests (2D and 3D, inviscid and viscous, 1 and 2 fluids, model_eqns 2 and 3), with a change of order rho u^2 (at most 1.5e-6 absolute at u = 1e-3), and the two Pitch Ramp tests. In the Pitch Ramp the plate starts to turn impulsively in a gas at p of order 0.1, where the wall-point acceleration is large against the pressure, so the state near the plate changes at order one: 168 of 15625 cells change by more than 0.1 in total energy (the largest by 0.87, density 0.50 to 0.84 at that cell), all in a band within about three cells of the plate. Cylinder, cuboid, airfoil, STL, no-slip and all other goldens are unchanged. Generated on NVIDIA and checked on AMD. The acceleration is also not evaluated at step 0 (there step_vel is the unset initial velocity, so a prescribed velocity that is nonzero at t = 0 would read as v/dt); the two 3D Pitch Ramp goldens include that guard, which moves the total energy by at most 1.2e-3 absolute (5e-4 of its largest value) against the version without it (the order-one numbers above are against master).
  • IBM label plus the eleven IB examples in the test list, without chemistry (--no-chemistry):
    • NVIDIA A100, nvfortran 25.11, OpenACC, 92 tests: 79 passed / 13 failed before regeneration (the twelve above and the new test without a golden); 92 passed, 0 failed after, on this branch.
    • AMD MI300A, amdflang (AMD flang 24.0.0git, AFAR), OpenMP offload, against the goldens generated on NVIDIA: 93 passed, 0 failed. That run is on the branch of the follow-up ib_slip_ccst PR (these 92 tests plus its one), so it exercises this PR's code with the option off; the branch of this PR alone was not built on AMD.
    • The four IB chemistry examples (ibm_reacting_surface, ibm_burning_grain, ibm_flameholder, ibm_lumped_wall) pass unchanged on NVIDIA.
    • The vortex errors above agree between NVIDIA and AMD to 1e-12.
  • Not run: GNU, Intel, Cray CCE (Frontier CCE runs in CI), multi-rank, speed.

Limitations

  • kappa is exact for circles, spheres and cylinder sides only; ellipses, airfoils, STL bodies and cuboids keep the zero-gradient condition at a slip wall (their goldens do not move).
  • The tangential velocity of a slip ghost is still a copy of the image point's, and the error of it grows with this PR on the vortex (above).
  • The acceleration of a prescribed body is a forward difference over one step (first order in dt, centred half a step ahead of the state it acts on). Its measured benefit is small. For an impulsive start at low pressure it can saturate the existing clamp; the two shipped prescribed-kinematics examples did not. Flapping plate at its shipped size (115 x 49 x 82, 1749 steps, A100): no NaN, forces bit-identical until the body moves, then Fz changes by 6.2% rms (mean +3%), Fx 1.7%, Ty 5.9%; no ghost pressure update (0 of 11.5e6) reached the clamp. Pitch-up on the shipped grid (MI300A, 150 steps from the ramp onset): no NaN, Fz changes by 0.94% rms; on a 4x coarser grid over the whole ramp (A100) Fz changes by 4.4% rms and no update (0 of 14.1e6) reached the clamp. The lift change has the sign of an added-mass reaction (more lift while the pitch accelerates, less while it decelerates); no analytic reference applies, so its size is not validated.
  • Chemistry with a slip or accelerating wall (h2o2 mechanism, 1200 K, 160 x 120 cells, 300 steps; a slip circle, an accelerating circle with and without slip, and an injecting grain with both): no NaN, mass fractions in [0, 1], and the ghost temperature equals the mirror temperature (ratio 1.0000, range 0.9986 to 1.0012) in master and this PR; the injecting grain gives the same ratio (0.874) in both. The wall-layer pressure changes by up to 11.5% in the accelerating cases. The A100 and MI300A agree. Fixed no-slip chemistry walls are unchanged.
  • The linearization behind the correction holds while the correction is order one; the existing [0.5, 2] clamp bounds it.

This PR was made with Claude Code (AI).

The wall-normal momentum balance dp/dn = -rho (a_w . n) + rho sum_i kappa_i u_i^2 gets its centripetal term from the image-point velocity relative to the wall, so a slip wall is no longer mirrored with the zero-gradient pressure. kappa is known for a circle, a sphere and the curved side of a cylinder, and zero for any other shape. The wall-velocity expression is shared by the ghost velocity and the pressure. The chemistry branch that keeps the image-point composition now takes the temperature from the ghost pressure it fills the density with.
…rescribed-motion body

patch_ib%accel and %angular_accel hold (vel - step_vel)/dt and (angular_vel - step_angular_vel)/dt at the first Runge-Kutta stage. The ghost pressure adds alpha x r_w and, for moving_ibm = 1, the translational acceleration (two-way bodies keep force/mass), so every moving body takes the wall-normal momentum balance, including prescribed-kinematics bodies, whose hinge rotation was carved out because the centroid acceleration was missing. A body of constant velocity has zero acceleration and an unchanged ghost state.
Adds a slip circle in a 100 m/s cross flow, whose ghost pressure carries the curvature term, and documents the wall-normal momentum balance of the ghost pressure at slip and moving walls.
Twelve goldens move: ten slip circle and sphere tests (2D and 3D, inviscid and viscous, one and two fluids, model_eqns = 2 and 3), whose ghost pressure now carries the curvature term (changes of order rho u^2, up to 1.5e-6 absolute at u = 1e-3), and the two prescribed-kinematics pitch ramps, whose ghost pressure now carries the acceleration of the wall point. One golden is new, for the slip circle in a cross flow. Generated on NVIDIA (nvfortran 25.11, OpenACC); cylinder, cuboid, airfoil, STL and no-slip goldens are unchanged.
@github-actions

Copy link
Copy Markdown

Lines of Code

File Lines Diff
src/simulation/m_ibm.fpp 2533 +22
src/simulation/m_time_steppers.fpp 942 +4
src/simulation/m_start_up.fpp 1463 +2
src/common/m_derived_types.fpp 538 +1
Directory Lines Diff
common 10801 +1
simulation 31201 +28
total 50534 +29

@github-actions

github-actions Bot commented Nov 12, 2026 •

Copy link
Copy Markdown

Claude Code Review

Head SHA: 22d8112

Files changed:

  • 6
  • docs/documentation/case.md
  • src/common/m_derived_types.fpp
  • src/simulation/m_ibm.fpp
  • src/simulation/m_start_up.fpp
  • src/simulation/m_time_steppers.fpp
  • toolchain/mfc/test/cases.py

Findings:

  • medium src/simulation/m_ibm.fpp: In the surface-reaction branch, T_IP = pres_IP*mw_IP/(rho*gas_constant) became pres_GP*..., but the neighbouring comment still says the heterogeneous solve retains the image-point state, and rho is still the image-point density. This line was not part of the stated change (curvature and acceleration ghost pressure), and it changes the temperature and energy of every reacting ghost point wherever pres_GP /= pres_IP. That now includes slip and moving walls. Either revert it or update the comment, and show with a reacting-IB test that the thermal and energy results stay consistent (confidence moderate; I could not see the surrounding code, since m_ibm.fpp is not in the context directory).
  • medium src/simulation/m_ibm.fpp: The old guard excluded moving_ibm == 1 with kin_model > 0. Its comment said hinged kinematics "rotate about a hinge, so the centroid's own acceleration would be needed". Those bodies now always take s_compute_ghost_point_pressure, and their accel and angular_accel come from differencing vel and angular_vel across RK stage 1. Prescribed-kinematics velocities are set by s_prescribed_kinematics at t_stage, which is also where the differenced values come from. That is plausible, but it is a behaviour change for flapping cases with no new test, and any golden files that use kin_model > 0 will shift. Please confirm that vel is the centroid velocity there, and add or verify a test (confidence moderate).
  • medium src/simulation/m_time_steppers.fpp: accel and angular_accel are computed only at s == 1 with (vel - step_vel)/dt. For RK3 stages 2 and 3 the ghost pressure reuses this first-stage value, and the body's velocity differs there. The docs describe this ("change over the first Runge-Kutta stage"), so it is a first-order approximation rather than a bug. Under cfl_dt, check that dt is the same step used for the stage update. For moving_ibm >= 2, accel is overridden by force/mass in the ghost-pressure routine, but angular_accel is still the stage difference rather than torque/moment, so the two are inconsistent (confidence low-moderate).
  • medium src/simulation/m_ibm.fpp: A fixed slip wall (moving_ibm == 0, slip) now goes through the ghost-pressure path. It no longer keeps pres_GP = pres_IP, so every existing slip-IB golden file will change. The PR adds only one new test (cross-flow circle). Please confirm that the expected golden-file updates are regenerated and verified, since the project requires verification for PRs that change CFD results.

@sbryngelson
sbryngelson marked this pull request as ready for review November 12, 2026 14:23
@sbryngelson
sbryngelson requested a review from a team as a code owner November 12, 2026 14:23
@sbryngelson
sbryngelson requested review from wilfonba and removed request for a team November 12, 2026 14:23
@codecov

codecov Bot commented Nov 12, 2026 •

Copy link
Copy Markdown

Codecov Report

❌ Patch coverage is 54.05405% with 17 lines in your changes missing coverage. Please review.
✅ Project coverage is 65.16%. Comparing base (9649d50) to head (47bdbbc).
⚠️ Report is 10 commits behind head on master.

Files with missing lines Patch % Lines
src/simulation/m_ibm.fpp 48.48% 11 Missing and 6 partials ⚠️
Additional details and impacted files
@@            Coverage Diff             @@
##           master    #2209      +/-   ##
==========================================
+ Coverage   63.71%   65.16%   +1.45%     
==========================================
  Files          90      110      +20     
  Lines       24569    29571    +5002     
  Branches     3584     4478     +894     
==========================================
+ Hits        15653    19270    +3617     
- Misses       6135     6743     +608     
- Partials     2781     3558     +777     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:
  • ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

@sbryngelson

Copy link
Copy Markdown
Member Author

Thanks. Taking the four findings in order.

  1. T_IP = pres_GP*...: in the branches where the ghost pressure changes, rho is recomputed from the ghost partial densities just above (s_compute_mixture_coefficients(alpha_rho_GP, ...)), so this is the image-point temperature again, not a mix of the ghost pressure with the image-point density; where the pressure is not changed the two are identical. I tested it with chemistry on a slip circle, an accelerating circle and an injecting grain (h2o2, 1200 K): the ghost temperature equals the mirror temperature (ratio 1.0000, range 0.9986 to 1.0012) in master and this PR, with the same ratio for the grain (0.874) in both. The PR body now says so.
  2. kin_model > 0: the carve-out is removed on purpose (its comment said the centroid acceleration was missing; the acceleration is now provided). Ran both shipped prescribed-kinematics examples: flapping plate at full size and pitch-up (shipped grid on the MI300A, coarse grid over the whole ramp on the A100), no NaN and no ghost pressure at the clamp; the added-mass sign is right and the magnitude is unvalidated. Added a step-0 guard (the acceleration is not evaluated at step 0, where step_vel is not yet a previous velocity), which moves the two 3D Pitch Ramp goldens by at most 1.2e-3 absolute in total energy (5e-4 of its largest value); those two goldens are regenerated.
  3. Stage-1 acceleration reused in stages 2 and 3: yes, first order in dt as documented; dt at stage 1 is the step's own dt under cfl_dt too. Left as is.
  4. Existing slip goldens change: yes, intended; the twelve regenerated goldens (and now the two Pitch Ramp ones) are listed in the body, generated on NVIDIA and checked on AMD.

This branch has not been deployed

No deployments
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Development

Successfully merging this pull request may close these issues.

1 participant