Repository navigation
Moving IBM: RK state weights and centroid/angle stage velocities are incorrect #2025
Description
The force/trajectory review in #2017 exposed a pre-existing error in the immersed-body Runge–Kutta update. With moving_ibm=2, kin_model=0, the body uses different RK3 state weights from the fluid, and centroid/angle increments use velocities after their stage update. This can introduce a net impulse in an otherwise equal-and-opposite force exchange and reduce trajectory convergence to first order.
@danieljvickers — opening this separately from the face-flux force change so the integrator correction can be reviewed and validated independently. Investigation and source-extracted tests were performed with Codex.
Source and mechanism
Checked against #2017 base a6add01e35239cf9e87de38972e968ad1799c04c; these errors predate that PR. See the RK coefficients, fluid update, and body update.
For step-start state y0, stage input y, and RHS f, the fluid uses:
stage 1: y = y0 + dt*f0
stage 2: y = (y + 3*y0 + dt*f1)/4
stage 3: y = (2*y + y0 + 2*dt*f2)/3
The body saves step_vel at stage 1, but places that step-start state under coefficient 1 and the current stage under coefficient 2. Its stage-2/3 velocity blends are consequently (v0 + 3*v + dt*a1)/4 and (2*v0 + v + 2*dt*a2)/3. The same reversed state weighting appears in angular velocity, centroid, and angles.
For equal-and-opposite stage forces, the original RK3 updates imply:
body impulse: dt*(F0/4 + F1/12 + 2*F2/3)
fluid impulse: -dt*(F0/6 + F1/6 + 2*F2/3)
net impulse: dt*(F0-F1)/12
This residual exists even before considering moving-mask effects or external forces. Constant force conceals the velocity-weight defect, but exposes the second problem: centroid and angle use the newly updated velocity/angular velocity, whereas their RHS must use the input-stage values. With constant acceleration from rest, the original first step gives displacement a*dt^2 instead of a*dt^2/2. Swapping weights alone leaves the position update first order; this kinematic problem also affects RK2.
Reproduction and verification
Standalone reproducers, candidate CPU diff, and results are available. With Python 3 and gfortran on PATH:
gh gist clone 7efff43a8c45c0137f0b53b56e50c0d6
cd 7efff43a8c45c0137f0b53b56e50c0d6
python3 check-rk.py
python3 minimal-rk-replay.pycheck-rk.py compiles the extracted body routine and tests original, weights-only, and complete corrections against an independent Butcher-tableau implementation. Geometry updates are stubbed. The 432 configurations cover RK1/2/3, six force laws, four resolutions, uniform/varying time steps, nonzero initial states, unequal mass/inertia, translation, and 2D constant-inertia rotation. They include a surrogate fluid exchanging opposite forces with the body.
- Complete correction agrees with the independent RK implementation within
1.31e-14and recovers the expected temporal orders. - Maximum corrected exchange momentum drift is
3.55e-15. The representative 20-step RK3 exchange gives original drift0.004016637versus corrected drift3.33e-16. - For
x'=v,v'=-x, initial(x,v)=(1,0), att=1:
| steps | original position error | weights-only error | complete correction |
|---|---|---|---|
| 20 | 1.498e-2 | 1.418e-2 | 2.987e-6 |
| 40 | 7.593e-3 | 7.144e-3 | 3.626e-7 |
| 80 | 3.822e-3 | 3.585e-3 | 4.465e-8 |
| 160 | 1.917e-3 | 1.796e-3 | 5.539e-9 |
The extracted checks passed runtime bounds/floating-point checks and optimization comparisons. Eighteen prescribed-branch/stage-time checks verify dispatch using an analytic stub; they do not validate the physical flapping/pitch formulas.
Proposed correction and completion criteria
Save each body's input-stage velocity/angular velocity before modifying them. In each blend, place the current stage under rk_coef(s,1) and the step-start state under rk_coef(s,2). Use saved input-stage velocities for the centroid/angle RHS. The published rk-candidate.diff demonstrates this in the extracted CPU routine; it is not a production/offload patch.
A production fix should preserve the prescribed-motion dispatch and declare any new per-body temporaries appropriately for GPU execution. kin_model=0 disables special prescribed models; it does not mean stationary. The force-coupled moving_ibm=2 branch uses it. Prescribed kin_model=1/2 branches bypass this recurrence.
Completion should include the small RK/order/exchange checks, a full coupled moving-body run with stage-resolved momentum accounting, and CPU/GPU plus relevant MPI regression checks. General 3D rotation/inertia requires separate validation; the extracted tests establish translation and 2D constant-inertia rotation only. A golden-file update alone would not verify the correction.
Moving-mask momentum accounting remains a separate concern, related to #1636. The surrogate exchange test does not establish conservation through cell reclassification or physical moving-body convergence. See the full force/RK discussion on #2017.
Activity
sbryngelson commentedon Oct 10, 2026
@danieljvickers — here is the production-source draft, with the patch, standalone tests, supporting literature, validation, and caveats.
Attachments:
This is a production-source draft of the correction tested during the PR #2017 force review, addressing issue #2025. Prepared with Codex. It is a draft, not a merged or fully coupled-validated fix.
moving-ib-rk.patchchanges onlysrc/simulation/m_time_steppers.fpp. It applies cleanly to both the reviewed basea6add01e35239cf9e87de38972e968ad1799c04cand PR heade85a5b0ff56fce105fdc5a9317f36ba4353e2179. It was prepared in an isolated worktree; the force diagnostic changes used in the PR review are excluded.What changes and why
For force-coupled translation and constant-inertia planar rotation, the body equations are a first-order system:
x'=v,v'=F/m,theta'=omega,omega'=T/I. All components need the intended RK stage state and RHS.The draft:
rk_coef(s,1)and the step-start state underrk_coef(s,2), matching the fluid's SSPRK recurrence. This corrects the blends for velocity, angular velocity, centroid, and angles.privateclause. Both OpenACC and OpenMP expansions were checked.mib_analytical()override for the kinematic RHS. The earlier extracted-routine candidate did not account for this nonempty macro. This draft does. Special prescribedkin_model=1/2motion continues to use its separate branch.For RK3 the intended recurrence is:
The original body velocity blends reverse the current/start weights. With equal-and-opposite stage forces this produces body impulse
h*(F0/4+F1/12+2*F2/3)while the fluid receives-h*(F0/6+F1/6+2*F2/3). Their sum ish*(F0-F1)/12. With the corrected weights, both use weights(1/6,1/6,2/3), so this exchange cancels, conditional on matching stage forces and a fixed bookkeeping domain.Using the just-updated velocity in the position RHS is a separate error. From rest under constant acceleration the original first step advances position by
a*h^2, rather thana*h^2/2for RK2/3. Swapping the weights alone leaves the trajectory first order. The input-stage snapshot corrects the ODE recurrence, including RK2. RK1 remains forward Euler and retains its expected first-order accuracy.Literature supporting this correction
Verification completed
GNU Fortran 15.2; successful configure, build, and install. Build log in the ZIP:
rk-production-build.log.1.31e-14and recovers the expected orders. Maximum corrected exchange momentum drift is3.55e-15.0.004016637; corrected drift is3.33e-16.4.22e-15. These hook tests explicitly supplymytimeat the RK RHS nodes; they do not validate the full solver's scheduling of that clock.-O0and-O2, bounds checking, and invalid/divide-by-zero/overflow traps. OpenACC/OpenMP preprocessing emits the new arrays as private.git diff --checkand all seven./mfc.sh precheck -j 8gates pass. Patch applicability was checked on both pinned base and head. The precheck bootstrap emitted a Fypp linemarker-patch reapplication warning; the generated module, compilation, and all checks still passed. Logs are included.The standalone checks require Python 3 and
gfortranon PATH:rk-reproducer.f90contains the original, weights-only, and production-corrected routines with small force/geometry stubs.production-routine.f90is the generated corrected routine.provenance.jsonrecords revision pins and the patch hash.Caveats and remaining completion criteria
mytime; the draft preserves that convention. Supplying the correct RK node time in a hook test is not proof that the coupled solver already does so. Special prescribed motion's dispatch tests also use stubs.To apply in a compatible clean MFC checkout: