research: karman reproduce uses committed pinball ids after sync_bodies

add_body after initialize returns a pending negative id; sync_bodies commits
the new bodies in order, so the pinball ids are n0..n0+2 (f2f88c phase3 did
the same). Fixes the KeyError(-1) in the GPU smoke.

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

docs/steady_pinball_fiv.md

@@ -1,104 +0,0 @@
-# Uncontrolled steady pinball: development protocol
-
-Frozen before results: 1280 x 512 D2Q9/FP32/MRT, unshifted double buffer;
-rho=1, mean U=.01, nu=.004, D=20, Re_D=50. Centres are (600,255.5),
-(626,270.5), (626,240.5), radius10. No upstream cylinder. Passive probes stay
-at (800,295.5), (800,255.5), (800,215.5), radius5. Every omega is zero.
-Parabolic inlet, regularized NEQ damping1; y bounce-back planes .5/510.5;
-NEQ-extrapolated outlet with backflow clamp (the configured blend is inactive in
-this MRT route). Initial fluid is
-parabolic equilibrium. Fixed warmup is 4 NX/U=512000 steps.
-
-Geometry and fluid parameters follow DynamisLab 609116ec26f19766f5e9ce2b0df427d4e24faae8,
-`legacy_test/core/legacy_env_builder.py:build_steady_cloak` and the retained
-`config_lbm_pinball_legacy_re50d_u001.json`. This is a modern Celeris reconstruction,
-not a bitwise Legacy reproduction: Legacy uses shifted incompressible MRT, a
-different curved interpolation and prescribed parabolic outlet velocity. Current
-Celeris uses standard MRT, two-branch BFL and extrapolated outlet velocity.
-The fixed-hold control uses the SAME modern settings as the free case.
-
-## Shared translation and ownership
-
-`Translation(..., shared=True)` is explicitly transverse only. Reusing one
-envelope INSTANCE declares common motion; equal but separate envelopes are not
-linked. Before initialization and every positive Simulation.run, group members
-must have identical applied offsets, effective target offsets and interval wall
-velocities. Stage each member with set_body, then call run(1) once. An incomplete
-or inconsistent target group is rejected before device uploads and can be fixed.
-Existing pending-target direct-stepper restrictions remain.
-
-Only members of that group may have overlapping whole swept bands. Their fixed
-surface clearance must exceed4 lattice cells. This conservative bound exceeds
-2 sqrt(2)+MAX_STEP = 3.329: M1 changed cells lie in old/new disks and read one-link
-neighbors; M2 reads at most two links outside its new disk. Old/new relative
-displacement is at most .5. Thus reads cannot reach another member's old/new
-solid or M1 write region. M1 writes are disjoint. M2 writes opposite population
-slots inside its own new disk, also disjoint. Separate first_slot ranges still
-allocate each body's logical band/debug storage. Shared band cells alone do not
-write anything: M1 returns unless membership changes; M2 returns unless a cut
-link exists. Unrelated moving/static body and probe overlaps still reject.
-
-This is a sufficient restricted-group condition, not support for near-contact
-or independent motion. The r4, distance8.3, dy=.5 counterexample remains rejected.
-CPU tests enumerate candidate donors/write footprints over subcell phases;
-GPU tests compare fixed vs hold and member-order permutations through crossings
-at wall speed .01 over the full +/-D envelope (8000 steps). The initial .5-speed
-fixture failed with NaNs; a single-circle control also failed at step69. That
-Mach~.866 fixture was unsuitable as a low-Mach numerical gate. Its failure is
-retained; the geometric MAX_STEP proof is unchanged, not a stability guarantee.
-CUDA kernels and arithmetic are unchanged.
-
-## New structural parameters (not inherited from the paper)
-
-One transverse state y,v; total assembly M=10 rho (3 pi D^2/4)=9424.77796,
-fn=U/(5D)=.0001, Tn=10000, zeta=.01, K=M(2 pi fn)^2,
-C=2 zeta M(2 pi fn). Initial y=v=0, no seed or later kick. Force is the sum
-of the three instantaneous CME_y+Wen_y values from the previous fluid step;
-node impulses are logged separately and never added. The damped oscillator is
-integrated exactly for that piecewise-constant force over dt=1. Wall velocity
-is y[n+1]-y[n], not the endpoint oscillator velocity. No filtering or retuning.
-
-Observe600000 steps=60Tn=300 D/U after warmup, in each of free and fixed-hold
-cases. Final20Tn is a descriptive comparison window; ten-Tn blocks expose drift.
-No settled/periodic claim follows merely from reaching the cap. Stop on |y|>D,
-nonfinite values, no_donor, rejected displacement or speed>.05. Report any quiet
-response honestly. No automatic extension. --smoke uses32 warmup+32 steps only.
-
-Native-step rows retain y,v, each body's CME/Wen/node terms and refill counters;
-each600000x15 FP64 file is72MB. No full-field output or policy training. Compile
-and run on one device serially, keep original outputs and dirty trees untouched.
-This stage makes no claim of physical grid convergence or vibration suppression.
-
-Run from checkout with PYTHONPATH=src and the existing Python environment:
-`python research/steady_pinball/run.py --smoke --out <new-external-directory>`
-then separately `--out <new-free-directory>` and `--fixed --out <new-fixed-directory>`.
-
-## Zero-damping development follow-up (predeclared 2026-10-31)
-
-The original d0bcefa run and its zeta=.01 results remain historical evidence.
-The current runner exposes three fixed `--case` choices, all with C=0 as requested:
-
-| Case | M / (rho 3 pi D^2/4) | fn | M | K | Purpose |
-|---|---:|---:|---:|---:|---|
-| Z0 | 10 | .0001 | 9424.77796 | .00372075320 | Remove damping only |
-| M2 | 2 | .0001 | 1884.95559 | .000744150640 | Reduce mass and stiffness together at fixed fn |
-| R2 | 2 | .00006 | 1884.95559 | .000267894231 | Change stiffness only relative to M2 |
-
-R2 deliberately places fn at the retained fixed-body lift peak, f=.00006
-(St=.12 on mean U). This is an informed development choice, not independent
-physical validation or a promise of lock-in. All geometries, fluid settings,
-512000-step fixed warmup, y0=v0=0, force time level and unfiltered RHS are unchanged.
-No run result may select a replacement parameter or an automatic extension.
-
-Each case has a 600000-step cap, the same 300 D/U physical duration. This equals
-60 natural periods for Z0/M2 and 36 for R2. Compare the common final200000 native
-steps, and six100000-step blocks, retaining actual natural/response frequencies
-and ongoing drift. Do not call those blocks ten R2 periods. No settled-VIV claim
-follows from a larger displacement. Report all three outcomes, including quiet
-responses or stopped cases, alongside the existing fixed and damped-free records.
-
-Unchanged stops: |y|>D, invalid move/speed, nonfinite values or no_donor. A 1800s
-per-case wall cap and GPU headroom check bound resource use. Preserve partial rows;
-do not add damping, clip trajectories, filter forces, or relax any stop. Each
-case adds72MB of scalar history, no full fields. The separate Legacy/modern
-baseline diagnosis does not change this FSI protocol.

research/karman_reproduce/run.py

@@ -0,0 +1,116 @@
+"""Replay the frozen Legacy karman_re100 PPO policy on current CelerisLab.
+
+Research only; nothing in src/ is changed. Plant: Legacy Karman cloak scene
+(1280x512, D=20, parabolic inlet, mean U0 .01, nu .004, Re_D 50; upstream
+cylinder r=20 at x=200; pinball front (600, 255.5), rear (626, 255.5+-15);
+probes r=5 at x=800, y=255.5+40/0/-40). Interface contract = DynamisLab
+legacy_test (the executor behind the Legacy rerun), not f2f88c phase3:
+
+- command per interval: A = (8 a + b) U0, b = (0, -4, 4), a = deterministic
+  policy output in [-1, 1]^3; wall speed command, omega = +A / R (both
+  solvers impose u_w = omega x r, CCW positive);
+- smoothing per LATTICE step, A_eff <- .9 A_eff + .1 A, with A_eff = 0 at the
+  start of every 800-step interval (legacy_test.runtime.run_historical_interval);
+- observation: interval means of the three pinball forces (front, top, bottom)
+  and the three probe (ux, uy) sums, normalized by the frozen norm.json with
+  legacy_test.runtime.policy_observation; the first observation is zero.
+
+Output (--out, outside the checkout): signals.npz with sensors (N,6), forces
+(N,6), actions (N,3) in the f2f88c phase3 layout, so its comparator applies.
+
+    PYTHONPATH=src python research/karman_reproduce/run.py --dynamis DYN --model .../d1a3o12_re100.zip \\
+        --norm .../karman_re100/norm.json --config .../config_lbm_pinball_legacy_compat.json \\
+        --out /tmp/new-dir [--device 1]
+"""
+from __future__ import annotations
+
+import argparse
+import json
+import os
+import subprocess
+import sys
+
+import numpy as np
+
+L0, U0, NX, NY = 20.0, 0.01, 1280, 512
+R, YC = L0 / 2, (NY - 1) / 2
+SI, INTERVALS, WARMUP, FIFO_LEN = 800, 200, int(4 * NX / U0), 150
+BIAS = np.array([0.0, -4.0, 4.0])
+
+
+def git_head(path):
+    return subprocess.run(["git", "-C", path, "rev-parse", "HEAD"], capture_output=True,
+                          text=True).stdout.strip()
+
+
+def main():
+    ap = argparse.ArgumentParser()
+    ap.add_argument("--dynamis", required=True, help="DynamisLab checkout (imported read-only)")
+    ap.add_argument("--model", required=True, help="frozen d1a3o12_re100.zip")
+    ap.add_argument("--norm", required=True, help="frozen karman_re100 norm.json")
+    ap.add_argument("--config", required=True)
+    ap.add_argument("--out", required=True)
+    ap.add_argument("--device", type=int, default=0)
+    ap.add_argument("--intervals", type=int, default=INTERVALS)
+    args = ap.parse_args()
+    here = os.path.realpath(os.path.join(os.path.dirname(__file__), "..", ".."))
+    if os.path.realpath(args.out).startswith(here):
+        raise SystemExit("--out must be outside the checkout")
+    os.makedirs(args.out)  # fresh directory only; never overwrite results
+
+    sys.dont_write_bytecode = True  # import only: write nothing into the DynamisLab tree
+    sys.path.insert(0, os.path.join(args.dynamis, "src"))
+    from stable_baselines3 import PPO
+    from drl_pinball.legacy_test.core.model_loader import _dummy_env
+    from drl_pinball.legacy_test.runtime import load_policy_norm, policy_observation
+    model = PPO.load(args.model, env=_dummy_env(12), device="cpu")
+    norm = load_policy_norm(args.norm)
+
+    from CelerisLab.simulation import Simulation
+    sim = Simulation(args.config, device_id=args.device)
+    try:
+        sim.add_body("circle", center=(10 * L0, YC, 0.0), radius=L0)
+        probes = [sim.add_body("sensor", center=(40 * L0, YC + dy, 0.0), radius=5.0)
+                  for dy in (40.0, 0.0, -40.0)]
+        sim.initialize()
+        # Legacy build_karman_cloak order: stabilize disturbance + probes, then the
+        # 150-interval zero-action target sampling (kept as flow time), then the pinball.
+        sim.run(WARMUP + FIFO_LEN * SI)
+        n0 = sim.bodies.count
+        for c in ((30 * L0, YC, 0.0), (626.0, YC + 15, 0.0), (626.0, YC - 15, 0.0)):
+            sim.add_body("circle", center=c, radius=R)  # pending (negative) id until sync
+        sim.sync_bodies()
+        pinball = list(range(n0, n0 + 3))  # committed ids, as in f2f88c phase3
+        sim.run(WARMUP)
+        cells = float(len(sim.bodies.get(probes[0]).get_sensor_list(NX, NY)[0]))  # probe sum = mean * cells
+
+        obs = np.zeros(12, np.float32)
+        sensors, forces, actions = (np.zeros((args.intervals, n), np.float32) for n in (6, 6, 3))
+        for n in range(args.intervals):
+            a = np.asarray(model.predict(obs, deterministic=True)[0], np.float32).reshape(3)
+            target, eff = (8.0 * a + BIAS) * U0, np.zeros(3)
+            for k in range(SI):
+                eff = 0.9 * eff + 0.1 * target
+                for pid, w in zip(pinball, eff / R):
+                    sim.set_body(pid, omega=float(w))
+                sim.run(1, zero_obs=(k == 0), sync_obs=(k == SI - 1))
+            s = [v * cells for pid in probes for v in sim.read_sensor(pid)[:2]]
+            f = [v for pid in pinball for v in sim.read_force(pid)[:2]]
+            raw = np.asarray(s + f, np.float32)
+            if not np.all(np.isfinite(raw)):
+                raise RuntimeError(f"non-finite observation at interval {n}")
+            sensors[n], forces[n], actions[n] = raw[:6], raw[6:], a
+            obs = policy_observation(raw, norm)
+        np.savez_compressed(os.path.join(args.out, "signals.npz"), sensors=sensors,
+                            forces=forces, actions=actions)
+        with open(os.path.join(args.out, "case.json"), "w") as fh:
+            json.dump(dict(celeris=git_head(here), dynamis=git_head(args.dynamis),
+                           config=os.path.realpath(args.config), model=os.path.realpath(args.model),
+                           norm=os.path.realpath(args.norm), warmup=WARMUP, target_sampling=FIFO_LEN * SI, interval=SI, intervals=args.intervals,
+                           probe_cells=cells), fh, indent=1)
+    finally:
+        sim.close()
+
+
+if __name__ == "__main__":
+    main()

research/steady_pinball/config.json

@@ -1,49 +0,0 @@
-{
-    "_doc": "Legacy geometry Re_D=50, U0=0.01, SI=800 baseline (kar_legacy_re50d_u001)",
-    "grid": {
-        "lattice_model": "D2Q9",
-        "nx": 1280,
-        "ny": 512,
-        "nz": 1
-    },
-    "physics": {
-        "data_type": "FP32",
-        "viscosity": 0.004,
-        "velocity": 0.01,
-        "rho": 1.0
-    },
-    "method": {
-        "collision": "MRT",
-        "streaming": "double_buffer",
-        "store_precision": "FP32",
-        "ddf_shifting": false,
-        "les": {
-            "enabled": false,
-            "cs": 0.16,
-            "closed_form": true
-        },
-        "trt": {
-            "magic_param": 0.1875
-        },
-        "inlet": {
-            "profile": "parabolic",
-            "scheme": "regularized",
-            "regularized_neq_damp": 1.0
-        },
-        "outlet": {
-            "mode": "neq_extrap",
-            "backflow_clamp": true,
-            "blend_alpha": 0.7,
-            "srt_neq_damp": 0.5
-        },
-        "y_wall_bc": "bounce_back",
-        "omega_guard": {
-            "min": 0.01,
-            "max": 1.99
-        }
-    },
-    "cuda": {
-        "threads_per_block": 256,
-        "compute_capability": "auto"
-    }
-}

research/steady_pinball/run.py

@@ -1,167 +0,0 @@
-"""Uncontrolled paper-geometry steady Re_D=50, one transverse assembly DOF.
-
-Three predeclared zero-damping structures; --smoke uses only 32 + 32 steps.
-See docs/steady_pinball_fiv.md for historical solver differences and new structure.
-"""
-from __future__ import annotations
-
-import argparse
-import json
-import math
-from pathlib import Path
-import subprocess
-import time
-
-import numpy as np
-
-D, U, NU = 20., .01, .004
-CENTRES = ((600., 255.5), (626., 270.5), (626., 240.5))
-PROBES = ((800., 295.5), (800., 255.5), (800., 215.5))
-WARMUP, STEPS, PERIOD = 512000, 600000, 10000
-STRUCTURES = {'Z0': (10., 1.e-4), 'M2': (2., 1.e-4), 'R2': (2., 6.e-5)}
-COLUMNS = ('step', 'y', 'v', 'cme_F', 'cme_U', 'cme_L', 'wen_F', 'wen_U', 'wen_L',
-           'node_F', 'node_U', 'node_L', 'eq_only', 'no_donor', 'crossings')
-ROOT = Path(__file__).resolve().parents[2]
-
-
-def config():
-    return json.loads(Path(__file__).with_name('config.json').read_text())
-
-
-class Structure:
-    """Exact damped oscillator update for a force held over dt=1; no force filter."""
-    def __init__(self, mass_ratio=10., frequency=1.e-4, zeta=0.):
-        if not (math.isfinite(mass_ratio) and mass_ratio > 0. and
-                math.isfinite(frequency) and frequency > 0. and
-                math.isfinite(zeta) and 0. <= zeta < 1.):
-            raise ValueError('positive finite mass ratio/frequency and 0 <= zeta < 1 required')
-        self.mass_ratio, self.frequency = mass_ratio, frequency
-        self.mass = mass_ratio * 3. * math.pi * D * D / 4.
-        self.omega = 2. * math.pi * frequency
-        self.zeta = zeta
-        self.stiffness = self.mass * self.omega**2
-        self.damping = 2. * self.zeta * self.mass * self.omega
-        self.y = self.v = 0.
-        a = self.zeta * self.omega
-        wd = self.omega * math.sqrt(1. - self.zeta**2)
-        self.ec = math.exp(-a) * math.cos(wd)
-        self.es = math.exp(-a) * math.sin(wd) / wd
-        self.a = a
-
-    def candidate(self, force):
-        equilibrium = force / self.stiffness
-        z = self.y - equilibrium
-        return (equilibrium + (self.ec + self.a*self.es)*z + self.es*self.v,
-                -self.omega**2*self.es*z + (self.ec - self.a*self.es)*self.v)
-
-
-def summarize(rows):
-    windows = []
-    for start in range(0, len(rows), 10*PERIOD):
-        block = rows[start:start + 10*PERIOD]
-        fy = block[:, 3:9].sum(axis=1)
-        windows.append(dict(start_step=int(block[0, 0]), end_step=int(block[-1, 0]),
-                            y_mean_D=float(block[:, 1].mean()/D),
-                            y_rms_D=float(np.sqrt(np.mean(block[:, 1]**2))/D),
-                            y_peak_D=float(np.abs(block[:, 1]).max()/D),
-                            fy_mean=float(fy.mean()), fy_rms=float(np.sqrt(np.mean(fy**2))),
-                            mean_power=float(np.mean(fy*block[:, 2]))))
-    return windows
-
-
-def main():
-    ap = argparse.ArgumentParser(description=__doc__)
-    ap.add_argument('--out', required=True, type=Path)
-    ap.add_argument('--fixed', action='store_true', help='same modern moving-path geometry held fixed')
-    ap.add_argument('--smoke', action='store_true', help='32+32 steps, no physical result')
-    ap.add_argument('--case', choices=STRUCTURES, default='Z0',
-                    help='predeclared structure; all cases have zero structural damping')
-    args = ap.parse_args()
-    out = args.out.resolve()
-    if out == ROOT or ROOT in out.parents:
-        ap.error('output must be outside the checkout')
-    out.mkdir(parents=True, exist_ok=False)
-    cfg = config()
-    (out/'config.json').write_text(json.dumps(cfg, indent=2)+'\n')
-    source = subprocess.check_output(['git', '-C', str(ROOT), 'rev-parse', 'HEAD'], text=True).strip()
-    s = Structure(*STRUCTURES[args.case], zeta=0.)
-    warmup, count = (32, 32) if args.smoke else (WARMUP, STEPS)
-    rows = np.lib.format.open_memmap(out/'steps.npy', mode='w+', dtype=np.float64,
-                                     shape=(count, len(COLUMNS)))
-    rows[:] = np.nan
-    from CelerisLab import Simulation
-    from CelerisLab.body.motion import MovingTables, Translation
-    names = MovingTables.DIAG
-    col = {k: names.index(k) for k in names}
-    status, completed, f_last = 'running', 0, 0.
-    started = time.monotonic()
-    report = dict(source=source, case=args.case, fixed=args.fixed, smoke=args.smoke, columns=COLUMNS,
-                  structure=dict(mass=s.mass, stiffness=s.stiffness, damping=s.damping,
-                                 mass_ratio_three_disks=s.mass_ratio,
-                                 zeta=s.zeta, fn=s.frequency, Tn=1./s.frequency, y0=0., v0=0.),
-                  warmup_steps=warmup, requested_steps=count, dt=1., D=D, U=U, nu=NU,
-                  duration_Tn=count*s.frequency, duration_tU_D=count*U/D,
-                  summary_block_steps=10*PERIOD,
-                  rhs='sum of three instantaneous CME_y + Wen_y, one-step lag, no filtering')
-    try:
-        with Simulation(str(out/'config.json')) as sim:
-            # Match Legacy object order; probes remain passive and fixed in the lab.
-            for c in PROBES:
-                sim.add_body('sensor', c, 5.)
-            envelope = Translation((0., 0., -D, D), .05, shared=True)
-            bodies = [sim.add_body('circle', c, D/2, translation=envelope) for c in CENTRES]
-            sim.initialize()
-            # At release obtain the last instantaneous force, not a warmup average.
-            remaining = warmup - 1
-            while remaining:
-                chunk = min(10000, remaining)
-                sim.run(chunk)
-                remaining -= chunk
-            sim.run(1)
-            d = sim.moving_diagnostics()
-            f_last = float(d[:, col['cme_y']].sum(dtype=np.float64) + d[:, col['wen_y']].sum(dtype=np.float64))
-            report['warmup_elapsed_s'] = time.monotonic()-started
-            print(json.dumps(dict(event='release', **report)), flush=True)
-            for n in range(count):
-                y1, v1 = (0., 0.) if args.fixed else s.candidate(f_last)
-                if not all(map(math.isfinite, (y1, v1))) or abs(y1) > D:
-                    status = 'stopped_nonfinite_or_travel_bound'
-                    break
-                if not args.fixed:
-                    for bid in bodies:
-                        sim.set_body(bid, offset=(0., y1), vy=y1-s.y)
-                sim.run(1)
-                d = sim.moving_diagnostics()
-                cme, wen = d[:, col['cme_y']], d[:, col['wen_y']]
-                node = d[:, col['cover_jy']] - d[:, col['uncover_jy']]
-                row = (n+1, y1, v1, *cme, *wen, *node,
-                       d[:, col['eq_only']].sum(), d[:, col['no_donor']].sum(), d[:, col['crossings']].sum())
-                if not np.isfinite(row).all() or row[-2] != 0:
-                    status = 'stopped_nonfinite_or_no_donor'
-                    break
-                rows[n] = row
-                completed = n+1
-                s.y, s.v = y1, v1
-                f_last = float(cme.sum(dtype=np.float64)+wen.sum(dtype=np.float64))
-                if completed % PERIOD == 0:
-                    rows.flush()
-                    print(json.dumps(dict(event='progress', step=completed, t_Tn=completed*s.frequency,
-                                          y_D=s.y/D, elapsed_s=time.monotonic()-started)), flush=True)
-            else:
-                status = 'smoke_completed' if args.smoke else 'duration_completed_not_qualified'
-    except BaseException as exc:
-        status = f'error: {type(exc).__name__}: {exc}'
-        raise
-    finally:
-        rows.flush()
-        report.update(status=status, completed_steps=completed, elapsed_s=time.monotonic()-started,
-                      windows=summarize(rows[:completed]),
-                      totals={key: float(rows[:completed, i].sum()) for key, i in
-                              [('eq_only', 12), ('no_donor', 13), ('crossings', 14)]},
-                      limitation='Uncontrolled development response only. No automatic settling, grid convergence, control or historical-solver equivalence claim.')
-        (out/'report.json').write_text(json.dumps(report, indent=2)+'\n')
-        print(json.dumps(report), flush=True)
-
-
-if __name__ == '__main__':
-    main()

src/CelerisLab/body/motion.py

@@ -14,7 +14,6 @@ from typing import Callable, Optional, Tuple
 
 MAX_STEP = 0.5     # Euclidean centre displacement per step (lattice cells)
 BAND_MARGIN = 2    # cells added around the swept circle AABB
-SHARED_CLEARANCE = 4.0  # > two diagonal links + MAX_STEP; docs/steady_pinball_fiv.md
 R_MIN = 4.0        # smallest radius covered by the donor enumeration
 
 Vec = Tuple[float, float]
@@ -32,7 +31,6 @@ class Translation:
     """Declared motion envelope: offset box (x_min, x_max, y_min, y_max) and speed bound."""
     offset_box: Tuple[float, float, float, float]
     max_speed: float
-    shared: bool = False  # same envelope instance: synchronous transverse motion
 
     def __post_init__(self):
         box = tuple(float(v) for v in self.offset_box)
@@ -43,29 +41,10 @@ class Translation:
             raise ValueError("offset_box must contain the zero offset")
         if not 0.0 < speed <= MAX_STEP:
             raise ValueError(f"max_speed must satisfy 0 < max_speed <= {MAX_STEP}")
-        if self.shared and box[:2] != (0.0, 0.0):
-            raise ValueError("shared translation is transverse only (x offsets must be zero)")
         object.__setattr__(self, "offset_box", box)
         object.__setattr__(self, "max_speed", speed)
 
 
-def validate_shared_translations(circles):
-    """Check common old/target offsets and wall velocities before device work.
-
-    Per-body setters may stage members in any order. Envelope identity (not
-    value equality) defines a group; a missing target means hold at applied.
-    """
-    groups = {}
-    for motion in circles:
-        if not motion.translation.shared:
-            continue
-        target = motion.pending if motion.pending is not None else (motion.applied, (0.0, 0.0))
-        state = (motion.applied, target)
-        key = id(motion.translation)
-        if groups.setdefault(key, state) != state:
-            raise ValueError("shared translation needs identical applied/target offsets and velocities")
-
-
 class MovingCircle:
     """Applied offset, pending target and the fail-closed commit of one circle."""
 
@@ -93,8 +72,6 @@ class MovingCircle:
         x0, x1, y0, y1 = self.translation.offset_box
         if not (x0 <= offset[0] <= x1 and y0 <= offset[1] <= y1):
             raise ValueError(f"offset {offset} outside declared box {self.translation.offset_box}")
-        if self.translation.shared and velocity[0] != 0.0:
-            raise ValueError("shared translation needs zero x velocity")
         step = math.hypot(offset[0] - self.applied[0], offset[1] - self.applied[1])
         if step > MAX_STEP:
             raise ValueError(f"per-step displacement {step} exceeds {MAX_STEP}")

src/CelerisLab/simulation.py

@@ -259,8 +259,6 @@ class Simulation:
         moving = self._moving_bodies()
         if not moving:
             return
-        from .body.motion import SHARED_CLEARANCE, validate_shared_translations
-        validate_shared_translations([obj.motion for obj in moving])
         cfg = self.lbm_cfg
         if (not cfg.is_d2q9 or cfg.streaming != "double_buffer" or cfg.ddf_shifting
                 or cfg.store_precision != "FP32"):
@@ -279,13 +277,6 @@ class Simulation:
                 if other is obj:
                     continue
                 (cx, cy), r = other.center[:2], float(other.radius)
-                if (other.motion is not None and obj.motion.translation.shared
-                        and obj.motion.translation is other.motion.translation):
-                    # Whole envelopes may overlap; simultaneous local stencils may not.
-                    gap = np.hypot(obj.center[0] - cx, obj.center[1] - cy) - obj.radius - r
-                    if gap <= SHARED_CLEARANCE:
-                        raise ValueError("shared translation needs surface clearance > 4 cells")
-                    continue
                 if other.motion is not None:
                     ox0, ox1, oy0, oy1 = other.motion.band()
                 else:
@@ -728,9 +719,6 @@ class Simulation:
         if steps < 0 or checkpoint_interval < 0:
             raise ValueError("steps and checkpoint_interval must be nonnegative")
         moving = self._moving_bodies()
-        if steps > 0 and moving:
-            from .body.motion import validate_shared_translations
-            validate_shared_translations([obj.motion for obj in moving])
         if steps > 0 and moving and self.stepper.moving is None:
             raise RuntimeError(MOVING_BLOCKED)
         if moving and checkpoint_interval > 0:

tests/integration/test_shared_translation_gpu.py

@@ -1,62 +0,0 @@
-"""Small serial GPU gates for shared transverse circles; no physical benchmark."""
-import json
-from pathlib import Path
-
-import numpy as np
-
-from CelerisLab import Simulation
-from CelerisLab.body.motion import MovingTables, Translation
-
-CENTRES = ((70., 79.5), (96., 94.5), (96., 64.5))
-
-
-def simulate(tmp_path, *, moving, offsets, reverse=False):
-    config = json.loads((Path(__file__).resolve().parents[2]/'research/steady_pinball/config.json').read_text())
-    config['grid'].update(nx=192, ny=160)
-    p = tmp_path/'config.json'
-    p.write_text(json.dumps(config))
-    with Simulation(str(p)) as sim:
-        envelope = Translation((0., 0., -20., 20.), .5, shared=True)
-        centres = list(reversed(CENTRES)) if reverse else list(CENTRES)
-        ids = [sim.add_body('circle', c, 10., translation=envelope if moving else None) for c in centres]
-        sim.initialize()
-        previous = 0.
-        crossings = 0.
-        for y in offsets:
-            if moving:
-                for bid in ids:
-                    sim.set_body(bid, offset=(0., float(y)), vy=float(y-previous))
-            sim.run(1)
-            if moving:
-                diag = sim.moving_diagnostics()
-                assert np.isfinite(diag).all()
-                assert np.all(diag[:, MovingTables.DIAG.index('no_donor')] == 0)
-                crossings += diag[:, MovingTables.DIAG.index('crossings')].sum()
-            previous = y
-        forces = np.array([sim.read_force(bid) for bid in ids])
-        if reverse:
-            forces = forces[::-1]
-        return sim.get_ddf(), sim.get_flags(), forces, crossings
-
-
-def test_fixed_vs_shared_hold(tmp_path):
-    fixed = simulate(tmp_path, moving=False, offsets=[0.]*32)
-    held = simulate(tmp_path, moving=True, offsets=[0.]*32)
-    assert np.array_equal(fixed[0], held[0], equal_nan=True)
-    assert np.array_equal(fixed[1], held[1])
-    # FP32 atomic force reduction order only; populations remain exactly equal.
-    np.testing.assert_allclose(fixed[2], held[2], rtol=0., atol=1e-4)
-    assert held[3] == 0
-
-
-def test_common_motion_member_order_has_no_population_effect(tmp_path):
-    # MAX_STEP=.5 is a geometric bound, not a low-Mach stability guarantee.
-    # Keep the complete +/-D excursion at wall speed .01, equal to mean inflow.
-    offsets = np.concatenate((np.arange(1, 2001), np.arange(1999, -2001, -1),
-                              np.arange(-1999, 1))) * .01
-    forward = simulate(tmp_path, moving=True, offsets=offsets)
-    reverse = simulate(tmp_path, moving=True, offsets=offsets, reverse=True)
-    assert forward[3] > 0 and forward[3] == reverse[3]
-    assert np.array_equal(forward[0], reverse[0], equal_nan=True)
-    assert np.array_equal(forward[1], reverse[1])
-    np.testing.assert_allclose(forward[2], reverse[2], rtol=0., atol=1e-4)

tests/unit/test_shared_translation.py

@@ -1,122 +0,0 @@
-"""Host admission and independent lattice-stencil checks; no CUDA context."""
-from types import SimpleNamespace
-
-import numpy as np
-import pytest
-
-from CelerisLab.body.motion import MovingCircle, MovingTables, Translation, validate_shared_translations
-from CelerisLab.lbm.descriptors import FLUID
-from tests.unit.test_moving_bodies import STREAM, body, moving_simulation, routed_simulation
-
-
-def scene(monkeypatch, envelope=None, gap=10.):
-    sim = routed_simulation(monkeypatch)
-    envelope = envelope or Translation((0., 0., -6., 6.), .5, shared=True)
-    for y in (40. - (20. + gap) / 2., 40. + (20. + gap) / 2.):
-        sim.bodies.add(body(60., y, 10., envelope))
-    sim.field = SimpleNamespace(nx=120, ny=80, build_channel_flags=lambda: np.full((80, 120), FLUID))
-    return sim
-
-
-def test_shared_bands_are_admitted_but_other_groups_are_not(monkeypatch):
-    sim = scene(monkeypatch)
-    sim._validate_moving_bodies()
-    a, b = [o.motion for o in sim.bodies.objects]
-    assert a.band()[3] >= b.band()[2]  # whole envelopes genuinely overlap
-    band, geom = MovingTables.layout([(a, 0), (b, 1)])
-    assert band[10] == band[2] * band[3] * 8
-    assert len(geom) == 14
-    b.translation = Translation((0., 0., -6., 6.), .5, shared=True)
-    with pytest.raises(ValueError, match='overlaps body'):
-        sim._validate_moving_bodies()  # equal values do not link independent motion
-    b.translation = a.translation
-    sim.bodies.add(body(60., 40., 2.))
-    with pytest.raises(ValueError, match='overlaps body'):
-        sim._validate_moving_bodies()
-
-
-@pytest.mark.parametrize('gap', (.3, 4.))
-def test_small_surface_gaps_rejected(monkeypatch, gap):
-    with pytest.raises(ValueError, match='clearance'):
-        scene(monkeypatch, gap=gap)._validate_moving_bodies()
-
-
-def test_shared_envelope_is_transverse_only():
-    with pytest.raises(ValueError, match='transverse'):
-        Translation((-.1, .1, -1., 1.), .5, shared=True)
-
-    envelope = Translation((0., 0., -1., 1.), .5, shared=True)
-    with pytest.raises(ValueError, match="zero x velocity"):
-        MovingCircle((0., 0.), 10., envelope).check_target((0., .1), (.1, .1))
-
-def test_group_rejection_is_before_writes_and_can_be_corrected(monkeypatch):
-    sim, writes = moving_simulation(monkeypatch)
-    env = Translation((0., 0., -2., 2.), .5, shared=True)
-    sim.bodies.objects = [body(40., y, 10., env) for y in (20., 50.)]
-    a, b = [o.motion for o in sim.bodies.objects]
-    sim.stepper.moving = object()
-    a.set_target((0., .2), (0., .2))
-    with pytest.raises(ValueError, match='identical'):
-        sim.run(1, stream=STREAM)
-    assert writes == [] and a.applied == b.applied == (0., 0.)
-    b.set_target((0., .2), (0., .1))
-    with pytest.raises(ValueError, match='identical'):
-        sim.run(1, stream=STREAM)
-    assert writes == []
-    b.set_target((0., .2), (0., .2))
-    validate_shared_translations([a, b])
-    b.applied = (0., -.1)
-    with pytest.raises(ValueError, match='identical'):
-        validate_shared_translations([a, b])
-
-
-C = np.array([(x, y) for x in (-1, 0, 1) for y in (-1, 0, 1) if (x, y) != (0, 0)])
-
-
-def footprint_conflicts(centres, radius, dy):
-    """Enumerate geometric write/read supports, not refill/BFL formulas.
-
-    Tests every candidate donor, conservatively including M2's second link.
-    Body-indexed slot numbers are separate from physical (cell,direction) writes.
-    """
-    centres = np.asarray(centres, float)
-    lo = np.floor(centres.min(0) - radius - 4).astype(int)
-    hi = np.ceil(centres.max(0) + radius + 4).astype(int)
-    yy, xx = np.mgrid[lo[1]:hi[1]+1, lo[0]:hi[0]+1]
-    points = np.column_stack((xx.ravel(), yy.ravel()))
-    def inside(p, c):
-        return ((p - c)**2).sum(axis=-1) < radius**2
-    old = np.array([inside(points, c) for c in centres])
-    new_centres = centres + (0., dy)
-    new = np.array([inside(points, c) for c in new_centres])
-    changed = old ^ new
-    conflicts = int(np.count_nonzero(changed.sum(0) > 1))
-    writes = set()
-    for i, (oc, nc) in enumerate(zip(centres, new_centres)):
-        others = [j for j in range(len(centres)) if j != i]
-        fresh = points[old[i] & ~new[i]]
-        for e in C:
-            donor = fresh + e
-            valid = ~inside(donor, oc) & ~inside(donor, nc)
-            for j in others:
-                conflicts += int(np.count_nonzero(valid & (inside(donor, centres[j]) | inside(donor, new_centres[j]))))
-            link = ~new[i] & inside(points + e, nc)
-            fluid = points[link]
-            solid = fluid + e
-            for p in solid:
-                slot = (int(p[0]), int(p[1]), int(e[0]), int(e[1]))
-                conflicts += int(slot in writes)
-                writes.add(slot)
-            for read in (fluid, fluid - e):
-                for j in others:
-                    conflicts += int(np.count_nonzero(inside(read, centres[j]) | inside(read, new_centres[j])))
-    return conflicts
-
-
-def test_counterexample_and_pinball_local_ownership():
-    # Close but disjoint circles are insufficient: parent/Grok independent example.
-    assert footprint_conflicts(((.25, .25), (.25, 8.55)), 4., .5) > 0
-    for phase in (0., .25, .5, .75):
-        for dy in (-.5, -.1, .1, .5):
-            centres = np.array(((0., 0.), (26., 15.), (26., -15.))) + phase
-            assert footprint_conflicts(centres, 10., dy) == 0

tests/unit/test_steady_pinball.py

@@ -1,84 +0,0 @@
-"""Real config loading and pure structural update for the frozen development case."""
-import importlib.util
-import math
-from pathlib import Path
-
-import numpy as np
-import pytest
-
-path = Path(__file__).resolve().parents[2] / 'research/steady_pinball/run.py'
-spec = importlib.util.spec_from_file_location('steady_pinball_run', path)
-runner = importlib.util.module_from_spec(spec)
-spec.loader.exec_module(runner)
-
-
-def test_actual_loader_and_paper_geometry(tmp_path):
-    import json
-    from CelerisLab.config import load_lbm_config
-    p = tmp_path/'config.json'
-    p.write_text(json.dumps(runner.config()))
-    cfg = load_lbm_config(str(p))
-    cfg.validate()
-    assert (cfg.nx, cfg.ny, cfg.viscosity, cfg.velocity) == (1280, 512, .004, .01)
-    assert cfg.collision == 'MRT' and cfg.y_wall_bc == 'bounce_back'
-    assert runner.U * runner.D / runner.NU == 50.
-    assert runner.CENTRES == ((600., 255.5), (626., 270.5), (626., 240.5))
-
-
-def test_structure_rest_equilibrium_and_energy_decay():
-    s = runner.Structure(zeta=.01)
-    assert s.candidate(0.) == (0., 0.)
-    assert s.mass == pytest.approx(10*3*math.pi*100)
-    s.y, s.v = .2, 0.
-    # One damped period: independent continuous solution at total elapsed time.
-    for _ in range(runner.PERIOD):
-        s.y, s.v = s.candidate(0.)
-    t = runner.PERIOD
-    a = s.zeta*s.omega
-    w = s.omega*math.sqrt(1-s.zeta**2)
-    expected = .2*math.exp(-a*t)*(math.cos(w*t)+a/w*math.sin(w*t))
-    assert s.y == pytest.approx(expected, abs=2e-12)
-    assert .5*s.mass*s.v*s.v+.5*s.stiffness*s.y*s.y < .5*s.stiffness*.2**2
-    s.y, s.v = .3, 0.
-    assert s.candidate(s.stiffness*.3) == pytest.approx((.3, 0.), abs=1e-15)
-
-
-def test_report_uses_completed_rows_and_net_force():
-    rows = np.zeros((4, 15))
-    rows[:, 0] = np.arange(1, 5)
-    rows[:, 1] = 2.
-    rows[:, 3:9] = 1.
-    got = runner.summarize(rows)
-    assert got[0]['end_step'] == 4
-    assert got[0]['y_rms_D'] == .1
-    assert got[0]['fy_rms'] == 6.
-
-
-@pytest.mark.parametrize('case', ['Z0', 'M2', 'R2'])
-def test_zero_damping_cases_conserve_unforced_energy_and_match_forced_solution(case):
-    s = runner.Structure(*runner.STRUCTURES[case])
-    assert s.zeta == s.damping == 0.
-    s.y, s.v = .2, 1.e-5
-    y0, v0 = s.y, s.v
-    energy = .5*s.mass*v0**2 + .5*s.stiffness*y0**2
-    for _ in range(12345):
-        s.y, s.v = s.candidate(0.)
-    assert .5*s.mass*s.v**2 + .5*s.stiffness*s.y**2 == pytest.approx(energy, rel=2e-11)
-    force, elapsed = 1.e-4, 3210
-    s.y, s.v = y0, v0
-    for _ in range(elapsed):
-        s.y, s.v = s.candidate(force)
-    phase = s.omega*elapsed
-    equilibrium = force/s.stiffness
-    expected_y = equilibrium+(y0-equilibrium)*math.cos(phase)+v0*math.sin(phase)/s.omega
-    expected_v = -(y0-equilibrium)*s.omega*math.sin(phase)+v0*math.cos(phase)
-    assert (s.y, s.v) == pytest.approx((expected_y, expected_v), abs=2e-11)
-
-
-def test_predeclared_comparisons_change_structure_only():
-    z, m, r = [runner.Structure(*runner.STRUCTURES[k]) for k in ['Z0', 'M2', 'R2']]
-    assert m.mass/z.mass == pytest.approx(.2)
-    assert m.stiffness/z.stiffness == pytest.approx(.2)
-    assert r.mass == m.mass
-    assert r.stiffness/m.stiffness == pytest.approx(.36)
-    assert [runner.STEPS*s.frequency for s in [z, m, r]] == pytest.approx([60., 60., 36.])