lbm: translating circles move with refill (S3a)
MovingCircleCells reclassifies swept-band cells from immutable old/new FP64 predicates: covered cells add their raw momentum and mass to diagnostics and become solid; uncovered cells are refilled in all nine slots from donors fluid in both geometries (mean density, NEQ from the best-aligned donor, wall velocity at the cell), counting eq_only and no_donor. MovingCircleLinks adds the translation velocity to the wall velocity. Simulation.run(1) uploads the motion record, runs M1 then M2, synchronizes, rejects a nonzero no_donor count and only then commits the target and uploads the hold record; any failure marks the Simulation failed. Targets block direct stepper calls and multi-step runs. get_flags reads device flags back when translating circles exist. Contract section 10b. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This commit is contained in:
co-authored by
Claude Opus 5.5
parent
b920a303a4
commit
d50b8ba1be
@@ -363,6 +363,34 @@ recorded as not runnable rather than replaced by another test.
|
||||
- Not validation: a pass of the zero-motion gates says nothing about moving
|
||||
geometry. GPU parity (G0a/G0b) is unverified until run on a free device.
|
||||
|
||||
## 10b. S3a — translation with refill (implemented on `claude/moving-bodies-s3a`)
|
||||
|
||||
- Targets are enabled again (4a). `set_body(offset, vx, vy)` validates against
|
||||
`applied`; a pending target blocks direct `LBMStepper.step(n > 0)` before any
|
||||
launch; `run(n > 1)` with a target raises before mutation; `run(0)` never
|
||||
consumes it; `run(1)` advances it.
|
||||
- Per advancing step, on one stream: upload the immutable FP64 records
|
||||
`{c_old, c_new, r, u}` (only when changed) → M1 `MovingCircleCells` → M2
|
||||
`MovingCircleLinks` (`u_w = u + omega x r` at the hit point) → static curved →
|
||||
`OneStep` → swap → force region → sensors → synchronize → read `no_donor` →
|
||||
commit (`applied := pending.offset`) → upload the hold record
|
||||
(`c_old = c_new = applied`, `u = 0`). Hold steps launch no M1 and are the
|
||||
S2a path. Any failure from the upload on marks the circles and the
|
||||
Simulation failed; nothing is committed.
|
||||
- M1 classifies with the old/new predicates only (no flag reads). Covered
|
||||
cells add raw `sum c G` / `sum G` of `G^{n-1}` to `cover_*` and become
|
||||
`SOLID|OBSTACLE|BC_CURVED`. Uncovered cells: donors are fluid in both
|
||||
geometries; `rho_bar` over all of them; NEQ from the donor with the largest
|
||||
unnormalized `c·n > 0` (first direction wins ties, `n` from `c_new`);
|
||||
`u_w` (translation + spin) at the cell position; all nine slots written;
|
||||
raw `sum c G`, `sum G` of the refill go to `uncover_*`; no donor →
|
||||
`f_eq(RHO, u_w)` plus `no_donor` (hard failure before commit); no NEQ donor
|
||||
→ `eq_only`.
|
||||
- `get_flags()` reads the device flags back (synchronized) when moving tables
|
||||
exist; otherwise it returns the host array as before.
|
||||
- Not selected: the structural force. CME+node vs CME+Wen stay OPEN until X1
|
||||
and G4 have run.
|
||||
|
||||
## 11. Remaining decisions
|
||||
|
||||
Structural RHS (section 6); motion-record layout; refill details; coupling
|
||||
|
||||
@@ -111,14 +111,21 @@ class MovingCircle:
|
||||
step()
|
||||
synchronize()
|
||||
except BaseException as exc:
|
||||
self.failed = True
|
||||
if self._on_failure is not None:
|
||||
self._on_failure(exc)
|
||||
self.fail(exc)
|
||||
raise
|
||||
self.commit()
|
||||
|
||||
def commit(self) -> None:
|
||||
"""Applied := pending offset (caller has synchronized a successful step)."""
|
||||
if self.pending is not None:
|
||||
self.applied = self.pending[0]
|
||||
self.pending = None
|
||||
|
||||
def fail(self, exc: BaseException) -> None:
|
||||
self.failed = True
|
||||
if self._on_failure is not None:
|
||||
self._on_failure(exc)
|
||||
|
||||
def band(self) -> Tuple[int, int, int, int]:
|
||||
"""Inclusive cell AABB (x0, x1, y0, y1) swept over the whole offset box."""
|
||||
x0, x1, y0, y1 = self.translation.offset_box
|
||||
@@ -169,19 +176,25 @@ def enumerate_donor_coverage(radius: float, positions: int = 8, directions: int
|
||||
|
||||
|
||||
class MovingTables:
|
||||
"""Device tables of the translating circles for ``MovingCircleLinks`` (fixed centre).
|
||||
"""Device tables of the translating circles for ``MovingCircleCells``/``Links``.
|
||||
|
||||
Immutable per run: band table (x0, y0, width, height, first_slot, body_id) and
|
||||
FP64 geometry (cx, cy, r) at the applied offset; plus the separate per-body
|
||||
diagnostics (cme_x, cme_y, wen_x, wen_y, link_dm, n_links), never part of obs.
|
||||
Band table (x0, y0, width, height, first_slot, body_id), immutable per run;
|
||||
FP64 motion records (cx_old, cy_old, cx_new, cy_new, r, ux, uy), immutable
|
||||
during a step; separate per-body diagnostics (see DIAG), never part of obs.
|
||||
"""
|
||||
BAND_INTS, DIAG_FLOATS = 6, 6
|
||||
BAND_INTS, GEOM_DOUBLES = 6, 7
|
||||
DIAG = ("cme_x", "cme_y", "wen_x", "wen_y", "link_dm", "n_links",
|
||||
"cover_jx", "cover_jy", "cover_mass", "uncover_jx", "uncover_jy", "uncover_mass",
|
||||
"eq_only", "no_donor", "crossings")
|
||||
DIAG_FLOATS = len(DIAG)
|
||||
|
||||
def __init__(self, band, geom):
|
||||
import numpy as np
|
||||
import pycuda.driver as cuda
|
||||
self.n_bodies = len(geom) // 3
|
||||
self.n_bodies = len(geom) // self.GEOM_DOUBLES
|
||||
self.n_slots = int(band[-2] + band[-4] * band[-3] * 8)
|
||||
self.n_cells = self.n_slots // 8
|
||||
self.geom = list(geom) # host copy of the uploaded records
|
||||
self.band_gpu = cuda.to_device(np.asarray(band, np.int32))
|
||||
self.geom_gpu = cuda.to_device(np.asarray(geom, np.float64))
|
||||
self.diag_gpu = cuda.to_device(np.zeros(self.n_bodies * self.DIAG_FLOATS, np.float32))
|
||||
@@ -201,9 +214,27 @@ class MovingTables:
|
||||
width, height = x1 - x0 + 1, y1 - y0 + 1
|
||||
band += [x0, y0, width, height, first, int(body_id)]
|
||||
first += width * height * 8
|
||||
geom += [*motion.centre(motion.applied), motion.radius]
|
||||
geom += MovingTables.record(motion)
|
||||
return band, geom
|
||||
|
||||
@staticmethod
|
||||
def record(motion):
|
||||
"""FP64 motion record of one circle for its next step (hold if nothing pending)."""
|
||||
c_old, c_new, u = motion.record()
|
||||
return [*c_old, *c_new, motion.radius, *u]
|
||||
|
||||
def upload_records(self, circles, stream) -> bool:
|
||||
"""Upload the records when they differ from the device copy; True if advancing."""
|
||||
import numpy as np
|
||||
import pycuda.driver as cuda
|
||||
geom = [v for motion, _ in circles for v in self.record(motion)]
|
||||
if geom != self.geom:
|
||||
self._staged = np.asarray(geom, np.float64) # keep alive until the copy completes
|
||||
cuda.memcpy_htod_async(self.geom_gpu, self._staged, stream)
|
||||
self.geom = geom
|
||||
g = self.GEOM_DOUBLES
|
||||
return any(geom[i:i + 2] != geom[i + 2:i + 4] for i in range(0, len(geom), g))
|
||||
|
||||
def enable_slot_debug(self):
|
||||
import numpy as np
|
||||
import pycuda.driver as cuda
|
||||
|
||||
@@ -1,6 +1,7 @@
|
||||
// Translating-circle cut links, fixed centre (S2a; docs/moving_bodies.md, M2).
|
||||
// One thread per (swept-band cell, direction) slot. Geometry comes from an FP64
|
||||
// {cx, cy, r} table per body; the operation order mirrors
|
||||
// Translating circles (docs/moving_bodies.md): M1 MovingCircleCells reclassifies
|
||||
// band cells and refills uncovered ones; M2 MovingCircleLinks applies the cut links.
|
||||
// Geometry comes from an immutable FP64 record per body
|
||||
// {cx_old, cy_old, cx_new, cy_new, r, ux, uy}; the M2 operation order mirrors
|
||||
// CircleGeometry.build_cut_links / find_circle_intersection with explicit _rn
|
||||
// operations (bit parity with the host builder is a test target, not a claim).
|
||||
// Moving circles are never in the packed static links, so this kernel is the
|
||||
@@ -10,7 +11,10 @@
|
||||
|
||||
#if DIM == 2
|
||||
constexpr int MOVING_BAND_INTS = 6; // x0, y0, width, height, first_slot, body_id
|
||||
constexpr int MOVING_DIAG_FLOATS = 6; // cme_x, cme_y, wen_x, wen_y, link_dm, n_links
|
||||
constexpr int MOVING_GEOM_DOUBLES = 7;
|
||||
// cme_x, cme_y, wen_x, wen_y, link_dm, n_links, cover_jx, cover_jy, cover_mass,
|
||||
// uncover_jx, uncover_jy, uncover_mass, eq_only, no_donor, crossings
|
||||
constexpr int MOVING_DIAG_FLOATS = 15;
|
||||
|
||||
__device__ __forceinline__ double mc_sq(double v) { return __dmul_rn(v, v); }
|
||||
|
||||
@@ -53,7 +57,8 @@ __global__ void MovingCircleLinks(
|
||||
const unsigned int dir = local % 8u + 1u;
|
||||
const int x = b[0] + (int)((local / 8u) % (unsigned int)b[2]);
|
||||
const int y = b[1] + (int)((local / 8u) / (unsigned int)b[2]);
|
||||
const double cx = geom[3 * m], cy = geom[3 * m + 1], r = geom[3 * m + 2];
|
||||
const double* g = geom + MOVING_GEOM_DOUBLES * m;
|
||||
const double cx = g[2], cy = g[3], r = g[4];
|
||||
const int ex = d_cx[dir], ey = d_cy[dir];
|
||||
float* debug = slot_debug == nullptr ? nullptr : slot_debug + 4 * tid;
|
||||
if (debug != nullptr) debug[0] = -1.0f;
|
||||
@@ -85,7 +90,7 @@ __global__ void MovingCircleLinks(
|
||||
|
||||
const int bid = b[5];
|
||||
const float omega = action_omega(action, bid);
|
||||
const float Uw = -omega * ry, Vw = omega * rx; // fixed centre: no translation term
|
||||
const float Uw = (float)g[5] - omega * ry, Vw = (float)g[6] + omega * rx; // u + omega x r
|
||||
float exchange[3], dm;
|
||||
apply_bouzidi_link(dir, q, linear_index((unsigned int)x, (unsigned int)y), fi,
|
||||
Uw, Vw, rx, ry, fallback, obs, bid, exchange, &dm);
|
||||
@@ -99,6 +104,80 @@ __global__ void MovingCircleLinks(
|
||||
atomicAdd(&d[5], 1.0f);
|
||||
}
|
||||
}
|
||||
// M1: one thread per band cell, launched only on steps with c_old != c_new.
|
||||
// Membership comes from the immutable old/new predicates, never from flags, so
|
||||
// classification, refill and the flag write of one cell cannot race with another.
|
||||
// ddf holds post-collision G^{n-1}; refilled cells get all nine slots at that level.
|
||||
__global__ void MovingCircleCells(
|
||||
fpxx* fi, uint16_t* flag, const int* band, const double* geom,
|
||||
unsigned int n_bodies, unsigned int n_cells, const float* action, float* diag)
|
||||
{
|
||||
const unsigned int tid = threadIdx.x + blockIdx.x * blockDim.x;
|
||||
if (tid >= n_cells) return;
|
||||
unsigned int m = 0;
|
||||
while (m + 1 < n_bodies && tid >= (unsigned int)band[(m + 1) * MOVING_BAND_INTS + 4] / 8u) ++m;
|
||||
const int* b = band + m * MOVING_BAND_INTS;
|
||||
const unsigned int local = tid - (unsigned int)b[4] / 8u;
|
||||
const int x = b[0] + (int)(local % (unsigned int)b[2]);
|
||||
const int y = b[1] + (int)(local / (unsigned int)b[2]);
|
||||
const double* g = geom + MOVING_GEOM_DOUBLES * m;
|
||||
const bool in_old = mc_inside(x, y, g[0], g[1], g[4]);
|
||||
const bool in_new = mc_inside(x, y, g[2], g[3], g[4]);
|
||||
if (in_old == in_new) return;
|
||||
const unsigned long k = linear_index((unsigned int)x, (unsigned int)y);
|
||||
float* d = diag + m * MOVING_DIAG_FLOATS;
|
||||
float rho = 0.0f, jx = 0.0f, jy = 0.0f;
|
||||
if (in_new) { // covered: its fluid content leaves the active set
|
||||
for (unsigned int i = 0; i < NQ; ++i) {
|
||||
const float f = load_ddf(fi, index_f(k, i));
|
||||
rho += f; jx += (float)d_cx[i] * f; jy += (float)d_cy[i] * f;
|
||||
}
|
||||
flag[k] = FLAG_SOLID | FLAG_OBSTACLE | FLAG_BC_CURVED;
|
||||
atomicAdd(&d[6], jx); atomicAdd(&d[7], jy); atomicAdd(&d[8], rho);
|
||||
atomicAdd(&d[14], 1.0f);
|
||||
return;
|
||||
}
|
||||
// Uncovered: donors are fluid in BOTH geometries; NEQ from the neighbour with the
|
||||
// largest unnormalized c.n > 0 (first direction wins ties), n from c_new.
|
||||
const float nx = (float)__dsub_rn((double)x, g[2]), ny = (float)__dsub_rn((double)y, g[3]);
|
||||
float rho_sum = 0.0f, best_dot = 0.0f;
|
||||
int n_valid = 0, best = -1;
|
||||
for (unsigned int i = 1; i < NQ; ++i) {
|
||||
const int xn = x + d_cx[i], yn = y + d_cy[i];
|
||||
if (mc_inside(xn, yn, g[0], g[1], g[4]) || mc_inside(xn, yn, g[2], g[3], g[4])) continue;
|
||||
const unsigned long kn = linear_index((unsigned int)xn, (unsigned int)yn);
|
||||
for (unsigned int q = 0; q < NQ; ++q) rho_sum += load_ddf(fi, index_f(kn, q));
|
||||
++n_valid;
|
||||
const float dot = (float)d_cx[i] * nx + (float)d_cy[i] * ny;
|
||||
if (dot > best_dot) { best_dot = dot; best = (int)i; }
|
||||
}
|
||||
const float omega = action_omega(action, b[5]);
|
||||
const float uwx = (float)g[5] - omega * ny, uwy = (float)g[6] + omega * nx; // at the cell
|
||||
float feq[NQ], neq[NQ];
|
||||
compute_feq(n_valid > 0 ? rho_sum / (float)n_valid : (float)RHO, uwx, uwy, feq);
|
||||
for (unsigned int i = 0; i < NQ; ++i) neq[i] = 0.0f;
|
||||
if (best >= 0) {
|
||||
const unsigned long kd = linear_index((unsigned int)(x + d_cx[best]), (unsigned int)(y + d_cy[best]));
|
||||
float gd[NQ], rd = 0.0f, mx = 0.0f, my = 0.0f, fd[NQ];
|
||||
for (unsigned int i = 0; i < NQ; ++i) {
|
||||
gd[i] = load_ddf(fi, index_f(kd, i));
|
||||
rd += gd[i]; mx += (float)d_cx[i] * gd[i]; my += (float)d_cy[i] * gd[i];
|
||||
}
|
||||
compute_feq(rd, mx / rd, my / rd, fd);
|
||||
for (unsigned int i = 0; i < NQ; ++i) neq[i] = gd[i] - fd[i];
|
||||
} else {
|
||||
atomicAdd(&d[12], 1.0f);
|
||||
}
|
||||
if (n_valid == 0) atomicAdd(&d[13], 1.0f); // hard failure, checked before commit
|
||||
for (unsigned int i = 0; i < NQ; ++i) {
|
||||
const float f = feq[i] + neq[i];
|
||||
store_ddf(fi, index_f(k, i), f);
|
||||
rho += f; jx += (float)d_cx[i] * f; jy += (float)d_cy[i] * f;
|
||||
}
|
||||
flag[k] = FLAG_FLUID;
|
||||
atomicAdd(&d[9], jx); atomicAdd(&d[10], jy); atomicAdd(&d[11], rho);
|
||||
atomicAdd(&d[14], 1.0f);
|
||||
}
|
||||
#endif // DIM == 2
|
||||
|
||||
#endif // CELERIS_MOVING_CIRCLE_CU
|
||||
|
||||
@@ -40,6 +40,7 @@ class LBMStepper:
|
||||
self.sensor_fn = module.get_function("SensorKernel")
|
||||
self.force_region_fn = module.get_function("ForceRegionKernel")
|
||||
self.moving_fn = None if self._esopull else module.get_function("MovingCircleLinks")
|
||||
self.moving_cells_fn = None if self._esopull else module.get_function("MovingCircleCells")
|
||||
|
||||
tpb = cfg.threads_per_block
|
||||
self.block = (tpb, 1, 1)
|
||||
@@ -77,6 +78,11 @@ class LBMStepper:
|
||||
"""Advance *n* time steps."""
|
||||
if n > 0 and self.blocked is not None:
|
||||
raise RuntimeError(self.blocked)
|
||||
self._step(n, action_gpu=action_gpu, obs_gpu=obs_gpu, stream=stream)
|
||||
|
||||
def _step(self, n, *, action_gpu, obs_gpu, stream=None, advance=False):
|
||||
"""Unguarded steps; ``advance`` runs one step that moves translating circles."""
|
||||
assert not advance or n == 1
|
||||
f = self.field
|
||||
launch_kw = {}
|
||||
if stream is not None:
|
||||
@@ -89,7 +95,7 @@ class LBMStepper:
|
||||
block=self.block, grid=self.boundary_grid, **launch_kw,
|
||||
)
|
||||
if self.moving is not None:
|
||||
self._launch_moving(action_gpu, obs_gpu, **launch_kw)
|
||||
self._launch_moving(action_gpu, obs_gpu, advance, **launch_kw)
|
||||
self._launch_curved(action_gpu, obs_gpu, **launch_kw)
|
||||
|
||||
if self._esopull:
|
||||
@@ -115,8 +121,14 @@ class LBMStepper:
|
||||
# Any completed GPU step invalidates host DDF mirror.
|
||||
f.invalidate_host_ddf_cache()
|
||||
|
||||
def _launch_moving(self, action_gpu, obs_gpu, **launch_kw):
|
||||
def _launch_moving(self, action_gpu, obs_gpu, advance, **launch_kw):
|
||||
m, tpb = self.moving, self.cfg.threads_per_block
|
||||
if advance: # M1 completes before M2 on the same stream
|
||||
self.moving_cells_fn(
|
||||
self.field.ddf_gpu, self.field.flag_gpu, m.band_gpu, m.geom_gpu,
|
||||
np.uint32(m.n_bodies), np.uint32(m.n_cells), action_gpu, m.diag_gpu,
|
||||
block=(tpb, 1, 1), grid=((m.n_cells + tpb - 1) // tpb, 1, 1), **launch_kw,
|
||||
)
|
||||
self.moving_fn(
|
||||
self.field.ddf_gpu, self.field.flag_gpu, m.band_gpu, m.geom_gpu,
|
||||
np.uint32(m.n_bodies), np.uint32(m.n_slots), action_gpu, obs_gpu,
|
||||
|
||||
@@ -29,6 +29,7 @@ from .cuda import compiler
|
||||
from .lbm.descriptors import FLUID
|
||||
|
||||
# Translating circles step only through their moving tables (fixed centre in S2a).
|
||||
PENDING_BLOCKED = "a pending translation target is advanced only by Simulation.run(1)"
|
||||
MOVING_BLOCKED = ("translating circles need the moving kernels and their tables (attached "
|
||||
"at initialize); "
|
||||
"positive-step runs are rejected")
|
||||
@@ -217,7 +218,13 @@ class Simulation:
|
||||
# Direct stepper calls stay blocked until the moving tables exist.
|
||||
tables = getattr(self, "_moving_tables", None)
|
||||
self.stepper.moving = tables
|
||||
self.stepper.blocked = MOVING_BLOCKED if self._moving_bodies() and tables is None else None
|
||||
moving = self._moving_bodies()
|
||||
if moving and tables is None:
|
||||
self.stepper.blocked = MOVING_BLOCKED
|
||||
elif any(obj.motion.pending is not None for obj in moving):
|
||||
self.stepper.blocked = PENDING_BLOCKED # direct steps never consume a target
|
||||
else:
|
||||
self.stepper.blocked = None
|
||||
|
||||
def _build_moving_tables(self):
|
||||
moving = self._moving_bodies()
|
||||
@@ -228,8 +235,8 @@ class Simulation:
|
||||
|
||||
@context_bound
|
||||
def moving_diagnostics(self) -> np.ndarray:
|
||||
"""Per translating circle (in body order): cme_x, cme_y, wen_x, wen_y, link_dm,
|
||||
n_links accumulated since obs were last zeroed (diagnostics only)."""
|
||||
"""Per translating circle (in body order), columns ``MovingTables.DIAG``,
|
||||
accumulated since obs were last zeroed (diagnostics only, not a force RHS)."""
|
||||
self._require_usable()
|
||||
tables = self.stepper.moving
|
||||
if tables is None:
|
||||
@@ -354,17 +361,19 @@ class Simulation:
|
||||
omega: rotation speed (lattice units). ``None`` = no change.
|
||||
vx: translation velocity x (future). ``None`` = no change.
|
||||
vy: translation velocity y (future). ``None`` = no change.
|
||||
offset: target centre offset of a translating circle. Displacement
|
||||
is disabled until the moving cell update exists (S3): any offset
|
||||
raises before mutation; translating circles hold at the applied
|
||||
offset and accept omega only.
|
||||
offset: target centre offset of a translating circle; with vx/vy
|
||||
(default 0) it replaces the pending target after validation
|
||||
against the applied offset (docs/moving_bodies.md, 4a). Only
|
||||
``run(1)`` advances it; direct stepper calls are blocked meanwhile.
|
||||
"""
|
||||
obj = self.bodies.get(id)
|
||||
motion = obj.motion
|
||||
if motion is not None:
|
||||
if offset is not None or vx is not None or vy is not None:
|
||||
raise ValueError("translating-circle displacement is disabled until S3; "
|
||||
"only omega may be set")
|
||||
if offset is None and (vx is not None or vy is not None):
|
||||
raise ValueError("a translating circle takes vx/vy only with an offset")
|
||||
target = None if offset is None else motion.check_target(
|
||||
offset, (0.0 if vx is None else vx, 0.0 if vy is None else vy))
|
||||
vx = vy = None # the motion target, not ObjectState, owns translation
|
||||
elif offset is not None:
|
||||
raise ValueError(f"body {id} was not declared with a translation")
|
||||
proposed = dict(vars(obj.state))
|
||||
@@ -382,6 +391,9 @@ class Simulation:
|
||||
if vy is not None:
|
||||
obj.state.vy = float(vy)
|
||||
changed = True
|
||||
if motion is not None and target is not None:
|
||||
motion.pending = target
|
||||
self._attach_moving()
|
||||
if changed:
|
||||
self.bodies._refresh_action_from_objects()
|
||||
|
||||
@@ -700,8 +712,13 @@ class Simulation:
|
||||
checkpoint_interval = operator.index(checkpoint_interval)
|
||||
if steps < 0 or checkpoint_interval < 0:
|
||||
raise ValueError("steps and checkpoint_interval must be nonnegative")
|
||||
if steps > 0 and self._moving_bodies() and self.stepper.moving is None:
|
||||
moving = self._moving_bodies()
|
||||
if steps > 0 and moving and self.stepper.moving is None:
|
||||
raise RuntimeError(MOVING_BLOCKED)
|
||||
if moving and checkpoint_interval > 0:
|
||||
raise ValueError("checkpoints are not supported with translating circles")
|
||||
if steps > 1 and any(obj.motion.pending is not None for obj in moving):
|
||||
raise ValueError("a pending translation target needs run(1)")
|
||||
self._assert_runtime_contracts()
|
||||
# Discard any uncommitted body edits before stepping.
|
||||
if self.bodies.has_pending_edit():
|
||||
@@ -736,6 +753,9 @@ class Simulation:
|
||||
if done < steps or done % checkpoint_interval == 0:
|
||||
stream.synchronize()
|
||||
self.save_checkpoint()
|
||||
elif self.stepper.moving is not None and steps > 0:
|
||||
self._run_moving(steps, stream, sync_obs)
|
||||
return
|
||||
else:
|
||||
self.stepper.step(
|
||||
steps,
|
||||
@@ -754,6 +774,34 @@ class Simulation:
|
||||
|
||||
stream.synchronize()
|
||||
|
||||
def _run_moving(self, steps, stream, sync_obs):
|
||||
"""Steps with translating circles: hold, or one advancing step that commits
|
||||
each pending target only after synchronization and a zero ``no_donor`` count
|
||||
(docs/moving_bodies.md 4a, 5). Any failure leaves the Simulation unusable."""
|
||||
tables = self.stepper.moving
|
||||
circles = [(obj.motion, obj.obj_id) for obj in self._moving_bodies()]
|
||||
try:
|
||||
advancing = tables.upload_records(circles, stream)
|
||||
self.stepper._step(steps, action_gpu=self.bodies.action_gpu,
|
||||
obs_gpu=self.bodies.obs_gpu, stream=stream, advance=advancing)
|
||||
self.bodies.increment_obs_steps(steps)
|
||||
if sync_obs:
|
||||
self.bodies.download_obs_full_async(stream)
|
||||
stream.synchronize()
|
||||
if advancing:
|
||||
no_donor = tables.download_diagnostics()[:, tables.DIAG.index("no_donor")]
|
||||
if np.any(no_donor > 0):
|
||||
raise RuntimeError("translating circle refill found no donor cell")
|
||||
for motion, _ in circles:
|
||||
motion.commit()
|
||||
tables.upload_records(circles, stream) # first hold: old = new = applied, u = 0
|
||||
stream.synchronize()
|
||||
self._attach_moving()
|
||||
except BaseException as exc:
|
||||
for motion, _ in circles: # device geometry/state unknown after any failure
|
||||
motion.fail(exc)
|
||||
raise
|
||||
|
||||
def step(self, n: int = 1):
|
||||
"""Advance *n* steps (convenience for interactive use)."""
|
||||
self.run(n)
|
||||
@@ -771,8 +819,16 @@ class Simulation:
|
||||
self.field.download_ddf(step_id=self.stepper.step_count)
|
||||
return self.field.ddf.copy()
|
||||
|
||||
@context_bound
|
||||
def get_flags(self) -> np.ndarray:
|
||||
return self.field.flag.copy()
|
||||
"""Flag array; read back from the device when translating circles may have moved."""
|
||||
if getattr(self.stepper, "moving", None) is None:
|
||||
return self.field.flag.copy()
|
||||
self._require_usable()
|
||||
self.stream.synchronize()
|
||||
flags = np.empty_like(self.field.flag)
|
||||
cuda.memcpy_dtoh(flags, self.field.flag_gpu)
|
||||
return flags
|
||||
|
||||
# -- Runtime parameter updates -------------------------------------------
|
||||
@context_bound
|
||||
|
||||
@@ -217,20 +217,21 @@ def test_add_body_and_set_body_route_targets(monkeypatch):
|
||||
bid = sim.add_body("circle", (50., 40.), 10., translation=Translation((-2., 2., -3., 3.), .1))
|
||||
still = sim.add_body("circle", (100., 40.), 5.)
|
||||
motion = sim.bodies.get(bid).motion
|
||||
assert sim.stepper.blocked is not None # declaring a moving circle blocks direct steps
|
||||
for initialized in (False, True): # S2a: displacement disabled before and after initialize
|
||||
sim._initialized = initialized
|
||||
for kwargs in (dict(offset=(0., .3), vy=.05), dict(offset=(0., 0.)), dict(vy=.05),
|
||||
dict(vx=.05), dict(offset=(0., .1), omega=.01)):
|
||||
with pytest.raises(ValueError, match="disabled until S3"):
|
||||
sim.set_body(bid, **kwargs)
|
||||
assert motion.pending is None and motion.applied == (0., 0.)
|
||||
sim._initialized = False
|
||||
sim._moving_tables = "tables"
|
||||
sim._attach_moving()
|
||||
assert sim.stepper.blocked is None
|
||||
sim.set_body(bid, offset=(0., .3), vy=.05)
|
||||
assert sim.stepper.blocked == module.PENDING_BLOCKED # direct steps never consume a target
|
||||
sim.set_body(bid, offset=(0., .4), vy=.05) # last wins
|
||||
for kwargs in (dict(offset=(0., .6)), dict(vy=.05), dict(vx=.05),
|
||||
dict(offset=(0., .1), vy=.2), dict(offset=(0., 0.), vy=.01)):
|
||||
with pytest.raises(ValueError):
|
||||
sim.set_body(bid, **kwargs)
|
||||
assert motion.pending == ((0., .4), (0., .05)) # invalid input keeps the target
|
||||
sim.set_body(bid, omega=.01) # omega-only keeps the pending target
|
||||
obj = sim.bodies.get(bid)
|
||||
assert obj.state.omega == 0.
|
||||
sim.set_body(bid, omega=.01) # spin on hold stays allowed
|
||||
assert obj.state.omega == .01 and (obj.state.vx, obj.state.vy) == (0., 0.)
|
||||
assert motion.pending is None
|
||||
assert motion.pending == ((0., .4), (0., .05)) and motion.applied == (0., 0.)
|
||||
with pytest.raises(ValueError, match="not declared"):
|
||||
sim.set_body(still, offset=(0., .1))
|
||||
sim.set_body(still, vx=.01) # stationary API unchanged
|
||||
@@ -243,6 +244,69 @@ def test_add_body_and_set_body_route_targets(monkeypatch):
|
||||
**kwargs)
|
||||
|
||||
|
||||
class Tables:
|
||||
"""MovingTables double: records uploads, steps and diagnostics reads."""
|
||||
DIAG = ("no_donor",)
|
||||
|
||||
def __init__(self, log, no_donor=0., fail=None):
|
||||
self.log, self.no_donor, self.fail = log, no_donor, fail
|
||||
|
||||
def upload_records(self, circles, stream):
|
||||
records = [m.record() for m, _ in circles]
|
||||
self.log.append(("upload", records))
|
||||
if self.fail == "upload":
|
||||
raise RuntimeError("upload failed")
|
||||
return any(r[0] != r[1] for r in records)
|
||||
|
||||
def zero_diagnostics(self, stream):
|
||||
self.log.append("zero_diag")
|
||||
|
||||
def download_diagnostics(self):
|
||||
return np.array([[self.no_donor]], np.float32)
|
||||
|
||||
|
||||
def advancing_simulation(monkeypatch, **tables):
|
||||
sim, writes = moving_simulation(monkeypatch)
|
||||
sim.bodies.objects[0].motion = circle()
|
||||
sim.stepper = SimpleNamespace(
|
||||
moving=Tables(writes, **tables), blocked=None,
|
||||
_step=lambda n, advance=False, **k: writes.append(("step", n, advance)))
|
||||
sim._moving_tables = sim.stepper.moving
|
||||
return sim, writes, sim.bodies.objects[0].motion
|
||||
|
||||
|
||||
def test_run_one_advances_and_commits_after_sync(monkeypatch):
|
||||
sim, writes, motion = advancing_simulation(monkeypatch)
|
||||
motion.set_target((0., .25), (0., .05))
|
||||
with pytest.raises(ValueError, match="run\\(1\\)"): # multi-step with a target
|
||||
sim.run(2, stream=STREAM)
|
||||
assert writes == [] and motion.pending is not None
|
||||
sim.run(0, stream=STREAM) # zero steps never consume the target
|
||||
assert motion.pending is not None
|
||||
sim.run(1, stream=STREAM)
|
||||
uploads = [w for w in writes if isinstance(w, tuple) and w[0] == "upload"]
|
||||
assert uploads[0][1] == [((50., 40.), (50., 40.25), (0., .05))]
|
||||
assert uploads[1][1] == [((50., 40.25), (50., 40.25), (0., 0.))] # first hold after move
|
||||
assert ("step", 1, True) in writes
|
||||
assert motion.applied == (0., .25) and motion.pending is None
|
||||
writes.clear()
|
||||
sim.run(3, stream=STREAM) # hold steps: no advance
|
||||
assert ("step", 3, False) in writes and motion.applied == (0., .25)
|
||||
|
||||
|
||||
@pytest.mark.parametrize("failure", ("upload", "no_donor"))
|
||||
def test_failed_advance_commits_nothing_and_fails_simulation(monkeypatch, failure):
|
||||
kwargs = dict(fail="upload") if failure == "upload" else dict(no_donor=1.)
|
||||
sim, writes, motion = advancing_simulation(monkeypatch, **kwargs)
|
||||
motion._on_failure = sim._fail_moving
|
||||
motion.set_target((0., .25), (0., .05))
|
||||
with pytest.raises(RuntimeError):
|
||||
sim.run(1, stream=STREAM)
|
||||
assert motion.applied == (0., 0.) and motion.failed and sim._failed is not None
|
||||
with pytest.raises(RuntimeError, match="Simulation failed"):
|
||||
sim.run(1, stream=STREAM)
|
||||
|
||||
|
||||
def test_moving_tables_layout_and_attachment(monkeypatch):
|
||||
from CelerisLab.body.motion import MovingTables
|
||||
a = MovingCircle((50., 40.), 10., Translation((-1., 1., -3., 3.), .1))
|
||||
@@ -250,7 +314,7 @@ def test_moving_tables_layout_and_attachment(monkeypatch):
|
||||
band, geom = MovingTables.layout([(a, 2), (b, 5)])
|
||||
wa, ha = 63 - 37 + 1, 55 - 25 + 1
|
||||
assert band == [37, 25, wa, ha, 0, 2, 83, 23, 16, 16, wa * ha * 8, 5]
|
||||
assert geom == [50., 40., 10., 90.25, 30.5, 5.]
|
||||
assert geom == [50., 40., 50., 40., 10., 0., 0., 90.25, 30.5, 90.25, 30.5, 5., 0., 0.] # hold
|
||||
sim = routed_simulation(monkeypatch)
|
||||
sim.add_body("circle", (50., 40.), 10., translation=Translation((0., 0., 0., 0.), .1))
|
||||
sim._moving_tables = "tables"
|
||||
|
||||
Reference in New Issue
Block a user