Compare commits
26
Commits
| Author | SHA1 | Date | |
|---|---|---|---|
|
|
8c1f3c4d38 | ||
|
|
8505a723a5 | ||
|
|
10cea29e8e | ||
|
|
cd73fcd689 | ||
|
|
ee6e68bbfd | ||
|
|
f68f5264b2 | ||
|
|
03b23cb511 | ||
|
|
8ce810cc99 | ||
|
|
653fd997a6 | ||
|
|
25cc5cfd68 | ||
|
|
6f9e40caa3 | ||
|
|
c21a365eca | ||
|
|
d50b8ba1be | ||
|
|
b920a303a4 | ||
|
|
e25828968b | ||
|
|
0f8a428f07 | ||
|
|
983020b625 | ||
|
|
e1ed27b438 | ||
|
|
71c12aedce | ||
|
|
36eac713f4 | ||
|
|
e682d16c3a | ||
|
|
b448bc3f71 | ||
|
|
8a8e8a03ff | ||
|
|
305ab2409b | ||
|
|
c986343a15 | ||
|
|
b2034f87e0 |
@@ -3,9 +3,12 @@
|
||||
Owner direction 2026-10-28: unified style over historical compatibility, inside `src`.
|
||||
No numerical change is intended. Stages after 5f69eba: 65339ab (one boundary/Jet record,
|
||||
host preparation helpers, pair accessors), c5a3cdb (`block_collide` families), 442979b
|
||||
(exchange metadata on the boundary record). Each was checked by parent/candidate
|
||||
config.h/PTX comparison and host/GPU equality runs; c5a3cdb changed the PTX of the two
|
||||
channel-root SRT modules only (omega load/clamp moved, outputs exact).
|
||||
(exchange metadata on the boundary record). All three received parent/candidate
|
||||
generated-config and host/GPU equality checks. Only the collision slice ran a new PTX
|
||||
matrix (a4fc46d, whose src is identical to c5a3cdb); 65339ab and 442979b reused prior
|
||||
PTX evidence after unchanged CUDA/compiler and generated-input checks. In the collision
|
||||
matrix, only the two channel-root SRT modules changed PTX (omega load/clamp moved);
|
||||
their parent/candidate outputs matched exactly in the exercised cases.
|
||||
|
||||
## Conventions (defined once, in `lbm.hierarchy.Hierarchy`)
|
||||
|
||||
|
||||
@@ -0,0 +1,403 @@
|
||||
# Moving circles on the uniform dense grid — S0 DRAFT contract
|
||||
|
||||
Status: **DRAFT, not frozen** for the force RHS, kernels and physical
|
||||
acceptance. Host-only S1 (state, validation, guards; no kernels) is on branch
|
||||
`claude/moving-bodies-s1`; positive-step runs with a moving body are rejected
|
||||
until S2. When S2 lifts that guard, `run(n > 1)` with a pending target must
|
||||
raise `ValueError` before mutation (4a), and `run(0)` must not upload a motion
|
||||
record or consume the pending target. Base: main
|
||||
`b2034f87e03656c5f67ab88562a70d6c338b9918`. Scope approved by Frank (via TRON,
|
||||
thread C0C78UFKGRM/1793335557.460509): uniform-grid prescribed translation and
|
||||
single-cylinder free transverse VIV. Out of scope: refinement/blocked path,
|
||||
pinball assembly, EsoPull, shifted storage, D3Q19, DynamisLab, manuscript.
|
||||
|
||||
Items marked **OPEN** are undecided; items marked **CANDIDATE** are proposals.
|
||||
|
||||
## 1. Working standards
|
||||
|
||||
- Lean, readable kernels in the existing idiom (`kernels/step/*.cu`, shared
|
||||
`apply_bouzidi_link`); Python manages state and validation. No framework,
|
||||
no compatibility layer, no new dependency.
|
||||
- The existing stationary/rotating path (`CurvedBoundaryKernel`, host link
|
||||
packing, `sync_bodies`) stays byte-identical. Moving circles are a separate,
|
||||
explicitly enabled path with its own numerical contract (this file).
|
||||
- One authority per state (section 4). No per-step full-domain download,
|
||||
rebuild or recompile; local swept-band device updates are in scope.
|
||||
- Small slices (S1 host spec, S2 kernels without motion, S3 motion + refill,
|
||||
S4 oscillator); each reviewed SHA is frozen before targeted runs.
|
||||
- Tolerances, benchmark parameters and the force composition used in the
|
||||
structural RHS are declared before the runs that judge them. Fresh vs carried
|
||||
evidence and existing baseline failures are recorded; smoke, convergence and
|
||||
physical validation are reported separately. No undeclared force filtering.
|
||||
- External code (waLBerla GPLv3, FluidX3D source-available, VIVSIM) is a
|
||||
concept reference only; nothing is copied.
|
||||
|
||||
## 2. Current update order (main b2034f87, dense double buffer)
|
||||
|
||||
`Simulation.run(n)`: upload host action (if `upload_act`) → zero obs (if
|
||||
`zero_obs`) → `DoubleBufferStepper.step(n)` → download obs → synchronize.
|
||||
|
||||
Per step (`lbm/stepper.py` `step`):
|
||||
|
||||
1. `CurvedBoundaryKernel` over the packed cut links: reads post-collision
|
||||
`G^{n-1}` at `k_f` (and `k_f - c` for q < 1/2), adds the moving-wall term
|
||||
(`Uw = -omega*ry, Vw = omega*rx` + optional per-link `surface_velocity`),
|
||||
writes the reflected value into slot `opp` of the solid node `k_s`, and
|
||||
accumulates `c(f_toward + f_reflected)` and torque into `obs`.
|
||||
2. `OneStep` pull/collide `ddf → temp` (collides FLUID cells; obstacle
|
||||
interiors flagged `SOLID|OBSTACLE|BC_CURVED`), then buffer swap.
|
||||
3. `ForceRegionKernel`, `SensorKernel`.
|
||||
|
||||
Raw `obs` force is the force **on the body** (link derivation in section 6;
|
||||
DFG docs report `Cd > 0` from it). Translation `vx, vy` in the action slots is
|
||||
currently stored but not read by any curved kernel.
|
||||
|
||||
## 3. Proposed update order (moving circles)
|
||||
|
||||
Per step n, before the existing step 1:
|
||||
|
||||
- **M1 `MovingCircleCells`** (one thread per swept-band cell): classify with the
|
||||
analytic predicates `in_old = inside(x, c^{n-1})`, `in_new = inside(x, c^n)`;
|
||||
write `flag = in_new ? SOLID|OBSTACLE|BC_CURVED : FLUID`; covered cells
|
||||
(`!in_old && in_new`) accumulate `node_cover`; fresh cells
|
||||
(`in_old && !in_new`) are refilled (section 5) and accumulate `node_uncover`.
|
||||
- **M2 `MovingCircleLinks`** (one thread per band cell × direction): compute
|
||||
q and hit point at `c^n` (FP64, same formula as `find_circle_intersection`,
|
||||
same clamp [0.01, 1] and fallback classes as `CircleGeometry.build_cut_links`),
|
||||
skip non-cut slots, call `apply_bouzidi_link` with
|
||||
`u_w = (vx^n, vy^n) + omega × (hit - c^n)`, accumulate the link terms
|
||||
(section 6).
|
||||
|
||||
Then the unchanged steps 1–3. M2 must follow M1 because a q < 1/2 link may
|
||||
read `k_f - c` at a fresh cell.
|
||||
|
||||
Exclusivity: moving bodies are never packed into the `CurvedBoundaryKernel`
|
||||
link arrays (`n_curved` covers stationary/rotating bodies only), so M2 is the
|
||||
sole writer of their solid slots and the sole source of their link force; no
|
||||
double write and no double `obs`. The stationary path's packing, launch and
|
||||
arithmetic are unchanged.
|
||||
|
||||
Completeness: after M1, every (fluid cell, direction) pair in the swept band
|
||||
whose neighbour is `in_new` is visited by M2. This holds because the band is
|
||||
the circle AABB at the declared extreme offsets plus 2 cells, and the host
|
||||
rejects per-step moves with `|c_new - c_old|_2 > 0.5` and bands touching the domain edge,
|
||||
other bodies, sensors or force regions. Each M2 write lands at
|
||||
`index_f(k_s, opp)`, the exact slot the adjacent fluid cell pulls, so every
|
||||
fluid pull from a solid slot of a moving body is BFL-written.
|
||||
|
||||
## 4. State, masks and authority
|
||||
|
||||
- Host `ObjectState` is the single authority for one displacement
|
||||
`offset = (dx, dy)` relative to the immutable base centre, plus
|
||||
`velocity = (vx, vy)` and `omega`. Centres are always derived
|
||||
(`c = base + offset`), never stored as independent state. Velocity is
|
||||
supplied, never differenced from offsets.
|
||||
- `c^{n-1}` is derived from the last *committed* offset, `c^n` from the
|
||||
offset being applied. Per body a small device **motion record**
|
||||
`{c^{n-1}, c^n, u^n}` is uploaded with the action buffer; action-slot meaning
|
||||
is not overloaded. (CANDIDATE layout; OPEN.)
|
||||
- Fail-closed commit: an enqueued upload is not a completed one. The new
|
||||
offset becomes the committed offset only after the stream synchronization
|
||||
that follows the upload and the step group succeeds. Any upload or later
|
||||
synchronization failure leaves host and device motion state unknown: the
|
||||
simulation is marked unusable (existing `_require_usable` lifecycle), and
|
||||
no API may report a host/device match afterwards.
|
||||
|
||||
### 4a. Applied geometry vs pending target (CANDIDATE, narrow API)
|
||||
|
||||
Two distinct facts per moving body, one authority each:
|
||||
`applied` = offset of the geometry last applied on the device, changed only
|
||||
by a successful commit; `pending` = an optional requested target
|
||||
`(offset, velocity)`, changed only by `set_body`. A setter returning proves
|
||||
nothing about the device. The uploaded motion record is derived each run:
|
||||
`{c_old = base + applied, c_new = base + (pending.offset if pending is not
|
||||
None else applied), u = pending.velocity if pending is not None else 0}`
|
||||
(explicit `is not None`, so a zero offset or zero velocity is a valid target). Translational wall velocity exists only on a step
|
||||
that advances to a pending target; a hold step (`c_old == c_new`) has
|
||||
translational `u_w = 0` (spin unchanged), so a non-advancing circle never
|
||||
acts as a sliding wall. One moving step consumes one target.
|
||||
|
||||
| situation | contract |
|
||||
|---|---|
|
||||
| `set_body` target | validated against `applied` (box, `|Δc|_2 ≤ 0.5` cell, speed bound) before mutation; on success replaces `pending` |
|
||||
| target with `offset == applied` and nonzero translational velocity | `ValueError` before mutation (would be a hold step with sliding wall velocity) |
|
||||
| two setters before a run | last wins; each validated against `applied`, not against the earlier pending target |
|
||||
| `run(1)` with pending | upload record, step, sync; on success `applied := pending.offset` (applied stores only the offset), `pending := None` |
|
||||
| `run(1)` without pending | hold: `c_old == c_new`, no cover/uncover, links at `applied`, translational `u_w = 0` |
|
||||
| `run(n > 1)` with pending | `ValueError` before mutation (single-step-only moving API) |
|
||||
| `run(n > 1)`, no pending | allowed: hold for n steps, translational `u_w = 0`, spin only |
|
||||
| `run(0)` | no motion upload, no commit; `pending` kept |
|
||||
| validation failure | `ValueError`, nothing mutated, retry allowed |
|
||||
| upload/launch/sync failure after enqueue | simulation unusable; `applied` and `pending` are not reported as device truth; no re-entry |
|
||||
|
||||
Displacement norm: Euclidean per-step `|c_new - c_old|_2 ≤ 0.5` lattice
|
||||
cells (declared bound; CANDIDATE value).
|
||||
- Membership is never read from a flag that M1 is overwriting. The host
|
||||
verifies at `add_body` that every swept-band cell has base flag FLUID, so the
|
||||
band's base state is known without reading `flag_gpu`.
|
||||
- `flag_gpu` is the device authority for band cells while a moving body
|
||||
exists; the host flag mirror is stale for those cells and is not used.
|
||||
|
||||
## 5. Refill of fresh cells (CANDIDATE)
|
||||
|
||||
- Stored time level: everything in `ddf` at M1 is post-collision `G^{n-1}`; the
|
||||
refill writes all 9 slots at that same level, so every slot of a fresh cell
|
||||
that the following pull or a q < 1/2 link may read is defined. Donor values
|
||||
are read at the same level `G^{n-1}` (old-fluid cells, untouched by M1).
|
||||
- `G = f_eq(rho_bar, u_w(x)) + G_neq(x_d)` where `rho_bar` = mean density of
|
||||
valid neighbours and `x_d` = the valid neighbour best aligned with the outward
|
||||
normal `(x - c^n)/|x - c^n|`; `G_neq(x_d) = G(x_d) - f_eq(rho(x_d), u(x_d))`.
|
||||
- Valid donor: inside the domain, base FLUID, `!in_old && !in_new`.
|
||||
- Sets: `V(x)` = the 8 lattice neighbours that are in the domain, base FLUID
|
||||
and `!in_old && !in_new`. `rho_bar` = mean density over `V(x)`. NEQ donor =
|
||||
the member of `V(x)` with the largest `c·n` among those with `c·n > 0`
|
||||
(`n` = outward normal at `c^n`); none → NEQ part 0 (counted as `eq_only`).
|
||||
- `V(x)` empty → no defined density: the cell is refilled with
|
||||
`f_eq(RHO, u_w(x))` so all 9 slots stay defined, the `no_donor` counter is
|
||||
incremented, and the host treats `no_donor > 0` after the run as a hard
|
||||
failure (simulation unusable). It is not a silent fallback.
|
||||
- Non-emptiness is NOT assumed from the displacement bound. Geometry
|
||||
assumptions: circle radius `R ≥ R_min` (CANDIDATE 4 cells), `|Δc|_2 ≤ 0.5`.
|
||||
S1 includes a host enumeration over sub-cell centre positions and
|
||||
displacement directions on a fine grid within these bounds that checks every
|
||||
fresh cell has nonempty `V(x)`; it is evidence for those bounds, not a
|
||||
proof, and runtime keeps the counter.
|
||||
- Fresh cells only write, old-fluid cells only read: no race.
|
||||
|
||||
## 6. Discrete momentum ledger (to review; not a selection)
|
||||
|
||||
Notation per step n: `F^n` = fluid set after M1; links L = (k_f, dir) at
|
||||
`c^n`; `f_t = G_dir(k_f)`; `f_r` = value M2 writes into `(k_s, opp)`.
|
||||
Fluid momentum before streaming:
|
||||
|
||||
P_pre = P(G^{n-1} on F^{n-1}) - sum_cover sum_q c_q G_q + sum_fresh sum_q c_q G_q^refill
|
||||
|
||||
During the pull, every population entering a fluid cell comes from a fluid
|
||||
cell except the BFL-written slots; `f_t` leaves into solid and is never pulled
|
||||
by a fluid cell. Collision conserves momentum for SRT/MRT without forcing.
|
||||
Hence, in exact arithmetic,
|
||||
|
||||
P_fluid^n - P_fluid^{n-1} = - sum_L c (f_t + f_r) - I_cover + I_uncover + B^n + S^n
|
||||
|
||||
with `B^n` external-boundary momentum (inlet/outlet/walls) and `S^n` forcing.
|
||||
So **`F_CME + node_cover - node_uncover`** is, by construction, the body
|
||||
momentum that balances the fluid ledger of *this* storage/update order.
|
||||
Branch details: halfway fallback links (`f_r = f_t` plus the moving-wall
|
||||
addend) and q < 1/2 links enter the link accumulators exactly like Bouzidi
|
||||
links, with `c(f_t + f_r)`. In the q < 1/2 branch `G_dir(k_f - c)` is an
|
||||
interpolation input only, not a fluid-momentum sink. Pull-slot completeness
|
||||
is the section 3 contract. The scope excludes `ddf_shifting` (rejected), so
|
||||
the unshifted force form applies. The balance is a fluid-side identity; it
|
||||
does not close body+fluid globally in a domain with inlet/outlet/forcing
|
||||
(`B^n`, `S^n`), and the solid-slot writes are boundary plumbing, not physical
|
||||
streaming packets.
|
||||
|
||||
The Wen 2014 relative-velocity link sum adds `-u_w (f_t - f_r) = u_w dm` per
|
||||
link (`dm = f_r - f_t`, the BFL mass defect). That term has no counterpart in
|
||||
the fluid ledger above. Peng et al. (PRE 2016 submitted, §II.C, Eqs. 20–21)
|
||||
show for a straight wall that the bounce-back error accumulated before a node
|
||||
switch cancels with the switching momentum. Whether `sum u_w dm` approximates
|
||||
`I_cover - I_uncover` here (implying double counting if both are used) is
|
||||
**NOT demonstrated** for this implementation: it needs a bounded experiment
|
||||
(prescribed slow translation; cumulative `sum u_w dm` vs cumulative node
|
||||
impulse, per crossing cycle). Momentum conservation is also distinct from
|
||||
Galilean invariance, which G4 tests.
|
||||
|
||||
External reference note: waLBerla's covered/refill momentum switches
|
||||
(`conserveMomentumWhenMapping`, `conserveMomentumWhenRestoring`) default to
|
||||
false (Astra), so their existence does not establish which combination its
|
||||
validation cases use. A bookkeeping identity alone does not select the
|
||||
physical force either.
|
||||
|
||||
Mass ledger (same proposition status). With `M^n` the active-fluid mass
|
||||
summed over `F^n` at the sample point,
|
||||
|
||||
M^n - M^{n-1} = sum_L dm - m_cover + m_uncover + B_m^n
|
||||
|
||||
(`dm = f_r - f_t`, `m_cover`/`m_uncover` = summed `sum_q G` of covered and
|
||||
refilled cells, `B_m` external-boundary mass flux; collision conserves mass).
|
||||
|
||||
Sample point (fixed): post-collision storage after the swap of step n, i.e.
|
||||
the same level `G` that M1 reads at step n+1; both ledgers are evaluated
|
||||
there. Ledger residual = left side minus every listed term. Zero mass flux
|
||||
does not imply zero momentum exchange: closed walls (bounce-back, free-slip)
|
||||
still exchange momentum, so "no inlet/outlet, no forcing" does NOT give
|
||||
`B = S = 0`. The rounding-level momentum check therefore uses either a
|
||||
momentum-isolated case (periodic in both directions, no forcing) or counts
|
||||
the wall impulse explicitly as part of `B`; the mass check may use a closed
|
||||
domain. Otherwise `B`, `S` are not separately measured, the residual is
|
||||
reported as their estimate and is not a pass criterion.
|
||||
Both ledger identities are **proposed, unreviewed**: they are not called
|
||||
proven until review and a host replay agree.
|
||||
|
||||
Accumulators (separate, per body, per step window; all reported raw):
|
||||
`link_cme`, `link_wen_correction` (= `sum u_w dm`), `node_cover`,
|
||||
`node_uncover` (momentum), `link_dm`, `mass_cover`, `mass_uncover` (mass),
|
||||
plus counters of crossings, `eq_only` and `no_donor`, so force spikes can be
|
||||
correlated with crossing events.
|
||||
|
||||
**OPEN — structural RHS candidates (not selected):**
|
||||
(a) `link_cme + node_cover - node_uncover`; (b) `link_cme + link_wen_correction`
|
||||
(no node impulse). Selection follows the ledger review and the bounded
|
||||
experiment above.
|
||||
|
||||
## 7. Rejection semantics
|
||||
|
||||
- Before any mutation (`ValueError`, no device write, host state unchanged):
|
||||
moving body with EsoPull, `ddf_shifting`, refinement/blocked stepper, D3Q19,
|
||||
Jets or surface actuation on that body, `sync_bodies`/checkpoint while a
|
||||
moving body exists; swept band touching domain boundaries, other bodies,
|
||||
sensors or force regions; `set_body` offset outside the declared box,
|
||||
per-step move `|c_new - c_old|_2 > 0.5` cell (the single bound of 4a), speed
|
||||
above the declared bound.
|
||||
- Host motion state is committed only after the upload *and* the following
|
||||
synchronization of the step group succeed (section 4); any failure in
|
||||
between makes the simulation unusable and no host/device match is reported.
|
||||
- The only soft device diagnostic is `eq_only` (NEQ part omitted, counted).
|
||||
`no_donor > 0` is a hard failure: it is checked after the step group and
|
||||
before the step is committed or any results are returned; on failure the
|
||||
simulation is unusable (section 5).
|
||||
|
||||
## 8. Coupling and units (CANDIDATE)
|
||||
|
||||
Lattice units, rho0 = 1, per unit span. Step n uses `c^n`, `u^n` and yields
|
||||
`F^n`. Weak staggered update
|
||||
`v^{n+1} = v^n + (F^n - c v^n - k y^n)/m`, `y^{n+1} = y^n + v^{n+1}`
|
||||
(Newmark is an alternative; OPEN). Mass, natural frequency, damping and `U_r`
|
||||
follow the chosen benchmark's own definitions. Leontini, Thompson & Hourigan,
|
||||
JFS 22:857–864 (2006), DOI 10.1016/j.jfluidstructs.2006.04.003, as read by
|
||||
Grok from the full text: Re = 200, `m* = m/m_f = 10` (m_f displaced fluid
|
||||
mass), `zeta = c/c_crit = 0.01`, `U* = U/(f_N D)` with `f_N` the natural
|
||||
frequency in fluid assuming an inviscid added-mass term; the paper also
|
||||
reports frequencies normalised by the in-vacuo value. Our runner derives `k`
|
||||
from that in-fluid `f_N` (added-mass coefficient as stated in the paper) and
|
||||
never from the vacuum formula; `c_crit` convention and the original
|
||||
comparison data are still OPEN. Stability/convergence is gated by fluid
|
||||
time/space refinement of the coupled run, not assumed from the mass ratio;
|
||||
structural sub-stepping does not substitute for fluid-coupling time
|
||||
refinement.
|
||||
|
||||
## 9. Gates (map; tolerances OPEN, declared before runs)
|
||||
|
||||
G0 zero-motion regression, split: G0a populations of the existing stationary
|
||||
and spin-only path bitwise unchanged vs main; G0b moving path at rest (with
|
||||
and without spin) reproduces the stationary link set and populations; force
|
||||
estimator outputs are compared separately from populations, since the moving
|
||||
path adds estimator terms · G1 raw force sign (`Cd > 0`) ·
|
||||
G2/G3 prescribed translation incl. ledger accumulators · G4 matched
|
||||
relative-flow (moving body in quiescent fluid vs fixed body in uniform
|
||||
stream), optional boost, co-moving sanity · G5 Dütsch et al. 1998 Re 100,
|
||||
KC 5 in-line (needs 2-component offset) · G6 Leontini free transverse VIV
|
||||
(m* = 10, ζ = 0.01, Re = 200, definitions pending) · G7 time-level contract
|
||||
frozen before the oscillator.
|
||||
|
||||
Candidate pre-run engineering tolerances (Astra; for discussion, not
|
||||
literature limits or measured outcomes; norms, averaging window, convergence
|
||||
protocol and physical reference uncertainty must be fixed before adoption, and
|
||||
none may be weakened after seeing runs):
|
||||
normalised momentum-ledger residual FP32 ≤ 5e-6 (FP64 oracle ≤ 2e-13);
|
||||
10-cycle relative mass drift ≤ 1e-4; prescribed-motion grid study D = 40/80/160
|
||||
with change in the force-response amplitude (observed quantity and norm to be
|
||||
named before adoption) ≤ 2 % and phase ≤ 2°; enlarged domain
|
||||
amplitude change ≤ 1 %; matched-frame force RMS ≤ 2 % and improving with
|
||||
refinement; VIV refinement amplitude ≤ 2 %, frequency ≤ 1 %, five consecutive
|
||||
cycle amplitudes within 1 %.
|
||||
|
||||
Conditions stated before any moving kernel (S2/S3), not only at the
|
||||
oscillator stage: (i) fluid refinement protocol for the prescribed-motion
|
||||
gates — diffusive scaling (dx halved, dt quartered, physical Re, KC and
|
||||
amplitude fixed), the observed quantities and norms named in advance;
|
||||
(ii) Galilean constant-boost G4-boost requires boundaries compatible with a
|
||||
uniform boost (periodic in the boost direction or boundaries that move with
|
||||
it); if the current dense configuration cannot represent that, G4-boost is
|
||||
recorded as not runnable rather than replaced by another test.
|
||||
|
||||
## 10. Proposed minimal S1 scope (host only; held until cross-review)
|
||||
|
||||
- `add_body` moving-circle declaration: base centre, declared offset box
|
||||
(x/y), speed bound; swept-band computation and band-placement validation.
|
||||
- `set_body(offset=..., velocity=..., omega=...)` validation before mutation
|
||||
(section 7), applied/pending state and the run-transition table of 4a,
|
||||
and the fail-closed commit of section 4.
|
||||
- Host enumeration check of refill-donor non-emptiness (section 5) within
|
||||
the declared `R_min` and displacement bound.
|
||||
- Motion-record host layout and upload path; exclusion of moving bodies from
|
||||
the packed `CurvedBoundaryKernel` links (stationary packing unchanged).
|
||||
- Explicit rejections for unsupported modes.
|
||||
- S1 execution guard: until the S2 moving kernels exist, any `run(n)` with
|
||||
`n > 0` (and `step(n)`) while a moving body is declared raises before any
|
||||
device mutation (no action/motion upload, no launch), including the
|
||||
nominal at-rest case, because the moving body is excluded from the static
|
||||
packed links and nothing else would apply its boundary. `run(0)` keeps its
|
||||
existing behaviour. Simulations without moving bodies are unaffected.
|
||||
Preparation/commit helpers are exercised by host tests with mocks; no
|
||||
operational unsupported stepping path is exposed.
|
||||
- Host tests (Astra runs them): rejection-before-mutation, band geometry,
|
||||
derived centres, fail-closed behaviour under an injected upload/sync
|
||||
failure, stationary packing byte-identical, execution guard (positive-step
|
||||
run with a moving body raises before upload/launch; stationary unaffected). No kernels, no GPU runs, no
|
||||
y-only VIV oscillator (S4), no force RHS.
|
||||
|
||||
## 10a. S2a — fixed-centre moving links (implemented on `claude/moving-bodies-s2a`)
|
||||
|
||||
- `MovingCircleLinks` (`kernels/step/moving_circle.cu`, dense double buffer)
|
||||
runs before the static curved kernel when moving tables are attached: one
|
||||
thread per (band cell, direction) slot, FP64 `{cx, cy, r}` at the applied
|
||||
offset, q/hit/clamp/fallback in the order of `build_cut_links` with explicit
|
||||
`_rn` operations (host/device bit parity is a test target, not a claim),
|
||||
then the unchanged `apply_bouzidi_link` with spin wall velocity only.
|
||||
- `obs` keeps its contract (CME and torque from `apply_bouzidi_link`). A
|
||||
separate per-body diagnostics buffer receives `cme`, the Wen term
|
||||
`u_w * dm`, `link_dm` and the link count from the same link outputs; it is
|
||||
zeroed with obs and only exists when moving circles do. Optional
|
||||
`slot_debug` (q, fallback, rx, ry per slot) is null in production.
|
||||
- Displacement is disabled: `set_body(offset=..., vx=..., vy=...)` on a
|
||||
translating circle raises before mutation, before and after initialize, so
|
||||
no pending target can exist and `run`/`step`/direct `LBMStepper.step` all
|
||||
hold at the applied geometry. `initialize()` never reads or consumes
|
||||
`pending`. `stepper.blocked` is cleared only when the moving tables are
|
||||
attached (initialize, stepper recreation).
|
||||
- S3 rule (not implemented): once targets are re-enabled, any pending target
|
||||
must raise in `run`/`step` and in `LBMStepper.step` alike before any launch,
|
||||
swap or step-count change; zero-step calls never consume it.
|
||||
- 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`) → synchronize. Hold steps launch no M1
|
||||
and are the S2a path. Commit point: the advancing step-group sync plus a
|
||||
zero `no_donor` count; a failure before it commits nothing. The hold-record
|
||||
upload stages the next step after the commit: if it (or its sync) fails,
|
||||
`applied` already matches the moved device geometry, but the Simulation is
|
||||
still marked failed and refuses reuse. Any failure from the first upload on
|
||||
marks the circles and the Simulation failed.
|
||||
- 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
|
||||
scheme; benchmark definitions and digitized reference data; tolerances;
|
||||
execution/resource feasibility (Astra, pending); whether a quiescent domain
|
||||
(zero inlet velocity) runs cleanly in the current config.
|
||||
@@ -0,0 +1,209 @@
|
||||
# Prescribed in-line oscillation: Dütsch Re 100 / KC 5 — PROTOCOL (review draft)
|
||||
|
||||
Production under test: `claude/moving-bodies-s3a` at `653fd997` (unchanged).
|
||||
Primary source: Dütsch, Durst, Becker & Lienhart, J. Fluid Mech. 360, 249–271
|
||||
(1998), doi:10.1017/S002211209800860X. Source facts below were verified from
|
||||
the paper (Slack C0C78UFKGRM/1793335557.460509, TRON 1793356317–1793356724,
|
||||
page numbers given). Everything marked ENGINEERING is our a-priori choice,
|
||||
fixed before any benchmark result; it is not a paper uncertainty.
|
||||
|
||||
## 0. Source facts (with provenance)
|
||||
|
||||
- Motion (p.251): `x_c = -A sin(2 pi f t)`; paper phase phi = 0 is the mean
|
||||
position moving left at `-U_max`. KC = 5 gives `A/D = 5/(2 pi) = .7957747`.
|
||||
- Start (p.254): the computation starts a quarter cycle before phi = 0, at
|
||||
`x = +A` with zero velocity (phi = -90°).
|
||||
- Experiment (pp.252, 254–256): D = 10 mm, L = .65 m, tank 1.2 × .7 × .7 m,
|
||||
rho = 998.2 kg/m³, mu = 1.00818 mPa·s, Re = 100, KC = 5; LDA with 1° phase
|
||||
bins and 200 samples per node and degree; no quantified experimental error
|
||||
recovered.
|
||||
- Profiles (Fig. 8, p.260): LAB-FIXED stations x1 = -6, 0, 6, 12 mm, i.e.
|
||||
x/D = -.6, 0, .6, 1.2, at phi = 180°, 210°, 330°; u1 and u2 versus x2.
|
||||
- Original numerics (§2): cylinder-fixed accelerated frame (`x_i = x~_i +
|
||||
x_i,s`, outer velocity `-u_c`, inertial source `-rho du_c/dt`), O-grid outer
|
||||
diameter 120 D, illustrated 192 circumferential nodes, 720 steps/cycle; tank
|
||||
influence neglected.
|
||||
- Force fit (verified): force ON THE BODY per unit span
|
||||
`F1 = -0.5 rho D c_d xdot|xdot| - (pi/4) rho D² c_i xddot` with `x = -A sin phi`;
|
||||
normalized by `0.5 rho D U_max²`: `c_d cos phi |cos phi| - (pi² c_i / KC) sin phi`.
|
||||
`c_i` is the ADDED-MASS coefficient; total inertia `c_m = c_i + 1` is a
|
||||
different convention and is not used.
|
||||
- Table 1 computed Morison coefficients (c_d, c_i) for 96 / 192 / 384 perimeter
|
||||
nodes at dt/T = 1/720: (2.15, 1.42), (2.10, 1.45), (2.09, 1.45). These are
|
||||
fitted computed coefficients, not a measured force and not an error bar. The
|
||||
paper states the coarse direct force is closer to the finest direct force
|
||||
than the Morison fit, but gives no waveform-error norm.
|
||||
|
||||
## 1. Kinematics, time levels and phase bookkeeping
|
||||
|
||||
- Our run is in the LAB frame: the circle moves, the fluid is initially at rest
|
||||
(rho = 1, u = 0); no fictitious inertial source is added. Lab-fixed profile
|
||||
stations are compared directly (no frame transform needed).
|
||||
- Geometry reference: the circle's base centre is the INITIAL position, the
|
||||
right turning point `x_mean + A`; the declared offset box is x in [-2A, 0],
|
||||
y = 0. The body is never teleported: offset 0 at s = 0.
|
||||
- Elapsed time s = n dt from the start. Analytic position relative to the mean
|
||||
`x(s) = A cos(2 pi s / T)`, velocity `u(s) = -U_max sin(2 pi s / T)`; paper
|
||||
phase `phi = (360 s / T - 90) mod 360`.
|
||||
- Step n moves from s_n to s_{n+1}: target offset `x(s_{n+1}) - A`; wall
|
||||
velocity = the single-step interval average `(x(s_{n+1}) - x(s_n)) / dt`
|
||||
(analytic, supplied, not differenced from stored offsets). The per-step
|
||||
force of `run(1)` / `read_force(normalize=True)` is attributed to that same
|
||||
interval.
|
||||
- Paper cycle k = 1, 2, ... spans `s in [T/4 + (k-1) T, T/4 + k T)`; the
|
||||
initial quarter cycle (s < T/4) is cycle 0 and never scored. Every logged
|
||||
sample carries (k, phi). Example: phi = 330° of cycle 1 is at s = 7T/6.
|
||||
- T is chosen as a multiple of 360 native steps, so phi = 180°, 210°, 330° fall
|
||||
exactly on step boundaries and 1° bins hold whole steps.
|
||||
|
||||
## 2. Resolution ladder (uniform grid, fixed Mach), ENGINEERING
|
||||
|
||||
| D (cells) | T (steps) | U_max = 5 D / T | nu_lat = U_max D / 100 | tau (3 nu + .5) |
|
||||
|---|---|---|---|---|
|
||||
| 20 | 2,520 | .039683 | .0079365 | .5238 |
|
||||
| 40 | 5,040 | .039683 | .015873 | .5476 |
|
||||
| 80 | 10,080 | .039683 | .031746 | .5952 |
|
||||
|
||||
- MRT collision. Mach fixed (Ma ≈ .0687); a separate Mach check at D = 40
|
||||
with T = 10,080 (U_max halved, tau .5238) is labelled MACH and not part of the
|
||||
ladder.
|
||||
- Per-step displacement ≤ U_max ≈ .04 cells (bound .5).
|
||||
|
||||
## 3. Domain and boundaries (finite-domain reproduction), ENGINEERING
|
||||
|
||||
The paper's 120 D O-grid is not reproduced; our rectangular box is a
|
||||
reproduction choice. The native dense channel (zero-velocity inlet,
|
||||
pressure-NEQ outlet, bounce-back y walls) is a finite-domain approximation,
|
||||
not infinite still fluid. The G4 x-periodic / y-slip override is a separately
|
||||
labelled sensitivity family and never mixed into the primary ladder. The
|
||||
runner owner confirms the zero-velocity inlet is supported as configured
|
||||
before any run.
|
||||
|
||||
Order (closed):
|
||||
1. Smoke + ledger at D = 40, native BC, 40 D × 40 D (finite box).
|
||||
2. Domain sensitivity before any literature statement: same BC, D = 40,
|
||||
L ∈ {20 D, 40 D, 60 D}; record the change in Fx waveform, c_d, c_i and the
|
||||
scored profiles.
|
||||
3. Optional labelled BC sensitivity: one x-periodic / y-slip case at D = 40,
|
||||
40 D; reported as a delta against native only.
|
||||
4. Only after step 2 passes: the D = 20 / 40 / 80 ladder on the chosen BC and
|
||||
domain (default 40 D if, against 60 D, `RMS(dFx)/F_ref` and |dc_d|, |dc_i|
|
||||
stay inside the section 6 bands).
|
||||
- Memory/runtime: measured by a smoke run before each new size; measurements
|
||||
inform scheduling only, never tolerances.
|
||||
|
||||
## 4. Observables (raw, per step, preserved)
|
||||
|
||||
- `obs` Fx, Fy (CME) and all `MovingTables.DIAG` columns per step: CME,
|
||||
CME + Wen and CME + cover − uncover reported separately; crossings, eq_only,
|
||||
no_donor (no_donor must stay 0).
|
||||
- Velocity u1, u2 on lab-fixed lines x/D = -.6, 0, .6, 1.2 (relative to the
|
||||
mean cylinder position) versus x2, at the scored phases. Lab frame
|
||||
throughout; if any body-frame intermediate is used, position and velocity
|
||||
are both transformed.
|
||||
- Sampling cost: raw force and diagnostics counters are read every step
|
||||
(small). Flags and velocity fields are read only at the steps inside the
|
||||
three prescribed 1° phase bins (φ ∈ [180, 181), [210, 211), [330, 331),
|
||||
T/360 steps each) of the cycles needed for settling and scoring, plus
|
||||
sparse declared audits; never every step.
|
||||
- Point set: fixed x2/D sample points per station, identical on every grid
|
||||
and declared before running; per point the interpolation weights are
|
||||
recorded; the scored set is the common subset valid on all compared grids
|
||||
at that phase, with its coverage recorded; point sets never change silently.
|
||||
- Interpolation: bilinear on cell centres using only cells that are fluid at
|
||||
that step (flags read back via `get_flags()`); weights renormalized over the
|
||||
fluid corners; a sample with no fluid corner, or inside the circle at that
|
||||
step, is MISSING and logged, never set to zero. One-degree bins keep valid
|
||||
fluid weights. The numerical average need not mimic the 200 LDA samples.
|
||||
|
||||
## 5. Periodic-regime gate (every case, before any statistic), ENGINEERING
|
||||
|
||||
Ledger rows, same schema for every case: dx, dt (in D and T units); U_max, nu,
|
||||
Re, KC (lattice and nondimensional); T (native steps, nondimensional);
|
||||
duration (native steps, s/T, paper cycles); discard; accepted window; cycle
|
||||
ledger; verdict PASS or FAIL-CLOSED.
|
||||
|
||||
Scales: `F_ref = 0.5 rho U_max² D`, `U_ref = U_max`; never a near-zero mean
|
||||
force as denominator.
|
||||
|
||||
- Settling metrics, cycle k−1 → k, phase-aligned by step within the cycle:
|
||||
- waveform `RMS(dFx) / F_ref ≤ .02`;
|
||||
- `|d mean|`, `|d RMS|`, `|d peak|` of Fx, each `/ F_ref ≤ .01`;
|
||||
- scored profiles (fluid-only) `|du| / U_max`, `|dv| / U_max ≤ .03` at the
|
||||
three phases × four stations;
|
||||
- Morison `|dc_d| ≤ .02`, `|dc_i| ≤ .02`.
|
||||
- Rule: all metrics for N = 3 consecutive whole paper cycles → settled.
|
||||
|
||||
Cycle indexing (worked example, binding). Paper cycle k spans
|
||||
`s/T in [0.25 + (k-1), 0.25 + k)`, i.e. φ0 → φ0; cycle 0 is the startup
|
||||
quarter `s/T < 0.25`. Pair p_k compares cycle k−1 with cycle k.
|
||||
|
||||
- Discarded: cycles 0, 1, 2 (s/T < 2.25). Pairs are formed only between
|
||||
cycles ≥ 3, so the first pair is p_4 (3 → 4).
|
||||
- Settling proof for a case: three consecutive pairs all pass; K_case = the
|
||||
first cycle AFTER the last of those pairs. Earliest possible: pairs 3→4,
|
||||
4→5, 5→6 pass ⇒ K_case = 7.
|
||||
- Common window: K = max over compared cases of K_case; every case scores the
|
||||
same three complete cycles K, K+1, K+2 and must stay settled throughout:
|
||||
pairs (K−1→K), (K→K+1), (K+1→K+2) must pass on every case, otherwise that
|
||||
case is FAIL-CLOSED (an earlier local settle does not waive this).
|
||||
- Cap: cycle k ends at s/T = k + 0.25, so the last admissible scored cycle is
|
||||
11 (ends 11.25 T) and K ≤ 9. Examples: K = 7 → cycles 7–9 (ends 9.25 T);
|
||||
K = 9 → cycles 9–11 (ends 11.25 T); any K_case > 9 → FAIL-CLOSED.
|
||||
- Every case is run to the end of cycle K + 2 (proportionally more native
|
||||
steps on finer grids).
|
||||
- Discard: the initial quarter cycle plus the first 2 full paper cycles
|
||||
(s/T < 2.25). Settling is searched from cycle 3; cap 12 periods from s = 0;
|
||||
unsettled → FAIL-CLOSED (no literature, cross-resolution or VIV statistic).
|
||||
- Common window: K = the largest first-settled cycle index over all compared
|
||||
cases; every case is scored on the same 3 cycles K .. K + 2 and must reach
|
||||
them within the cap (proportionally more native steps on finer grids), or
|
||||
the comparison is FAIL-CLOSED.
|
||||
- The earlier G4 timing check (100 coarse vs 400 diffusive-fine steps) was a
|
||||
matched short transient, not periodic evidence.
|
||||
|
||||
## 6. Metrics and per-observable criteria (ENGINEERING, not paper uncertainty)
|
||||
|
||||
Three criteria stay separate; the G4 same-grid 2% target is not reused.
|
||||
|
||||
- Morison fit: least squares of the normalized per-step force over the
|
||||
accepted window against the two basis functions, each averaged over its
|
||||
single step interval, fixed paper phase, no free phase shift; c_d and
|
||||
added-mass c_i per force candidate.
|
||||
- Cross-resolution: every scored observable's error against its reference must
|
||||
decrease monotonically D20 → D40 → D80 (necessary, not sufficient).
|
||||
- Against Table 1 (computed fit, not measured force): at D ≥ 40,
|
||||
`|c_d - 2.09| ≤ .15` and `|c_i - 1.45| ≤ .15`, and monotone.
|
||||
- Force waveform against its own Morison reconstruction: residual reported
|
||||
only, no pass/fail (the paper reports fit error above its coarse numerical
|
||||
error and gives no percentage).
|
||||
- Fig. 8 LDA profiles: QUALITATIVE ONLY. The original figure could not be
|
||||
accessed or rendered; no digitized ordinates and no measurement uncertainty
|
||||
exist. After a digitization with figure number, tool and uncertainty note:
|
||||
`RMS(|du|, |dv|) / U_max ≤ .15` at D ≥ 40 and monotone; LDA error stays
|
||||
UNQUANTIFIED and is stated with every comparison.
|
||||
- Fy sanity (not literature): `|mean Fy| / F_ref ≤ .05` in the accepted window.
|
||||
- No force candidate is selected for an RHS by these results; each is
|
||||
reported.
|
||||
|
||||
Superseded (earlier drafts): the velocity-jump start, guessed stations and
|
||||
phases, "C_M" for c_i, and the uncertainty wording; cycles 4–6, the 30-cycle
|
||||
cap and the earlier tolerance numbers.
|
||||
|
||||
## 6a. Driver unit tests (CPU, before any run)
|
||||
|
||||
- Phase/cycle arithmetic: s ↔ (k, φ), cycle boundaries, bin membership,
|
||||
window/cap selection for K = 7, 8, 9, 10.
|
||||
- Solid-aware interpolation: all-fluid, partial-solid (renormalized weights),
|
||||
no-fluid (MISSING) and inside-circle cases.
|
||||
- Interval-averaged Morison bases and the two-parameter fit on a synthetic
|
||||
waveform with known c_d, c_i.
|
||||
|
||||
## 7. Order and scope
|
||||
|
||||
Protocol review → one-shot Dütsch driver (analytic `set_body` + `run(1)`,
|
||||
cycle ledger, lab-fixed fluid-only station sampler, two-basis Morison fit; no
|
||||
reusable runner until the first smoke shows shared helpers pay off; no
|
||||
production change) → smoke (memory/runtime) → domain sensitivity →
|
||||
ladder → physical-validation verdict → only then single-cylinder VIV. No AMR,
|
||||
no pinball, no manuscript work.
|
||||
@@ -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
|
||||
|
||||
@@ -0,0 +1,259 @@
|
||||
# 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,
|
||||
on_failure: Optional[Callable[[BaseException], None]] = None):
|
||||
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
|
||||
self._on_failure = on_failure # owner lifecycle hook (Simulation marks itself failed)
|
||||
|
||||
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, 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:
|
||||
"""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
|
||||
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
|
||||
|
||||
|
||||
class MovingTables:
|
||||
"""Device tables of the translating circles for ``MovingCircleCells``/``Links``.
|
||||
|
||||
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, 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) // 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))
|
||||
# Owned page-locked staging for async record uploads (allocated once).
|
||||
self._geom_pinned = cuda.pagelocked_empty(len(geom), np.float64)
|
||||
self.diag_gpu = cuda.to_device(np.zeros(self.n_bodies * self.DIAG_FLOATS, np.float32))
|
||||
self.slot_debug_gpu = None # optional (q, fallback, rx, ry) per slot for parity tests
|
||||
|
||||
@classmethod
|
||||
def upload(cls, circles):
|
||||
"""``circles``: [(MovingCircle, body_id)] in body order."""
|
||||
return cls(*cls.layout(circles))
|
||||
|
||||
@staticmethod
|
||||
def layout(circles):
|
||||
"""Host band and geometry tables (no device work)."""
|
||||
band, geom, first = [], [], 0
|
||||
for motion, body_id in circles:
|
||||
x0, x1, y0, y1 = motion.band()
|
||||
width, height = x1 - x0 + 1, y1 - y0 + 1
|
||||
band += [x0, y0, width, height, first, int(body_id)]
|
||||
first += width * height * 8
|
||||
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:
|
||||
"""Enqueue the records when they differ from the device copy; True if advancing.
|
||||
|
||||
The pinned staging buffer is reused: callers synchronize ``stream`` after every
|
||||
upload before the next one (``Simulation._run_moving``), so no copy is in flight
|
||||
when it is rewritten. ``geom`` is the enqueued state; after a failed copy the
|
||||
owner is failed and the cache is never trusted again.
|
||||
"""
|
||||
import pycuda.driver as cuda
|
||||
geom = [v for motion, _ in circles for v in self.record(motion)]
|
||||
if geom != self.geom:
|
||||
self._geom_pinned[:] = geom
|
||||
cuda.memcpy_htod_async(self.geom_gpu, self._geom_pinned, 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
|
||||
self.slot_debug_gpu = cuda.to_device(np.zeros(self.n_slots * 4, np.float32))
|
||||
|
||||
def zero_diagnostics(self, stream):
|
||||
import pycuda.driver as cuda
|
||||
cuda.memset_d32_async(self.diag_gpu, 0, self.n_bodies * self.DIAG_FLOATS, stream)
|
||||
|
||||
def download_diagnostics(self):
|
||||
import numpy as np
|
||||
import pycuda.driver as cuda
|
||||
out = np.empty((self.n_bodies, self.DIAG_FLOATS), np.float32)
|
||||
cuda.memcpy_dtoh(out, self.diag_gpu)
|
||||
return out
|
||||
@@ -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):
|
||||
|
||||
@@ -0,0 +1,183 @@
|
||||
// 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
|
||||
// only writer of their solid slots.
|
||||
#ifndef CELERIS_MOVING_CIRCLE_CU
|
||||
#define CELERIS_MOVING_CIRCLE_CU
|
||||
|
||||
#if DIM == 2
|
||||
constexpr int MOVING_BAND_INTS = 6; // x0, y0, width, height, first_slot, body_id
|
||||
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); }
|
||||
|
||||
__device__ __forceinline__ bool mc_inside(int x, int y, double cx, double cy, double r)
|
||||
{
|
||||
return __dadd_rn(mc_sq(__dsub_rn((double)x, cx)), mc_sq(__dsub_rn((double)y, cy))) < mc_sq(r);
|
||||
}
|
||||
|
||||
// Smallest root t in [0, 1] of |p + t c - centre| = r (find_circle_intersection order).
|
||||
__device__ __forceinline__ bool mc_hit(int x, int y, int ex, int ey,
|
||||
double cx, double cy, double r, double& t)
|
||||
{
|
||||
const double dx = (double)ex, dy = (double)ey;
|
||||
const double px = __dsub_rn((double)x, cx), py = __dsub_rn((double)y, cy);
|
||||
const double a = __dadd_rn(mc_sq(dx), mc_sq(dy));
|
||||
const double b = __dmul_rn(2.0, __dadd_rn(__dmul_rn(dx, px), __dmul_rn(dy, py)));
|
||||
const double c = __dsub_rn(__dadd_rn(mc_sq(px), mc_sq(py)), mc_sq(r));
|
||||
const double det = __dsub_rn(mc_sq(b), __dmul_rn(__dmul_rn(4.0, a), c));
|
||||
if (det < 0.0) return false;
|
||||
const double s = __dsqrt_rn(det), a2 = __dmul_rn(2.0, a);
|
||||
const double t1 = __ddiv_rn(__dsub_rn(-b, s), a2);
|
||||
const double t2 = __ddiv_rn(__dadd_rn(-b, s), a2);
|
||||
if (0.0 <= t1 && t1 <= 1.0) { t = t1; return true; }
|
||||
if (0.0 <= t2 && t2 <= 1.0) { t = t2; return true; }
|
||||
return false;
|
||||
}
|
||||
|
||||
__global__ void MovingCircleLinks(
|
||||
fpxx* fi, const uint16_t* flag,
|
||||
const int* band, const double* geom,
|
||||
unsigned int n_bodies, unsigned int n_slots,
|
||||
const float* action, float* obs, float* diag, float* slot_debug)
|
||||
{
|
||||
const unsigned int tid = threadIdx.x + blockIdx.x * blockDim.x;
|
||||
if (tid >= n_slots) return;
|
||||
unsigned int m = 0;
|
||||
while (m + 1 < n_bodies && tid >= (unsigned int)band[(m + 1) * MOVING_BAND_INTS + 4]) ++m;
|
||||
const int* b = band + m * MOVING_BAND_INTS;
|
||||
const unsigned int local = tid - (unsigned int)b[4];
|
||||
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* 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;
|
||||
|
||||
double t;
|
||||
if (mc_inside(x, y, cx, cy, r) || !mc_inside(x + ex, y + ey, cx, cy, r)
|
||||
|| !mc_hit(x, y, ex, ey, cx, cy, r, t))
|
||||
return;
|
||||
const double hx = __dadd_rn((double)x, __dmul_rn(t, (double)ex));
|
||||
const double hy = __dadd_rn((double)y, __dmul_rn(t, (double)ey));
|
||||
const double hdx = __dsub_rn(hx, (double)x), hdy = __dsub_rn(hy, (double)y);
|
||||
const double norm = __dsqrt_rn((double)(ex * ex + ey * ey));
|
||||
double qd = __ddiv_rn(__dsqrt_rn(__dadd_rn(mc_sq(hdx), mc_sq(hdy))), norm);
|
||||
qd = fmax(0.01, fmin(qd, 1.0));
|
||||
const float q = (float)qd;
|
||||
const float rx = (float)__dsub_rn(hx, cx), ry = (float)__dsub_rn(hy, cy);
|
||||
|
||||
// Same two-stage order as build_cut_links + pack_cut_links_to_soa: the fallback
|
||||
// class is chosen from the clamped FP64 q (q < 1/2 needs a fluid donor at x - c);
|
||||
// only then are q, rx, ry rounded to FP32 for the BFL branch and spin velocity.
|
||||
unsigned char fallback = CURVED_FALLBACK_BOUZIDI;
|
||||
if (qd < 0.5) {
|
||||
const int xd = x - ex, yd = y - ey;
|
||||
if (xd < 0 || xd >= (int)NX || yd < 0 || yd >= (int)NY || mc_inside(xd, yd, cx, cy, r)
|
||||
|| (flag[linear_index((unsigned int)xd, (unsigned int)yd)] & FLAG_FLUID) == 0)
|
||||
fallback = CURVED_FALLBACK_HALFWAY;
|
||||
}
|
||||
if (debug != nullptr) { debug[0] = q; debug[1] = (float)fallback; debug[2] = rx; debug[3] = ry; }
|
||||
|
||||
const int bid = b[5];
|
||||
const float omega = action_omega(action, bid);
|
||||
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);
|
||||
if (diag != nullptr) {
|
||||
float* d = diag + m * MOVING_DIAG_FLOATS;
|
||||
atomicAdd(&d[0], exchange[0]);
|
||||
atomicAdd(&d[1], exchange[1]);
|
||||
atomicAdd(&d[2], Uw * dm); // Wen 2014 correction -u_w (f_t - f_r)
|
||||
atomicAdd(&d[3], Vw * dm);
|
||||
atomicAdd(&d[4], dm);
|
||||
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
|
||||
@@ -71,6 +71,7 @@ extern "C"
|
||||
#include "step/one_step_double.cu"
|
||||
#include "step/one_step_esopull.cu"
|
||||
#include "step/aux_kernels.cu"
|
||||
#include "step/moving_circle.cu"
|
||||
#include "step/esopull_macro.cu"
|
||||
#include "step/esopull_sync.cu"
|
||||
|
||||
|
||||
@@ -20,6 +20,10 @@ class LBMStepper:
|
||||
self.cfg = cfg
|
||||
|
||||
self._esopull = (cfg.streaming == "esopull")
|
||||
# Set by the owner when a step would be wrong (translating circles
|
||||
# excluded from the packed links before the moving kernels exist).
|
||||
self.blocked: str | None = None
|
||||
self.moving = None # MovingTables attached by the owner (translating circles)
|
||||
|
||||
if self._esopull:
|
||||
self.step_fn = module.get_function("EsoPullStep")
|
||||
@@ -35,6 +39,8 @@ class LBMStepper:
|
||||
self.apply_esopull_curved_fn = module.get_function("ApplyEsoPullCurvedKernel")
|
||||
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)
|
||||
@@ -70,6 +76,13 @@ class LBMStepper:
|
||||
stream: cuda.Stream | None = None,
|
||||
):
|
||||
"""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:
|
||||
@@ -81,6 +94,8 @@ class LBMStepper:
|
||||
f.ddf_gpu, self.boundary_values_gpu, np.uint64(self._step_count),
|
||||
block=self.block, grid=self.boundary_grid, **launch_kw,
|
||||
)
|
||||
if self.moving is not None:
|
||||
self._launch_moving(action_gpu, obs_gpu, advance, **launch_kw)
|
||||
self._launch_curved(action_gpu, obs_gpu, **launch_kw)
|
||||
|
||||
if self._esopull:
|
||||
@@ -106,6 +121,21 @@ class LBMStepper:
|
||||
# Any completed GPU step invalidates host DDF mirror.
|
||||
f.invalidate_host_ddf_cache()
|
||||
|
||||
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,
|
||||
m.diag_gpu, m.slot_debug_gpu if m.slot_debug_gpu is not None else np.uintp(0),
|
||||
block=(tpb, 1, 1), grid=((m.n_slots + tpb - 1) // tpb, 1, 1), **launch_kw,
|
||||
)
|
||||
|
||||
def _launch_curved(self, action_gpu, obs_gpu, **launch_kw):
|
||||
f = self.field
|
||||
if f.n_curved == 0:
|
||||
|
||||
@@ -26,6 +26,13 @@ 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
|
||||
|
||||
# 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")
|
||||
from .lbm.field import LBMField
|
||||
from .lbm.stepper import LBMStepper
|
||||
from .body.objects import SimObject
|
||||
@@ -97,7 +104,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 +121,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 +134,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 +146,11 @@ 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,
|
||||
on_failure=self._fail_moving)
|
||||
self.stepper.blocked = MOVING_BLOCKED
|
||||
elif type_lower == "sensor":
|
||||
from .body.geometry.circle import CircleGeometry
|
||||
cx, cy = float(center[0]), float(center[1])
|
||||
@@ -190,6 +210,78 @@ class Simulation:
|
||||
raise
|
||||
self.bodies.jets[index] = jet # later geometry syncs keep the new amplitude
|
||||
|
||||
def _fail_moving(self, exc: BaseException) -> None:
|
||||
# Device geometry is unknown after a failed moving step: same lifecycle as Jets.
|
||||
self._failed = f"moving circle step failed: {exc!r}"
|
||||
self.stepper.blocked = self._failed # direct stepper calls are refused as well
|
||||
|
||||
def _attach_moving(self) -> None:
|
||||
# Direct stepper calls stay blocked until the moving tables exist.
|
||||
tables = getattr(self, "_moving_tables", None)
|
||||
self.stepper.moving = tables
|
||||
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()
|
||||
if not moving:
|
||||
return None
|
||||
from .body.motion import MovingTables
|
||||
return MovingTables.upload([(obj.motion, obj.obj_id) for obj in moving])
|
||||
|
||||
@context_bound
|
||||
def moving_diagnostics(self) -> np.ndarray:
|
||||
"""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:
|
||||
return np.zeros((0, 6), np.float32)
|
||||
self.stream.synchronize()
|
||||
return tables.download_diagnostics()
|
||||
|
||||
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 +352,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 +362,21 @@ 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). 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 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 +392,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()
|
||||
|
||||
@@ -443,6 +551,7 @@ class Simulation:
|
||||
self.field, self._module, self.lbm_cfg,
|
||||
)
|
||||
self.stepper._step_count = _prev_step_count
|
||||
self._attach_moving()
|
||||
|
||||
self._assert_object_count_contract(expected_count=expected_n_objects)
|
||||
|
||||
@@ -468,6 +577,8 @@ 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()
|
||||
self._attach_moving()
|
||||
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)
|
||||
@@ -484,6 +595,8 @@ class Simulation:
|
||||
self.bodies.sync_to_gpu(self.field, rebuild_flags=False)
|
||||
else:
|
||||
self.bodies.ensure_packed_buffers()
|
||||
self._moving_tables = self._build_moving_tables()
|
||||
self._attach_moving()
|
||||
self._assert_runtime_contracts()
|
||||
self._initialized = True
|
||||
|
||||
@@ -499,6 +612,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 +713,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")
|
||||
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():
|
||||
@@ -615,6 +736,8 @@ class Simulation:
|
||||
# Zero obs segments before step group
|
||||
if zero_obs:
|
||||
self.bodies.zero_obs_async(stream)
|
||||
if self.stepper.moving is not None:
|
||||
self.stepper.moving.zero_diagnostics(stream)
|
||||
|
||||
if checkpoint_interval > 0:
|
||||
done = 0
|
||||
@@ -631,6 +754,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,
|
||||
@@ -649,6 +775,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)
|
||||
@@ -666,8 +820,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
|
||||
@@ -698,6 +860,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 +871,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,
|
||||
@@ -726,6 +890,7 @@ class Simulation:
|
||||
cuda.Context.synchronize()
|
||||
# Release owned GPU resources while their context is current.
|
||||
self.stepper = None
|
||||
self._moving_tables = None
|
||||
self.bodies = None
|
||||
self.field = None
|
||||
self._module = None
|
||||
|
||||
@@ -0,0 +1,536 @@
|
||||
"""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 Owner, make_simulation # noqa: E402
|
||||
|
||||
|
||||
class ClosableOwner(Owner):
|
||||
"""Owner double that refuses activation after close, like CudaContext.activate."""
|
||||
def activate(self):
|
||||
if self.is_closed:
|
||||
raise RuntimeError("CUDA context is closed")
|
||||
return super().activate()
|
||||
|
||||
|
||||
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=[], action_gpu=None, obs_gpu=None,
|
||||
has_pending_edit=lambda: False,
|
||||
_upload_action_async=lambda s: writes.append("action"),
|
||||
zero_obs_async=lambda s: writes.append("obs"),
|
||||
increment_obs_steps=lambda n: writes.append("count"),
|
||||
download_obs_full_async=lambda s: writes.append("download"))
|
||||
sim.stepper = SimpleNamespace(step=lambda n, **k: writes.append(("step", n)), blocked=None,
|
||||
moving=None)
|
||||
sim._assert_runtime_contracts = lambda: None
|
||||
return sim, writes
|
||||
|
||||
|
||||
STREAM = SimpleNamespace(synchronize=lambda: None)
|
||||
|
||||
|
||||
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))
|
||||
sim.bodies.objects[0].motion.applied = (0., 0.)
|
||||
for call in (lambda: sim.run(1, stream=STREAM), lambda: sim.run(5, stream=STREAM),
|
||||
lambda: sim.step(), lambda: sim.run(3, stream=STREAM, 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)
|
||||
sim.bodies.objects[0].state = SimpleNamespace(omega=.02) # spin does not lift the guard
|
||||
with pytest.raises(RuntimeError, match="moving kernels"):
|
||||
sim.run(2, stream=STREAM)
|
||||
assert writes == [] and sim.bodies.objects[0].motion.applied == (0., 0.)
|
||||
|
||||
|
||||
def test_run_zero_keeps_pending_and_static_runs_proceed(monkeypatch):
|
||||
sim, writes = moving_simulation(monkeypatch)
|
||||
motion = sim.bodies.objects[0].motion
|
||||
motion.set_target((0., .1), (0., .01))
|
||||
sim.run(0, stream=STREAM) # existing action/obs work only; no motion upload or commit
|
||||
assert motion.pending == ((0., .1), (0., .01)) and motion.applied == (0., 0.)
|
||||
assert ("step", 0) in writes
|
||||
static, writes = moving_simulation(monkeypatch, moving=False)
|
||||
static.save_checkpoint = lambda *a: writes.append("checkpoint")
|
||||
static.run(3, stream=STREAM)
|
||||
static.run(4, stream=STREAM, checkpoint_interval=2)
|
||||
assert writes.count(("step", 3)) == 1 and writes.count(("step", 2)) == 2
|
||||
|
||||
|
||||
def test_direct_stepper_entry_is_guarded():
|
||||
from CelerisLab.lbm.stepper import LBMStepper
|
||||
launches = []
|
||||
module = SimpleNamespace(get_function=lambda name: lambda *a, **k: launches.append(name))
|
||||
cfg = SimpleNamespace(streaming="double_buffer", threads_per_block=128, is_d2q9=True,
|
||||
nx=64, ny=32, nz=1)
|
||||
field = SimpleNamespace(ddf_gpu="ddf", temp_gpu="temp", n_curved=0)
|
||||
stepper = LBMStepper(field, module, cfg)
|
||||
assert stepper.blocked is None
|
||||
stepper.blocked = "translating circles need the moving kernels"
|
||||
for n in (1, 3):
|
||||
with pytest.raises(RuntimeError, match="moving kernels"):
|
||||
stepper.step(n, action_gpu=None, obs_gpu=None)
|
||||
assert launches == [] and (field.ddf_gpu, field.temp_gpu) == ("ddf", "temp")
|
||||
assert stepper.step_count == 0
|
||||
stepper.step(0, action_gpu=None, obs_gpu=None) # zero steps stay a no-op
|
||||
assert launches == [] and stepper.step_count == 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()
|
||||
|
||||
|
||||
# -- Real-object routing, initialize rejections, packing, lifecycle -------------
|
||||
import numpy as np # noqa: E402
|
||||
|
||||
from CelerisLab.body.geometry.circle import CircleGeometry # noqa: E402
|
||||
from CelerisLab.body.manager import ObjectManager # noqa: E402
|
||||
from CelerisLab.body.objects import SimObject # noqa: E402
|
||||
from CelerisLab.lbm.descriptors import FLUID, SOLID # noqa: E402
|
||||
|
||||
|
||||
class Cfg:
|
||||
dim, is_d3q19, is_d2q9 = 2, False, True
|
||||
streaming, ddf_shifting, store_precision, y_wall_bc = "double_buffer", False, "FP32", "bounce_back"
|
||||
|
||||
|
||||
def body(x, y, r, translation=None, **kw):
|
||||
obj = SimObject(obj_id=-1, geometry=CircleGeometry(x, y, r), center=(x, y), radius=r, **kw)
|
||||
if translation is not None:
|
||||
obj.motion = MovingCircle((x, y), r, translation)
|
||||
return obj
|
||||
|
||||
|
||||
def routed_simulation(monkeypatch):
|
||||
sim = make_simulation(monkeypatch)
|
||||
sim._failed, sim._initialized, sim.lbm_cfg = None, False, Cfg()
|
||||
sim.stepper = SimpleNamespace(blocked=None)
|
||||
sim.bodies = ObjectManager(nx=120, ny=80, nz=1, nq=9, cfg=Cfg())
|
||||
return sim
|
||||
|
||||
|
||||
def test_add_body_and_set_body_route_targets(monkeypatch):
|
||||
sim = routed_simulation(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
|
||||
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 == .01 and (obj.state.vx, obj.state.vy) == (0., 0.)
|
||||
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
|
||||
assert sim.bodies.get(still).state.vx == .01 and sim.bodies.get(still).motion is None
|
||||
for kwargs in (dict(type="sensor"), dict(surface_actuation=dict(angle=0., half_width=.2,
|
||||
normal_velocity=.01))):
|
||||
kind = kwargs.pop("type", "circle")
|
||||
with pytest.raises(ValueError, match="translation"):
|
||||
sim.add_body(kind, (20., 20.), 5., translation=Translation((0., 0., 0., 0.), .1),
|
||||
**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))
|
||||
uploads = sum(1 for w in self.log if isinstance(w, tuple) and w[0] == "upload")
|
||||
if self.fail == "upload" or (self.fail == "hold_upload" and uploads == 2):
|
||||
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)))
|
||||
from CelerisLab.lbm.stepper import LBMStepper
|
||||
sim.stepper.step = LBMStepper.step.__get__(sim.stepper) # real public entry
|
||||
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", "first_sync", "no_donor", "hold_upload",
|
||||
"hold_sync"))
|
||||
def test_failed_advance_fails_simulation(monkeypatch, failure):
|
||||
kwargs = dict(no_donor=1.) if failure == "no_donor" else dict(fail=failure)
|
||||
sim, writes, motion = advancing_simulation(monkeypatch, **kwargs)
|
||||
motion._on_failure = sim._fail_moving
|
||||
motion.set_target((0., .25), (0., .05))
|
||||
syncs = []
|
||||
|
||||
def synchronize(): # stands in for an asynchronous failure surfacing at sync
|
||||
syncs.append(1)
|
||||
if (failure, len(syncs)) in (("first_sync", 1), ("hold_sync", 2)):
|
||||
raise RuntimeError("sync failed")
|
||||
with pytest.raises(RuntimeError):
|
||||
sim.run(1, stream=SimpleNamespace(synchronize=synchronize))
|
||||
# Before the step-group sync and no_donor check nothing is committed; the hold
|
||||
# record is staged after the commit (applied already matches the moved device).
|
||||
committed = failure in ("hold_upload", "hold_sync")
|
||||
assert motion.applied == ((0., .25) if committed else (0., 0.))
|
||||
assert (motion.pending is None) == committed
|
||||
assert motion.failed and sim._failed is not None
|
||||
assert sim.stepper.blocked == sim._failed # direct stepper reuse refused too
|
||||
with pytest.raises(RuntimeError, match="Simulation failed"):
|
||||
sim.run(1, stream=STREAM)
|
||||
|
||||
|
||||
def test_failed_hold_step_blocks_the_direct_stepper(monkeypatch):
|
||||
# No target is pending, so only the failure latch can refuse direct steps.
|
||||
from CelerisLab.lbm.stepper import LBMStepper
|
||||
sim, writes, motion = advancing_simulation(monkeypatch, fail="upload")
|
||||
motion._on_failure = sim._fail_moving
|
||||
launches = []
|
||||
kernels = SimpleNamespace(get_function=lambda name: lambda *a, **k: launches.append(name))
|
||||
cfg = SimpleNamespace(streaming="double_buffer", threads_per_block=128, is_d2q9=True,
|
||||
nx=64, ny=32, nz=1)
|
||||
field = SimpleNamespace(ddf_gpu="ddf", temp_gpu="temp", flag_gpu="flag", n_curved=0)
|
||||
stepper = LBMStepper(field, kernels, cfg)
|
||||
stepper.moving = sim.stepper.moving
|
||||
sim.stepper = stepper
|
||||
assert motion.pending is None and stepper.blocked is None
|
||||
with pytest.raises(RuntimeError, match="upload failed"):
|
||||
sim.run(2, stream=STREAM) # hold step
|
||||
assert sim._failed is not None and stepper.blocked == sim._failed
|
||||
with pytest.raises(RuntimeError, match="moving circle step failed"):
|
||||
stepper.step(1, action_gpu=None, obs_gpu=None)
|
||||
assert launches == [] and stepper.step_count == 0
|
||||
assert (field.ddf_gpu, field.temp_gpu) == ("ddf", "temp")
|
||||
|
||||
|
||||
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))
|
||||
b = MovingCircle((90.25, 30.5), 5., Translation((0., 0., 0., 0.), .1))
|
||||
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., 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"
|
||||
sim._attach_moving()
|
||||
assert sim.stepper.moving == "tables" and sim.stepper.blocked is None
|
||||
sim._moving_tables = None
|
||||
sim._attach_moving()
|
||||
assert sim.stepper.moving is None and sim.stepper.blocked is not None
|
||||
|
||||
|
||||
class CountingOwner(Owner):
|
||||
"""Owner double that counts nested activations like CudaContext."""
|
||||
def __init__(self):
|
||||
super().__init__()
|
||||
self.active = 0
|
||||
|
||||
def activate(self):
|
||||
from contextlib import contextmanager
|
||||
|
||||
@contextmanager
|
||||
def scope():
|
||||
self.active += 1
|
||||
try:
|
||||
yield self
|
||||
finally:
|
||||
self.active -= 1
|
||||
return scope()
|
||||
|
||||
|
||||
def test_public_diagnostics_keep_owner_active_through_sync_and_download(monkeypatch):
|
||||
sim = make_simulation(monkeypatch)
|
||||
sim.ctx, sim._failed = CountingOwner(), None
|
||||
seen = []
|
||||
sim._stream = SimpleNamespace(synchronize=lambda: seen.append(("sync", sim.ctx.active)))
|
||||
expected = np.arange(6, dtype=np.float32).reshape(1, 6)
|
||||
|
||||
def download():
|
||||
seen.append(("download", sim.ctx.active))
|
||||
return expected
|
||||
sim.stepper = SimpleNamespace(moving=SimpleNamespace(download_diagnostics=download))
|
||||
assert sim.moving_diagnostics() is expected
|
||||
assert [name for name, _ in seen] == ["sync", "download"]
|
||||
assert all(active > 0 for _, active in seen) and sim.ctx.active == 0
|
||||
|
||||
|
||||
def test_stepper_launches_moving_links_before_static_curved():
|
||||
from CelerisLab.lbm.stepper import LBMStepper
|
||||
launches = []
|
||||
module = SimpleNamespace(get_function=lambda name: lambda *a, **k: launches.append((name, a, k)))
|
||||
cfg = SimpleNamespace(streaming="double_buffer", threads_per_block=128, is_d2q9=True,
|
||||
nx=64, ny=32, nz=1)
|
||||
curved = SimpleNamespace(**{f"{name}_gpu": name for name in (
|
||||
"fluid_idx", "dir", "q", "rx", "ry", "rz", "fallback_class", "body_id")},
|
||||
surface_velocity=None, surface_velocity_gpu=None, link_exchange_gpu=None)
|
||||
field = SimpleNamespace(ddf_gpu="ddf", temp_gpu="temp", flag_gpu="flag", n_curved=5,
|
||||
curved=curved, invalidate_host_ddf_cache=lambda: None)
|
||||
stepper = LBMStepper(field, module, cfg)
|
||||
stepper._launch_force_region = stepper._launch_sensor = lambda *a, **k: None
|
||||
stepper.moving = SimpleNamespace(band_gpu="band", geom_gpu="geom", diag_gpu="diag",
|
||||
slot_debug_gpu=None, n_bodies=1, n_slots=300)
|
||||
stepper.step(2, action_gpu="act", obs_gpu="obs")
|
||||
names = [name for name, _, _ in launches]
|
||||
assert names == ["MovingCircleLinks", "CurvedBoundaryKernel", "OneStep"] * 2
|
||||
assert stepper.step_count == 2
|
||||
assert launches[1][1][0] == "ddf" and launches[3][1][0] == "temp" # swap between steps
|
||||
args, kw = launches[0][1], launches[0][2]
|
||||
assert args[:4] == ("ddf", "flag", "band", "geom") and int(args[5]) == 300
|
||||
assert kw["grid"] == (3, 1, 1) and kw["block"] == (128, 1, 1)
|
||||
|
||||
|
||||
def test_recompile_installs_a_guarded_stepper(monkeypatch):
|
||||
from CelerisLab.lbm.stepper import LBMStepper
|
||||
sim = routed_simulation(monkeypatch)
|
||||
sim.add_body("circle", (50., 40.), 10., translation=Translation((0., 0., -1., 1.), .1))
|
||||
launches = []
|
||||
kernel_module = SimpleNamespace(get_function=lambda name: lambda *a, **k: launches.append(name))
|
||||
monkeypatch.setattr(module.compiler, "build_module", lambda *a, **k: kernel_module)
|
||||
sim.lbm_cfg = SimpleNamespace(validate_runtime=lambda: None, viscosity=.01, omega=1.,
|
||||
streaming="double_buffer", threads_per_block=128,
|
||||
is_d2q9=True, nx=120, ny=80, nz=1)
|
||||
sim._compiled_viscosity = .01
|
||||
sim._assert_allocation_contract = lambda: None
|
||||
sim._resolve_compile_arch = lambda: "sm_70"
|
||||
sim._assert_object_count_contract = lambda **k: None
|
||||
sim.field = SimpleNamespace(module=None, invalidate_module_constant_cache=lambda: None,
|
||||
_upload_params=lambda: None, _runtime_params={},
|
||||
ddf_gpu="ddf", temp_gpu="temp", n_curved=0)
|
||||
sim.stepper = SimpleNamespace(_step_count=7, blocked=None)
|
||||
sim._recompile()
|
||||
assert isinstance(sim.stepper, LBMStepper) and sim.stepper.step_count == 7
|
||||
with pytest.raises(RuntimeError, match="moving kernels"):
|
||||
sim.stepper.step(1, action_gpu=None, obs_gpu=None)
|
||||
assert launches == [] and (sim.field.ddf_gpu, sim.field.temp_gpu) == ("ddf", "temp")
|
||||
assert sim.stepper.step_count == 7
|
||||
|
||||
|
||||
def validate(monkeypatch, objects, flags=None, **cfg):
|
||||
sim = routed_simulation(monkeypatch)
|
||||
for obj in objects:
|
||||
sim.bodies.add(obj)
|
||||
for name, value in cfg.items():
|
||||
setattr(sim.lbm_cfg, name, value)
|
||||
base = np.full(120 * 80, FLUID, np.uint16) if flags is None else flags
|
||||
sim.field = SimpleNamespace(nx=120, ny=80, build_channel_flags=lambda: base.copy())
|
||||
sim._validate_moving_bodies()
|
||||
|
||||
|
||||
def test_initialize_rejects_band_and_modes(monkeypatch):
|
||||
t = Translation((-1., 1., -3., 3.), .1)
|
||||
validate(monkeypatch, [body(50., 40., 10., t), body(100., 40., 5.)]) # accepted
|
||||
walls = np.full(120 * 80, FLUID, np.uint16)
|
||||
walls[60 + 30 * 120] = SOLID
|
||||
cases = [([body(10., 40., 10., t)], None, {}, "domain edge"),
|
||||
([body(50., 40., 10., t)], walls, {}, "non-fluid"),
|
||||
([body(50., 40., 10., t), body(68., 40., 4.)], None, {}, "overlaps"),
|
||||
([body(50., 40., 10., t), body(50., 60., 3., is_sensor=True)], None, {}, "overlaps"),
|
||||
([body(50., 40., 10., t), body(64., 40., 3., is_force_region=True)], None, {}, "overlaps"),
|
||||
([body(50., 40., 10., t)], None, {"streaming": "esopull"}, "double_buffer"),
|
||||
([body(50., 40., 10., t)], None, {"ddf_shifting": True}, "unshifted"),
|
||||
([body(50., 40., 10., t)], None, {"is_d2q9": False}, "D2Q9")]
|
||||
for objects, flags, cfg, match in cases:
|
||||
with pytest.raises(ValueError, match=match):
|
||||
validate(monkeypatch, objects, flags, **cfg)
|
||||
|
||||
|
||||
def test_initialize_validates_before_mutation(monkeypatch):
|
||||
sim = routed_simulation(monkeypatch)
|
||||
sim.bodies.add(body(10., 40., 10., Translation((0., 0., 0., 0.), .1)))
|
||||
uploads = []
|
||||
sim.field = SimpleNamespace(nx=120, ny=80, build_channel_flags=lambda: np.full(9600, FLUID, np.uint16),
|
||||
upload_flags=lambda: uploads.append("flags"))
|
||||
sim._assert_allocation_contract = lambda: None
|
||||
with pytest.raises(ValueError, match="domain edge"):
|
||||
sim.initialize()
|
||||
assert uploads == [] and not sim._initialized
|
||||
|
||||
|
||||
def test_moving_circle_excluded_from_static_links():
|
||||
def manager(*objects):
|
||||
mgr = ObjectManager(nx=120, ny=80, nz=1, nq=9, cfg=Cfg())
|
||||
ids = [mgr.add(obj) for obj in objects]
|
||||
return mgr, ids
|
||||
reference_mgr, _ = manager(body(100., 40., 5.), body(20., 40., 6.))
|
||||
mixed_mgr, ids = manager(body(100., 40., 5.), body(20., 40., 6.),
|
||||
body(55., 40., 10., Translation((0., 0., -1., 1.), .1)))
|
||||
reference = reference_mgr.build_compact_lists_for(reference_mgr.objects)
|
||||
mixed = mixed_mgr.build_compact_lists_for(mixed_mgr.objects)
|
||||
for a, b in zip(reference[:8], mixed[:8]): # static link columns byte-identical
|
||||
assert np.asarray(a).tobytes() == np.asarray(b).tobytes()
|
||||
assert len(reference[0]) > 0 and ids[2] not in set(np.asarray(mixed[3]).tolist())
|
||||
|
||||
|
||||
def test_failed_moving_step_marks_simulation_failed(monkeypatch):
|
||||
sim = routed_simulation(monkeypatch)
|
||||
sim.ctx = ClosableOwner()
|
||||
bid = sim.add_body("circle", (50., 40.), 10., translation=Translation((0., 0., -1., 1.), .1))
|
||||
motion = sim.bodies.get(bid).motion
|
||||
motion.set_target((0., .2), (0., .05))
|
||||
|
||||
def fail(*_):
|
||||
raise RuntimeError("sync failed")
|
||||
with pytest.raises(RuntimeError, match="sync failed"):
|
||||
motion.advance(lambda *_: None, lambda: None, fail)
|
||||
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,
|
||||
sim.save_checkpoint, lambda: sim.load_checkpoint("x.h5")):
|
||||
with pytest.raises(RuntimeError, match="Simulation failed"):
|
||||
call()
|
||||
with pytest.raises(RuntimeError, match="unknown"): # no retry on the failed circle
|
||||
motion.advance(lambda *_: None, lambda: None, lambda: None)
|
||||
with sim.ctx.activate(): # a refused close keeps the failed state and resources
|
||||
with pytest.raises(RuntimeError, match="active CUDA"):
|
||||
sim.close()
|
||||
assert sim._failed is not None and sim.bodies is not None and not sim._closed
|
||||
sim.close() # retry outside the activation succeeds
|
||||
assert sim._closed and sim._failed is not None
|
||||
with pytest.raises(RuntimeError, match="CUDA context is closed"): # no re-entry after close
|
||||
sim.run(1)
|
||||
with pytest.raises(RuntimeError, match="unknown"): # the failed circle stays refused
|
||||
motion.record()
|
||||
@@ -294,8 +294,11 @@ def test_passive_overlap_failure_keeps_live_state(device, entry, second_x):
|
||||
sim = SimpleNamespace(bodies=bodies, field=field, lbm_cfg=field.cfg,
|
||||
_initialized=entry == "runtime_edit",
|
||||
_assert_allocation_contract=lambda: None,
|
||||
stream=SimpleNamespace(synchronize=device.synchronize))
|
||||
sim._require_usable = Simulation._require_usable.__get__(sim)
|
||||
stream=SimpleNamespace(synchronize=device.synchronize),
|
||||
stepper=SimpleNamespace(blocked=None))
|
||||
for name in ("_require_usable", "_moving_bodies", "_reject_moving_bodies",
|
||||
"_validate_moving_bodies", "_attach_moving"):
|
||||
setattr(sim, name, getattr(Simulation, name).__get__(sim))
|
||||
with pytest.raises(ValueError, match="active fluid receivers|unique cut-link writers"):
|
||||
if entry == "sync":
|
||||
bodies.sync_to_gpu(field)
|
||||
|
||||
Reference in New Issue
Block a user