CelerisLab/docs/moving_bodies.md

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

Status: **DRAFT for review, not frozen.** No implementation exists. 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 of 1 cell or more 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.
- 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`.
- No valid donor for the NEQ part → equilibrium only (counted). No valid
  neighbour at all → cannot happen when the per-step move is < 1 cell
  (validated on host); counted anyway.
- 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.

Accumulators (separate, per body, per step window; all reported raw):
`link_cme`, `link_wen_correction` (= `sum u_w dm`), `node_cover`,
`node_uncover`, plus counters of crossings and refill fallbacks, 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 ≥ 1 cell, speed above the declared bound.
- Host state is committed only after a successful upload (same contract as Jet
  amplitude updates); an uncertain device-write failure makes the simulation
  unusable.
- Device-side soft diagnostics (refill fallbacks) are counters, not errors.

## 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 (stationary path unchanged; moving path at rest
reproduces stationary link set and DDF) · 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 final-step amplitude change ≤ 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 %.

## 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) and the fail-closed commit of section 4.
- Motion-record host layout and upload path; exclusion of moving bodies from
  the packed `CurvedBoundaryKernel` links (stationary packing unchanged).
- Explicit rejections for unsupported modes.
- 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. 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.