Compare commits

...
Author SHA1 Message Date
Astra (AI agent) d0968e6e91 research: reuse quadrature and clarify snapshot and VIV coordinates 2026-10-31 12:04:56 +08:00
Frank14f 4f17c3c586 motion: exercise real step ownership and move host test helpers out of src 2026-10-31 11:24:09 +08:00
3f7f303e17 research: VIV driver reads diagnostics by MovingTables.DIAG name; coupling-contract README
Index values are unchanged (1, 3, 7, 10, 13); the run recorded at fb507aa1
used the same columns. Research-only.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
2026-10-31 11:20:57 +08:00
ClaudeandClaude Opus 5.5 d5b13d7f31 motion: mark MovingCircle.advance as a test helper, not the production path
Docstring only.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
2026-10-31 01:24:35 +00:00
ClaudeandClaude Opus 5.5 353a7b03c4 stepper: bind moving-circle kernels only for D2Q9 double buffer
Both kernels are compiled under DIM == 2, but LBMStepper requested them for
every double-buffer module, so stationary D3Q19 double-buffer construction
would fail at get_function. Moving circles already reject D3Q19 and EsoPull
at initialize. Stale latch comment updated. Parametrized CPU test with a
module double that lacks the symbols.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
2026-10-31 01:23:54 +00:00
ClaudeandClaude Opus 5.5 0d6df6e6b8 simulation: consume a zero-motion translation target; empty diagnostics shape
A pending target equal to the applied offset uploaded the same record, so
_run_moving saw no advance and never committed it: pending stayed latched,
blocking run(n > 1) and direct steps. Commit whenever a target was pending
(no_donor check and hold-record re-upload still only on advancing steps).
moving_diagnostics() without moving circles now returns (0, DIAG_FLOATS)
instead of the stale (0, 6). Tests for both.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
2026-10-31 01:23:04 +00:00
11 changed files with 243 additions and 112 deletions
+6
View File
@@ -394,6 +394,12 @@ before `run()`, the force will be reset to zero. For the common usage pattern
| `sim.load_checkpoint(path)` | Restore from HDF5 (config must match) | | `sim.load_checkpoint(path)` | Restore from HDF5 (config must match) |
| `sim.snapshot()` / `sim.restore()` | In-memory field snapshot | | `sim.snapshot()` / `sim.restore()` | In-memory field snapshot |
In-memory snapshots save populations only. Restore keeps the current clock,
controls, observations, geometry and flags; use it only while geometry is
unchanged. This is not a restart or an FSI rollback. In particular, do not
restore a snapshot taken before a translating circle changed position.
This geometry restriction is a caller precondition, not a runtime check.
#### Low-level access #### Low-level access
| Attribute | Description | | Attribute | Description |
+4
View File
@@ -30,6 +30,10 @@ cells (M1 `MovingCircleCells`), moving-wall Bouzidi with wall velocity
mass, eq_only, no_donor, crossings) accumulated since obs were zeroed; mass, eq_only, no_donor, crossings) accumulated since obs were zeroed;
`obs`/`read_force` keep their contract (CME on the body). `obs`/`read_force` keep their contract (CME on the body).
- `get_flags()` reads device flags back when moving circles exist. - `get_flags()` reads device flags back when moving circles exist.
- `snapshot()` / `restore()` save and restore populations only, keeping the
current clock, controls, observations, geometry and flags. Use them only
with unchanged geometry; restoring across a translation is unsupported
and is not rejected at runtime. They do not provide FSI rollback.
*Supported:* dense double buffer, D2Q9, FP32 unshifted, SRT/MRT, circles, *Supported:* dense double buffer, D2Q9, FP32 unshifted, SRT/MRT, circles,
spin. *Rejected before mutation:* EsoPull, ddf shifting, D3Q19, refinement / spin. *Rejected before mutation:* EsoPull, ddf shifting, D3Q19, refinement /
+11 -4
View File
@@ -29,10 +29,17 @@ Claude, Grok and Astra at 1abfe954; implemented in research/viv_re150/.
- D = 40, MRT, dense double buffer; U_lat = .04, nu = U D / 150 = .010667, - D = 40, MRT, dense double buffer; U_lat = .04, nu = U D / 150 = .010667,
tau = .532. tau = .532.
- Box [0, 75 D] × [0, 60 D] (3000 × 2400) with the cylinder centre at - Nominal 75 D × 60 D lattice allocation (3000 × 2400), based on the
(30 D, 30 D): 30 D from the inlet, 45 D from the outlet, centred in y Kolahdouz §4.1.1 layout. Coordinate clarification of the existing run:
(Kolahdouz §4.1.1). Qin's own box ([0, 40 D] × [0, 20 D], centre (10 D, 10 D)) the driver uses `(30 D, (ny - 1)/2)` = `(1200, 1199.5)` in lattice
is the target paper's, not ours. coordinates, or `(30 D, 29.9875 D)`. The earlier `(30 D, 30 D)` wording
described the nominal box, not the exact executed centre.
Free-slip reflection planes are at `y = .5` and `y = ny - 1.5 = 2398.5`;
their separation is 2398 cells = 59.95 D, with the centre exactly midway,
29.975 D from each wall. Node counts are not physical wall-to-wall lengths.
This corrects the description, without changing the historical geometry
or results. Qin's own box ([0, 40 D] × [0, 20 D], centre (10 D, 10 D)) is
the target paper's, not ours.
- Boundaries: uniform Zou–He inlet, NEQ outlet (rho = 1, backflow clamp), - Boundaries: uniform Zou–He inlet, NEQ outlet (rho = 1, backflow clamp),
free-slip y walls. Kolahdouz uses a traction outlet: declared mismatch. free-slip y walls. Kolahdouz uses a traction outlet: declared mismatch.
- Qin's own VIV box (40 D × 20 D) is not used; a smaller box would be a - Qin's own VIV box (40 D × 20 D) is not used; a smaller box would be a
+8 -1
View File
@@ -8,6 +8,7 @@ s/T in [0.25 + (k-1), 0.25 + k); cycle 0 is the startup quarter.
from __future__ import annotations from __future__ import annotations
import math import math
from functools import lru_cache
import numpy as np import numpy as np
@@ -53,10 +54,16 @@ def bin_phase(s, T):
return None return None
@lru_cache(maxsize=8)
def _quadrature_rule(points):
"""Reuse the nodes and weights, which depend only on quadrature order."""
return np.polynomial.legendre.leggauss(points)
def interval_bases(n, T, points=16): def interval_bases(n, T, points=16):
"""Single-step interval averages over [n, n + 1] of cos phi |cos phi| and sin phi """Single-step interval averages over [n, n + 1] of cos phi |cos phi| and sin phi
(phi in radians = 2 pi s / T - pi / 2), Gauss-Legendre quadrature.""" (phi in radians = 2 pi s / T - pi / 2), Gauss-Legendre quadrature."""
x, w = np.polynomial.legendre.leggauss(points) x, w = _quadrature_rule(points)
s = int(n) % T + 0.5 * (x + 1.0) # periodic: identical bases in every cycle s = int(n) % T + 0.5 * (x + 1.0) # periodic: identical bases in every cycle
phi = 2.0 * np.pi * s / T - 0.5 * np.pi phi = 2.0 * np.pi * s / T - 0.5 * np.pi
c = np.cos(phi) c = np.cos(phi)
+23
View File
@@ -0,0 +1,23 @@
# Research VIV driver — coupling contract (not production)
`run.py` couples a 1-DOF undamped spring–mass (`viv.py`) to the production
moving-circle API. Everything specific to fluid–structure coupling lives here,
not in `src/`.
- Per step: known (y_n, v_n) and the last completed force F_last → exact
harmonic candidate (y_{n+1}, v_{n+1}) → `set_body(offset=(0, y_{n+1}),
vy=y_{n+1} − y_n)` (interval wall velocity, lattice dt = 1) → `run(1)` →
commit y, v and F_last only after the step succeeded. One step of explicit
force lag (weak coupling); no sub-iteration.
- Structural RHS = CME + Wen (y), read from `moving_diagnostics()` by the
names in `MovingTables.DIAG`; raw CME, the `obs` force and node terms are
logged as diagnostics. This is a stage choice, not a validated force law.
- Production guarantees used: fail-closed target commit (synchronized step,
`no_donor == 0`), input rejection before mutation, a failed device step
marks the Simulation failed.
- Not provided, by design: device-resident structure state or batched
(multi-step) coupling, strong/sub-iterated coupling (would need state
rollback), more than one DOF. These are future work and would be a separate
execution path, not an abstraction over this driver.
- Recorded result and limits: `docs/validation/moving_bodies_stage1.md`;
predeclaration: `docs/validation/viv_re150_predeclare.md`.
+9 -7
View File
@@ -1,8 +1,9 @@
"""One-shot predeclared VIV run (docs/validation/viv_re150_predeclare.md). """One-shot predeclared VIV run (docs/validation/viv_re150_predeclare.md).
Re 150, U* = 5, m/(rho D^2) = 2, zeta = 0, D = 40 on a 75 D x 60 D native Re 150, U* = 5, m/(rho D^2) = 2, zeta = 0, D = 40 on a nominal 75 D x 60 D
channel (uniform inlet, NEQ outlet, free-slip y walls), cylinder at lattice (uniform inlet, NEQ outlet, free-slip y walls). The cylinder is at
(30 D, 30 D). Structural RHS = CME + Wen (diagnostics: raw CME, node terms). (30 D, (ny - 1)/2); y-wall planes are .5 and ny - 1.5 in lattice coordinates.
Structural RHS = CME + Wen (diagnostics: raw CME, node terms).
Research driver only: Simulation API, no production change. Research driver only: Simulation API, no production change.
PYTHONPATH=src python research/viv_re150/run.py --out /path/outside/checkout PYTHONPATH=src python research/viv_re150/run.py --out /path/outside/checkout
@@ -65,7 +66,8 @@ def main():
step_cap = args.steps step_cap = args.steps
centre = (30.0 * Dn, 0.5 * (ny - 1)) centre = (30.0 * Dn, 0.5 * (ny - 1))
from CelerisLab.body.motion import Translation from CelerisLab.body.motion import MovingTables, Translation
col = {k: MovingTables.DIAG.index(k) for k in ("cme_y", "wen_y", "cover_jy", "uncover_jy", "no_donor")}
from CelerisLab.simulation import Simulation from CelerisLab.simulation import Simulation
with tempfile.NamedTemporaryFile("w", suffix=".json", delete=False) as fh: with tempfile.NamedTemporaryFile("w", suffix=".json", delete=False) as fh:
json.dump(cfg, fh) json.dump(cfg, fh)
@@ -95,9 +97,9 @@ def main():
sim.run(1) # device failures raise and mark the Simulation failed sim.run(1) # device failures raise and mark the Simulation failed
d = sim.moving_diagnostics()[0] d = sim.moving_diagnostics()[0]
f_obs = sim.read_force(body)[1] f_obs = sim.read_force(body)[1]
f_rhs = float(d[1] + d[3]) # CME + Wen (y) f_rhs = float(d[col["cme_y"]] + d[col["wen_y"]]) # CME + Wen (y)
row = (n, y1, v1, f_rhs, float(d[1]), float(d[3]), float(f_obs), row = (n, y1, v1, f_rhs, float(d[col["cme_y"]]), float(d[col["wen_y"]]), float(f_obs),
float(d[7] - d[10]), float(d[13])) float(d[col["cover_jy"]] - d[col["uncover_jy"]]), float(d[col["no_donor"]]))
if not viv.finite(row): # every recorded output, not only the RHS if not viv.finite(row): # every recorded output, not only the RHS
status = "fail_nonfinite_output" status = "fail_nonfinite_output"
break break
+3 -61
View File
@@ -1,10 +1,10 @@
# CelerisLab/body/motion.py # CelerisLab/body/motion.py
"""Host state of a translating circle (dense double buffer, D2Q9 only). """Translation targets and device tables (dense double buffer, D2Q9 only).
Contract: docs/moving_bodies.md, sections 4a, 5 and 7. The applied geometry Contract: docs/moving_bodies.md, sections 4a, 5 and 7. The applied geometry
(last successfully committed offset) and an optional pending target are two (last successfully committed offset) and an optional pending target are two
separate facts; centres are always ``base + offset``. Nothing here touches separate facts; centres are always ``base + offset``. MovingCircle owns host
the device: the caller supplies upload/step/synchronize callables. targets, MovingTables owns device records, and Simulation commits completed steps.
""" """
from __future__ import annotations from __future__ import annotations
@@ -18,10 +18,6 @@ R_MIN = 4.0 # smallest radius covered by the donor enumeration
Vec = Tuple[float, float] 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: def _vec(value, name: str) -> Vec:
x, y = (float(v) for v in value) x, y = (float(v) for v in value)
@@ -98,23 +94,6 @@ class MovingCircle:
offset, velocity = self.applied, (0.0, 0.0) offset, velocity = self.applied, (0.0, 0.0)
return self.centre(self.applied), self.centre(offset), velocity 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, never reports applied state again,
and the owner's ``on_failure`` hook is called.
"""
record = self.record()
try:
upload(record)
step()
synchronize()
except BaseException as exc:
self.fail(exc)
raise
self.commit()
def commit(self) -> None: def commit(self) -> None:
"""Applied := pending offset (caller has synchronized a successful step).""" """Applied := pending offset (caller has synchronized a successful step)."""
if self.pending is not None: if self.pending is not None:
@@ -138,43 +117,6 @@ class MovingCircle:
raise RuntimeError("moving circle state is unknown after a failed step") 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
class MovingTables: class MovingTables:
"""Device tables of the translating circles for ``MovingCircleCells``/``Links``. """Device tables of the translating circles for ``MovingCircleCells``/``Links``.
+5 -3
View File
@@ -21,7 +21,7 @@ class LBMStepper:
self._esopull = (cfg.streaming == "esopull") self._esopull = (cfg.streaming == "esopull")
# Set by the owner when a step would be wrong (translating circles # Set by the owner when a step would be wrong (translating circles
# excluded from the packed links before the moving kernels exist). # without attached moving tables, a pending target, or a failed step).
self.blocked: str | None = None self.blocked: str | None = None
self.moving = None # MovingTables attached by the owner (translating circles) self.moving = None # MovingTables attached by the owner (translating circles)
@@ -39,8 +39,10 @@ class LBMStepper:
self.apply_esopull_curved_fn = module.get_function("ApplyEsoPullCurvedKernel") self.apply_esopull_curved_fn = module.get_function("ApplyEsoPullCurvedKernel")
self.sensor_fn = module.get_function("SensorKernel") self.sensor_fn = module.get_function("SensorKernel")
self.force_region_fn = module.get_function("ForceRegionKernel") self.force_region_fn = module.get_function("ForceRegionKernel")
self.moving_fn = None if self._esopull else module.get_function("MovingCircleLinks") # Moving-circle kernels exist only in D2Q9 double-buffer modules.
self.moving_cells_fn = None if self._esopull else module.get_function("MovingCircleCells") moving = not self._esopull and cfg.is_d2q9
self.moving_fn = module.get_function("MovingCircleLinks") if moving else None
self.moving_cells_fn = module.get_function("MovingCircleCells") if moving else None
tpb = cfg.threads_per_block tpb = cfg.threads_per_block
self.block = (tpb, 1, 1) self.block = (tpb, 1, 1)
+7 -3
View File
@@ -242,7 +242,8 @@ class Simulation:
self._require_usable() self._require_usable()
tables = self.stepper.moving tables = self.stepper.moving
if tables is None: if tables is None:
return np.zeros((0, 6), np.float32) from .body.motion import MovingTables
return np.zeros((0, MovingTables.DIAG_FLOATS), np.float32)
self.stream.synchronize() self.stream.synchronize()
return tables.download_diagnostics() return tables.download_diagnostics()
@@ -785,6 +786,7 @@ class Simulation:
(docs/moving_bodies.md 4a, 5). Any failure leaves the Simulation unusable.""" (docs/moving_bodies.md 4a, 5). Any failure leaves the Simulation unusable."""
tables = self.stepper.moving tables = self.stepper.moving
circles = [(obj.motion, obj.obj_id) for obj in self._moving_bodies()] circles = [(obj.motion, obj.obj_id) for obj in self._moving_bodies()]
pending = any(motion.pending is not None for motion, _ in circles)
try: try:
advancing = tables.upload_records(circles, stream) advancing = tables.upload_records(circles, stream)
self.stepper._step(steps, action_gpu=self.bodies.action_gpu, self.stepper._step(steps, action_gpu=self.bodies.action_gpu,
@@ -797,10 +799,12 @@ class Simulation:
no_donor = tables.download_diagnostics()[:, tables.DIAG.index("no_donor")] no_donor = tables.download_diagnostics()[:, tables.DIAG.index("no_donor")]
if np.any(no_donor > 0): if np.any(no_donor > 0):
raise RuntimeError("translating circle refill found no donor cell") raise RuntimeError("translating circle refill found no donor cell")
if pending: # a target equal to the applied offset is consumed as a hold step
for motion, _ in circles: for motion, _ in circles:
motion.commit() motion.commit()
tables.upload_records(circles, stream) # first hold: old = new = applied, u = 0 if advancing:
stream.synchronize() tables.upload_records(circles, stream) # first hold: old = new = applied, u = 0
stream.synchronize()
self._attach_moving() self._attach_moving()
except BaseException as exc: except BaseException as exc:
for motion, _ in circles: # device geometry/state unknown after any failure for motion, _ in circles: # device geometry/state unknown after any failure
+29
View File
@@ -82,6 +82,35 @@ def test_interval_bases_and_morison_fit_recover_known_coefficients():
assert p.morison_fit(f, b[:, 0], b[:, 1]) == pytest.approx((2.09, 1.45), abs=1e-9) assert p.morison_fit(f, b[:, 0], b[:, 1]) == pytest.approx((2.09, 1.45), abs=1e-9)
def test_quadrature_rule_reused_across_steps_and_periods(monkeypatch):
original = np.polynomial.legendre.leggauss
orders = []
def counted(points):
orders.append(points)
return original(points)
p._quadrature_rule.cache_clear()
monkeypatch.setattr(np.polynomial.legendre, "leggauss", counted)
try:
for points in (16, 8, 16):
for period in (T, 2 * T):
for n in (0, 123, period + 123):
p.interval_bases(n, period, points)
assert orders == [16, 8]
finally:
p._quadrature_rule.cache_clear()
def test_interval_sine_basis_matches_exact_integral():
for period in (T, 2 * T):
half_width = math.pi / period
for n in (0, period // 4, period // 2, period - 1, 7 * period + 123):
midpoint = 2 * math.pi * ((n % period) + .5) / period - math.pi / 2
expected = math.sin(midpoint) * math.sin(half_width) / half_width
assert p.interval_bases(n, period)[1] == pytest.approx(expected, rel=0, abs=2e-15)
def test_geometry_repeats_exactly_across_cycles(): def test_geometry_repeats_exactly_across_cycles():
A = p.amplitude(20) A = p.amplitude(20)
for m in (0, 1, T // 4, T // 2 - 1, 7 * T // 12): for m in (0, 1, T // 4, T // 2 - 1, 7 * T // 12):
+138 -33
View File
@@ -1,14 +1,13 @@
"""Translating-circle host contract, S1 (docs/moving_bodies.md 4a, 5, 7). """Translating-circle host contract (docs/moving_bodies.md 4a, 5, 7).
Host doubles only: no kernels exist yet, so every positive-step run with a Real Simulation entrypoints with device doubles; no CUDA context or GPU work.
moving body must be rejected before any device write.
""" """
import math
from types import SimpleNamespace from types import SimpleNamespace
import pytest import pytest
from CelerisLab.body.motion import (MAX_STEP, R_MIN, MovingCircle, Translation, from CelerisLab.body.motion import MAX_STEP, R_MIN, MovingCircle, Translation
enumerate_donor_coverage)
BAD = (float("nan"), float("inf")) BAD = (float("nan"), float("inf"))
@@ -46,46 +45,96 @@ def test_targets_validate_against_applied_and_last_wins():
assert c.applied == (0., 0.) assert c.applied == (0., 0.)
def test_record_commit_and_hold(): def test_record_commit_and_hold(monkeypatch):
c = circle() sim, writes, c = routed_advancing_simulation(monkeypatch)
assert c.record() == ((50., 40.), (50., 40.), (0., 0.)) # hold: no translation velocity assert c.record() == ((50., 40.), (50., 40.), (0., 0.))
c.set_target((0., .25), (0., .05)) sim.set_body(0, offset=(0., .25), vy=.05)
assert c.record() == ((50., 40.), (50., 40.25), (0., .05)) assert c.record() == ((50., 40.), (50., 40.25), (0., .05))
uploads = [] seen_at_sync = []
c.advance(uploads.append, lambda: None, lambda: None) sim.run(1, stream=SimpleNamespace(synchronize=lambda: seen_at_sync.append(c.applied)))
assert uploads == [((50., 40.), (50., 40.25), (0., .05))] assert seen_at_sync == [(0., 0.), (0., .25)] # commit between step sync and hold sync
uploads = [w[1][0] for w in writes if isinstance(w, tuple) and w[0] == "upload"]
assert uploads == [((50., 40.), (50., 40.25), (0., .05)),
((50., 40.25), (50., 40.25), (0., 0.))]
assert c.applied == (0., .25) and c.pending is None # offset-only commit 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 sim.run(1, stream=STREAM) # no new target: hold step
uploads = [w[1][0] for w in writes if isinstance(w, tuple) and w[0] == "upload"]
assert uploads[-1] == ((50., 40.25), (50., 40.25), (0., 0.)) and c.applied == (0., .25) 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 sim.set_body(0, offset=(0., 0.), vy=-.05) # zero offset is a real target
assert c.pending == ((0., 0.), (0., -.05)) assert c.pending == ((0., 0.), (0., -.05))
c.advance(uploads.append, lambda: None, lambda: None) sim.run(1, stream=STREAM)
assert c.applied == (0., 0.) and c.pending is None assert c.applied == (0., 0.) and c.pending is None
@pytest.mark.parametrize("stage", ("upload", "step", "synchronize")) @pytest.mark.parametrize("stage", ("upload", "step", "synchronize"))
def test_failure_after_upload_attempt_is_fail_closed(stage): def test_failure_after_upload_attempt_is_fail_closed(monkeypatch, stage):
c = circle() sim, writes, c = routed_advancing_simulation(monkeypatch, fail="upload" if stage == "upload" else None)
c.set_target((0., .25), (0., .05)) sim.set_body(0, offset=(0., .25), vy=.05)
def call(name): def fail(*a, **k):
def run(*_): raise RuntimeError(f"{stage} failed")
if name == stage: if stage == "step":
raise RuntimeError(f"{name} failed") sim.stepper._step = fail
return run stream = SimpleNamespace(synchronize=fail) if stage == "synchronize" else STREAM
with pytest.raises(RuntimeError, match=f"{stage} failed"): with pytest.raises(RuntimeError, match=f"{stage} failed"):
c.advance(call("upload"), call("step"), call("synchronize")) sim.run(1, stream=stream)
assert c.failed and c.applied == (0., 0.) assert c.failed and c.applied == (0., 0.) and c.pending == ((0., .25), (0., .05))
for later in (c.record, lambda: c.set_target((0., .1), (0., 0.)), assert sim._failed is not None and sim.stepper.blocked == sim._failed
lambda: c.advance(*(lambda *_: None,) * 3)): for later in (c.record, lambda: sim.set_body(0, offset=(0., .1), vy=0.)):
with pytest.raises(RuntimeError, match="unknown"): with pytest.raises(RuntimeError, match="unknown"):
later() later()
with pytest.raises(RuntimeError, match="Simulation failed"):
sim.run(1, stream=STREAM)
with pytest.raises(RuntimeError, match="moving circle step failed"):
sim.stepper.step(1, action_gpu=None, obs_gpu=None)
def test_band_covers_offset_box_with_margin(): 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) assert circle().band() == (50 - 2 - 10 - 2, 50 + 2 + 10 + 2, 40 - 3 - 10 - 2, 40 + 3 + 10 + 2)
# 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 _inside(x: int, y: int, c, r: float) -> bool:
return (x - c[0]) ** 2 + (y - c[1]) ** 2 < r * r
def empty_donor_cells(c_old, c_new, 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
@pytest.mark.parametrize("radius", (R_MIN, 4.3, 7.5, 12.0)) @pytest.mark.parametrize("radius", (R_MIN, 4.3, 7.5, 12.0))
def test_refill_donor_sets_nonempty_within_bounds(radius): 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. # Evidence for R >= R_MIN and |dc|_2 <= MAX_STEP on a fixed sample, not a proof.
@@ -278,6 +327,26 @@ def advancing_simulation(monkeypatch, **tables):
return sim, writes, sim.bodies.objects[0].motion return sim, writes, sim.bodies.objects[0].motion
def routed_advancing_simulation(monkeypatch, **tables):
"""Real body registration, target routing, run and owner hook; double device I/O."""
sim = routed_simulation(monkeypatch)
bid = sim.add_body("circle", (50., 40.), 10., translation=Translation((-2., 2., -3., 3.), .1))
sim._initialized = True
sim._assert_runtime_contracts = lambda: None
writes = []
monkeypatch.setattr(sim.bodies, "_upload_action_async", lambda s: writes.append("action"))
monkeypatch.setattr(sim.bodies, "zero_obs_async", lambda s: writes.append("obs"))
monkeypatch.setattr(sim.bodies, "download_obs_full_async", lambda s: writes.append("download"))
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)
sim._moving_tables = sim.stepper.moving
sim._attach_moving()
return sim, writes, sim.bodies.get(bid).motion
def test_run_one_advances_and_commits_after_sync(monkeypatch): def test_run_one_advances_and_commits_after_sync(monkeypatch):
sim, writes, motion = advancing_simulation(monkeypatch) sim, writes, motion = advancing_simulation(monkeypatch)
motion.set_target((0., .25), (0., .05)) motion.set_target((0., .25), (0., .05))
@@ -297,6 +366,44 @@ def test_run_one_advances_and_commits_after_sync(monkeypatch):
assert ("step", 3, False) in writes and motion.applied == (0., .25) assert ("step", 3, False) in writes and motion.applied == (0., .25)
def test_zero_motion_target_is_consumed_by_run_one(monkeypatch):
sim, writes, motion = advancing_simulation(monkeypatch)
motion.set_target((0., 0.), (0., 0.)) # equal to applied: a hold target
sim._attach_moving()
assert sim.stepper.blocked == module.PENDING_BLOCKED
sim.run(1, stream=STREAM)
assert motion.pending is None and motion.applied == (0., 0.)
assert sim.stepper.blocked is None and ("step", 1, False) in writes
writes.clear()
sim.run(3, stream=STREAM) # no stale latch
assert ("step", 3, False) in writes
@pytest.mark.parametrize("streaming,d2q9,bound", [("double_buffer", True, True),
("double_buffer", False, False),
("esopull", True, False)])
def test_moving_kernels_bound_only_for_d2q9_double_buffer(monkeypatch, streaming, d2q9, bound):
from CelerisLab.lbm import stepper as stepper_module
monkeypatch.setattr(stepper_module.cuda, "mem_alloc", lambda n: object())
moving = {"MovingCircleLinks", "MovingCircleCells"}
def get_function(name):
if name in moving and not (d2q9 and streaming == "double_buffer"):
raise KeyError(name) # absent from D3Q19 / EsoPull modules
return lambda *a, **k: None
cfg = SimpleNamespace(streaming=streaming, threads_per_block=128, is_d2q9=d2q9,
nx=16, ny=8, nz=4, nq=9 if d2q9 else 19)
stepper = stepper_module.LBMStepper(SimpleNamespace(), SimpleNamespace(get_function=get_function), cfg)
assert (stepper.moving_fn is not None) == bound and (stepper.moving_cells_fn is not None) == bound
def test_moving_diagnostics_empty_shape_matches_layout(monkeypatch):
from CelerisLab.body.motion import MovingTables
sim = make_simulation(monkeypatch)
sim._failed, sim.stepper = None, SimpleNamespace(moving=None)
assert sim.moving_diagnostics().shape == (0, MovingTables.DIAG_FLOATS)
@pytest.mark.parametrize("failure", ("upload", "first_sync", "no_donor", "hold_upload", @pytest.mark.parametrize("failure", ("upload", "first_sync", "no_donor", "hold_upload",
"hold_sync")) "hold_sync"))
def test_failed_advance_fails_simulation(monkeypatch, failure): def test_failed_advance_fails_simulation(monkeypatch, failure):
@@ -507,23 +614,21 @@ def test_moving_circle_excluded_from_static_links():
def test_failed_moving_step_marks_simulation_failed(monkeypatch): def test_failed_moving_step_marks_simulation_failed(monkeypatch):
sim = routed_simulation(monkeypatch) sim, writes, motion = routed_advancing_simulation(monkeypatch)
sim.ctx = ClosableOwner() sim.ctx = ClosableOwner()
bid = sim.add_body("circle", (50., 40.), 10., translation=Translation((0., 0., -1., 1.), .1)) sim.set_body(0, offset=(0., .2), vy=.05)
motion = sim.bodies.get(bid).motion
motion.set_target((0., .2), (0., .05))
def fail(*_): def fail(*_):
raise RuntimeError("sync failed") raise RuntimeError("sync failed")
with pytest.raises(RuntimeError, match="sync failed"): with pytest.raises(RuntimeError, match="sync failed"):
motion.advance(lambda *_: None, lambda: None, fail) sim.run(1, stream=SimpleNamespace(synchronize=fail))
assert sim._failed is not None and motion.applied == (0., 0.) assert sim._failed is not None and motion.applied == (0., 0.)
for call in (lambda: sim.run(0), lambda: sim.run(1), sim.step, sim.sync_bodies, for call in (lambda: sim.run(0), lambda: sim.run(1), sim.step, sim.sync_bodies,
sim.save_checkpoint, lambda: sim.load_checkpoint("x.h5")): sim.save_checkpoint, lambda: sim.load_checkpoint("x.h5")):
with pytest.raises(RuntimeError, match="Simulation failed"): with pytest.raises(RuntimeError, match="Simulation failed"):
call() call()
with pytest.raises(RuntimeError, match="unknown"): # no retry on the failed circle with pytest.raises(RuntimeError, match="unknown"): # no retry on the failed circle
motion.advance(lambda *_: None, lambda: None, lambda: None) sim.set_body(0, offset=(0., .1), vy=.05)
with sim.ctx.activate(): # a refused close keeps the failed state and resources with sim.ctx.activate(): # a refused close keeps the failed state and resources
with pytest.raises(RuntimeError, match="active CUDA"): with pytest.raises(RuntimeError, match="active CUDA"):
sim.close() sim.close()