CelerisLab/docs/moving_bodies.md

Moving circles on the uniform dense grid

Status and user guide (current)

Development-stage feature, integrated from the reviewed slices S1–S3a (production code identical to 653fd997). Physical evidence and its limits: docs/validation/moving_bodies_stage1.md. The force used as a structural RHS is a declared stage choice, not an established result.

What it does. Translating circles (with optional spin) on the dense double-buffer D2Q9 path: per-step cut links from an FP64 geometry record (M2 MovingCircleLinks), reclassification and refill of covered/uncovered cells (M1 MovingCircleCells), moving-wall Bouzidi with wall velocity u + omega x r, and separate per-body diagnostics.

API.

• Simulation.add_body("circle", center=(x, y), radius=r, translation=Translation(offset_box, max_speed)) before initialize(); offset_box = (x_min, x_max, y_min, y_max) in lattice cells around the base centre (must contain 0), 0 < max_speed ≤ .5.
• set_body(id, offset=(dx, dy), vx=..., vy=..., omega=...) sets a pending target (validated against the applied offset: box, Euclidean step ≤ .5 cell, speed, no sliding hold); vx/vy are the interval wall velocity.
• run(1) advances a pending target and commits it only after the step synchronizes with no_donor == 0; run(n > 1) with a pending target raises; without a target run(n) holds the geometry (spin allowed). Direct LBMStepper.step refuses while a target is pending.
• moving_diagnostics() returns per body the columns of MovingTables.DIAG (CME, Wen term, link mass, cover/uncover momentum and mass, eq_only, no_donor, crossings) accumulated since obs were zeroed; obs/read_force keep their contract (CME on the body).
• get_flags() reads device flags back when moving circles exist.

Supported: dense double buffer, D2Q9, FP32 unshifted, SRT/MRT, circles, spin. Rejected before mutation: EsoPull, ddf shifting, D3Q19, refinement / blocked steppers, Jets or surface actuation on a moving circle, sync_bodies and checkpoints while a moving circle exists, swept bands that touch the domain edge, other bodies, sensors or force regions.

Failure lifecycle. Device-side failures (upload, step, sync, no_donor) mark the Simulation failed; later run/step/sync/checkpoint calls and direct stepper calls refuse until close().

Known limitations. No mass correction at curved walls (static or moving); cover/uncover and Wen terms are diagnostics, not a selected force law; the host loop costs a fixed per-step synchronization (≈ 1 ms/step measured on a 7.2 M-cell grid, run(1) class); outer boundaries are the single dense channel package (see configs/CONFIG.md).

Design record (staged contract, kept as written)

The sections below are the contract the slices were reviewed against (S0 → S3a). They are kept unchanged for provenance; where they say "draft", "proposed" or "until S2/S3", the status above supersedes them.

Original scope (Frank via TRON, thread C0C78UFKGRM/1793335557.460509): uniform-grid prescribed translation and single-cylinder free transverse VIV; base main b2034f87e03656c5f67ab88562a70d6c338b9918.

1. Working standards

• Lean, readable kernels in the existing idiom (kernels/step/*.cu, shared apply_bouzidi_link); Python manages state and validation. No framework, no compatibility layer, no new dependency.
• The existing stationary/rotating path (CurvedBoundaryKernel, host link packing, sync_bodies) stays byte-identical. Moving circles are a separate, explicitly enabled path with its own numerical contract (this file).
• One authority per state (section 4). No per-step full-domain download, rebuild or recompile; local swept-band device updates are in scope.
• Small slices (S1 host spec, S2 kernels without motion, S3 motion + refill, S4 oscillator); each reviewed SHA is frozen before targeted runs.
• Tolerances, benchmark parameters and the force composition used in the structural RHS are declared before the runs that judge them. Fresh vs carried evidence and existing baseline failures are recorded; smoke, convergence and physical validation are reported separately. No undeclared force filtering.
• External code (waLBerla GPLv3, FluidX3D source-available, VIVSIM) is a concept reference only; nothing is copied.

2. Current update order (main b2034f87, dense double buffer)

Simulation.run(n): upload host action (if upload_act) → zero obs (if zero_obs) → DoubleBufferStepper.step(n) → download obs → synchronize.

Per step (lbm/stepper.py step):

1. CurvedBoundaryKernel over the packed cut links: reads post-collision G^{n-1} at k_f (and k_f - c for q < 1/2), adds the moving-wall term (Uw = -omega*ry, Vw = omega*rx + optional per-link surface_velocity), writes the reflected value into slot opp of the solid node k_s, and accumulates c(f_toward + f_reflected) and torque into obs.
2. OneStep pull/collide ddf → temp (collides FLUID cells; obstacle interiors flagged SOLID|OBSTACLE|BC_CURVED), then buffer swap.
3. ForceRegionKernel, SensorKernel.

Raw obs force is the force on the body (link derivation in section 6; DFG docs report Cd > 0 from it). Translation vx, vy in the action slots is currently stored but not read by any curved kernel.

3. Proposed update order (moving circles)

Per step n, before the existing step 1:

• M1 MovingCircleCells (one thread per swept-band cell): classify with the analytic predicates in_old = inside(x, c^{n-1}), in_new = inside(x, c^n); write flag = in_new ? SOLID|OBSTACLE|BC_CURVED : FLUID; covered cells (!in_old && in_new) accumulate node_cover; fresh cells (in_old && !in_new) are refilled (section 5) and accumulate node_uncover.
• M2 MovingCircleLinks (one thread per band cell × direction): compute q and hit point at c^n (FP64, same formula as find_circle_intersection, same clamp [0.01, 1] and fallback classes as CircleGeometry.build_cut_links), skip non-cut slots, call apply_bouzidi_link with u_w = (vx^n, vy^n) + omega × (hit - c^n), accumulate the link terms (section 6).

Then the unchanged steps 1–3. M2 must follow M1 because a q < 1/2 link may read k_f - c at a fresh cell.

Exclusivity: moving bodies are never packed into the CurvedBoundaryKernel link arrays (n_curved covers stationary/rotating bodies only), so M2 is the sole writer of their solid slots and the sole source of their link force; no double write and no double obs. The stationary path's packing, launch and arithmetic are unchanged.

Completeness: after M1, every (fluid cell, direction) pair in the swept band whose neighbour is in_new is visited by M2. This holds because the band is the circle AABB at the declared extreme offsets plus 2 cells, and the host rejects per-step moves with |c_new - c_old|_2 > 0.5 and bands touching the domain edge, other bodies, sensors or force regions. Each M2 write lands at index_f(k_s, opp), the exact slot the adjacent fluid cell pulls, so every fluid pull from a solid slot of a moving body is BFL-written.

4. State, masks and authority

• Host ObjectState is the single authority for one displacement offset = (dx, dy) relative to the immutable base centre, plus velocity = (vx, vy) and omega. Centres are always derived (c = base + offset), never stored as independent state. Velocity is supplied, never differenced from offsets.
• c^{n-1} is derived from the last committed offset, c^n from the offset being applied. Per body a small device motion record {c^{n-1}, c^n, u^n} is uploaded with the action buffer; action-slot meaning is not overloaded. (CANDIDATE layout; OPEN.)
• Fail-closed commit: an enqueued upload is not a completed one. The new offset becomes the committed offset only after the stream synchronization that follows the upload and the step group succeeds. Any upload or later synchronization failure leaves host and device motion state unknown: the simulation is marked unusable (existing _require_usable lifecycle), and no API may report a host/device match afterwards.

4a. Applied geometry vs pending target (CANDIDATE, narrow API)

Two distinct facts per moving body, one authority each: applied = offset of the geometry last applied on the device, changed only by a successful commit; pending = an optional requested target (offset, velocity), changed only by set_body. A setter returning proves nothing about the device. The uploaded motion record is derived each run: {c_old = base + applied, c_new = base + (pending.offset if pending is not None else applied), u = pending.velocity if pending is not None else 0} (explicit is not None, so a zero offset or zero velocity is a valid target). Translational wall velocity exists only on a step that advances to a pending target; a hold step (c_old == c_new) has translational u_w = 0 (spin unchanged), so a non-advancing circle never acts as a sliding wall. One moving step consumes one target.

situation	contract
set_body target	validated against applied (box, `
target with offset == applied and nonzero translational velocity	ValueError before mutation (would be a hold step with sliding wall velocity)
two setters before a run	last wins; each validated against applied, not against the earlier pending target
run(1) with pending	upload record, step, sync; on success applied := pending.offset (applied stores only the offset), pending := None
run(1) without pending	hold: c_old == c_new, no cover/uncover, links at applied, translational u_w = 0
run(n > 1) with pending	ValueError before mutation (single-step-only moving API)
run(n > 1), no pending	allowed: hold for n steps, translational u_w = 0, spin only
run(0)	no motion upload, no commit; pending kept
validation failure	ValueError, nothing mutated, retry allowed
upload/launch/sync failure after enqueue	simulation unusable; applied and pending are not reported as device truth; no re-entry

Displacement norm: Euclidean per-step |c_new - c_old|_2 ≤ 0.5 lattice cells (declared bound; CANDIDATE value).

• Membership is never read from a flag that M1 is overwriting. The host verifies at add_body that every swept-band cell has base flag FLUID, so the band's base state is known without reading flag_gpu.
• flag_gpu is the device authority for band cells while a moving body exists; the host flag mirror is stale for those cells and is not used.

5. Refill of fresh cells (CANDIDATE)

• Stored time level: everything in ddf at M1 is post-collision G^{n-1}; the refill writes all 9 slots at that same level, so every slot of a fresh cell that the following pull or a q < 1/2 link may read is defined. Donor values are read at the same level G^{n-1} (old-fluid cells, untouched by M1).
• G = f_eq(rho_bar, u_w(x)) + G_neq(x_d) where rho_bar = mean density of valid neighbours and x_d = the valid neighbour best aligned with the outward normal (x - c^n)/|x - c^n|; G_neq(x_d) = G(x_d) - f_eq(rho(x_d), u(x_d)).
• Valid donor: inside the domain, base FLUID, !in_old && !in_new.
• Sets: V(x) = the 8 lattice neighbours that are in the domain, base FLUID and !in_old && !in_new. rho_bar = mean density over V(x). NEQ donor = the member of V(x) with the largest c·n among those with c·n > 0 (n = outward normal at c^n); none → NEQ part 0 (counted as eq_only).
• V(x) empty → no defined density: the cell is refilled with f_eq(RHO, u_w(x)) so all 9 slots stay defined, the no_donor counter is incremented, and the host treats no_donor > 0 after the run as a hard failure (simulation unusable). It is not a silent fallback.
• Non-emptiness is NOT assumed from the displacement bound. Geometry assumptions: circle radius R ≥ R_min (CANDIDATE 4 cells), |Δc|_2 ≤ 0.5. S1 includes a host enumeration over sub-cell centre positions and displacement directions on a fine grid within these bounds that checks every fresh cell has nonempty V(x); it is evidence for those bounds, not a proof, and runtime keeps the counter.
• Fresh cells only write, old-fluid cells only read: no race.

6. Discrete momentum ledger (to review; not a selection)

Notation per step n: F^n = fluid set after M1; links L = (k_f, dir) at c^n; f_t = G_dir(k_f); f_r = value M2 writes into (k_s, opp). Fluid momentum before streaming:

P_pre = P(G^{n-1} on F^{n-1}) - sum_cover sum_q c_q G_q + sum_fresh sum_q c_q G_q^refill

During the pull, every population entering a fluid cell comes from a fluid cell except the BFL-written slots; f_t leaves into solid and is never pulled by a fluid cell. Collision conserves momentum for SRT/MRT without forcing. Hence, in exact arithmetic,

P_fluid^n - P_fluid^{n-1} = - sum_L c (f_t + f_r) - I_cover + I_uncover + B^n + S^n

with B^n external-boundary momentum (inlet/outlet/walls) and S^n forcing. So F_CME + node_cover - node_uncover is, by construction, the body momentum that balances the fluid ledger of this storage/update order. Branch details: halfway fallback links (f_r = f_t plus the moving-wall addend) and q < 1/2 links enter the link accumulators exactly like Bouzidi links, with c(f_t + f_r). In the q < 1/2 branch G_dir(k_f - c) is an interpolation input only, not a fluid-momentum sink. Pull-slot completeness is the section 3 contract. The scope excludes ddf_shifting (rejected), so the unshifted force form applies. The balance is a fluid-side identity; it does not close body+fluid globally in a domain with inlet/outlet/forcing (B^n, S^n), and the solid-slot writes are boundary plumbing, not physical streaming packets.

The Wen 2014 relative-velocity link sum adds -u_w (f_t - f_r) = u_w dm per link (dm = f_r - f_t, the BFL mass defect). That term has no counterpart in the fluid ledger above. Peng et al. (PRE 2016 submitted, §II.C, Eqs. 20–21) show for a straight wall that the bounce-back error accumulated before a node switch cancels with the switching momentum. Whether sum u_w dm approximates I_cover - I_uncover here (implying double counting if both are used) is NOT demonstrated for this implementation: it needs a bounded experiment (prescribed slow translation; cumulative sum u_w dm vs cumulative node impulse, per crossing cycle). Momentum conservation is also distinct from Galilean invariance, which G4 tests.

External reference note: waLBerla's covered/refill momentum switches (conserveMomentumWhenMapping, conserveMomentumWhenRestoring) default to false (Astra), so their existence does not establish which combination its validation cases use. A bookkeeping identity alone does not select the physical force either.

Mass ledger (same proposition status). With M^n the active-fluid mass summed over F^n at the sample point,

M^n - M^{n-1} = sum_L dm - m_cover + m_uncover + B_m^n

(dm = f_r - f_t, m_cover/m_uncover = summed sum_q G of covered and refilled cells, B_m external-boundary mass flux; collision conserves mass).

Sample point (fixed): post-collision storage after the swap of step n, i.e. the same level G that M1 reads at step n+1; both ledgers are evaluated there. Ledger residual = left side minus every listed term. Zero mass flux does not imply zero momentum exchange: closed walls (bounce-back, free-slip) still exchange momentum, so "no inlet/outlet, no forcing" does NOT give B = S = 0. The rounding-level momentum check therefore uses either a momentum-isolated case (periodic in both directions, no forcing) or counts the wall impulse explicitly as part of B; the mass check may use a closed domain. Otherwise B, S are not separately measured, the residual is reported as their estimate and is not a pass criterion. Both ledger identities are proposed, unreviewed: they are not called proven until review and a host replay agree.

Accumulators (separate, per body, per step window; all reported raw): link_cme, link_wen_correction (= sum u_w dm), node_cover, node_uncover (momentum), link_dm, mass_cover, mass_uncover (mass), plus counters of crossings, eq_only and no_donor, so force spikes can be correlated with crossing events.

OPEN — structural RHS candidates (not selected): (a) link_cme + node_cover - node_uncover; (b) link_cme + link_wen_correction (no node impulse). Selection follows the ledger review and the bounded experiment above.

7. Rejection semantics

• Before any mutation (ValueError, no device write, host state unchanged): moving body with EsoPull, ddf_shifting, refinement/blocked stepper, D3Q19, Jets or surface actuation on that body, sync_bodies/checkpoint while a moving body exists; swept band touching domain boundaries, other bodies, sensors or force regions; set_body offset outside the declared box, per-step move |c_new - c_old|_2 > 0.5 cell (the single bound of 4a), speed above the declared bound.
• Host motion state is committed only after the upload and the following synchronization of the step group succeed (section 4); any failure in between makes the simulation unusable and no host/device match is reported.
• The only soft device diagnostic is eq_only (NEQ part omitted, counted). no_donor > 0 is a hard failure: it is checked after the step group and before the step is committed or any results are returned; on failure the simulation is unusable (section 5).

8. Coupling and units (CANDIDATE)

Lattice units, rho0 = 1, per unit span. Step n uses c^n, u^n and yields F^n. Weak staggered update v^{n+1} = v^n + (F^n - c v^n - k y^n)/m, y^{n+1} = y^n + v^{n+1} (Newmark is an alternative; OPEN). Mass, natural frequency, damping and U_r follow the chosen benchmark's own definitions. Leontini, Thompson & Hourigan, JFS 22:857–864 (2006), DOI 10.1016/j.jfluidstructs.2006.04.003, as read by Grok from the full text: Re = 200, m* = m/m_f = 10 (m_f displaced fluid mass), zeta = c/c_crit = 0.01, U* = U/(f_N D) with f_N the natural frequency in fluid assuming an inviscid added-mass term; the paper also reports frequencies normalised by the in-vacuo value. Our runner derives k from that in-fluid f_N (added-mass coefficient as stated in the paper) and never from the vacuum formula; c_crit convention and the original comparison data are still OPEN. Stability/convergence is gated by fluid time/space refinement of the coupled run, not assumed from the mass ratio; structural sub-stepping does not substitute for fluid-coupling time refinement.

9. Gates (map; tolerances OPEN, declared before runs)

G0 zero-motion regression, split: G0a populations of the existing stationary and spin-only path bitwise unchanged vs main; G0b moving path at rest (with and without spin) reproduces the stationary link set and populations; force estimator outputs are compared separately from populations, since the moving path adds estimator terms · G1 raw force sign (Cd > 0) · G2/G3 prescribed translation incl. ledger accumulators · G4 matched relative-flow (moving body in quiescent fluid vs fixed body in uniform stream), optional boost, co-moving sanity · G5 Dütsch et al. 1998 Re 100, KC 5 in-line (needs 2-component offset) · G6 Leontini free transverse VIV (m* = 10, ζ = 0.01, Re = 200, definitions pending) · G7 time-level contract frozen before the oscillator.

Candidate pre-run engineering tolerances (Astra; for discussion, not literature limits or measured outcomes; norms, averaging window, convergence protocol and physical reference uncertainty must be fixed before adoption, and none may be weakened after seeing runs): normalised momentum-ledger residual FP32 ≤ 5e-6 (FP64 oracle ≤ 2e-13); 10-cycle relative mass drift ≤ 1e-4; prescribed-motion grid study D = 40/80/160 with change in the force-response amplitude (observed quantity and norm to be named before adoption) ≤ 2 % and phase ≤ 2°; enlarged domain amplitude change ≤ 1 %; matched-frame force RMS ≤ 2 % and improving with refinement; VIV refinement amplitude ≤ 2 %, frequency ≤ 1 %, five consecutive cycle amplitudes within 1 %.

Conditions stated before any moving kernel (S2/S3), not only at the oscillator stage: (i) fluid refinement protocol for the prescribed-motion gates — diffusive scaling (dx halved, dt quartered, physical Re, KC and amplitude fixed), the observed quantities and norms named in advance; (ii) Galilean constant-boost G4-boost requires boundaries compatible with a uniform boost (periodic in the boost direction or boundaries that move with it); if the current dense configuration cannot represent that, G4-boost is recorded as not runnable rather than replaced by another test.

10. Proposed minimal S1 scope (host only; held until cross-review)

• add_body moving-circle declaration: base centre, declared offset box (x/y), speed bound; swept-band computation and band-placement validation.
• set_body(offset=..., velocity=..., omega=...) validation before mutation (section 7), applied/pending state and the run-transition table of 4a, and the fail-closed commit of section 4.
• Host enumeration check of refill-donor non-emptiness (section 5) within the declared R_min and displacement bound.
• Motion-record host layout and upload path; exclusion of moving bodies from the packed CurvedBoundaryKernel links (stationary packing unchanged).
• Explicit rejections for unsupported modes.
• S1 execution guard: until the S2 moving kernels exist, any run(n) with n > 0 (and step(n)) while a moving body is declared raises before any device mutation (no action/motion upload, no launch), including the nominal at-rest case, because the moving body is excluded from the static packed links and nothing else would apply its boundary. run(0) keeps its existing behaviour. Simulations without moving bodies are unaffected. Preparation/commit helpers are exercised by host tests with mocks; no operational unsupported stepping path is exposed.
• Host tests (Astra runs them): rejection-before-mutation, band geometry, derived centres, fail-closed behaviour under an injected upload/sync failure, stationary packing byte-identical, execution guard (positive-step run with a moving body raises before upload/launch; stationary unaffected). No kernels, no GPU runs, no y-only VIV oscillator (S4), no force RHS.

10a. S2a — fixed-centre moving links (implemented on claude/moving-bodies-s2a)

• MovingCircleLinks (kernels/step/moving_circle.cu, dense double buffer) runs before the static curved kernel when moving tables are attached: one thread per (band cell, direction) slot, FP64 {cx, cy, r} at the applied offset, q/hit/clamp/fallback in the order of build_cut_links with explicit _rn operations (host/device bit parity is a test target, not a claim), then the unchanged apply_bouzidi_link with spin wall velocity only.
• obs keeps its contract (CME and torque from apply_bouzidi_link). A separate per-body diagnostics buffer receives cme, the Wen term u_w * dm, link_dm and the link count from the same link outputs; it is zeroed with obs and only exists when moving circles do. Optional slot_debug (q, fallback, rx, ry per slot) is null in production.
• Displacement is disabled: set_body(offset=..., vx=..., vy=...) on a translating circle raises before mutation, before and after initialize, so no pending target can exist and run/step/direct LBMStepper.step all hold at the applied geometry. initialize() never reads or consumes pending. stepper.blocked is cleared only when the moving tables are attached (initialize, stepper recreation).
• S3 rule (not implemented): once targets are re-enabled, any pending target must raise in run/step and in LBMStepper.step alike before any launch, swap or step-count change; zero-step calls never consume it.
• Not validation: a pass of the zero-motion gates says nothing about moving geometry. GPU parity (G0a/G0b) is unverified until run on a free device.

10b. S3a — translation with refill (implemented on claude/moving-bodies-s3a)

• Targets are enabled again (4a). set_body(offset, vx, vy) validates against applied; a pending target blocks direct LBMStepper.step(n > 0) before any launch; run(n > 1) with a target raises before mutation; run(0) never consumes it; run(1) advances it.
• Per advancing step, on one stream: upload the immutable FP64 records {c_old, c_new, r, u} (only when changed) → M1 MovingCircleCells → M2 MovingCircleLinks (u_w = u + omega x r at the hit point) → static curved → OneStep → swap → force region → sensors → synchronize → read no_donor → commit (applied := pending.offset) → upload the hold record (c_old = c_new = applied, u = 0) → synchronize. Hold steps launch no M1 and are the S2a path. Commit point: the advancing step-group sync plus a zero no_donor count; a failure before it commits nothing. The hold-record upload stages the next step after the commit: if it (or its sync) fails, applied already matches the moved device geometry, but the Simulation is still marked failed and refuses reuse. Any failure from the first upload on marks the circles and the Simulation failed.
• M1 classifies with the old/new predicates only (no flag reads). Covered cells add raw sum c G / sum G of G^{n-1} to cover_* and become SOLID|OBSTACLE|BC_CURVED. Uncovered cells: donors are fluid in both geometries; rho_bar over all of them; NEQ from the donor with the largest unnormalized c·n > 0 (first direction wins ties, n from c_new); u_w (translation + spin) at the cell position; all nine slots written; raw sum c G, sum G of the refill go to uncover_*; no donor → f_eq(RHO, u_w) plus no_donor (hard failure before commit); no NEQ donor → eq_only.
• get_flags() reads the device flags back (synchronized) when moving tables exist; otherwise it returns the host array as before.
• Not selected: the structural force. CME+node vs CME+Wen stay OPEN until X1 and G4 have run.

11. Remaining decisions

Structural RHS (section 6); motion-record layout; refill details; coupling scheme; benchmark definitions and digitized reference data; tolerances; execution/resource feasibility (Astra, pending); whether a quiescent domain (zero inlet velocity) runs cleanly in the current config.