docs: Dütsch protocol — adopt Grok's cycle-index labels (K = first scored cycle)

Same cycles and pairs as 8505a723, relabelled: proof pairs 3→4..5→6 give
K_case = 7, window K..K+2 with remain-settled pairs, last scored cycle 11.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

docs/moving_bodies.md

@@ -363,6 +363,38 @@ recorded as not runnable rather than replaced by another test.
 - 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

docs/validation/duetsch_re100_kc5_protocol.md

@@ -0,0 +1,209 @@
+# Prescribed in-line oscillation: Dütsch Re 100 / KC 5 — PROTOCOL (review draft)
+
+Production under test: `claude/moving-bodies-s3a` at `653fd997` (unchanged).
+Primary source: Dütsch, Durst, Becker & Lienhart, J. Fluid Mech. 360, 249–271
+(1998), doi:10.1017/S002211209800860X. Source facts below were verified from
+the paper (Slack C0C78UFKGRM/1793335557.460509, TRON 1793356317–1793356724,
+page numbers given). Everything marked ENGINEERING is our a-priori choice,
+fixed before any benchmark result; it is not a paper uncertainty.
+
+## 0. Source facts (with provenance)
+
+- Motion (p.251): `x_c = -A sin(2 pi f t)`; paper phase phi = 0 is the mean
+  position moving left at `-U_max`. KC = 5 gives `A/D = 5/(2 pi) = .7957747`.
+- Start (p.254): the computation starts a quarter cycle before phi = 0, at
+  `x = +A` with zero velocity (phi = -90°).
+- Experiment (pp.252, 254–256): D = 10 mm, L = .65 m, tank 1.2 × .7 × .7 m,
+  rho = 998.2 kg/m³, mu = 1.00818 mPa·s, Re = 100, KC = 5; LDA with 1° phase
+  bins and 200 samples per node and degree; no quantified experimental error
+  recovered.
+- Profiles (Fig. 8, p.260): LAB-FIXED stations x1 = -6, 0, 6, 12 mm, i.e.
+  x/D = -.6, 0, .6, 1.2, at phi = 180°, 210°, 330°; u1 and u2 versus x2.
+- Original numerics (§2): cylinder-fixed accelerated frame (`x_i = x~_i +
+  x_i,s`, outer velocity `-u_c`, inertial source `-rho du_c/dt`), O-grid outer
+  diameter 120 D, illustrated 192 circumferential nodes, 720 steps/cycle; tank
+  influence neglected.
+- Force fit (verified): force ON THE BODY per unit span
+  `F1 = -0.5 rho D c_d xdot|xdot| - (pi/4) rho D² c_i xddot` with `x = -A sin phi`;
+  normalized by `0.5 rho D U_max²`: `c_d cos phi |cos phi| - (pi² c_i / KC) sin phi`.
+  `c_i` is the ADDED-MASS coefficient; total inertia `c_m = c_i + 1` is a
+  different convention and is not used.
+- Table 1 computed Morison coefficients (c_d, c_i) for 96 / 192 / 384 perimeter
+  nodes at dt/T = 1/720: (2.15, 1.42), (2.10, 1.45), (2.09, 1.45). These are
+  fitted computed coefficients, not a measured force and not an error bar. The
+  paper states the coarse direct force is closer to the finest direct force
+  than the Morison fit, but gives no waveform-error norm.
+
+## 1. Kinematics, time levels and phase bookkeeping
+
+- Our run is in the LAB frame: the circle moves, the fluid is initially at rest
+  (rho = 1, u = 0); no fictitious inertial source is added. Lab-fixed profile
+  stations are compared directly (no frame transform needed).
+- Geometry reference: the circle's base centre is the INITIAL position, the
+  right turning point `x_mean + A`; the declared offset box is x in [-2A, 0],
+  y = 0. The body is never teleported: offset 0 at s = 0.
+- Elapsed time s = n dt from the start. Analytic position relative to the mean
+  `x(s) = A cos(2 pi s / T)`, velocity `u(s) = -U_max sin(2 pi s / T)`; paper
+  phase `phi = (360 s / T - 90) mod 360`.
+- Step n moves from s_n to s_{n+1}: target offset `x(s_{n+1}) - A`; wall
+  velocity = the single-step interval average `(x(s_{n+1}) - x(s_n)) / dt`
+  (analytic, supplied, not differenced from stored offsets). The per-step
+  force of `run(1)` / `read_force(normalize=True)` is attributed to that same
+  interval.
+- Paper cycle k = 1, 2, ... spans `s in [T/4 + (k-1) T, T/4 + k T)`; the
+  initial quarter cycle (s < T/4) is cycle 0 and never scored. Every logged
+  sample carries (k, phi). Example: phi = 330° of cycle 1 is at s = 7T/6.
+- T is chosen as a multiple of 360 native steps, so phi = 180°, 210°, 330° fall
+  exactly on step boundaries and 1° bins hold whole steps.
+
+## 2. Resolution ladder (uniform grid, fixed Mach), ENGINEERING
+
+| D (cells) | T (steps) | U_max = 5 D / T | nu_lat = U_max D / 100 | tau (3 nu + .5) |
+|---|---|---|---|---|
+| 20 | 2,520 | .039683 | .0079365 | .5238 |
+| 40 | 5,040 | .039683 | .015873 | .5476 |
+| 80 | 10,080 | .039683 | .031746 | .5952 |
+
+- MRT collision. Mach fixed (Ma ≈ .0687); a separate Mach check at D = 40
+  with T = 10,080 (U_max halved, tau .5238) is labelled MACH and not part of the
+  ladder.
+- Per-step displacement ≤ U_max ≈ .04 cells (bound .5).
+
+## 3. Domain and boundaries (finite-domain reproduction), ENGINEERING
+
+The paper's 120 D O-grid is not reproduced; our rectangular box is a
+reproduction choice. The native dense channel (zero-velocity inlet,
+pressure-NEQ outlet, bounce-back y walls) is a finite-domain approximation,
+not infinite still fluid. The G4 x-periodic / y-slip override is a separately
+labelled sensitivity family and never mixed into the primary ladder. The
+runner owner confirms the zero-velocity inlet is supported as configured
+before any run.
+
+Order (closed):
+1. Smoke + ledger at D = 40, native BC, 40 D × 40 D (finite box).
+2. Domain sensitivity before any literature statement: same BC, D = 40,
+   L ∈ {20 D, 40 D, 60 D}; record the change in Fx waveform, c_d, c_i and the
+   scored profiles.
+3. Optional labelled BC sensitivity: one x-periodic / y-slip case at D = 40,
+   40 D; reported as a delta against native only.
+4. Only after step 2 passes: the D = 20 / 40 / 80 ladder on the chosen BC and
+   domain (default 40 D if, against 60 D, `RMS(dFx)/F_ref` and |dc_d|, |dc_i|
+   stay inside the section 6 bands).
+- Memory/runtime: measured by a smoke run before each new size; measurements
+  inform scheduling only, never tolerances.
+
+## 4. Observables (raw, per step, preserved)
+
+- `obs` Fx, Fy (CME) and all `MovingTables.DIAG` columns per step: CME,
+  CME + Wen and CME + cover − uncover reported separately; crossings, eq_only,
+  no_donor (no_donor must stay 0).
+- Velocity u1, u2 on lab-fixed lines x/D = -.6, 0, .6, 1.2 (relative to the
+  mean cylinder position) versus x2, at the scored phases. Lab frame
+  throughout; if any body-frame intermediate is used, position and velocity
+  are both transformed.
+- Sampling cost: raw force and diagnostics counters are read every step
+  (small). Flags and velocity fields are read only at the steps inside the
+  three prescribed 1° phase bins (φ ∈ [180, 181), [210, 211), [330, 331),
+  T/360 steps each) of the cycles needed for settling and scoring, plus
+  sparse declared audits; never every step.
+- Point set: fixed x2/D sample points per station, identical on every grid
+  and declared before running; per point the interpolation weights are
+  recorded; the scored set is the common subset valid on all compared grids
+  at that phase, with its coverage recorded; point sets never change silently.
+- Interpolation: bilinear on cell centres using only cells that are fluid at
+  that step (flags read back via `get_flags()`); weights renormalized over the
+  fluid corners; a sample with no fluid corner, or inside the circle at that
+  step, is MISSING and logged, never set to zero. One-degree bins keep valid
+  fluid weights. The numerical average need not mimic the 200 LDA samples.
+
+## 5. Periodic-regime gate (every case, before any statistic), ENGINEERING
+
+Ledger rows, same schema for every case: dx, dt (in D and T units); U_max, nu,
+Re, KC (lattice and nondimensional); T (native steps, nondimensional);
+duration (native steps, s/T, paper cycles); discard; accepted window; cycle
+ledger; verdict PASS or FAIL-CLOSED.
+
+Scales: `F_ref = 0.5 rho U_max² D`, `U_ref = U_max`; never a near-zero mean
+force as denominator.
+
+- Settling metrics, cycle k−1 → k, phase-aligned by step within the cycle:
+  - waveform `RMS(dFx) / F_ref ≤ .02`;
+  - `|d mean|`, `|d RMS|`, `|d peak|` of Fx, each `/ F_ref ≤ .01`;
+  - scored profiles (fluid-only) `|du| / U_max`, `|dv| / U_max ≤ .03` at the
+    three phases × four stations;
+  - Morison `|dc_d| ≤ .02`, `|dc_i| ≤ .02`.
+- Rule: all metrics for N = 3 consecutive whole paper cycles → settled.
+
+Cycle indexing (worked example, binding). Paper cycle k spans
+`s/T in [0.25 + (k-1), 0.25 + k)`, i.e. φ0 → φ0; cycle 0 is the startup
+quarter `s/T < 0.25`. Pair p_k compares cycle k−1 with cycle k.
+
+- Discarded: cycles 0, 1, 2 (s/T < 2.25). Pairs are formed only between
+  cycles ≥ 3, so the first pair is p_4 (3 → 4).
+- Settling proof for a case: three consecutive pairs all pass; K_case = the
+  first cycle AFTER the last of those pairs. Earliest possible: pairs 3→4,
+  4→5, 5→6 pass ⇒ K_case = 7.
+- Common window: K = max over compared cases of K_case; every case scores the
+  same three complete cycles K, K+1, K+2 and must stay settled throughout:
+  pairs (K−1→K), (K→K+1), (K+1→K+2) must pass on every case, otherwise that
+  case is FAIL-CLOSED (an earlier local settle does not waive this).
+- Cap: cycle k ends at s/T = k + 0.25, so the last admissible scored cycle is
+  11 (ends 11.25 T) and K ≤ 9. Examples: K = 7 → cycles 7–9 (ends 9.25 T);
+  K = 9 → cycles 9–11 (ends 11.25 T); any K_case > 9 → FAIL-CLOSED.
+- Every case is run to the end of cycle K + 2 (proportionally more native
+  steps on finer grids).
+- Discard: the initial quarter cycle plus the first 2 full paper cycles
+  (s/T < 2.25). Settling is searched from cycle 3; cap 12 periods from s = 0;
+  unsettled → FAIL-CLOSED (no literature, cross-resolution or VIV statistic).
+- Common window: K = the largest first-settled cycle index over all compared
+  cases; every case is scored on the same 3 cycles K .. K + 2 and must reach
+  them within the cap (proportionally more native steps on finer grids), or
+  the comparison is FAIL-CLOSED.
+- The earlier G4 timing check (100 coarse vs 400 diffusive-fine steps) was a
+  matched short transient, not periodic evidence.
+
+## 6. Metrics and per-observable criteria (ENGINEERING, not paper uncertainty)
+
+Three criteria stay separate; the G4 same-grid 2% target is not reused.
+
+- Morison fit: least squares of the normalized per-step force over the
+  accepted window against the two basis functions, each averaged over its
+  single step interval, fixed paper phase, no free phase shift; c_d and
+  added-mass c_i per force candidate.
+- Cross-resolution: every scored observable's error against its reference must
+  decrease monotonically D20 → D40 → D80 (necessary, not sufficient).
+- Against Table 1 (computed fit, not measured force): at D ≥ 40,
+  `|c_d - 2.09| ≤ .15` and `|c_i - 1.45| ≤ .15`, and monotone.
+- Force waveform against its own Morison reconstruction: residual reported
+  only, no pass/fail (the paper reports fit error above its coarse numerical
+  error and gives no percentage).
+- Fig. 8 LDA profiles: QUALITATIVE ONLY. The original figure could not be
+  accessed or rendered; no digitized ordinates and no measurement uncertainty
+  exist. After a digitization with figure number, tool and uncertainty note:
+  `RMS(|du|, |dv|) / U_max ≤ .15` at D ≥ 40 and monotone; LDA error stays
+  UNQUANTIFIED and is stated with every comparison.
+- Fy sanity (not literature): `|mean Fy| / F_ref ≤ .05` in the accepted window.
+- No force candidate is selected for an RHS by these results; each is
+  reported.
+
+Superseded (earlier drafts): the velocity-jump start, guessed stations and
+phases, "C_M" for c_i, and the uncertainty wording; cycles 4–6, the 30-cycle
+cap and the earlier tolerance numbers.
+
+## 6a. Driver unit tests (CPU, before any run)
+
+- Phase/cycle arithmetic: s ↔ (k, φ), cycle boundaries, bin membership,
+  window/cap selection for K = 7, 8, 9, 10.
+- Solid-aware interpolation: all-fluid, partial-solid (renormalized weights),
+  no-fluid (MISSING) and inside-circle cases.
+- Interval-averaged Morison bases and the two-parameter fit on a synthetic
+  waveform with known c_d, c_i.
+
+## 7. Order and scope
+
+Protocol review → one-shot Dütsch driver (analytic `set_body` + `run(1)`,
+cycle ledger, lab-fixed fluid-only station sampler, two-basis Morison fit; no
+reusable runner until the first smoke shows shared helpers pay off; no
+production change) → smoke (memory/runtime) → domain sensitivity →
+ladder → physical-validation verdict → only then single-cylinder VIV. No AMR,
+no pinball, no manuscript work.

src/CelerisLab/body/motion.py

@@ -111,14 +111,21 @@ class MovingCircle:
             step()
             synchronize()
         except BaseException as exc:
-            self.failed = True
-            if self._on_failure is not None:
-                self._on_failure(exc)
+            self.fail(exc)
             raise
+        self.commit()
+
+    def commit(self) -> None:
+        """Applied := pending offset (caller has synchronized a successful step)."""
         if self.pending is not None:
             self.applied = self.pending[0]
             self.pending = None
 
+    def fail(self, exc: BaseException) -> None:
+        self.failed = True
+        if self._on_failure is not None:
+            self._on_failure(exc)
+
     def band(self) -> Tuple[int, int, int, int]:
         """Inclusive cell AABB (x0, x1, y0, y1) swept over the whole offset box."""
         x0, x1, y0, y1 = self.translation.offset_box
@@ -169,21 +176,29 @@ def enumerate_donor_coverage(radius: float, positions: int = 8, directions: int
 
 
 class MovingTables:
-    """Device tables of the translating circles for ``MovingCircleLinks`` (fixed centre).
+    """Device tables of the translating circles for ``MovingCircleCells``/``Links``.
 
-    Immutable per run: band table (x0, y0, width, height, first_slot, body_id) and
-    FP64 geometry (cx, cy, r) at the applied offset; plus the separate per-body
-    diagnostics (cme_x, cme_y, wen_x, wen_y, link_dm, n_links), never part of obs.
+    Band table (x0, y0, width, height, first_slot, body_id), immutable per run;
+    FP64 motion records (cx_old, cy_old, cx_new, cy_new, r, ux, uy), immutable
+    during a step; separate per-body diagnostics (see DIAG), never part of obs.
     """
-    BAND_INTS, DIAG_FLOATS = 6, 6
+    BAND_INTS, GEOM_DOUBLES = 6, 7
+    DIAG = ("cme_x", "cme_y", "wen_x", "wen_y", "link_dm", "n_links",
+            "cover_jx", "cover_jy", "cover_mass", "uncover_jx", "uncover_jy", "uncover_mass",
+            "eq_only", "no_donor", "crossings")
+    DIAG_FLOATS = len(DIAG)
 
     def __init__(self, band, geom):
         import numpy as np
         import pycuda.driver as cuda
-        self.n_bodies = len(geom) // 3
+        self.n_bodies = len(geom) // self.GEOM_DOUBLES
         self.n_slots = int(band[-2] + band[-4] * band[-3] * 8)
+        self.n_cells = self.n_slots // 8
+        self.geom = list(geom)  # host copy of the uploaded records
         self.band_gpu = cuda.to_device(np.asarray(band, np.int32))
         self.geom_gpu = cuda.to_device(np.asarray(geom, np.float64))
+        # Owned page-locked staging for async record uploads (allocated once).
+        self._geom_pinned = cuda.pagelocked_empty(len(geom), np.float64)
         self.diag_gpu = cuda.to_device(np.zeros(self.n_bodies * self.DIAG_FLOATS, np.float32))
         self.slot_debug_gpu = None  # optional (q, fallback, rx, ry) per slot for parity tests
 
@@ -201,9 +216,32 @@ class MovingTables:
             width, height = x1 - x0 + 1, y1 - y0 + 1
             band += [x0, y0, width, height, first, int(body_id)]
             first += width * height * 8
-            geom += [*motion.centre(motion.applied), motion.radius]
+            geom += MovingTables.record(motion)
         return band, geom
 
+    @staticmethod
+    def record(motion):
+        """FP64 motion record of one circle for its next step (hold if nothing pending)."""
+        c_old, c_new, u = motion.record()
+        return [*c_old, *c_new, motion.radius, *u]
+
+    def upload_records(self, circles, stream) -> bool:
+        """Enqueue the records when they differ from the device copy; True if advancing.
+
+        The pinned staging buffer is reused: callers synchronize ``stream`` after every
+        upload before the next one (``Simulation._run_moving``), so no copy is in flight
+        when it is rewritten. ``geom`` is the enqueued state; after a failed copy the
+        owner is failed and the cache is never trusted again.
+        """
+        import pycuda.driver as cuda
+        geom = [v for motion, _ in circles for v in self.record(motion)]
+        if geom != self.geom:
+            self._geom_pinned[:] = geom
+            cuda.memcpy_htod_async(self.geom_gpu, self._geom_pinned, stream)
+            self.geom = geom
+        g = self.GEOM_DOUBLES
+        return any(geom[i:i + 2] != geom[i + 2:i + 4] for i in range(0, len(geom), g))
+
     def enable_slot_debug(self):
         import numpy as np
         import pycuda.driver as cuda

src/CelerisLab/lbm/kernels/step/moving_circle.cu

@@ -1,6 +1,7 @@
-// Translating-circle cut links, fixed centre (S2a; docs/moving_bodies.md, M2).
-// One thread per (swept-band cell, direction) slot. Geometry comes from an FP64
-// {cx, cy, r} table per body; the operation order mirrors
+// Translating circles (docs/moving_bodies.md): M1 MovingCircleCells reclassifies
+// band cells and refills uncovered ones; M2 MovingCircleLinks applies the cut links.
+// Geometry comes from an immutable FP64 record per body
+// {cx_old, cy_old, cx_new, cy_new, r, ux, uy}; the M2 operation order mirrors
 // CircleGeometry.build_cut_links / find_circle_intersection with explicit _rn
 // operations (bit parity with the host builder is a test target, not a claim).
 // Moving circles are never in the packed static links, so this kernel is the
@@ -10,7 +11,10 @@
 
 #if DIM == 2
 constexpr int MOVING_BAND_INTS = 6;   // x0, y0, width, height, first_slot, body_id
-constexpr int MOVING_DIAG_FLOATS = 6; // cme_x, cme_y, wen_x, wen_y, link_dm, n_links
+constexpr int MOVING_GEOM_DOUBLES = 7;
+// cme_x, cme_y, wen_x, wen_y, link_dm, n_links, cover_jx, cover_jy, cover_mass,
+// uncover_jx, uncover_jy, uncover_mass, eq_only, no_donor, crossings
+constexpr int MOVING_DIAG_FLOATS = 15;
 
 __device__ __forceinline__ double mc_sq(double v) { return __dmul_rn(v, v); }
 
@@ -53,7 +57,8 @@ __global__ void MovingCircleLinks(
     const unsigned int dir = local % 8u + 1u;
     const int x = b[0] + (int)((local / 8u) % (unsigned int)b[2]);
     const int y = b[1] + (int)((local / 8u) / (unsigned int)b[2]);
-    const double cx = geom[3 * m], cy = geom[3 * m + 1], r = geom[3 * m + 2];
+    const double* g = geom + MOVING_GEOM_DOUBLES * m;
+    const double cx = g[2], cy = g[3], r = g[4];
     const int ex = d_cx[dir], ey = d_cy[dir];
     float* debug = slot_debug == nullptr ? nullptr : slot_debug + 4 * tid;
     if (debug != nullptr) debug[0] = -1.0f;
@@ -85,7 +90,7 @@ __global__ void MovingCircleLinks(
 
     const int bid = b[5];
     const float omega = action_omega(action, bid);
-    const float Uw = -omega * ry, Vw = omega * rx;  // fixed centre: no translation term
+    const float Uw = (float)g[5] - omega * ry, Vw = (float)g[6] + omega * rx;  // u + omega x r
     float exchange[3], dm;
     apply_bouzidi_link(dir, q, linear_index((unsigned int)x, (unsigned int)y), fi,
                        Uw, Vw, rx, ry, fallback, obs, bid, exchange, &dm);
@@ -99,6 +104,80 @@ __global__ void MovingCircleLinks(
         atomicAdd(&d[5], 1.0f);
     }
 }
+// M1: one thread per band cell, launched only on steps with c_old != c_new.
+// Membership comes from the immutable old/new predicates, never from flags, so
+// classification, refill and the flag write of one cell cannot race with another.
+// ddf holds post-collision G^{n-1}; refilled cells get all nine slots at that level.
+__global__ void MovingCircleCells(
+    fpxx* fi, uint16_t* flag, const int* band, const double* geom,
+    unsigned int n_bodies, unsigned int n_cells, const float* action, float* diag)
+{
+    const unsigned int tid = threadIdx.x + blockIdx.x * blockDim.x;
+    if (tid >= n_cells) return;
+    unsigned int m = 0;
+    while (m + 1 < n_bodies && tid >= (unsigned int)band[(m + 1) * MOVING_BAND_INTS + 4] / 8u) ++m;
+    const int* b = band + m * MOVING_BAND_INTS;
+    const unsigned int local = tid - (unsigned int)b[4] / 8u;
+    const int x = b[0] + (int)(local % (unsigned int)b[2]);
+    const int y = b[1] + (int)(local / (unsigned int)b[2]);
+    const double* g = geom + MOVING_GEOM_DOUBLES * m;
+    const bool in_old = mc_inside(x, y, g[0], g[1], g[4]);
+    const bool in_new = mc_inside(x, y, g[2], g[3], g[4]);
+    if (in_old == in_new) return;
+    const unsigned long k = linear_index((unsigned int)x, (unsigned int)y);
+    float* d = diag + m * MOVING_DIAG_FLOATS;
+    float rho = 0.0f, jx = 0.0f, jy = 0.0f;
+    if (in_new) {  // covered: its fluid content leaves the active set
+        for (unsigned int i = 0; i < NQ; ++i) {
+            const float f = load_ddf(fi, index_f(k, i));
+            rho += f; jx += (float)d_cx[i] * f; jy += (float)d_cy[i] * f;
+        }
+        flag[k] = FLAG_SOLID | FLAG_OBSTACLE | FLAG_BC_CURVED;
+        atomicAdd(&d[6], jx); atomicAdd(&d[7], jy); atomicAdd(&d[8], rho);
+        atomicAdd(&d[14], 1.0f);
+        return;
+    }
+    // Uncovered: donors are fluid in BOTH geometries; NEQ from the neighbour with the
+    // largest unnormalized c.n > 0 (first direction wins ties), n from c_new.
+    const float nx = (float)__dsub_rn((double)x, g[2]), ny = (float)__dsub_rn((double)y, g[3]);
+    float rho_sum = 0.0f, best_dot = 0.0f;
+    int n_valid = 0, best = -1;
+    for (unsigned int i = 1; i < NQ; ++i) {
+        const int xn = x + d_cx[i], yn = y + d_cy[i];
+        if (mc_inside(xn, yn, g[0], g[1], g[4]) || mc_inside(xn, yn, g[2], g[3], g[4])) continue;
+        const unsigned long kn = linear_index((unsigned int)xn, (unsigned int)yn);
+        for (unsigned int q = 0; q < NQ; ++q) rho_sum += load_ddf(fi, index_f(kn, q));
+        ++n_valid;
+        const float dot = (float)d_cx[i] * nx + (float)d_cy[i] * ny;
+        if (dot > best_dot) { best_dot = dot; best = (int)i; }
+    }
+    const float omega = action_omega(action, b[5]);
+    const float uwx = (float)g[5] - omega * ny, uwy = (float)g[6] + omega * nx;  // at the cell
+    float feq[NQ], neq[NQ];
+    compute_feq(n_valid > 0 ? rho_sum / (float)n_valid : (float)RHO, uwx, uwy, feq);
+    for (unsigned int i = 0; i < NQ; ++i) neq[i] = 0.0f;
+    if (best >= 0) {
+        const unsigned long kd = linear_index((unsigned int)(x + d_cx[best]), (unsigned int)(y + d_cy[best]));
+        float gd[NQ], rd = 0.0f, mx = 0.0f, my = 0.0f, fd[NQ];
+        for (unsigned int i = 0; i < NQ; ++i) {
+            gd[i] = load_ddf(fi, index_f(kd, i));
+            rd += gd[i]; mx += (float)d_cx[i] * gd[i]; my += (float)d_cy[i] * gd[i];
+        }
+        compute_feq(rd, mx / rd, my / rd, fd);
+        for (unsigned int i = 0; i < NQ; ++i) neq[i] = gd[i] - fd[i];
+    } else {
+        atomicAdd(&d[12], 1.0f);
+    }
+    if (n_valid == 0) atomicAdd(&d[13], 1.0f);  // hard failure, checked before commit
+    for (unsigned int i = 0; i < NQ; ++i) {
+        const float f = feq[i] + neq[i];
+        store_ddf(fi, index_f(k, i), f);
+        rho += f; jx += (float)d_cx[i] * f; jy += (float)d_cy[i] * f;
+    }
+    flag[k] = FLAG_FLUID;
+    atomicAdd(&d[9], jx); atomicAdd(&d[10], jy); atomicAdd(&d[11], rho);
+    atomicAdd(&d[14], 1.0f);
+}
 #endif  // DIM == 2
 
 #endif  // CELERIS_MOVING_CIRCLE_CU

src/CelerisLab/lbm/stepper.py

@@ -40,6 +40,7 @@ class LBMStepper:
         self.sensor_fn = module.get_function("SensorKernel")
         self.force_region_fn = module.get_function("ForceRegionKernel")
         self.moving_fn = None if self._esopull else module.get_function("MovingCircleLinks")
+        self.moving_cells_fn = None if self._esopull else module.get_function("MovingCircleCells")
 
         tpb = cfg.threads_per_block
         self.block = (tpb, 1, 1)
@@ -77,6 +78,11 @@ class LBMStepper:
         """Advance *n* time steps."""
         if n > 0 and self.blocked is not None:
             raise RuntimeError(self.blocked)
+        self._step(n, action_gpu=action_gpu, obs_gpu=obs_gpu, stream=stream)
+
+    def _step(self, n, *, action_gpu, obs_gpu, stream=None, advance=False):
+        """Unguarded steps; ``advance`` runs one step that moves translating circles."""
+        assert not advance or n == 1
         f = self.field
         launch_kw = {}
         if stream is not None:
@@ -89,7 +95,7 @@ class LBMStepper:
                     block=self.block, grid=self.boundary_grid, **launch_kw,
                 )
             if self.moving is not None:
-                self._launch_moving(action_gpu, obs_gpu, **launch_kw)
+                self._launch_moving(action_gpu, obs_gpu, advance, **launch_kw)
             self._launch_curved(action_gpu, obs_gpu, **launch_kw)
 
             if self._esopull:
@@ -115,8 +121,14 @@ class LBMStepper:
             # Any completed GPU step invalidates host DDF mirror.
             f.invalidate_host_ddf_cache()
 
-    def _launch_moving(self, action_gpu, obs_gpu, **launch_kw):
+    def _launch_moving(self, action_gpu, obs_gpu, advance, **launch_kw):
         m, tpb = self.moving, self.cfg.threads_per_block
+        if advance:  # M1 completes before M2 on the same stream
+            self.moving_cells_fn(
+                self.field.ddf_gpu, self.field.flag_gpu, m.band_gpu, m.geom_gpu,
+                np.uint32(m.n_bodies), np.uint32(m.n_cells), action_gpu, m.diag_gpu,
+                block=(tpb, 1, 1), grid=((m.n_cells + tpb - 1) // tpb, 1, 1), **launch_kw,
+            )
         self.moving_fn(
             self.field.ddf_gpu, self.field.flag_gpu, m.band_gpu, m.geom_gpu,
             np.uint32(m.n_bodies), np.uint32(m.n_slots), action_gpu, obs_gpu,

src/CelerisLab/simulation.py

@@ -29,6 +29,7 @@ from .cuda import compiler
 from .lbm.descriptors import FLUID
 
 # Translating circles step only through their moving tables (fixed centre in S2a).
+PENDING_BLOCKED = "a pending translation target is advanced only by Simulation.run(1)"
 MOVING_BLOCKED = ("translating circles need the moving kernels and their tables (attached "
                   "at initialize); "
                   "positive-step runs are rejected")
@@ -212,12 +213,19 @@ class Simulation:
     def _fail_moving(self, exc: BaseException) -> None:
         # Device geometry is unknown after a failed moving step: same lifecycle as Jets.
         self._failed = f"moving circle step failed: {exc!r}"
+        self.stepper.blocked = self._failed  # direct stepper calls are refused as well
 
     def _attach_moving(self) -> None:
         # Direct stepper calls stay blocked until the moving tables exist.
         tables = getattr(self, "_moving_tables", None)
         self.stepper.moving = tables
-        self.stepper.blocked = MOVING_BLOCKED if self._moving_bodies() and tables is None else None
+        moving = self._moving_bodies()
+        if moving and tables is None:
+            self.stepper.blocked = MOVING_BLOCKED
+        elif any(obj.motion.pending is not None for obj in moving):
+            self.stepper.blocked = PENDING_BLOCKED  # direct steps never consume a target
+        else:
+            self.stepper.blocked = None
 
     def _build_moving_tables(self):
         moving = self._moving_bodies()
@@ -228,8 +236,8 @@ class Simulation:
 
     @context_bound
     def moving_diagnostics(self) -> np.ndarray:
-        """Per translating circle (in body order): cme_x, cme_y, wen_x, wen_y, link_dm,
-        n_links accumulated since obs were last zeroed (diagnostics only)."""
+        """Per translating circle (in body order), columns ``MovingTables.DIAG``,
+        accumulated since obs were last zeroed (diagnostics only, not a force RHS)."""
         self._require_usable()
         tables = self.stepper.moving
         if tables is None:
@@ -354,17 +362,19 @@ class Simulation:
             omega: rotation speed (lattice units).  ``None`` = no change.
             vx: translation velocity x (future).  ``None`` = no change.
             vy: translation velocity y (future).  ``None`` = no change.
-            offset: target centre offset of a translating circle. Displacement
-                is disabled until the moving cell update exists (S3): any offset
-                raises before mutation; translating circles hold at the applied
-                offset and accept omega only.
+            offset: target centre offset of a translating circle; with vx/vy
+                (default 0) it replaces the pending target after validation
+                against the applied offset (docs/moving_bodies.md, 4a). Only
+                ``run(1)`` advances it; direct stepper calls are blocked meanwhile.
         """
         obj = self.bodies.get(id)
         motion = obj.motion
         if motion is not None:
-            if offset is not None or vx is not None or vy is not None:
-                raise ValueError("translating-circle displacement is disabled until S3; "
-                                 "only omega may be set")
+            if offset is None and (vx is not None or vy is not None):
+                raise ValueError("a translating circle takes vx/vy only with an offset")
+            target = None if offset is None else motion.check_target(
+                offset, (0.0 if vx is None else vx, 0.0 if vy is None else vy))
+            vx = vy = None  # the motion target, not ObjectState, owns translation
         elif offset is not None:
             raise ValueError(f"body {id} was not declared with a translation")
         proposed = dict(vars(obj.state))
@@ -382,6 +392,9 @@ class Simulation:
         if vy is not None:
             obj.state.vy = float(vy)
             changed = True
+        if motion is not None and target is not None:
+            motion.pending = target
+            self._attach_moving()
         if changed:
             self.bodies._refresh_action_from_objects()
 
@@ -700,8 +713,13 @@ class Simulation:
         checkpoint_interval = operator.index(checkpoint_interval)
         if steps < 0 or checkpoint_interval < 0:
             raise ValueError("steps and checkpoint_interval must be nonnegative")
-        if steps > 0 and self._moving_bodies() and self.stepper.moving is None:
+        moving = self._moving_bodies()
+        if steps > 0 and moving and self.stepper.moving is None:
             raise RuntimeError(MOVING_BLOCKED)
+        if moving and checkpoint_interval > 0:
+            raise ValueError("checkpoints are not supported with translating circles")
+        if steps > 1 and any(obj.motion.pending is not None for obj in moving):
+            raise ValueError("a pending translation target needs run(1)")
         self._assert_runtime_contracts()
         # Discard any uncommitted body edits before stepping.
         if self.bodies.has_pending_edit():
@@ -736,6 +754,9 @@ class Simulation:
                 if done < steps or done % checkpoint_interval == 0:
                     stream.synchronize()
                     self.save_checkpoint()
+        elif self.stepper.moving is not None and steps > 0:
+            self._run_moving(steps, stream, sync_obs)
+            return
         else:
             self.stepper.step(
                 steps,
@@ -754,6 +775,34 @@ class Simulation:
 
         stream.synchronize()
 
+    def _run_moving(self, steps, stream, sync_obs):
+        """Steps with translating circles: hold, or one advancing step that commits
+        each pending target only after synchronization and a zero ``no_donor`` count
+        (docs/moving_bodies.md 4a, 5). Any failure leaves the Simulation unusable."""
+        tables = self.stepper.moving
+        circles = [(obj.motion, obj.obj_id) for obj in self._moving_bodies()]
+        try:
+            advancing = tables.upload_records(circles, stream)
+            self.stepper._step(steps, action_gpu=self.bodies.action_gpu,
+                               obs_gpu=self.bodies.obs_gpu, stream=stream, advance=advancing)
+            self.bodies.increment_obs_steps(steps)
+            if sync_obs:
+                self.bodies.download_obs_full_async(stream)
+            stream.synchronize()
+            if advancing:
+                no_donor = tables.download_diagnostics()[:, tables.DIAG.index("no_donor")]
+                if np.any(no_donor > 0):
+                    raise RuntimeError("translating circle refill found no donor cell")
+                for motion, _ in circles:
+                    motion.commit()
+                tables.upload_records(circles, stream)  # first hold: old = new = applied, u = 0
+                stream.synchronize()
+                self._attach_moving()
+        except BaseException as exc:
+            for motion, _ in circles:  # device geometry/state unknown after any failure
+                motion.fail(exc)
+            raise
+
     def step(self, n: int = 1):
         """Advance *n* steps (convenience for interactive use)."""
         self.run(n)
@@ -771,8 +820,16 @@ class Simulation:
         self.field.download_ddf(step_id=self.stepper.step_count)
         return self.field.ddf.copy()
 
+    @context_bound
     def get_flags(self) -> np.ndarray:
-        return self.field.flag.copy()
+        """Flag array; read back from the device when translating circles may have moved."""
+        if getattr(self.stepper, "moving", None) is None:
+            return self.field.flag.copy()
+        self._require_usable()
+        self.stream.synchronize()
+        flags = np.empty_like(self.field.flag)
+        cuda.memcpy_dtoh(flags, self.field.flag_gpu)
+        return flags
 
     # -- Runtime parameter updates -------------------------------------------
     @context_bound

tests/unit/test_moving_bodies.py

@@ -217,20 +217,21 @@ def test_add_body_and_set_body_route_targets(monkeypatch):
     bid = sim.add_body("circle", (50., 40.), 10., translation=Translation((-2., 2., -3., 3.), .1))
     still = sim.add_body("circle", (100., 40.), 5.)
     motion = sim.bodies.get(bid).motion
-    assert sim.stepper.blocked is not None  # declaring a moving circle blocks direct steps
-    for initialized in (False, True):  # S2a: displacement disabled before and after initialize
-        sim._initialized = initialized
-        for kwargs in (dict(offset=(0., .3), vy=.05), dict(offset=(0., 0.)), dict(vy=.05),
-                       dict(vx=.05), dict(offset=(0., .1), omega=.01)):
-            with pytest.raises(ValueError, match="disabled until S3"):
-                sim.set_body(bid, **kwargs)
-            assert motion.pending is None and motion.applied == (0., 0.)
-    sim._initialized = False
+    sim._moving_tables = "tables"
+    sim._attach_moving()
+    assert sim.stepper.blocked is None
+    sim.set_body(bid, offset=(0., .3), vy=.05)
+    assert sim.stepper.blocked == module.PENDING_BLOCKED  # direct steps never consume a target
+    sim.set_body(bid, offset=(0., .4), vy=.05)  # last wins
+    for kwargs in (dict(offset=(0., .6)), dict(vy=.05), dict(vx=.05),
+                   dict(offset=(0., .1), vy=.2), dict(offset=(0., 0.), vy=.01)):
+        with pytest.raises(ValueError):
+            sim.set_body(bid, **kwargs)
+        assert motion.pending == ((0., .4), (0., .05))  # invalid input keeps the target
+    sim.set_body(bid, omega=.01)  # omega-only keeps the pending target
     obj = sim.bodies.get(bid)
-    assert obj.state.omega == 0.
-    sim.set_body(bid, omega=.01)  # spin on hold stays allowed
     assert obj.state.omega == .01 and (obj.state.vx, obj.state.vy) == (0., 0.)
-    assert motion.pending is None
+    assert motion.pending == ((0., .4), (0., .05)) and motion.applied == (0., 0.)
     with pytest.raises(ValueError, match="not declared"):
         sim.set_body(still, offset=(0., .1))
     sim.set_body(still, vx=.01)  # stationary API unchanged
@@ -243,6 +244,108 @@ def test_add_body_and_set_body_route_targets(monkeypatch):
                          **kwargs)
 
 
+class Tables:
+    """MovingTables double: records uploads, steps and diagnostics reads."""
+    DIAG = ("no_donor",)
+
+    def __init__(self, log, no_donor=0., fail=None):
+        self.log, self.no_donor, self.fail = log, no_donor, fail
+
+    def upload_records(self, circles, stream):
+        records = [m.record() for m, _ in circles]
+        self.log.append(("upload", records))
+        uploads = sum(1 for w in self.log if isinstance(w, tuple) and w[0] == "upload")
+        if self.fail == "upload" or (self.fail == "hold_upload" and uploads == 2):
+            raise RuntimeError("upload failed")
+        return any(r[0] != r[1] for r in records)
+
+    def zero_diagnostics(self, stream):
+        self.log.append("zero_diag")
+
+    def download_diagnostics(self):
+        return np.array([[self.no_donor]], np.float32)
+
+
+def advancing_simulation(monkeypatch, **tables):
+    sim, writes = moving_simulation(monkeypatch)
+    sim.bodies.objects[0].motion = circle()
+    sim.stepper = SimpleNamespace(
+        moving=Tables(writes, **tables), blocked=None,
+        _step=lambda n, advance=False, **k: writes.append(("step", n, advance)))
+    from CelerisLab.lbm.stepper import LBMStepper
+    sim.stepper.step = LBMStepper.step.__get__(sim.stepper)  # real public entry
+    sim._moving_tables = sim.stepper.moving
+    return sim, writes, sim.bodies.objects[0].motion
+
+
+def test_run_one_advances_and_commits_after_sync(monkeypatch):
+    sim, writes, motion = advancing_simulation(monkeypatch)
+    motion.set_target((0., .25), (0., .05))
+    with pytest.raises(ValueError, match="run\\(1\\)"):  # multi-step with a target
+        sim.run(2, stream=STREAM)
+    assert writes == [] and motion.pending is not None
+    sim.run(0, stream=STREAM)  # zero steps never consume the target
+    assert motion.pending is not None
+    sim.run(1, stream=STREAM)
+    uploads = [w for w in writes if isinstance(w, tuple) and w[0] == "upload"]
+    assert uploads[0][1] == [((50., 40.), (50., 40.25), (0., .05))]
+    assert uploads[1][1] == [((50., 40.25), (50., 40.25), (0., 0.))]  # first hold after move
+    assert ("step", 1, True) in writes
+    assert motion.applied == (0., .25) and motion.pending is None
+    writes.clear()
+    sim.run(3, stream=STREAM)  # hold steps: no advance
+    assert ("step", 3, False) in writes and motion.applied == (0., .25)
+
+
+@pytest.mark.parametrize("failure", ("upload", "first_sync", "no_donor", "hold_upload",
+                                     "hold_sync"))
+def test_failed_advance_fails_simulation(monkeypatch, failure):
+    kwargs = dict(no_donor=1.) if failure == "no_donor" else dict(fail=failure)
+    sim, writes, motion = advancing_simulation(monkeypatch, **kwargs)
+    motion._on_failure = sim._fail_moving
+    motion.set_target((0., .25), (0., .05))
+    syncs = []
+
+    def synchronize():  # stands in for an asynchronous failure surfacing at sync
+        syncs.append(1)
+        if (failure, len(syncs)) in (("first_sync", 1), ("hold_sync", 2)):
+            raise RuntimeError("sync failed")
+    with pytest.raises(RuntimeError):
+        sim.run(1, stream=SimpleNamespace(synchronize=synchronize))
+    # Before the step-group sync and no_donor check nothing is committed; the hold
+    # record is staged after the commit (applied already matches the moved device).
+    committed = failure in ("hold_upload", "hold_sync")
+    assert motion.applied == ((0., .25) if committed else (0., 0.))
+    assert (motion.pending is None) == committed
+    assert motion.failed and sim._failed is not None
+    assert sim.stepper.blocked == sim._failed  # direct stepper reuse refused too
+    with pytest.raises(RuntimeError, match="Simulation failed"):
+        sim.run(1, stream=STREAM)
+
+
+def test_failed_hold_step_blocks_the_direct_stepper(monkeypatch):
+    # No target is pending, so only the failure latch can refuse direct steps.
+    from CelerisLab.lbm.stepper import LBMStepper
+    sim, writes, motion = advancing_simulation(monkeypatch, fail="upload")
+    motion._on_failure = sim._fail_moving
+    launches = []
+    kernels = SimpleNamespace(get_function=lambda name: lambda *a, **k: launches.append(name))
+    cfg = SimpleNamespace(streaming="double_buffer", threads_per_block=128, is_d2q9=True,
+                          nx=64, ny=32, nz=1)
+    field = SimpleNamespace(ddf_gpu="ddf", temp_gpu="temp", flag_gpu="flag", n_curved=0)
+    stepper = LBMStepper(field, kernels, cfg)
+    stepper.moving = sim.stepper.moving
+    sim.stepper = stepper
+    assert motion.pending is None and stepper.blocked is None
+    with pytest.raises(RuntimeError, match="upload failed"):
+        sim.run(2, stream=STREAM)  # hold step
+    assert sim._failed is not None and stepper.blocked == sim._failed
+    with pytest.raises(RuntimeError, match="moving circle step failed"):
+        stepper.step(1, action_gpu=None, obs_gpu=None)
+    assert launches == [] and stepper.step_count == 0
+    assert (field.ddf_gpu, field.temp_gpu) == ("ddf", "temp")
+
+
 def test_moving_tables_layout_and_attachment(monkeypatch):
     from CelerisLab.body.motion import MovingTables
     a = MovingCircle((50., 40.), 10., Translation((-1., 1., -3., 3.), .1))
@@ -250,7 +353,7 @@ def test_moving_tables_layout_and_attachment(monkeypatch):
     band, geom = MovingTables.layout([(a, 2), (b, 5)])
     wa, ha = 63 - 37 + 1, 55 - 25 + 1
     assert band == [37, 25, wa, ha, 0, 2, 83, 23, 16, 16, wa * ha * 8, 5]
-    assert geom == [50., 40., 10., 90.25, 30.5, 5.]
+    assert geom == [50., 40., 50., 40., 10., 0., 0., 90.25, 30.5, 90.25, 30.5, 5., 0., 0.]  # hold
     sim = routed_simulation(monkeypatch)
     sim.add_body("circle", (50., 40.), 10., translation=Translation((0., 0., 0., 0.), .1))
     sim._moving_tables = "tables"