diff --git a/README.md b/README.md index 6e48280..cec8cd4 100644 --- a/README.md +++ b/README.md @@ -394,6 +394,12 @@ before `run()`, the force will be reset to zero. For the common usage pattern | `sim.load_checkpoint(path)` | Restore from HDF5 (config must match) | | `sim.snapshot()` / `sim.restore()` | In-memory field snapshot | +In-memory snapshots save populations only. Restore keeps the current clock, +controls, observations, geometry and flags; use it only while geometry is +unchanged. This is not a restart or an FSI rollback. In particular, do not +restore a snapshot taken before a translating circle changed position. +This geometry restriction is a caller precondition, not a runtime check. + #### Low-level access | Attribute | Description | diff --git a/docs/moving_bodies.md b/docs/moving_bodies.md index ffab1f0..590f15f 100644 --- a/docs/moving_bodies.md +++ b/docs/moving_bodies.md @@ -30,6 +30,10 @@ cells (M1 `MovingCircleCells`), moving-wall Bouzidi with wall velocity 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. +- `snapshot()` / `restore()` save and restore populations only, keeping the + current clock, controls, observations, geometry and flags. Use them only + with unchanged geometry; restoring across a translation is unsupported + and is not rejected at runtime. They do not provide FSI rollback. *Supported:* dense double buffer, D2Q9, FP32 unshifted, SRT/MRT, circles, spin. *Rejected before mutation:* EsoPull, ddf shifting, D3Q19, refinement / diff --git a/docs/validation/viv_re150_predeclare.md b/docs/validation/viv_re150_predeclare.md index 59e47cf..9a9813d 100644 --- a/docs/validation/viv_re150_predeclare.md +++ b/docs/validation/viv_re150_predeclare.md @@ -29,10 +29,17 @@ Claude, Grok and Astra at 1abfe954; implemented in research/viv_re150/. - D = 40, MRT, dense double buffer; U_lat = .04, nu = U D / 150 = .010667, 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. +- Nominal 75 D × 60 D lattice allocation (3000 × 2400), based on the + Kolahdouz §4.1.1 layout. Coordinate clarification of the existing run: + the driver uses `(30 D, (ny - 1)/2)` = `(1200, 1199.5)` in lattice + coordinates, or `(30 D, 29.9875 D)`. The earlier `(30 D, 30 D)` wording + described the nominal box, not the exact executed centre. + Free-slip reflection planes are at `y = .5` and `y = ny - 1.5 = 2398.5`; + their separation is 2398 cells = 59.95 D, with the centre exactly midway, + 29.975 D from each wall. Node counts are not physical wall-to-wall lengths. + This corrects the description, without changing the historical geometry + or results. Qin's own box ([0, 40 D] × [0, 20 D], centre (10 D, 10 D)) is + the target paper's, not ours. - Boundaries: uniform Zou–He inlet, NEQ outlet (rho = 1, backflow clamp), 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 diff --git a/research/duetsch_re100_kc5/protocol.py b/research/duetsch_re100_kc5/protocol.py index 1d65e49..ca0db20 100644 --- a/research/duetsch_re100_kc5/protocol.py +++ b/research/duetsch_re100_kc5/protocol.py @@ -8,6 +8,7 @@ s/T in [0.25 + (k-1), 0.25 + k); cycle 0 is the startup quarter. from __future__ import annotations import math +from functools import lru_cache import numpy as np @@ -53,10 +54,16 @@ def bin_phase(s, T): return None +@lru_cache(maxsize=8) +def _quadrature_rule(points): + """Reuse the nodes and weights, which depend only on quadrature order.""" + return np.polynomial.legendre.leggauss(points) + + def interval_bases(n, T, points=16): """Single-step interval averages over [n, n + 1] of cos phi |cos phi| and sin phi (phi in radians = 2 pi s / T - pi / 2), Gauss-Legendre quadrature.""" - x, w = np.polynomial.legendre.leggauss(points) + x, w = _quadrature_rule(points) s = int(n) % T + 0.5 * (x + 1.0) # periodic: identical bases in every cycle phi = 2.0 * np.pi * s / T - 0.5 * np.pi c = np.cos(phi) diff --git a/research/viv_re150/run.py b/research/viv_re150/run.py index 31b9586..38445a0 100644 --- a/research/viv_re150/run.py +++ b/research/viv_re150/run.py @@ -1,8 +1,9 @@ """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). +Re 150, U* = 5, m/(rho D^2) = 2, zeta = 0, D = 40 on a nominal 75 D x 60 D +lattice (uniform inlet, NEQ outlet, free-slip y walls). The cylinder is at +(30 D, (ny - 1)/2); y-wall planes are .5 and ny - 1.5 in lattice coordinates. +Structural RHS = CME + Wen (diagnostics: raw CME, node terms). Research driver only: Simulation API, no production change. PYTHONPATH=src python research/viv_re150/run.py --out /path/outside/checkout diff --git a/tests/unit/test_duetsch_protocol.py b/tests/unit/test_duetsch_protocol.py index d70b924..a498b96 100644 --- a/tests/unit/test_duetsch_protocol.py +++ b/tests/unit/test_duetsch_protocol.py @@ -82,6 +82,35 @@ def test_interval_bases_and_morison_fit_recover_known_coefficients(): assert p.morison_fit(f, b[:, 0], b[:, 1]) == pytest.approx((2.09, 1.45), abs=1e-9) +def test_quadrature_rule_reused_across_steps_and_periods(monkeypatch): + original = np.polynomial.legendre.leggauss + orders = [] + + def counted(points): + orders.append(points) + return original(points) + + p._quadrature_rule.cache_clear() + monkeypatch.setattr(np.polynomial.legendre, "leggauss", counted) + try: + for points in (16, 8, 16): + for period in (T, 2 * T): + for n in (0, 123, period + 123): + p.interval_bases(n, period, points) + assert orders == [16, 8] + finally: + p._quadrature_rule.cache_clear() + + +def test_interval_sine_basis_matches_exact_integral(): + for period in (T, 2 * T): + half_width = math.pi / period + for n in (0, period // 4, period // 2, period - 1, 7 * period + 123): + midpoint = 2 * math.pi * ((n % period) + .5) / period - math.pi / 2 + expected = math.sin(midpoint) * math.sin(half_width) / half_width + assert p.interval_bases(n, period)[1] == pytest.approx(expected, rel=0, abs=2e-15) + + def test_geometry_repeats_exactly_across_cycles(): A = p.amplitude(20) for m in (0, 1, T // 4, T // 2 - 1, 7 * T // 12):