Activity

1. sbryngelson

   added 2 commits that reference this issue on Jul 10, 2026

   test: soften moving-IB circle case to avoid the un-GCL'd CFL blow-up (M…

   ...

   4c877aa

   test: keep two-body moving-IB (43AF9F25) at its original wall speed

   ...

   5e51d3f
2. sbryngelson

   added

   enhancementNew feature or requestNew feature or request

   on Jul 12, 2026
3. sbryngelson commented on Oct 10, 2026

   

   sbryngelson

   on Oct 10, 2026

   MemberAuthor

   More actions

   @danieljvickers — posting the full moving-body investigation here, with links to #2017 (force), #1958 (second-order wall treatment), #2025 (body RK defect), and #2026 (bilinear/trilinear interpolation).

   These are separate corrections. The tested version of #1958 and our body RK correction do not fix the moving-body force/transport inconsistency below: the counterexample still fails with both applied. This issue tracks the geometric-conservation work; it is not itself an implemented fix. Swept-volume/state transport addresses changing masks, while a consistent relative-interface flux is also needed to recover the physical moving-wall load. Bilinear/trilinear interpolation was not included in these tests, and improving interpolation alone cannot remove the continuum advective term derived below.

   The existing RK draft, supporting notes, and native attachments are available in the #2025 comment: RK ZIP and RK patch.

   The report below contains the actual solver tests, exact solution, stage-budget derivation, literature, and qualifications. The new 55-case bundle and diagnostic build patches are prepared locally; their native attachment upload is still pending. The diagnostic patches are reproducibility/instrumentation patches, not a general production fix.

   Moving IBM: exact-solution tests and stage budgets

   42 additional full serial CPU solver cases completed, following the 13 uniform-state cases. The principal finding is a moving-body force inconsistency in PR #2017: its raw fixed-grid momentum flux includes advective transport that is not the hydrodynamic load on a moving impermeable body. An exact zero-force solution reproduces the error, a uniform boost changes the reported force, and two-way coupling turns it into spurious body motion.

   Source and test scope

   • Base: a6add01e35239cf9e87de38972e968ad1799c04c.
   • PR Take the IB force from the momentum fluxes across the body's faces #2017: e85a5b0ff56fce105fdc5a9317f36ba4353e2179.
   • PR 2nd order IBM #1958: 82bd69ffabe040a8efde9e9985c527b0605a445c; its diff applies cleanly on PR2017 plus our RK correction and diagnostic hooks.
   • All moving tests use the corrected body RK draft. The original production draft worktree retains that source change only. Diagnostics and PR1958 are in the separate moving-budget-worktree.
   • Double precision, serial CPU, fastmath off. No GPU/MPI/AMR or general rotating-body validation.
   • 15 cases with PR2017 + RK + stage diagnostics; 15 corresponding cases with PR1958 enabled and minimum IP distance 1.5h; three mirror-IP controls with PR1958; five angular-density moving/two-way cases; two old-force controls; two resting-frame angular-density controls (one static flag and one using the moving-IBM path with zero velocity).

   Exact zero-force counterexample

   The periodic unit square contains a slip circle of radius R=0.12 centered at (0.35+Ut, 0.5), with U=0.25. The ideal gas has gamma=1.4, pressure p=2, and uniform velocity (U,0). Define X and Y as minimum-image distances from the moving center and r²=X²+Y². The density is

   $$ \rho=1+A\frac{X}{R}\exp\left[-\frac{r^2-R^2}{2R^2}-\left(\frac{r^2-R^2}{\ell^2}\right)^2\right], \qquad A=0.1,\quad \ell^2=0.025. $$

   Its periodic-seam tails are below double precision. Density stays between 0.9 and 1.1. On the body, rho=1+A cos(theta), and its radial derivative vanishes. The zero-gradient density condition is therefore compatible at the wall. This smooth density translates with the body while p and velocity remain constant. Continuity is rho_t+U rho_x=0; momentum and energy reduce to the same advection identity. The fluid and body have equal normal velocities, so impermeability holds. Uniform pressure gives exactly zero net force and torque, for any body mass.

   Initial conservative states use five-point Gaussian cell averages: rho, rho U, zero y-momentum, p/(gamma-1)+rho U²/2, and alpha=1. Comparing five- and nine-point quadrature gives a maximum angular-profile discrepancy of 1.4e-12 at 32² and below 2e-15 at 64²/128². The full solver then runs to t=0.5, translating by 4, 8, or 16 cells.

   With PR2017 + RK + PR1958 (IP minimum 1.5h):

   Grid	First-stage reported Fx	RK-weighted time-mean Fx	Exact Fx
   32²	-0.00236224	-0.00160655	0
   64²	-0.00235878	-0.00212972	0
   128²	-0.00234585	-0.00227795	0

   Pressure and velocity remain at roundoff in these prescribed-motion cases. The false force is already present at the first RHS evaluation of the supplied exact initial state; it cannot be attributed to accumulated body RK error.

   For a normal pointing from body to fluid, the raw Euler momentum flux contains rho u(u dot n)+p n. Here the advective contribution predicts

   $$ F_{x,\mathrm{raw}}=-U^2\int_{\partial B}\rho n_x,ds=-U^2 A\pi R=-0.002356194490192345. $$

   The measured mean approaches this nonzero value. Moreover, at every measured stage in the prescribed-motion cases, removing the advective contribution diagnostically (Fx + U * integrated_fluid_mass_RHS) reduces the force to roundoff. This subtraction is a mechanism check for this exact constant-velocity solution, not a validated general correction.

   In the resting frame, both body and fluid velocity are zero and the same angular-density profile is stationary. The PR2017 force is then zero to roundoff, including a control that keeps moving_ibm=1 and exercises the same moving-IBM code path at zero speed. Its force therefore changes under a uniform Galilean boost even though the physical load is unchanged.

   Two-way coupling exposes the consequence. At 64², the exact body velocity remains 0.25 and its final x is 0.475. With PR2017 + RK + PR1958, mass 1 gives final velocity 0.248985434 and x=0.47474277; mass 0.01 gives velocity 0.231208154 and x=0.47013750. The light-body pressure error reaches 0.00485469. Old pressure-gradient-force controls, using the corrected RK build and the original wall treatment, retain body velocity and position to roundoff; their force is below 3e-16. The exact uniform pressure/velocity makes the force mechanism comparison meaningful despite the wall-treatment difference. The diagnostic subtraction in the combined-PR build independently isolates the advective contribution.

   What the revised force does verify

   On the fixed mask used for each RHS, PR2017 balances integrated fluid momentum RHS and body force to approximately 1e-15. This is a useful discrete compatibility property. It does not include the transport associated with subsequently moving the domain represented by that mask.

   Let Q0 be a component of the fluid total on the stage-input mask, Qstar_old the pre-reconstruction state integrated on that old mask, Qstar_new the same state on the rebuilt mask, and Qfinal_new the corrected state on the new mask. The stage difference decomposes into:

   1. Fixed-mask state evolution, Qstar_old-Q0.
   2. Reclassification, J=Qstar_new-Qstar_old.
   3. Correction on the new fluid set, C=Qfinal_new-Qstar_new.

   Here C=0: the correction writes solid/ghost cells, not cells classified as fluid. It can still alter the states later exposed by motion.

   RK also blends the step-start conservative state using the stage's old mask. Define T_s as the step-start state integrated on that stage mask minus the step-start state integrated on the original step mask. For RK3 the exact full-step identity is

   $$ Q_{n+1}-Q_n=\Delta t\left(\frac{R_1}{6}+\frac{R_2}{6}+\frac{2R_3}{3}\right) +\frac{J_1}{6}+\frac{2J_2}{3}+J_3 +\frac{T_2}{2}+\frac{T_3}{3} +\frac{C_1}{6}+\frac{2C_2}{3}+C_3. $$

   This identity was checked on measured stage totals for continuity, both momentum components, energy, and alpha, with closures below 6e-14. Thus the extra terms are measured consequences of reclassification and RK state reuse, not an unexplained residual.

   Nonzero-exchange case

   A prescribed U=0.25 body translates through initially still, uniform fluid with p=2, rho=1, in the periodic box. This produces a physical response; pressure changes in this case are not automatically spurious.

   At 64² with both PRs, over t=0.5:

   • Discrete fluid x-momentum change: +0.0105903600.
   • RK-integrated reported body-force impulse: -0.0109229005.
   • Sum: -0.000332540524. The measured geometry/RK terms reproduce this to roundoff.
   • With half dt, the sum is -0.000339411164.
   • On steps with no mask changes, the momentum residual is below 9e-17.
   • The whole-cell fluid mass increases by 0.00450311973; halving dt gives 0.00450473346. Initial and final binary represented volumes agree because the body translates by an integer number of cells.

   These are exact diagnostics of MFC's whole-cell masked totals. Such totals approximate the physical moving-domain integral; the test is not an independent exact cut-volume quadrature. Nevertheless, the results establish that the implemented fixed-mask flux balance is not the full moving-mask momentum balance. A changed subcell offset changes the residual and creates a y component; time refinement alone does not remove it.

   Energy is also recorded. In this case the fluid energy change plus integrated force work is 0.0309720836. The diagnostic splits this into the energy-RHS/force-work difference (0.0361837793) and the geometry/RK contribution (-0.00521169578). The raw numerical momentum-flux force is not yet a physical traction/work pair for a moving wall; these numbers should not be promoted to a standalone exact physical-energy error without that qualification.

   Density accuracy and image-point controls

   A separate radially symmetric profile, rho=1+0.1 exp(-((r²-R²)/0.025)²), also translates exactly with uniform pressure and velocity. At 32²/64²/128² the existing mirror-IP wall gives maximum density errors 0.0159444/0.00506278/0.00170822, corresponding to observed rates 1.66 and 1.57. The matching no-body errors are 0.00223610/0.000326823/0.0000223032.

   PR1958 with mirror IP gives essentially the same density errors, as expected for an already exact constant velocity. Its minimum IP distance of 1.5h changes density sampling and gives errors 0.0298666/0.00879898/0.00279631 (rates 1.76 and 1.65). This test does not evaluate the intended improvement of velocity extrapolation in a sheared flow. Larger image-point distance improves stencil availability but changes the density reconstruction error here. Halving dt at 64² barely changes the density error. Three grids do not establish asymptotic order, and Dan's bilinear/trilinear work is not applied.

   Changing prescribed-body mass from 1 to 1000 produces 216 byte-identical snapshot files in both radial-density and nonzero-exchange controls. The pressure acceleration correction is guarded by moving_ibm>1; prescribed motion does not incorrectly use hydrodynamic force/mass.

   Fix direction supported by the tests

   For a moving interface, the relative flux is F_n-(v_wall dot n)Q. Impermeability then gives zero mass flux, pressure traction in the inviscid momentum exchange, and pressure work p(v_wall dot n) in energy. This removes the advective force in the exact counterexample. A general implementation must combine that exchange with swept-volume/state transport and RK geometry/state bookkeeping; rotation makes wall velocity vary over the surface. Simply giving all overwritten ghost momentum to the body is not a justified correction.

   Monasse et al. explicitly separate pressure-force exchange from swept conserved quantities and construct a conservative compressible rigid-body coupling: paper, sections 4.2–4.4 and 5. This is the relevant direction for #1636 and moving-body support in #2017. Corrected RK and higher-order interpolation remain useful, but neither removes the raw advective transport term by itself.

   Artifacts and reproduction

   run_nonuniform.py runs each 15-case suite. run_extra_tests.py reproduces the mirror, angular-density, old-force, and resting-frame controls. analyze_nonuniform.py closes stage/full-step budgets. verify_asymmetric.py checks the force signature, body motion, controls, and quadrature. plot_nonuniform.py and plot_asymmetric.py produce PNG/PDF figures.

   The preprocessor is the base+RK build. Base+RK simulation supplies old-force controls. For PR2017, apply budget-build.diff on the pinned PR head, build serial simulation through ./mfc.sh, and preserve that executable as binaries/pr2017-rk-budget-simulation. Then apply pr-1958.diff and rebuild for combined-PR tests. Alternatively, combined-budget-build.diff applies the entire combined source change directly to the pinned PR2017 head. These diagnostic patches are serial CPU test instruments, not proposed production changes.

   The original rk-production-draft and the isolated moving-budget-worktree are siblings of this artifact directory. From their parent:

   rk-production-draft/build/venv/bin/python moving-ib-robustness/run_nonuniform.py
   rk-production-draft/build/venv/bin/python moving-ib-robustness/analyze_nonuniform.py
   IB_SECOND_ORDER=1 rk-production-draft/build/venv/bin/python moving-ib-robustness/run_nonuniform.py
   IB_SECOND_ORDER=1 rk-production-draft/build/venv/bin/python moving-ib-robustness/analyze_nonuniform.py
   rk-production-draft/build/venv/bin/python moving-ib-robustness/run_extra_tests.py
   rk-production-draft/build/venv/bin/python moving-ib-robustness/verify_asymmetric.py

   Per-case inputs, logs, force histories, full stage CSVs, and analysis are retained. Provenance JSON records source diffs and binary hashes. The compact ZIP includes initial/final snapshots; intermediate snapshots remain in the local case directories. Findings are from actual solver runs, not the earlier interpolation arithmetic example. No new findings have been posted publicly in this continuation.