diff --git a/src/CelerisLab/body/manager.py b/src/CelerisLab/body/manager.py index e1ab25a..188f6da 100644 --- a/src/CelerisLab/body/manager.py +++ b/src/CelerisLab/body/manager.py @@ -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 diff --git a/src/CelerisLab/body/motion.py b/src/CelerisLab/body/motion.py new file mode 100644 index 0000000..32ae02d --- /dev/null +++ b/src/CelerisLab/body/motion.py @@ -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 diff --git a/src/CelerisLab/body/objects.py b/src/CelerisLab/body/objects.py index 19bb74c..75a6bb0 100644 --- a/src/CelerisLab/body/objects.py +++ b/src/CelerisLab/body/objects.py @@ -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): diff --git a/src/CelerisLab/simulation.py b/src/CelerisLab/simulation.py index a9ad027..7c74be2 100644 --- a/src/CelerisLab/simulation.py +++ b/src/CelerisLab/simulation.py @@ -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, diff --git a/tests/unit/test_moving_bodies.py b/tests/unit/test_moving_bodies.py new file mode 100644 index 0000000..e5261ca --- /dev/null +++ b/tests/unit/test_moving_bodies.py @@ -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()