Compare commits
20
Commits
| Author | SHA1 | Date | |
|---|---|---|---|
|
|
4bcf5ef10b | ||
|
|
d5b13d7f31 | ||
|
|
353a7b03c4 | ||
|
|
0d6df6e6b8 | ||
|
|
7886906517 | ||
|
|
48ee0982c9 | ||
|
|
307468b671 | ||
|
|
fb507aa1f2 | ||
|
|
a9ca6a9a16 | ||
|
|
bb83b3476d | ||
|
|
f08b5e671f | ||
|
|
707fa8b5fc | ||
|
|
1abfe954cd | ||
|
|
0c9e041cd1 | ||
|
|
2564329e68 | ||
|
|
78637faffc | ||
|
|
778d3185e5 | ||
|
|
0899799e00 | ||
|
|
5c2b497a43 | ||
|
|
731629f5a3 |
+2
-2
@@ -6,8 +6,8 @@ it is; it adds no framework. Line references are to candidate 0f8fdbc3.
|
||||
| Class | Contents | How it changes | Owner rule |
|
||||
|---|---|---|---|
|
||||
| 1. Module and its globals | Compile-time (visible to the compiler, fixed per module): `config.h` macros incl. `NX`/`NY`, `COLLISION_MODEL`, `BLOCK_ROOT_BC`, `WALL_*`. Runtime per-module globals (not visible to the compiler; written by host code): `__constant__` `d_params`, `d_w`, `d_cx`, `d_cy`, `d_block_channel`, `d_block_wall`; plus bound kernel handles | Globals are written once by the owner right after `load_module`. Uniform grid: `Simulation._recompile` (construction, body-topology sync) loads a new module and re-uploads; `update_runtime_params` rewrites `d_params` (omega, force, n_objects) between runs. Steppers: never rewritten after construction, so they are immutable for the instance's lifetime although they are runtime values. The blocked slot/reservoir root reads its slot from the `d_block_wall` constant (the former dense route used `WALL_*` macros; both were bit-identical in the parity gates). | Every owner loads its own module (`compiler.load_module` = `module_from_file`, no cache). A loaded module is never shared between live owners: globals are per module, so sharing would share mutable `d_params`. |
|
||||
| 2. Geometry, layout, tables | flags, curved-link tables, block layouts and halo pairs, transfer (explode/coalesce) tables, block kinds | Immutable for the life of an instance. Uniform grid: rebuilt by body sync. `FixedPatchStepper`: a geometry change is a new instance (`remap_finest`), committed only after a successful transfer. | No in-place table edits. |
|
||||
| 3. Controls and dependent boundary columns | Jet amplitude -> per-link `surface_velocity`; body actions (uniform grid) | Between steps, without a table or module rebuild (`set_jet_amplitude`). | Input (index, level, amplitude) is validated before any device write; invalid input changes nothing and the instance stays usable. The host copy commits only after the upload returned. Any error after validation leaves the device column undefined and fails the instance: `FixedPatchStepper` status `failed`, `Simulation._failed`; every later device-touching public method (`Simulation`: `stream`, `set_jet_amplitude`, `initialize`, `sync_bodies`, `run`/`step`, `read_surface_exchange`, `read_body(s)`, `get_macroscopic`, `get_ddf`, `update_runtime_params`, `snapshot`, `restore`, save/load checkpoint) raises; cached host readers and host-side staging (e.g. `get_flags`, `add_body`/`remove_body` before a sync) do not touch the device and stay callable. `close()` still releases resources and stays retryable. |
|
||||
| 2. Geometry, layout, tables | flags, curved-link tables, block layouts and halo pairs, transfer (explode/coalesce) tables, block kinds | Immutable for the life of an instance, except translating circles on the uniform grid: their swept-band flags and cut links are rewritten on the device every advancing step (`MovingCircleCells`/`MovingCircleLinks`, docs/moving_bodies.md); the host flag mirror is then stale for those cells. Uniform grid otherwise: rebuilt by body sync. `FixedPatchStepper`: a geometry change is a new instance (`remap_finest`), committed only after a successful transfer. | No in-place table edits. |
|
||||
| 3. Controls and dependent boundary columns | Jet amplitude -> per-link `surface_velocity`; body actions (uniform grid) | Between steps, without a table or module rebuild (`set_jet_amplitude`). | Input (index, level, amplitude) is validated before any device write; invalid input changes nothing and the instance stays usable. The host copy commits only after the upload returned. Any error after validation leaves the device column undefined and fails the instance: `FixedPatchStepper` status `failed`, `Simulation._failed`; every later device-touching public method (`Simulation`: `stream`, `set_jet_amplitude`, `initialize`, `sync_bodies`, `run`/`step`, `read_surface_exchange`, `read_body(s)`, `get_macroscopic`, `get_ddf`, `update_runtime_params`, `snapshot`, `restore`, save/load checkpoint) raises; cached host readers and host-side staging (e.g. `get_flags` without translating circles, `add_body`/`remove_body` before a sync) do not touch the device and stay callable; with translating circles `get_flags` reads the device flags back and refuses after a failure. Translating-circle targets follow the same rule: `set_body` validates before mutation, the applied offset commits only after a successful synchronized step, and a failed moving step fails the instance (docs/moving_bodies.md). `close()` still releases resources and stays retryable. |
|
||||
| 4. Per-step state | population double buffers and their parity, raw-exchange ledgers, observations | Device-only during `step(n)` / `run(n)`; buffer swaps are host-handle swaps. | Valid on the host only after the single wait at the end of the call. |
|
||||
|
||||
Storage routes are separate concerns:
|
||||
|
||||
+59
-15
@@ -1,17 +1,61 @@
|
||||
# Moving circles on the uniform dense grid — S0 DRAFT contract
|
||||
# Moving circles on the uniform dense grid
|
||||
|
||||
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.
|
||||
## Status and user guide (current)
|
||||
|
||||
Items marked **OPEN** are undecided; items marked **CANDIDATE** are proposals.
|
||||
Development-stage feature, integrated from the reviewed slices S1–S3a
|
||||
(production code identical to `653fd997`). Physical evidence and its limits:
|
||||
`docs/validation/moving_bodies_stage1.md`. The force used as a structural RHS
|
||||
is a declared stage choice, not an established result.
|
||||
|
||||
*What it does.* Translating circles (with optional spin) on the dense
|
||||
double-buffer D2Q9 path: per-step cut links from an FP64 geometry record
|
||||
(M2 `MovingCircleLinks`), reclassification and refill of covered/uncovered
|
||||
cells (M1 `MovingCircleCells`), moving-wall Bouzidi with wall velocity
|
||||
`u + omega x r`, and separate per-body diagnostics.
|
||||
|
||||
*API.*
|
||||
- `Simulation.add_body("circle", center=(x, y), radius=r,
|
||||
translation=Translation(offset_box, max_speed))` before `initialize()`;
|
||||
`offset_box = (x_min, x_max, y_min, y_max)` in lattice cells around the
|
||||
base centre (must contain 0), `0 < max_speed ≤ .5`.
|
||||
- `set_body(id, offset=(dx, dy), vx=..., vy=..., omega=...)` sets a pending
|
||||
target (validated against the applied offset: box, Euclidean step ≤ .5
|
||||
cell, speed, no sliding hold); `vx`/`vy` are the interval wall velocity.
|
||||
- `run(1)` advances a pending target and commits it only after the step
|
||||
synchronizes with `no_donor == 0`; `run(n > 1)` with a pending target
|
||||
raises; without a target `run(n)` holds the geometry (spin allowed).
|
||||
Direct `LBMStepper.step` refuses while a target is pending.
|
||||
- `moving_diagnostics()` returns per body the columns of
|
||||
`MovingTables.DIAG` (CME, Wen term, link mass, cover/uncover momentum and
|
||||
mass, eq_only, no_donor, crossings) accumulated since obs were zeroed;
|
||||
`obs`/`read_force` keep their contract (CME on the body).
|
||||
- `get_flags()` reads device flags back when moving circles exist.
|
||||
|
||||
*Supported:* dense double buffer, D2Q9, FP32 unshifted, SRT/MRT, circles,
|
||||
spin. *Rejected before mutation:* EsoPull, ddf shifting, D3Q19, refinement /
|
||||
blocked steppers, Jets or surface actuation on a moving circle,
|
||||
`sync_bodies` and checkpoints while a moving circle exists, swept bands that
|
||||
touch the domain edge, other bodies, sensors or force regions.
|
||||
|
||||
*Failure lifecycle.* Device-side failures (upload, step, sync, no_donor)
|
||||
mark the Simulation failed; later run/step/sync/checkpoint calls and direct
|
||||
stepper calls refuse until `close()`.
|
||||
|
||||
*Known limitations.* No mass correction at curved walls (static or moving);
|
||||
cover/uncover and Wen terms are diagnostics, not a selected force law; the
|
||||
host loop costs a fixed per-step synchronization (≈ 1 ms/step measured on a
|
||||
7.2 M-cell grid, run(1) class); outer boundaries are the single dense
|
||||
channel package (see `configs/CONFIG.md`).
|
||||
|
||||
## Design record (staged contract, kept as written)
|
||||
|
||||
The sections below are the contract the slices were reviewed against
|
||||
(S0 → S3a). They are kept unchanged for provenance; where they say "draft",
|
||||
"proposed" or "until S2/S3", the status above supersedes them.
|
||||
|
||||
Original scope (Frank via TRON, thread C0C78UFKGRM/1793335557.460509):
|
||||
uniform-grid prescribed translation and single-cylinder free transverse VIV;
|
||||
base main `b2034f87e03656c5f67ab88562a70d6c338b9918`.
|
||||
|
||||
## 1. Working standards
|
||||
|
||||
@@ -126,9 +170,9 @@ acts as a sliding wall. One moving step consumes one target.
|
||||
| `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 |
|
||||
| `run(0)` | no motion upload, no commit; `pending` kept; action upload and obs/diagnostic zeroing still run |
|
||||
| validation failure (incl. action validation and pinned staging in `run`) | `ValueError`/error before any device work, nothing mutated, retry allowed |
|
||||
| failure of any device work in `run` (action upload, obs/diagnostic zeroing, record upload, launch, readback, sync), also for `run(0)` | simulation unusable, every translating circle failed; `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).
|
||||
|
||||
@@ -82,6 +82,13 @@ Tests: `tests/unit/test_jets.py`, `test_jets_review.py`, `test_shifted_jets.py`,
|
||||
- DFG 2D-1 Re 20 (V4, SRT, t* 60, pressure-NEQ outlet and local Zou-He inlet): refined and uniform
|
||||
runs were NOT_STEADY under the frozen steady rule, so no interval verdict was given.
|
||||
|
||||
## Moving circles (dense `Simulation`, uniform grid only)
|
||||
|
||||
Translating circles with optional spin are supported on the dense double
|
||||
buffer D2Q9 FP32 path only (not on `FixedPatchStepper` or refined levels).
|
||||
API, supported/rejected modes and limits: `docs/moving_bodies.md`; stage-1
|
||||
evidence: `docs/validation/moving_bodies_stage1.md`.
|
||||
|
||||
## Not supported in this release
|
||||
|
||||
3D, dynamic/automatic AMR (tagging, time-varying hierarchies beyond `remap_finest`),
|
||||
|
||||
@@ -0,0 +1,46 @@
|
||||
# Moving circles — stage-1 evidence (development level)
|
||||
|
||||
Owner criterion (Yanqi via TRON): prescribed motion and one free VIV case both
|
||||
"reasonably close" to a predeclared reference = development-stage success.
|
||||
This is practical closeness, not a formal validation claim. Every number
|
||||
below is from Astra's receipts as posted in Slack C0C78UFKGRM thread
|
||||
1793335557.460509; raw outputs stay on Horsea under `/tmp/celeris-results-*`
|
||||
and the Mac evidence directory (not in git).
|
||||
|
||||
## Code
|
||||
|
||||
- Production: `src/` identical to `653fd997` (slices S1 `983020b6`,
|
||||
S2a `b920a303`, S3a `653fd997`). CPU suites and GPU hold/crossing gates
|
||||
passed on those slices (G0a/G0b population-bitwise; X1 crossing ledger).
|
||||
- Research drivers: `research/duetsch_re100_kc5/` (`dd37883d`, capture
|
||||
`a0aea3b9`) and `research/viv_re150/` (`fb507aa1`), run as committed.
|
||||
|
||||
## Prescribed in-line oscillation (Dütsch et al. 1998, Re 100, KC 5)
|
||||
|
||||
- Protocol: `docs/validation/duetsch_re100_kc5_protocol.md`.
|
||||
- Frozen settling gate (.002 F_ref waveform): **FAILED** and kept on record:
|
||||
D40 domain stage (20/40/60 D, 11.25 T) and the 60 D duration diagnostic
|
||||
(31.25 T, 0/28 pairs; last waveform drift .00613 F_ref).
|
||||
- Unscored late-cycle CME fits c_d ≈ 2.02–2.05, c_i ≈ 1.46–1.48 vs the
|
||||
paper's computed (2.09, 1.45); fields clean and symmetric (owner view).
|
||||
The owner accepts ~0.61 % cycle drift at D40 for this stage; the gate is
|
||||
not re-scored.
|
||||
- Fig. 8 LDA profiles: qualitative only (no digitized ordinates).
|
||||
|
||||
## Free transverse VIV (Kolahdouz §4.1.1 setup, Qin §4.5 target)
|
||||
|
||||
- Predeclaration: `docs/validation/viv_re150_predeclare.md` (agreed at
|
||||
`1abfe954` under its former name `leontini_viv_predeclare.md`); run at `fb507aa1`, single case, D = 40, 75 D × 60 D.
|
||||
- Result: status `scored`, 336,384 steps (t U/D = 336.4); A_max = .5464 D vs
|
||||
.520 D (+5.1 %, inside the declared ±10 % band); mean half peak-to-peak
|
||||
.5448 D; f/f_n = .955; block drifts .78/.29/.18 %, last-five spreads
|
||||
.15/.11/.10 %; finite, no_donor = 0, no seed; 337.7 s, ~1.0 ms/step.
|
||||
|
||||
## Limitations (binding)
|
||||
|
||||
- One resolution (D40) per case; no grid or time convergence.
|
||||
- Structural RHS = CME + Wen is a theory-motivated stage choice; CME-only
|
||||
and node-impulse variants were not run as RHS; no force law is validated.
|
||||
- Boundary mismatches: NEQ outlet with backflow clamp vs traction outlet;
|
||||
finite box; the strict Dütsch settling failure is unexplained.
|
||||
- No mass correction at curved walls; ~1 ms/step host-coupled overhead.
|
||||
@@ -0,0 +1,148 @@
|
||||
# Single-cylinder free transverse VIV, Re 150 (Kolahdouz/Qin) — PREDECLARATION (agreed at 1abfe954)
|
||||
|
||||
(File renamed from `leontini_viv_predeclare.md`; the Leontini Re 200 case is only the queued second case below.)
|
||||
|
||||
Stage goal (Yanqi via TRON): prescribed Dütsch case + one VIV case both
|
||||
"reasonably close" to a predeclared reference = development-stage success, then
|
||||
cleanup. The frozen Dütsch settling failure (0/28 pairs at 31.25 T) stays on
|
||||
record. Practical closeness is not a formal validation claim. Agreed by
|
||||
Claude, Grok and Astra at 1abfe954; implemented in research/viv_re150/.
|
||||
|
||||
## Primary case (one case)
|
||||
|
||||
- Reference: Qin et al., arXiv:2003.12186, §4.5 (maximum amplitudes at
|
||||
U_red = 3…8: .075, .551, .520, .459, .370, .084 D in the text, as quoted by
|
||||
Grok/Astra) for the target; setup from Kolahdouz et al., arXiv:2003.12046,
|
||||
§4.1.1 (read 2026-10-30 via ar5iv): U* sweep 3–8 with amplitudes only in
|
||||
Fig. 3(b) (no text values), domain [−30, 45] D × [−30, 30] D with the
|
||||
cylinder "initially at rest and centered at the origin", uniform inflow on
|
||||
the left, zero normal traction and tangential velocity at the outlet, zero
|
||||
normal velocity and tangential traction on y = ±30 D; no initial
|
||||
perturbation and no total simulated time stated.
|
||||
- Re = 150, U* = U/(f_n D) = 5, m/(rho D²) = 2 (m* = rho_s/rho_f = 8/pi),
|
||||
zeta = 0, f_n = sqrt(K/M)/(2 pi) in vacuum. Transverse DOF only.
|
||||
- Target: A_max = .520 D (peak amplitude, Qin §4.5). Multi-code reference
|
||||
(Kolahdouz ILE-NS, Qin II-LBM). Reference uncertainty: not stated; the
|
||||
spread between codes is not quantified here.
|
||||
|
||||
## Our reproduction (declared differences)
|
||||
|
||||
- D = 40, MRT, dense double buffer; U_lat = .04, nu = U D / 150 = .010667,
|
||||
tau = .532.
|
||||
- Box [0, 75 D] × [0, 60 D] (3000 × 2400) with the cylinder centre at
|
||||
(30 D, 30 D): 30 D from the inlet, 45 D from the outlet, centred in y
|
||||
(Kolahdouz §4.1.1). Qin's own box ([0, 40 D] × [0, 20 D], centre (10 D, 10 D))
|
||||
is the target paper's, not ours.
|
||||
- Boundaries: uniform Zou–He inlet, NEQ outlet (rho = 1, backflow clamp),
|
||||
free-slip y walls. Kolahdouz uses a traction outlet: declared mismatch.
|
||||
- Qin's own VIV box (40 D × 20 D) is not used; a smaller box would be a
|
||||
declared cost cut with more blockage, not "the same case".
|
||||
- Allowed y offset ±1.5 D; peak transverse speed ≈ .026 cells/step (bound .5).
|
||||
|
||||
## Structure, coupling and time levels
|
||||
|
||||
- M = 2 rho D² (lattice rho = 1), f_n = U/(5 D) = 2e-4 per step, K = M (2 pi f_n)²,
|
||||
zeta = 0. Equation of motion uses the structural mass only.
|
||||
- Per-step sequence (binding): known (y_n, v_n) and the last completed force
|
||||
F_last → candidate (y_{n+1}, v_{n+1}) by the update below → `set_body` with
|
||||
offset y_{n+1} and wall velocity (y_{n+1} − y_n)/dt with lattice dt = 1 (interval average; a
|
||||
different quantity from the endpoint v_{n+1}) → `run(1)` → only on success:
|
||||
commit y_n ← y_{n+1}, v_n ← v_{n+1}, F_last ← force read for that step (failure
|
||||
semantics: see the window section). Whether the
|
||||
measured force represents the interval's start, midpoint or end is not
|
||||
asserted; the one-step lag is declared, not corrected.
|
||||
- Initial state: fluid initialized by `Simulation.initialize()` to the uniform
|
||||
inlet velocity U everywhere (the code's actual initial condition; the
|
||||
source gives no fluid initial condition), body at rest at the declared
|
||||
centre, y_0 = v_0 = 0, F_last = 0.
|
||||
- Structural update over one step with constant force F = F_last (exact
|
||||
harmonic solution, zeta = 0): y_e = F/K, w = 2 pi f_n,
|
||||
y^{n+1} = y_e + (y^n − y_e) cos w + (v^n / w) sin w,
|
||||
v^{n+1} = −(y^n − y_e) w sin w + v^n cos w.
|
||||
- Primary RHS (single coupled run): CME + Wen (theory-motivated stage choice:
|
||||
Galilean-invariant link correction; G4 could not separate it from CME).
|
||||
Raw CME and the node terms are logged as diagnostics of that trajectory;
|
||||
the node impulse is not added. A run with CME-only as RHS is a separate,
|
||||
labelled sensitivity case, decided separately. Neither RHS is claimed proven.
|
||||
|
||||
## Window, failure semantics, acceptance (engineering, development stage)
|
||||
|
||||
Clocks: lattice dt = 1; T_n = 1/f_n = 5,000 steps. Response cycles are
|
||||
counted from the measured signal; T_n is used only for the seed check.
|
||||
|
||||
- Start: body at rest at the declared centre (the only initial condition the
|
||||
source states), y_0 = v_0 = 0, F_last = 0; fluid uniform at U (as
|
||||
initialized by the code; our choice, not the source's).
|
||||
- Seed (conditional, one time, declared): at step n_seed = 200,000
|
||||
(t = 40 T_n), if max |y_t| over states t in [175,000, 200,000) is below
|
||||
.01 D, then
|
||||
in that step, after the check and before the candidate update, set
|
||||
v_n += 2 pi f_n (.01 D) (F_last unchanged) and proceed normally; the
|
||||
geometry moves continuously through the update. `seed_applied`, step and
|
||||
size are logged; never repeated or changed.
|
||||
- Cycle boundary (discrete, fixed threshold): step n closes a cycle iff
|
||||
y[n−1] < 0 and y[n] ≥ 0 (upward crossing of the equilibrium y = 0);
|
||||
y_0 = 0 is not a crossing; no boundary fewer than 2 samples after the
|
||||
previous one.
|
||||
- Per cycle: positive peak max(y), negative peak min(y), mean offset
|
||||
(time-mean of y), amplitude a = (max − min)/2 (the only amplitude used for
|
||||
settling).
|
||||
- Eligibility: a completed cycle is counted (toward the settling proof and
|
||||
the 150 budget) only if a ≥ .01 D; ineligible cycles are logged, never
|
||||
counted, so near-zero motion cannot settle or pass.
|
||||
- Seed reset: if the seed is applied, the incomplete cycle at n_seed is
|
||||
discarded and the settling proof starts empty; the 150 budget and the
|
||||
native cap keep running from t = 0.
|
||||
- Windows: adjacent non-overlapping blocks of 10 counted cycles;
|
||||
W_k = mean amplitude of block k (W_k ≥ .01 D by eligibility, so no zero
|
||||
denominator).
|
||||
- Block rule R(k), k ≥ 2: |W_k − W_{k−1}| / W_{k−1} ≤ .02 and, for the last 5
|
||||
cycles of block k, max_i |a_i − ā_5| / ā_5 ≤ .01 with ā_5 their mean.
|
||||
- Settled at the end of block k when R(k) holds. Scoring = the next two
|
||||
blocks k+1, k+2 (20 counted cycles); both must satisfy R. If one fails,
|
||||
that block becomes the new reference: the settling proof restarts from it
|
||||
and continues within the caps (deterministic; no other retry). A pass
|
||||
requires one complete stable 20-cycle scoring window before a cap; success
|
||||
stops after its 20th cycle.
|
||||
- Caps (both always armed; the first hit wins and is reported): 150 counted
|
||||
response cycles including the 20 scoring cycles (so settling must complete
|
||||
by counted cycle 130) → `timeout_response_cycles`; 1,000,000 native steps
|
||||
(200 T_n) → `timeout_native_steps`.
|
||||
- Durations reported: total steps, t U/D, counted and uncounted cycles,
|
||||
measured T_resp (mean boundary interval); N_resp / f_n only as a labelled
|
||||
secondary clock.
|
||||
- Growth is not a failure: initial VIV growth is physical. Checked every
|
||||
step (FAIL-CLOSED): non-finite y, v or force; |y| > 1.0 D (the tighter
|
||||
gate; about twice the largest reference peak, .551 D) or the ±1.5 D offset
|
||||
box reached; no_donor > 0; a failed upload/step/sync.
|
||||
- On any failure the driver stops and reports; it never retries,
|
||||
sub-iterates, retunes the RHS or parameters, or starts a second case.
|
||||
Device-side failures (upload/step/sync, no_donor) mark the Simulation
|
||||
failed through the existing lifecycle; input rejection by `set_body`
|
||||
(offset box, speed, step bound) and timeouts end the driver run without
|
||||
that mark.
|
||||
- Target metric over the 20 scoring cycles: A_max/D = max over those cycles
|
||||
of max(|positive peak|, |negative peak|) about y = 0 (Qin reports "maximum
|
||||
amplitude"); also reported: mean amplitude, positive and negative peaks,
|
||||
mean offset, f/f_n (no target).
|
||||
- Acceptance: |A_max − .520 D| ≤ .052 D (±10 %), stated together with the
|
||||
boundary mismatch and the absence of a quoted reference uncertainty. A pass
|
||||
is "reasonably close" for the stage: not grid convergence, not force-RHS
|
||||
validation, and it does not touch the Dütsch settling record.
|
||||
- Driver reads force and the small diagnostics per step only; field reads
|
||||
only at a few predeclared snapshots. No runtime claim: there is no measured
|
||||
timing on this grid yet (the D40/60D receipt is a different, sampled case).
|
||||
|
||||
## Plumbing smoke (not the benchmark)
|
||||
|
||||
`run.py --smoke --D 20 --steps 32` runs exactly 32 steps of the same BC
|
||||
package and coupling loop (no seed, settle or score; never scored);
|
||||
status `smoke_ok` iff every step is finite, no_donor stays 0 and the body
|
||||
stays inside its box. The physical case is only the default invocation
|
||||
(`--D 40`, no `--smoke`).
|
||||
|
||||
## Queued, not in this stage
|
||||
|
||||
Leontini et al. JFS 22 (2006) Re = 200, m* = 10, zeta = 0.01, U* = 5.0 (needs
|
||||
the owners' c_crit convention and a Fig. 1 digitization). D80, literature
|
||||
claims, pinball, AMR, BC redesign, threshold changes.
|
||||
@@ -0,0 +1,134 @@
|
||||
"""One-shot predeclared VIV run (docs/validation/viv_re150_predeclare.md).
|
||||
|
||||
Re 150, U* = 5, m/(rho D^2) = 2, zeta = 0, D = 40 on a 75 D x 60 D native
|
||||
channel (uniform inlet, NEQ outlet, free-slip y walls), cylinder at
|
||||
(30 D, 30 D). Structural RHS = CME + Wen (diagnostics: raw CME, node terms).
|
||||
Research driver only: Simulation API, no production change.
|
||||
|
||||
PYTHONPATH=src python research/viv_re150/run.py --out /path/outside/checkout
|
||||
"""
|
||||
from __future__ import annotations
|
||||
|
||||
import argparse
|
||||
import json
|
||||
import math
|
||||
import os
|
||||
import tempfile
|
||||
|
||||
import numpy as np
|
||||
|
||||
import viv
|
||||
|
||||
D, U = 40, 0.04
|
||||
TEMPLATE = os.path.join(os.path.dirname(os.path.abspath(__file__)), "..", "..", "src",
|
||||
"CelerisLab", "configs", "config_lbm.json")
|
||||
|
||||
|
||||
def config(D=D, U=U):
|
||||
nx, ny = 75 * D, 60 * D
|
||||
with open(TEMPLATE) as fh:
|
||||
cfg = json.load(fh)
|
||||
cfg["grid"].update(lattice_model="D2Q9", nx=nx, ny=ny, nz=1)
|
||||
cfg["physics"].update(data_type="FP32", viscosity=U * D / viv.RE, velocity=U, rho=1.0)
|
||||
cfg["method"].update(collision="MRT", streaming="double_buffer", store_precision="FP32",
|
||||
ddf_shifting=False, y_wall_bc="free_slip")
|
||||
cfg["method"]["les"]["enabled"] = False
|
||||
cfg["method"]["inlet"]["profile"] = "uniform"
|
||||
return cfg, nx, ny
|
||||
|
||||
|
||||
def main():
|
||||
ap = argparse.ArgumentParser()
|
||||
ap.add_argument("--out", required=True)
|
||||
ap.add_argument("--D", type=int, default=D)
|
||||
ap.add_argument("--smoke", action="store_true",
|
||||
help="plumbing smoke only: exactly --steps steps (<= 64), never scored")
|
||||
ap.add_argument("--steps", type=int, default=0)
|
||||
args = ap.parse_args()
|
||||
if args.smoke != bool(args.steps) or args.steps > 64 or args.steps < 0:
|
||||
raise SystemExit("--smoke requires --steps N with 1 <= N <= 64 (and only with --smoke)")
|
||||
if args.D != D and not args.smoke:
|
||||
raise SystemExit("the physical case is predeclared at D = 40")
|
||||
if os.path.realpath(args.out).startswith(os.path.realpath(os.path.join(os.path.dirname(__file__), "..", ".."))):
|
||||
raise SystemExit("--out must be outside the checkout")
|
||||
os.makedirs(args.out, exist_ok=True)
|
||||
for name in ("report.json", "steps.npy"): # never overwrite earlier evidence
|
||||
if os.path.exists(os.path.join(args.out, name)):
|
||||
raise SystemExit(f"{args.out} already holds {name}; use a fresh output directory")
|
||||
Dn = args.D
|
||||
cfg, nx, ny = config(Dn)
|
||||
s, cyc = viv.Structure(Dn, U), viv.Cycles(Dn)
|
||||
T_n = 1.0 / s.f_n
|
||||
step_cap, n_seed = int(round(viv.STEP_CAP_TN * T_n)), int(round(viv.SEED_TN * T_n))
|
||||
quiet = int(round(5 * T_n))
|
||||
if args.smoke:
|
||||
step_cap = args.steps
|
||||
centre = (30.0 * Dn, 0.5 * (ny - 1))
|
||||
|
||||
from CelerisLab.body.motion import Translation
|
||||
from CelerisLab.simulation import Simulation
|
||||
with tempfile.NamedTemporaryFile("w", suffix=".json", delete=False) as fh:
|
||||
json.dump(cfg, fh)
|
||||
sim = Simulation(fh.name)
|
||||
rows, status, seed_step, f_last = [], "running", None, 0.0
|
||||
try:
|
||||
body = sim.add_body("circle", centre, 0.5 * Dn,
|
||||
translation=Translation((0.0, 0.0, -1.5 * Dn, 1.5 * Dn), 0.5))
|
||||
sim.initialize()
|
||||
for n in range(step_cap):
|
||||
# rows[k] holds the state y_{k+1}; the declared window is y_t, t in [n_seed - 5 T_n, n_seed)
|
||||
if not args.smoke and n == n_seed and \
|
||||
max(abs(r[1]) for r in rows[n_seed - quiet - 1:n_seed - 1]) < viv.ELIGIBLE * Dn:
|
||||
s.v += s.seed_velocity(Dn) # after the quiet check, before the update
|
||||
seed_step = n
|
||||
cyc.discard_current()
|
||||
cyc.reset_proof()
|
||||
y1, v1 = s.candidate(f_last)
|
||||
if not (math.isfinite(y1) and math.isfinite(v1)) or abs(y1) > viv.Y_LIMIT * Dn:
|
||||
status = "fail_bound_or_nonfinite"
|
||||
break
|
||||
try:
|
||||
sim.set_body(body, offset=(0.0, y1), vy=y1 - s.y)
|
||||
except ValueError as exc: # input rejection: stop, no retry
|
||||
status = f"rejected: {exc}"
|
||||
break
|
||||
sim.run(1) # device failures raise and mark the Simulation failed
|
||||
d = sim.moving_diagnostics()[0]
|
||||
f_obs = sim.read_force(body)[1]
|
||||
f_rhs = float(d[1] + d[3]) # CME + Wen (y)
|
||||
row = (n, y1, v1, f_rhs, float(d[1]), float(d[3]), float(f_obs),
|
||||
float(d[7] - d[10]), float(d[13]))
|
||||
if not viv.finite(row): # every recorded output, not only the RHS
|
||||
status = "fail_nonfinite_output"
|
||||
break
|
||||
s.y, s.v, f_last = y1, v1, f_rhs # commit only after a successful step
|
||||
rows.append(row)
|
||||
if args.smoke:
|
||||
continue # no settle / score in the plumbing smoke
|
||||
if cyc.update(s.y) == "scored":
|
||||
status = "scored"
|
||||
break
|
||||
if cyc.counted >= viv.CYCLE_CAP:
|
||||
status = "timeout_response_cycles"
|
||||
break
|
||||
else:
|
||||
status = ("smoke_ok" if all(r[8] == 0.0 for r in rows) else "smoke_no_donor") \
|
||||
if args.smoke else "timeout_native_steps"
|
||||
except BaseException as exc:
|
||||
status = f"error: {exc!r}"
|
||||
raise
|
||||
finally:
|
||||
np.save(os.path.join(args.out, "steps.npy"), np.asarray(rows, float))
|
||||
report = dict(status=status, smoke=args.smoke, D=Dn, steps=len(rows), t_U_over_D=len(rows) * U / Dn, T_n=T_n,
|
||||
seed_step=seed_step, counted_cycles=cyc.counted, cycles=cyc.log,
|
||||
columns=["n", "y", "v", "F_rhs_cme_wen", "F_cme", "F_wen", "F_obs",
|
||||
"node_cover_minus_uncover_y", "no_donor"])
|
||||
if status == "scored" and not args.smoke:
|
||||
report["score"] = viv.score(cyc.window)
|
||||
with open(os.path.join(args.out, "report.json"), "w") as fh:
|
||||
json.dump(report, fh, indent=1)
|
||||
sim.close()
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
@@ -0,0 +1,126 @@
|
||||
"""Pure logic of the predeclared Re 150 free transverse VIV case (stdlib only).
|
||||
|
||||
Contract: docs/validation/viv_re150_predeclare.md. Oscillator update with a
|
||||
constant force over one lattice step, response-cycle tracking with fixed y = 0
|
||||
crossings, eligibility, 10-cycle blocks, settling/scoring with deterministic
|
||||
restart, and the conditional seed.
|
||||
"""
|
||||
from __future__ import annotations
|
||||
|
||||
import math
|
||||
|
||||
RE, U_STAR, MASS_NUMBER = 150.0, 5.0, 2.0 # m/(rho D^2) = 2 -> m* = 8/pi
|
||||
TARGET_A, BAND = 0.520, 0.10 # Qin §4.5 maximum amplitude, ±10 %
|
||||
ELIGIBLE, BLOCK, LAST5 = 0.01, 10, 5 # in D; cycles; cycles
|
||||
WINDOW_TOL, LAST5_TOL = 0.02, 0.01
|
||||
CYCLE_CAP, STEP_CAP_TN, SEED_TN, Y_LIMIT = 150, 200, 40, 1.0
|
||||
|
||||
|
||||
def finite(values):
|
||||
"""True iff every recorded value is finite (fail-closed output check)."""
|
||||
return all(math.isfinite(float(v)) for v in values)
|
||||
|
||||
|
||||
class Structure:
|
||||
"""Undamped spring-mass, exact harmonic update over dt = 1 with constant force."""
|
||||
|
||||
def __init__(self, D, u):
|
||||
self.m = MASS_NUMBER * D * D # lattice rho = 1
|
||||
self.f_n = u / (U_STAR * D)
|
||||
self.w = 2.0 * math.pi * self.f_n
|
||||
self.k = self.m * self.w * self.w
|
||||
self.y = self.v = 0.0
|
||||
|
||||
def candidate(self, force):
|
||||
ye = force / self.k
|
||||
c, s = math.cos(self.w), math.sin(self.w)
|
||||
return (ye + (self.y - ye) * c + self.v / self.w * s,
|
||||
-(self.y - ye) * self.w * s + self.v * c)
|
||||
|
||||
def seed_velocity(self, D):
|
||||
return self.w * ELIGIBLE * D # amplitude ~ .01 D
|
||||
|
||||
|
||||
class Cycles:
|
||||
"""Upward y = 0 crossings, eligibility, blocks, settle/score with restart."""
|
||||
|
||||
def __init__(self, D):
|
||||
self.D = D
|
||||
self.reset_proof()
|
||||
self.counted = 0
|
||||
self.log = [] # one dict per completed cycle
|
||||
self._cur = None
|
||||
self._prev = 0.0
|
||||
self._since = 0
|
||||
|
||||
def reset_proof(self):
|
||||
self.amps = [] # eligible amplitudes since proof start
|
||||
self.blocks = [] # block means
|
||||
self.settled_at = None # block index where R held
|
||||
self.window = None # completed scoring cycles (list of dicts)
|
||||
|
||||
def discard_current(self):
|
||||
self._cur, self._since = None, 0
|
||||
|
||||
def rule(self, k):
|
||||
"""R(k): block-mean change ≤ 2 % and last-five spread ≤ 1 %."""
|
||||
if k < 1:
|
||||
return False
|
||||
w0, w1 = self.blocks[k - 1], self.blocks[k]
|
||||
last = self.amps[(k + 1) * BLOCK - LAST5:(k + 1) * BLOCK]
|
||||
a5 = sum(last) / LAST5
|
||||
return abs(w1 - w0) / w0 <= WINDOW_TOL and max(abs(a - a5) for a in last) / a5 <= LAST5_TOL
|
||||
|
||||
def update(self, y):
|
||||
"""Feed y after a committed step; returns 'scored' when a stable window completes."""
|
||||
out = None
|
||||
crossing = self._prev < 0.0 <= y and self._since >= 2
|
||||
if crossing:
|
||||
if self._cur is not None:
|
||||
out = self._close()
|
||||
self._cur, self._since = dict(ymax=y, ymin=y, ysum=0.0, n=0), 0
|
||||
if self._cur is not None:
|
||||
c = self._cur
|
||||
c["ymax"], c["ymin"] = max(c["ymax"], y), min(c["ymin"], y)
|
||||
c["ysum"] += y
|
||||
c["n"] += 1
|
||||
self._since += 1
|
||||
self._prev = y
|
||||
return out
|
||||
|
||||
def _close(self):
|
||||
c = self._cur
|
||||
a = 0.5 * (c["ymax"] - c["ymin"]) / self.D
|
||||
cyc = dict(ymax=c["ymax"] / self.D, ymin=c["ymin"] / self.D,
|
||||
offset=c["ysum"] / c["n"] / self.D, amp=a, steps=c["n"], eligible=a >= ELIGIBLE)
|
||||
self.log.append(cyc)
|
||||
if not cyc["eligible"]:
|
||||
return None
|
||||
self.counted += 1
|
||||
self.amps.append(a)
|
||||
if self.settled_at is not None:
|
||||
self.window.append(cyc)
|
||||
if len(self.amps) % BLOCK:
|
||||
return None
|
||||
self.blocks.append(sum(self.amps[-BLOCK:]) / BLOCK)
|
||||
k = len(self.blocks) - 1
|
||||
if self.settled_at is None:
|
||||
if self.rule(k):
|
||||
self.settled_at, self.window = k, []
|
||||
return None
|
||||
if not self.rule(k): # scoring block failed: restart from it
|
||||
amps, blocks = self.amps[-BLOCK:], self.blocks[-1:]
|
||||
self.reset_proof()
|
||||
self.amps, self.blocks = amps, blocks
|
||||
return None
|
||||
return "scored" if k == self.settled_at + 2 else None
|
||||
|
||||
|
||||
def score(window):
|
||||
"""Target metric and reported quantities over the 20 scoring cycles."""
|
||||
a_max = max(max(abs(c["ymax"]), abs(c["ymin"])) for c in window)
|
||||
return dict(A_max=a_max, mean_amp=sum(c["amp"] for c in window) / len(window),
|
||||
ymax=max(c["ymax"] for c in window), ymin=min(c["ymin"] for c in window),
|
||||
offset=sum(c["offset"] for c in window) / len(window),
|
||||
T_resp=sum(c["steps"] for c in window) / len(window),
|
||||
within_band=abs(a_max - TARGET_A) <= BAND * TARGET_A)
|
||||
@@ -629,10 +629,17 @@ class ObjectManager:
|
||||
cuda.memcpy_dtoh_async(self.obs_pinned, self.obs_gpu, stream)
|
||||
|
||||
def _upload_action_async(self, stream: cuda.Stream) -> None:
|
||||
"""Async H2D of host action array to action_gpu.
|
||||
"""Checked async H2D of host action array to action_gpu."""
|
||||
if self._prepare_action_upload():
|
||||
self._upload_prepared_action_async(stream)
|
||||
|
||||
def _prepare_action_upload(self) -> bool:
|
||||
"""Validate the action array and stage it in pinned memory; no device work.
|
||||
|
||||
Single gate for every action path (setters or direct ``action`` edits):
|
||||
surface-actuated and Jet-driven bodies must keep zero motion slots.
|
||||
Returns False when there is nothing to upload. Restaging assumes no
|
||||
copy from ``_action_pinned`` is in flight (``run`` ends with a sync).
|
||||
"""
|
||||
slot = 3*(self.cfg.dim if self.cfg else 2)
|
||||
for obj in self.objects:
|
||||
@@ -643,10 +650,14 @@ class ObjectManager:
|
||||
self.validate_jet_motion(obj.obj_id, {"vx": self.action[base], "vy": self.action[base+1],
|
||||
"omega": self.action[base+slot-1]})
|
||||
if self.action_gpu is None or self.action.size == 0:
|
||||
return
|
||||
return False
|
||||
if self._action_pinned is None or self._action_pinned.size != self.action.size:
|
||||
self._action_pinned = cuda.pagelocked_empty(self.action.size, np.float32)
|
||||
self._action_pinned[:] = self.action
|
||||
return True
|
||||
|
||||
def _upload_prepared_action_async(self, stream: cuda.Stream) -> None:
|
||||
"""Enqueue H2D of the staged ``_action_pinned`` (see ``_prepare_action_upload``)."""
|
||||
cuda.memcpy_htod_async(self.action_gpu, self._action_pinned, stream)
|
||||
|
||||
@property
|
||||
|
||||
@@ -99,7 +99,12 @@ class MovingCircle:
|
||||
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.
|
||||
"""Single-circle transaction helper (unit tests): upload, step, sync, commit.
|
||||
|
||||
Not the production path: ``Simulation._run_moving`` batches all circles,
|
||||
adds the ``no_donor`` gate and uses ``commit``/``fail`` directly.
|
||||
|
||||
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,
|
||||
|
||||
@@ -86,6 +86,20 @@ step/one_step_*.cu → kernel 编排
|
||||
- 转速 `omega` 属于运行时状态,可在步进期间更新;该更新不改变 `N_OBJS` 或 `OBS_*` 内存布局,因此**不需要**因 `omega` 变化触发重新编译。
|
||||
- 力矩观测从 `obs` 的 torque segment 读取;force / torque / sensor 段统一由 `config_obs.h` 宏描述。
|
||||
|
||||
### 外边界能力(Outer boundaries, dense `Simulation`)
|
||||
|
||||
The dense path supports one channel package, not four independently
|
||||
selectable sides:
|
||||
- left x = 0: velocity inlet (`method.inlet.scheme`, `profile`, `physics.velocity`);
|
||||
- right x = nx−1: outlet (`method.outlet.mode`; `neq_extrap` holds rho = 1 with
|
||||
a `max(u, 0)` backflow clamp when `backflow_clamp` is set);
|
||||
- bottom and top: one shared `y_wall_bc` (`bounce_back` halfway or
|
||||
`free_slip`); corners are walls.
|
||||
`y_wall_bc = "on_site_velocity"` switches to x-periodic with on-site y planes
|
||||
and requires an object-free domain. There is no x/y-periodic or symmetry
|
||||
option for domains with bodies; adding outer BC types is a separate design
|
||||
question, not covered by these keys.
|
||||
|
||||
### 稳定性默认窗口(Stability defaults)
|
||||
|
||||
- 默认 `method.omega_guard = {min: 0.01, max: 1.99}`,用于约束碰撞频率与 LES 有效粘度路径。
|
||||
|
||||
@@ -21,7 +21,7 @@ class LBMStepper:
|
||||
|
||||
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).
|
||||
# without attached moving tables, a pending target, or a failed step).
|
||||
self.blocked: str | None = None
|
||||
self.moving = None # MovingTables attached by the owner (translating circles)
|
||||
|
||||
@@ -39,8 +39,10 @@ class LBMStepper:
|
||||
self.apply_esopull_curved_fn = module.get_function("ApplyEsoPullCurvedKernel")
|
||||
self.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")
|
||||
# Moving-circle kernels exist only in D2Q9 double-buffer modules.
|
||||
moving = not self._esopull and cfg.is_d2q9
|
||||
self.moving_fn = module.get_function("MovingCircleLinks") if moving else None
|
||||
self.moving_cells_fn = module.get_function("MovingCircleCells") if moving else None
|
||||
|
||||
tpb = cfg.threads_per_block
|
||||
self.block = (tpb, 1, 1)
|
||||
|
||||
@@ -28,7 +28,7 @@ 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).
|
||||
# Translating circles step only through their moving tables (attached at initialize).
|
||||
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); "
|
||||
@@ -122,8 +122,9 @@ class Simulation:
|
||||
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.
|
||||
circle (before initialize only; docs/moving_bodies.md).
|
||||
Positive-step runs are rejected until ``initialize()`` has
|
||||
attached its moving tables.
|
||||
**kwargs: passed to the geometry constructor.
|
||||
|
||||
Returns:
|
||||
@@ -215,6 +216,12 @@ class Simulation:
|
||||
self._failed = f"moving circle step failed: {exc!r}"
|
||||
self.stepper.blocked = self._failed # direct stepper calls are refused as well
|
||||
|
||||
def _fail_moving_run(self, exc: BaseException) -> None:
|
||||
# Device work was attempted with translating circles (or their tables) present.
|
||||
for obj in self._moving_bodies():
|
||||
obj.motion.fail(exc)
|
||||
self._fail_moving(exc) # also when tables are attached but no circle is registered
|
||||
|
||||
def _attach_moving(self) -> None:
|
||||
# Direct stepper calls stay blocked until the moving tables exist.
|
||||
tables = getattr(self, "_moving_tables", None)
|
||||
@@ -241,7 +248,8 @@ class Simulation:
|
||||
self._require_usable()
|
||||
tables = self.stepper.moving
|
||||
if tables is None:
|
||||
return np.zeros((0, 6), np.float32)
|
||||
from .body.motion import MovingTables
|
||||
return np.zeros((0, MovingTables.DIAG_FLOATS), np.float32)
|
||||
self.stream.synchronize()
|
||||
return tables.download_diagnostics()
|
||||
|
||||
@@ -360,8 +368,11 @@ class Simulation:
|
||||
Args:
|
||||
id: body_id from add_body().
|
||||
omega: rotation speed (lattice units). ``None`` = no change.
|
||||
vx: translation velocity x (future). ``None`` = no change.
|
||||
vy: translation velocity y (future). ``None`` = no change.
|
||||
vx: translation velocity x. For a translating circle: the
|
||||
interval wall velocity of the target, only together with
|
||||
``offset``; otherwise stored in the body state (not applied
|
||||
to stationary walls). ``None`` = no change.
|
||||
vy: translation velocity y; same rules as ``vx``.
|
||||
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
|
||||
@@ -729,9 +740,20 @@ class Simulation:
|
||||
if stream is None:
|
||||
stream = self.stream
|
||||
|
||||
# Async upload action
|
||||
if upload_act and self.bodies.count > 0:
|
||||
self.bodies._upload_action_async(stream)
|
||||
# Host-only action validation and staging: a rejection here stays retryable.
|
||||
staged = upload_act and self.bodies.count > 0 and self.bodies._prepare_action_upload()
|
||||
watched = bool(moving) or self.stepper.moving is not None
|
||||
try:
|
||||
self._run_steps(steps, stream, staged, sync_obs, zero_obs, checkpoint_interval)
|
||||
except BaseException as exc:
|
||||
if watched: # device work was attempted: motion state is unknown
|
||||
self._fail_moving_run(exc)
|
||||
raise
|
||||
|
||||
def _run_steps(self, steps, stream, staged, sync_obs, zero_obs, checkpoint_interval):
|
||||
"""Device part of ``run``: action upload, obs zeroing, step group, sync."""
|
||||
if staged:
|
||||
self.bodies._upload_prepared_action_async(stream)
|
||||
|
||||
# Zero obs segments before step group
|
||||
if zero_obs:
|
||||
@@ -778,30 +800,28 @@ class Simulation:
|
||||
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."""
|
||||
(docs/moving_bodies.md 4a, 5). ``run`` fails the Simulation on any exception."""
|
||||
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()
|
||||
pending = any(motion.pending is not None for motion, _ in circles)
|
||||
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")
|
||||
if pending: # a target equal to the applied offset is consumed as a hold step
|
||||
for motion, _ in circles:
|
||||
motion.commit()
|
||||
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
|
||||
self._attach_moving()
|
||||
|
||||
def step(self, n: int = 1):
|
||||
"""Advance *n* steps (convenience for interactive use)."""
|
||||
|
||||
@@ -11,6 +11,20 @@ _spec = importlib.util.spec_from_file_location(
|
||||
p = importlib.util.module_from_spec(_spec)
|
||||
_spec.loader.exec_module(p)
|
||||
|
||||
def _load(rel, name):
|
||||
"""Load a research driver module under a unique name (no bare ``import run``)."""
|
||||
import sys
|
||||
path = Path(__file__).resolve().parents[2] / rel
|
||||
sys.path.insert(0, str(path.parent)) # for the driver's own sibling imports
|
||||
try:
|
||||
spec = importlib.util.spec_from_file_location(name, path)
|
||||
mod = importlib.util.module_from_spec(spec)
|
||||
spec.loader.exec_module(mod)
|
||||
finally:
|
||||
sys.path.remove(str(path.parent))
|
||||
return mod
|
||||
|
||||
|
||||
T = 2520
|
||||
|
||||
|
||||
@@ -93,12 +107,7 @@ def test_exact_point_inside_body_is_missing_even_with_fluid_corners():
|
||||
|
||||
def test_generated_config_passes_the_real_loader(tmp_path):
|
||||
import json
|
||||
import sys
|
||||
sys.path.insert(0, str(Path(__file__).resolve().parents[2] / "research/duetsch_re100_kc5"))
|
||||
try:
|
||||
import run
|
||||
finally:
|
||||
sys.path.pop(0)
|
||||
run = _load("research/duetsch_re100_kc5/run.py", "duetsch_run")
|
||||
from CelerisLab.config import load_lbm_config
|
||||
cfg_dict, u_max, n = run.config(20, 10.0, 2520)
|
||||
path = tmp_path / "config_lbm.json"
|
||||
@@ -111,13 +120,7 @@ def test_generated_config_passes_the_real_loader(tmp_path):
|
||||
|
||||
|
||||
def _analyze():
|
||||
import sys
|
||||
sys.path.insert(0, str(Path(__file__).resolve().parents[2] / "research/duetsch_re100_kc5"))
|
||||
try:
|
||||
import analyze
|
||||
finally:
|
||||
sys.path.pop(0)
|
||||
return analyze
|
||||
return _load("research/duetsch_re100_kc5/analyze.py", "duetsch_analyze")
|
||||
|
||||
|
||||
def test_hf_energy_is_parseval_weighted():
|
||||
|
||||
@@ -113,7 +113,8 @@ def moving_simulation(monkeypatch, *, moving=True):
|
||||
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"),
|
||||
_prepare_action_upload=lambda: True,
|
||||
_upload_prepared_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"))
|
||||
@@ -297,6 +298,44 @@ def test_run_one_advances_and_commits_after_sync(monkeypatch):
|
||||
assert ("step", 3, False) in writes and motion.applied == (0., .25)
|
||||
|
||||
|
||||
def test_zero_motion_target_is_consumed_by_run_one(monkeypatch):
|
||||
sim, writes, motion = advancing_simulation(monkeypatch)
|
||||
motion.set_target((0., 0.), (0., 0.)) # equal to applied: a hold target
|
||||
sim._attach_moving()
|
||||
assert sim.stepper.blocked == module.PENDING_BLOCKED
|
||||
sim.run(1, stream=STREAM)
|
||||
assert motion.pending is None and motion.applied == (0., 0.)
|
||||
assert sim.stepper.blocked is None and ("step", 1, False) in writes
|
||||
writes.clear()
|
||||
sim.run(3, stream=STREAM) # no stale latch
|
||||
assert ("step", 3, False) in writes
|
||||
|
||||
|
||||
@pytest.mark.parametrize("streaming,d2q9,bound", [("double_buffer", True, True),
|
||||
("double_buffer", False, False),
|
||||
("esopull", True, False)])
|
||||
def test_moving_kernels_bound_only_for_d2q9_double_buffer(monkeypatch, streaming, d2q9, bound):
|
||||
from CelerisLab.lbm import stepper as stepper_module
|
||||
monkeypatch.setattr(stepper_module.cuda, "mem_alloc", lambda n: object())
|
||||
moving = {"MovingCircleLinks", "MovingCircleCells"}
|
||||
|
||||
def get_function(name):
|
||||
if name in moving and not (d2q9 and streaming == "double_buffer"):
|
||||
raise KeyError(name) # absent from D3Q19 / EsoPull modules
|
||||
return lambda *a, **k: None
|
||||
cfg = SimpleNamespace(streaming=streaming, threads_per_block=128, is_d2q9=d2q9,
|
||||
nx=16, ny=8, nz=4, nq=9 if d2q9 else 19)
|
||||
stepper = stepper_module.LBMStepper(SimpleNamespace(), SimpleNamespace(get_function=get_function), cfg)
|
||||
assert (stepper.moving_fn is not None) == bound and (stepper.moving_cells_fn is not None) == bound
|
||||
|
||||
|
||||
def test_moving_diagnostics_empty_shape_matches_layout(monkeypatch):
|
||||
from CelerisLab.body.motion import MovingTables
|
||||
sim = make_simulation(monkeypatch)
|
||||
sim._failed, sim.stepper = None, SimpleNamespace(moving=None)
|
||||
assert sim.moving_diagnostics().shape == (0, MovingTables.DIAG_FLOATS)
|
||||
|
||||
|
||||
@pytest.mark.parametrize("failure", ("upload", "first_sync", "no_donor", "hold_upload",
|
||||
"hold_sync"))
|
||||
def test_failed_advance_fails_simulation(monkeypatch, failure):
|
||||
@@ -323,6 +362,78 @@ def test_failed_advance_fails_simulation(monkeypatch, failure):
|
||||
sim.run(1, stream=STREAM)
|
||||
|
||||
|
||||
def raising(log, name, site):
|
||||
def call(*a, **k):
|
||||
log.append(name)
|
||||
if name == site:
|
||||
raise RuntimeError(f"{name} failed")
|
||||
return call
|
||||
|
||||
|
||||
def two_circle_simulation(monkeypatch, site):
|
||||
# No owner hooks: run itself must fail the Simulation and every circle.
|
||||
sim, writes, first = advancing_simulation(monkeypatch)
|
||||
second = SimpleNamespace(motion=circle(), obj_id=1, center=(80., 40.), radius=10.)
|
||||
sim.bodies.objects.append(second)
|
||||
sim.bodies._upload_prepared_action_async = raising(writes, "action", site)
|
||||
sim.bodies.zero_obs_async = raising(writes, "obs", site)
|
||||
sim.bodies.download_obs_full_async = raising(writes, "download", site)
|
||||
sim.stepper.moving.zero_diagnostics = raising(writes, "zero_diag", site)
|
||||
step = raising(writes, "step", site)
|
||||
sim.stepper._step = lambda n, advance=False, **k: step()
|
||||
return sim, (first, second.motion), SimpleNamespace(synchronize=raising(writes, "sync", site))
|
||||
|
||||
|
||||
@pytest.mark.parametrize("steps,site", [(n, s) for n in (0, 1, 3)
|
||||
for s in ("action", "obs", "zero_diag", "step", "download", "sync")
|
||||
if not (n == 0 and s == "step")]) # run(0) launches nothing
|
||||
def test_device_failure_in_run_fails_every_circle(monkeypatch, steps, site):
|
||||
sim, motions, stream = two_circle_simulation(monkeypatch, site)
|
||||
if steps < 2:
|
||||
motions[0].set_target((0., .25), (0., .05))
|
||||
with pytest.raises(RuntimeError, match=f"{site} failed"):
|
||||
sim.run(steps, stream=stream)
|
||||
assert all(m.failed for m in motions) and sim._failed is not None
|
||||
assert sim.stepper.blocked == sim._failed
|
||||
assert motions[0].applied == (0., 0.) # nothing committed before the first sync
|
||||
with pytest.raises(RuntimeError, match="Simulation failed"):
|
||||
sim.run(1, stream=STREAM)
|
||||
|
||||
|
||||
def test_action_rejection_stays_retryable_before_device_work(monkeypatch):
|
||||
sim, writes, motion = advancing_simulation(monkeypatch)
|
||||
motion.set_target((0., .25), (0., .05))
|
||||
|
||||
def reject():
|
||||
raise ValueError("bad action")
|
||||
sim.bodies._prepare_action_upload = reject
|
||||
for steps in (0, 1):
|
||||
with pytest.raises(ValueError, match="bad action"):
|
||||
sim.run(steps, stream=STREAM)
|
||||
assert writes == [] and sim._failed is None and sim.stepper.blocked is None
|
||||
assert not motion.failed and motion.pending is not None
|
||||
sim.bodies._prepare_action_upload = lambda: True
|
||||
sim.run(1, stream=STREAM)
|
||||
assert motion.applied == (0., .25) and motion.pending is None
|
||||
|
||||
|
||||
def test_attached_tables_without_circles_fail_closed(monkeypatch):
|
||||
sim, writes, _ = advancing_simulation(monkeypatch)
|
||||
sim.bodies.objects = []
|
||||
sim.stepper.moving.zero_diagnostics = raising(writes, "zero_diag", "zero_diag")
|
||||
with pytest.raises(RuntimeError, match="zero_diag failed"):
|
||||
sim.run(0, stream=STREAM)
|
||||
assert sim._failed is not None and sim.stepper.blocked == sim._failed
|
||||
|
||||
|
||||
def test_stationary_run_failure_policy_unchanged(monkeypatch):
|
||||
sim, writes = moving_simulation(monkeypatch, moving=False)
|
||||
sim.bodies.zero_obs_async = raising(writes, "obs", "obs")
|
||||
with pytest.raises(RuntimeError, match="obs failed"):
|
||||
sim.run(2, stream=STREAM)
|
||||
assert sim._failed is None and sim.stepper.blocked is None
|
||||
|
||||
|
||||
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
|
||||
|
||||
@@ -0,0 +1,105 @@
|
||||
"""CPU tests of the predeclared VIV logic (research/viv_re150/viv.py; stdlib only)."""
|
||||
import importlib.util
|
||||
import math
|
||||
from pathlib import Path
|
||||
|
||||
_spec = importlib.util.spec_from_file_location(
|
||||
"viv_logic", Path(__file__).resolve().parents[2] / "research/viv_re150/viv.py")
|
||||
viv = importlib.util.module_from_spec(_spec)
|
||||
_spec.loader.exec_module(viv)
|
||||
|
||||
def _load(rel, name):
|
||||
"""Load a research driver module under a unique name (no bare ``import run``)."""
|
||||
import sys
|
||||
path = Path(__file__).resolve().parents[2] / rel
|
||||
sys.path.insert(0, str(path.parent)) # for the driver's own sibling imports
|
||||
try:
|
||||
spec = importlib.util.spec_from_file_location(name, path)
|
||||
mod = importlib.util.module_from_spec(spec)
|
||||
spec.loader.exec_module(mod)
|
||||
finally:
|
||||
sys.path.remove(str(path.parent))
|
||||
return mod
|
||||
|
||||
|
||||
D, U = 40, 0.04
|
||||
|
||||
|
||||
def _feed(cyc, amps, period=100):
|
||||
"""Feed sine cycles of the given amplitudes (in D); return the first non-None result."""
|
||||
n = 0
|
||||
for a in amps:
|
||||
for i in range(period):
|
||||
out = cyc.update(a * D * math.sin(2 * math.pi * (i + 1) / period))
|
||||
n += 1
|
||||
if out:
|
||||
return out, n
|
||||
return None, n
|
||||
|
||||
|
||||
def test_structure_exact_update_conserves_free_oscillation_and_equilibrium():
|
||||
s = viv.Structure(D, U)
|
||||
assert math.isclose(1 / s.f_n, 5000.0) and math.isclose(s.m, 2 * D * D)
|
||||
s.y, s.v = 0.3 * D, 0.0
|
||||
e0 = 0.5 * s.k * s.y ** 2
|
||||
for _ in range(12345):
|
||||
s.y, s.v = s.candidate(0.0)
|
||||
assert math.isclose(0.5 * s.k * s.y ** 2 + 0.5 * s.m * s.v ** 2, e0, rel_tol=1e-9)
|
||||
s.y, s.v = 2.0 / s.k, 0.0
|
||||
assert all(math.isclose(a, b, abs_tol=1e-12) for a, b in zip(s.candidate(2.0), (2.0 / s.k, 0.0)))
|
||||
assert math.isclose(s.seed_velocity(D) / s.w, viv.ELIGIBLE * D)
|
||||
|
||||
|
||||
def test_constant_amplitude_settles_then_scores_twenty_cycles():
|
||||
cyc = viv.Cycles(D)
|
||||
out, _ = _feed(cyc, [0.5] * 60)
|
||||
assert out == "scored" and len(cyc.window) == 20 and cyc.settled_at == 1
|
||||
r = viv.score(cyc.window)
|
||||
assert math.isclose(r["A_max"], 0.5, rel_tol=1e-3) and abs(r["offset"]) < 1e-3
|
||||
assert r["T_resp"] == 100 and r["within_band"] # |.500 - .520| <= .052
|
||||
|
||||
|
||||
def test_small_cycles_are_not_counted_and_y0_is_not_a_crossing():
|
||||
cyc = viv.Cycles(D)
|
||||
assert cyc.update(0.0) is None and cyc.log == []
|
||||
out, _ = _feed(cyc, [0.005] * 50)
|
||||
assert out is None and cyc.counted == 0 and cyc.log and not any(c["eligible"] for c in cyc.log)
|
||||
|
||||
|
||||
def test_scoring_failure_restarts_from_the_failing_block():
|
||||
cyc = viv.Cycles(D)
|
||||
amps = [0.5] * 20 + [0.55] * 10 + [0.55] * 40 # settle at block 1, block 2 jumps 10 %
|
||||
out, _ = _feed(cyc, amps)
|
||||
assert out == "scored"
|
||||
assert all(math.isclose(c["amp"], 0.55, rel_tol=1e-3) for c in cyc.window)
|
||||
|
||||
|
||||
def test_seed_reset_discards_incomplete_cycle_and_proof():
|
||||
cyc = viv.Cycles(D)
|
||||
_feed(cyc, [0.5] * 15)
|
||||
cyc.update(0.2 * D) # inside an incomplete cycle
|
||||
counted = cyc.counted
|
||||
cyc.discard_current()
|
||||
cyc.reset_proof()
|
||||
assert cyc.amps == [] and cyc.blocks == [] and cyc.counted == counted
|
||||
|
||||
|
||||
def test_viv_config_passes_the_real_loader(tmp_path):
|
||||
import json
|
||||
run = _load("research/viv_re150/run.py", "viv_run")
|
||||
from CelerisLab.config import load_lbm_config
|
||||
cfg_dict, nx, ny = run.config()
|
||||
path = tmp_path / "config_lbm.json"
|
||||
path.write_text(json.dumps(cfg_dict))
|
||||
cfg = load_lbm_config(str(path))
|
||||
assert (cfg.nx, cfg.ny) == (nx, ny) == (3000, 2400)
|
||||
assert cfg.collision == "MRT" and cfg.y_wall_bc == "free_slip" and cfg.velocity == .04
|
||||
assert math.isclose(cfg.viscosity, .04 * 40 / 150)
|
||||
|
||||
|
||||
def test_finite_rejects_any_nonfinite_recorded_output():
|
||||
row = (3, .1, .2, .3, .4, .5, .6, .7, 0.0)
|
||||
assert viv.finite(row)
|
||||
for i in range(len(row)):
|
||||
for bad in (float("nan"), float("inf"), -float("inf")):
|
||||
assert not viv.finite(row[:i] + (bad,) + row[i + 1:])
|
||||
Reference in New Issue
Block a user