CelerisLab/docs/moving_bodies.md

# Moving circles on the uniform dense grid — S0 DRAFT contract

Status: **DRAFT, not frozen** for the force RHS, kernels and physical
acceptance. Host-only S1 (state, validation, guards; no kernels) is on branch
`claude/moving-bodies-s1`; positive-step runs with a moving body are rejected
until S2. When S2 lifts that guard, `run(n > 1)` with a pending target must
raise `ValueError` before mutation (4a), and `run(0)` must not upload a motion
record or consume the pending target. Base: main
`b2034f87e03656c5f67ab88562a70d6c338b9918`. Scope approved by Frank (via TRON,
thread C0C78UFKGRM/1793335557.460509): uniform-grid prescribed translation and
single-cylinder free transverse VIV. Out of scope: refinement/blocked path,
pinball assembly, EsoPull, shifted storage, D3Q19, DynamisLab, manuscript.

Items marked **OPEN** are undecided; items marked **CANDIDATE** are proposals.

## 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, `|Δc|_2 ≤ 0.5` cell, speed bound) before mutation; on success replaces `pending` |
| 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.

## 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.