body: host-only translating-circle state (S1, no kernels)

MovingCircle keeps the applied offset and an optional pending target apart,
validates targets against the applied offset (box, Euclidean 0.5-cell step,
speed, no sliding hold) and commits only after upload, step and sync succeed;
any failure after the upload attempt is fail-closed. Simulation.add_body
(translation=...) declares one; moving circles are excluded from the packed
static links, unsupported modes and swept-band placement are rejected at
initialize, and every positive-step run raises before any device write until
the moving kernels exist. Host tests include a refill-donor enumeration.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This commit is contained in:
Claude
2026-10-30 06:34:06 +00:00
co-authored by Claude Opus 5.5
parent c986343a15
commit 305ab2409b
5 changed files with 370 additions and 3 deletions
+3 -1
View File
@@ -405,7 +405,9 @@ class ObjectManager:
for obj in objects:
obj.validate_surface_actuation(self.cfg)
if (not obj.is_sensor) and (not obj.is_force_region) and hasattr(obj, 'get_curved_list'):
# Moving circles never enter the packed static links (docs/moving_bodies.md, section 3).
if (not obj.is_sensor) and (not obj.is_force_region) and obj.motion is None \
and hasattr(obj, 'get_curved_list'):
(
fluid_idx, dirs, q_vals, body_ids,
rx_vals, ry_vals, rz_vals, fallback_vals
+163
View File
@@ -0,0 +1,163 @@
# CelerisLab/body/motion.py
"""Host state of a translating circle (dense double buffer, D2Q9 only).
Contract: docs/moving_bodies.md, sections 4a, 5 and 7. The applied geometry
(last successfully committed offset) and an optional pending target are two
separate facts; centres are always ``base + offset``. Nothing here touches
the device: the caller supplies upload/step/synchronize callables.
"""
from __future__ import annotations
import math
from dataclasses import dataclass
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
R_MIN = 4.0 # smallest radius covered by the donor enumeration
Vec = Tuple[float, float]
# D2Q9 non-rest directions, paired order of lbm/descriptors.
_CX = (1, -1, 0, 0, 1, -1, 1, -1)
_CY = (0, 0, 1, -1, 1, -1, -1, 1)
def _vec(value, name: str) -> Vec:
x, y = (float(v) for v in value)
if not (math.isfinite(x) and math.isfinite(y)):
raise ValueError(f"{name} must be finite, got {value!r}")
return x, y
@dataclass(frozen=True)
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
def __post_init__(self):
box = tuple(float(v) for v in self.offset_box)
speed = float(self.max_speed)
if len(box) != 4 or not all(math.isfinite(v) for v in box + (speed,)):
raise ValueError("offset_box needs 4 finite values and max_speed must be finite")
if not (box[0] <= 0.0 <= box[1] and box[2] <= 0.0 <= box[3]):
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}")
object.__setattr__(self, "offset_box", box)
object.__setattr__(self, "max_speed", speed)
class MovingCircle:
"""Applied offset, pending target and the fail-closed commit of one circle."""
def __init__(self, base, radius: float, translation: Translation):
self.base = _vec(base, "base")
self.radius = float(radius)
if not (math.isfinite(self.radius) and self.radius >= R_MIN):
raise ValueError(f"moving circle radius must be >= {R_MIN}")
if not isinstance(translation, Translation):
raise TypeError("translation must be a Translation")
self.translation = translation
self.applied: Vec = (0.0, 0.0)
self.pending: Optional[Tuple[Vec, Vec]] = None
self.failed = False
def centre(self, offset: Vec) -> Vec:
return self.base[0] + offset[0], self.base[1] + offset[1]
def check_target(self, offset, velocity) -> Tuple[Vec, Vec]:
"""Validate a target against ``applied``; returns it, mutates nothing."""
self._require_usable()
offset, velocity = _vec(offset, "offset"), _vec(velocity, "velocity")
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}")
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}")
if math.hypot(*velocity) > self.translation.max_speed:
raise ValueError(f"velocity {velocity} exceeds max_speed {self.translation.max_speed}")
if offset == self.applied and velocity != (0.0, 0.0):
raise ValueError("a hold target (offset == applied) needs zero translational velocity")
return offset, velocity
def set_target(self, offset, velocity) -> None:
"""Replace the pending target (last wins); invalid input keeps the old one."""
self.pending = self.check_target(offset, velocity)
def record(self) -> Tuple[Vec, Vec, Vec]:
"""Motion record (c_old, c_new, u) for the next single moving step."""
self._require_usable()
if self.pending is not None:
offset, velocity = self.pending
else:
offset, velocity = self.applied, (0.0, 0.0)
return self.centre(self.applied), self.centre(offset), velocity
def advance(self, upload: Callable, step: Callable, synchronize: Callable) -> None:
"""Upload the record, run one step and commit only after a successful sync.
Any exception after the upload is attempted leaves device geometry
unknown: the circle is marked failed and never reports applied state again.
"""
record = self.record()
try:
upload(record)
step()
synchronize()
except BaseException:
self.failed = True
raise
if self.pending is not None:
self.applied = self.pending[0]
self.pending = None
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
r, m = self.radius, BAND_MARGIN
return (math.floor(self.base[0] + x0 - r) - m, math.ceil(self.base[0] + x1 + r) + m,
math.floor(self.base[1] + y0 - r) - m, math.ceil(self.base[1] + y1 + r) + m)
def _require_usable(self) -> None:
if self.failed:
raise RuntimeError("moving circle state is unknown after a failed step")
def _inside(x: int, y: int, c: Vec, r: float) -> bool:
return (x - c[0]) ** 2 + (y - c[1]) ** 2 < r * r
def empty_donor_cells(c_old: Vec, c_new: Vec, radius: float):
"""Fresh cells (inside c_old, outside c_new) whose donor set V(x) is empty."""
r = radius
xs = range(math.floor(min(c_old[0], c_new[0]) - r) - 1, math.ceil(max(c_old[0], c_new[0]) + r) + 2)
ys = range(math.floor(min(c_old[1], c_new[1]) - r) - 1, math.ceil(max(c_old[1], c_new[1]) + r) + 2)
empty = []
for x in xs:
for y in ys:
if not _inside(x, y, c_old, r) or _inside(x, y, c_new, r):
continue
if not any(not _inside(x + cx, y + cy, c_old, r) and not _inside(x + cx, y + cy, c_new, r)
for cx, cy in zip(_CX, _CY)):
empty.append((x, y))
return empty
def enumerate_donor_coverage(radius: float, positions: int = 8, directions: int = 16) -> int:
"""Count empty-donor fresh cells over sub-cell centres and moves of MAX_STEP.
Evidence for the declared bounds (radius >= R_MIN, |dc|_2 <= MAX_STEP),
not a proof; the device keeps the ``no_donor`` counter.
"""
count = 0
for i in range(positions):
for j in range(positions):
c_old = (i / positions, j / positions)
for k in range(directions):
a = 2.0 * math.pi * k / directions
c_new = (c_old[0] + MAX_STEP * math.cos(a), c_old[1] + MAX_STEP * math.sin(a))
count += len(empty_donor_cells(c_old, c_new, radius))
return count
+1
View File
@@ -140,6 +140,7 @@ class SimObject:
x=float(c[0]), y=float(c[1]), z=z0,
)
self.control = ObjectControl()
self.motion = None # MovingCircle for a translating circle (body/motion.py)
@property
def surface_actuation(self):
+72 -2
View File
@@ -26,6 +26,7 @@ import pycuda.driver as cuda
from .config import LBMConfig, BodyConfig, load_lbm_config, load_body_config
from .cuda.context import CudaContext, context_bound
from .cuda import compiler
from .lbm.descriptors import FLUID
from .lbm.field import LBMField
from .lbm.stepper import LBMStepper
from .body.objects import SimObject
@@ -97,7 +98,8 @@ class Simulation:
# -- Object management ---------------------------------------------------
def add_body(self, type: str = "circle", center=None,
radius: float = 0.0, *, surface_actuation=None, **kwargs) -> int:
radius: float = 0.0, *, surface_actuation=None, translation=None,
**kwargs) -> int:
"""Add a simulation body. Returns body_id.
Before ``initialize()`` the object is added to the formal registry
@@ -113,6 +115,9 @@ class Simulation:
angle, half_width (radians), normal_velocity (positive outward).
D2Q9/FP32/unshifted/DB only; set before initialization or stage
a new body. No runtime motion or amplitude setter.
translation: ``body.motion.Translation`` declaring a translating
circle (before initialize only; docs/moving_bodies.md). Until
the moving kernels exist, positive-step runs are rejected.
**kwargs: passed to the geometry constructor.
Returns:
@@ -123,6 +128,10 @@ class Simulation:
type_lower = str(type).strip().lower()
if surface_actuation is not None and type_lower != "circle":
raise ValueError("Surface actuation is only supported on circles")
if translation is not None and (type_lower != "circle" or surface_actuation is not None
or self._initialized):
raise ValueError("translation needs a circle without surface actuation, "
"added before initialize()")
if type_lower == "circle":
from .body.geometry.circle import CircleGeometry
@@ -131,6 +140,9 @@ class Simulation:
obj = SimObject(obj_id=-1, geometry=geom,
center=center, radius=radius,
is_sensor=False, surface_actuation=surface_actuation)
if translation is not None:
from .body.motion import MovingCircle
obj.motion = MovingCircle((cx, cy), radius, translation)
elif type_lower == "sensor":
from .body.geometry.circle import CircleGeometry
cx, cy = float(center[0]), float(center[1])
@@ -190,6 +202,43 @@ class Simulation:
raise
self.bodies.jets[index] = jet # later geometry syncs keep the new amplitude
def _moving_bodies(self):
return [obj for obj in self.bodies.objects if obj.motion is not None]
def _reject_moving_bodies(self, operation: str) -> None:
if self._moving_bodies():
raise ValueError(f"{operation} is not supported with translating circles")
def _validate_moving_bodies(self) -> None:
"""Reject unsupported modes and swept-band placement before any mutation."""
moving = self._moving_bodies()
if not moving:
return
cfg = self.lbm_cfg
if (not cfg.is_d2q9 or cfg.streaming != "double_buffer" or cfg.ddf_shifting
or cfg.store_precision != "FP32"):
raise ValueError("translating circles need D2Q9/FP32/unshifted/double_buffer")
if self.bodies.jets:
raise ValueError("Jets are not supported together with translating circles")
nx, ny = self.field.nx, self.field.ny
base = self.field.build_channel_flags().reshape(ny, nx)
for obj in moving:
x0, x1, y0, y1 = obj.motion.band()
if x0 < 1 or y0 < 1 or x1 > nx - 2 or y1 > ny - 2:
raise ValueError(f"swept band {(x0, x1, y0, y1)} touches the domain edge")
if not np.all((base[y0:y1 + 1, x0:x1 + 1] & FLUID) != 0):
raise ValueError(f"swept band {(x0, x1, y0, y1)} contains non-fluid base cells")
for other in self.bodies.objects:
if other is obj:
continue
(cx, cy), r = other.center[:2], float(other.radius)
if other.motion is not None:
ox0, ox1, oy0, oy1 = other.motion.band()
else:
ox0, ox1, oy0, oy1 = cx - r - 2, cx + r + 2, cy - r - 2, cy + r + 2
if ox0 <= x1 and x0 <= ox1 and oy0 <= y1 and y0 <= oy1:
raise ValueError(f"swept band {(x0, x1, y0, y1)} overlaps body {other.obj_id}")
def _require_usable(self):
failed = getattr(self, "_failed", None)
if failed is not None:
@@ -260,7 +309,7 @@ class Simulation:
# -- Runtime control -----------------------------------------------------
def set_body(self, id: int, omega: float = None,
vx: float = None, vy: float = None) -> None:
vx: float = None, vy: float = None, offset=None) -> None:
"""Set runtime body state (host action array only, no H2D).
The GPU action buffer is updated at the next ``run()`` call.
@@ -270,8 +319,20 @@ 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; with vx/vy
(default 0) it replaces the pending target after validation
against the applied offset (docs/moving_bodies.md, 4a).
"""
obj = self.bodies.get(id)
motion = obj.motion
if motion is not None:
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))
proposed.update({k: float(v) for k, v in (("omega", omega), ("vx", vx), ("vy", vy))
if v is not None})
@@ -287,6 +348,8 @@ 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
if changed:
self.bodies._refresh_action_from_objects()
@@ -468,6 +531,7 @@ class Simulation:
# Validate the complete geometry before flags, populations or modules
# change. Compact-list construction publishes no outputs on failure.
next_flags = self.bodies.build_flags(self.field.build_channel_flags())
self._validate_moving_bodies()
if self.bodies.jets and self.bodies.count == 0:
raise ValueError("Jets need bodies with cut links")
self.bodies.build_compact_lists_for(self.bodies.objects, next_flags)
@@ -499,6 +563,7 @@ class Simulation:
self._require_usable()
if not self._initialized:
raise RuntimeError("Call initialize() before sync_bodies()")
self._reject_moving_bodies("sync_bodies()")
self._assert_allocation_contract()
# 1. Synchronize GPU
@@ -599,6 +664,9 @@ 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():
raise RuntimeError("translating circles need the moving kernels (not implemented); "
"positive-step runs are rejected")
self._assert_runtime_contracts()
# Discard any uncommitted body edits before stepping.
if self.bodies.has_pending_edit():
@@ -698,6 +766,7 @@ class Simulation:
"""Save full simulation state to an HDF5 file."""
self._require_usable()
from .common.checkpoint import save_checkpoint
self._reject_moving_bodies("save_checkpoint()")
if self.bodies.jets:
raise NotImplementedError("Jets are not checkpointed yet")
return save_checkpoint(self.field, self.stepper, self.lbm_cfg,
@@ -708,6 +777,7 @@ class Simulation:
"""Restore simulation state from an HDF5 checkpoint."""
self._require_usable()
from .common.checkpoint import load_checkpoint
self._reject_moving_bodies("load_checkpoint()")
if not self._initialized:
self.initialize()
load_checkpoint(path, self.field, self.stepper, self.lbm_cfg,
+131
View File
@@ -0,0 +1,131 @@
"""Translating-circle host contract, S1 (docs/moving_bodies.md 4a, 5, 7).
Host doubles only: no kernels exist yet, so every positive-step run with a
moving body must be rejected before any device write.
"""
from types import SimpleNamespace
import pytest
from CelerisLab.body.motion import (MAX_STEP, R_MIN, MovingCircle, Translation,
enumerate_donor_coverage)
BAD = (float("nan"), float("inf"))
def circle(box=(-2., 2., -3., 3.), speed=.1):
return MovingCircle((50., 40.), 10., Translation(box, speed))
def test_translation_and_radius_validation():
for box in ((0., 1., -1., 1.)[:3], (.5, 1., -1., 1.), (-1., 1., BAD[0], 1.)):
with pytest.raises(ValueError):
Translation(box, .1)
for speed in (0., MAX_STEP + 1e-9) + BAD:
with pytest.raises(ValueError):
Translation((0., 0., -1., 1.), speed)
with pytest.raises(ValueError):
MovingCircle((0., 0.), R_MIN - .1, Translation((0., 0., -1., 1.), .1))
def test_targets_validate_against_applied_and_last_wins():
c = circle()
c.set_target((0., .3), (0., .05))
c.set_target((0., .5), (0., .05)) # second setter replaces, checked vs applied (0, 0)
assert c.pending == ((0., .5), (0., .05))
for offset, velocity in (((0., .51), (0., 0.)), # Euclidean bound
((.36, .36), (0., 0.)), # |dc|_2 = .509
((0., 3.1), (0., 0.)), # outside box
((0., .1), (0., .11)), # speed bound
((0., 0.), (0., .01)), # hold with sliding velocity
((BAD[0], 0.), (0., 0.)), ((0., 0.), (0., BAD[1]))):
with pytest.raises(ValueError):
c.set_target(offset, velocity)
assert c.pending == ((0., .5), (0., .05)) # invalid replacement keeps the old target
c.set_target((.35, .35), (0., 0.)) # |dc|_2 = .495 inside the bound
assert c.applied == (0., 0.)
def test_record_commit_and_hold():
c = circle()
assert c.record() == ((50., 40.), (50., 40.), (0., 0.)) # hold: no translation velocity
c.set_target((0., .25), (0., .05))
assert c.record() == ((50., 40.), (50., 40.25), (0., .05))
uploads = []
c.advance(uploads.append, lambda: None, lambda: None)
assert uploads == [((50., 40.), (50., 40.25), (0., .05))]
assert c.applied == (0., .25) and c.pending is None # offset-only commit
c.advance(uploads.append, lambda: None, lambda: None) # no new target: hold step
assert uploads[-1] == ((50., 40.25), (50., 40.25), (0., 0.)) and c.applied == (0., .25)
c.set_target((0., 0.), (0., -.05)) # return to zero offset is a real target
assert c.pending == ((0., 0.), (0., -.05))
c.advance(uploads.append, lambda: None, lambda: None)
assert c.applied == (0., 0.) and c.pending is None
@pytest.mark.parametrize("stage", ("upload", "step", "synchronize"))
def test_failure_after_upload_attempt_is_fail_closed(stage):
c = circle()
c.set_target((0., .25), (0., .05))
def call(name):
def run(*_):
if name == stage:
raise RuntimeError(f"{name} failed")
return run
with pytest.raises(RuntimeError, match=f"{stage} failed"):
c.advance(call("upload"), call("step"), call("synchronize"))
assert c.failed and c.applied == (0., 0.)
for later in (c.record, lambda: c.set_target((0., .1), (0., 0.)),
lambda: c.advance(*(lambda *_: None,) * 3)):
with pytest.raises(RuntimeError, match="unknown"):
later()
def test_band_covers_offset_box_with_margin():
assert circle().band() == (50 - 2 - 10 - 2, 50 + 2 + 10 + 2, 40 - 3 - 10 - 2, 40 + 3 + 10 + 2)
@pytest.mark.parametrize("radius", (R_MIN, 4.3, 7.5, 12.0))
def test_refill_donor_sets_nonempty_within_bounds(radius):
# Evidence for R >= R_MIN and |dc|_2 <= MAX_STEP on a fixed sample, not a proof.
assert enumerate_donor_coverage(radius) == 0
pytest.importorskip("pycuda.driver")
from CelerisLab import simulation as module # noqa: E402
from tests.unit.test_simulation_lifecycle import make_simulation # noqa: E402
def moving_simulation(monkeypatch, *, moving=True):
sim = make_simulation(monkeypatch)
obj = SimpleNamespace(motion=circle() if moving else None, obj_id=0,
center=(50., 40.), radius=10.)
writes = []
sim._failed, sim._initialized = None, True
sim.bodies = SimpleNamespace(objects=[obj], count=1, jets=[],
has_pending_edit=lambda: False,
_upload_action_async=lambda s: writes.append("action"),
zero_obs_async=lambda s: writes.append("obs"))
sim.stepper = SimpleNamespace(step=lambda *a, **k: writes.append("step"))
return sim, writes
def test_positive_steps_rejected_before_device_writes(monkeypatch):
sim, writes = moving_simulation(monkeypatch)
sim.bodies.objects[0].motion.set_target((0., .1), (0., .01))
for call in (lambda: sim.run(1), lambda: sim.run(5), lambda: sim.step(),
lambda: sim.run(3, checkpoint_interval=1)):
with pytest.raises(RuntimeError, match="moving kernels"):
call()
with pytest.raises(RuntimeError, match="moving kernels"): # nominal at rest as well
sim.bodies.objects[0].motion.pending = None
sim.run(1)
assert writes == [] and sim.bodies.objects[0].motion.applied == (0., 0.)
def test_moving_body_rejects_topology_and_checkpoint(monkeypatch):
sim, _ = moving_simulation(monkeypatch)
for call in (sim.sync_bodies, sim.save_checkpoint, lambda: sim.load_checkpoint("x.h5")):
with pytest.raises(ValueError, match="translating circles"):
call()