diff --git a/docs/architecture/engine.md b/docs/architecture/engine.md index bc1dc8c..b2b8500 100644 --- a/docs/architecture/engine.md +++ b/docs/architecture/engine.md @@ -255,16 +255,25 @@ consistent with it. **Implementation:** `src/pyflow/engine/numerics/pressure_coupling.py` (`docs/planning/roadmap.md` TASK-021 Pressure Coupling Interface, Stage 3; TASK-027 PISO Pressure Coupling, Stage 4; TASK-033 Pressure -Correction Loop, Stage 5). The interface and its MVP scheme, `PISO`, -both live there -- the fifth of the six `adr/ADR-003` components whose -registered name resolves to a real implementation. **Genuinely -multi-pass since TASK-033 (Stage 5, 2026-08-29)**, not only the single, -real dt-scaled correction pass TASK-027 shipped -- `docs/architecture/ -icds.md`'s own Pressure-Velocity Coupling entry records the full -resolution (the momentum coefficient Rhie-Chow needs, `a_P = V/dt`, -paired with the same compact Laplacian the Poisson matrix already uses). -TASK-021 also builds `src/pyflow/engine/numerics/assembly.py`, the -registry all six of these layers resolve a configured name through. +Correction Loop, TASK-034 Navier-Stokes Timestep, Stage 5). The +interface and its MVP scheme, `PISO`, both live there -- the fifth of +the six `adr/ADR-003` components whose registered name resolves to a +real implementation. **Genuinely multi-pass since TASK-033 (Stage 5, +2026-08-29)**, not only the single, real dt-scaled correction pass +TASK-027 shipped -- `docs/architecture/icds.md`'s own Pressure-Velocity +Coupling entry records the full resolution (the momentum coefficient +Rhie-Chow needs, `a_P = V/dt`, paired with the same compact Laplacian +the Poisson matrix already uses). **Periodic-boundary-aware since +TASK-034 (Stage 5, 2026-08-29)** -- `GreenGaussGradient`/ +`GreenGaussDivergence` had no periodic case at all before this, which +blocked reaching this layer at all on a periodic domain; `icds.md`'s own +entry has the full finding. TASK-021 also builds `src/pyflow/engine/ +numerics/assembly.py`, the registry all six of these layers resolve a +configured name through; `src/pyflow/engine/simulation.py`'s own +`navier_stokes_step` (TASK-034) is what actually assembles the +predictor/corrector/corrected-state sequence this layer's own Contract +implies, calling whichever `PressureCoupling` was configured rather than +this layer's MVP scheme by name. **Upgrade path:** PISO → SIMPLE / SIMPLEC / other strategies depending on transient-vs-steady-state regime (`upgrade-paths.md` diff --git a/docs/architecture/icds.md b/docs/architecture/icds.md index 0a7e8df..4ff5d6b 100644 --- a/docs/architecture/icds.md +++ b/docs/architecture/icds.md @@ -284,6 +284,30 @@ own tunables" shape `ConjugateGradientSolver` already established, not a widening of the shared interface. `docs/planning/roadmap.md` TASK-033's own Design decisions record the full numerical investigation. +**Done, TASK-034, 2026-08-29: `PISO` gained periodic-boundary support, +and this stage's own timestep assembles the whole pipeline for real.** +Found while building Stage 5 Completion Criterion 4's own "uniform flow +on a fully periodic domain" null test: `PISO`'s pressure treatment had +no periodic case at all before this -- `GreenGaussGradient`/ +`GreenGaussDivergence` raised `UnconfiguredBoundaryFaceError` +unconditionally for any periodic boundary face, which blocked even +*measuring* an already divergence-free field's divergence, let alone +correcting it. Both gained a `periodic_pairs` constructor parameter, the +same shape `CentralDifferenceDiffusion` already had since TASK-030; +`PISO` threads it to both plus its own `_rhie_chow_divergence` correction +loop, and to the Poisson matrix's own diffusion scheme, which had been +silently passed a hardcoded empty mapping regardless of what `PISO` +itself was told. Also new: `pyflow.engine.simulation.navier_stokes_step`, +the predictor/corrector/corrected-state assembly this ICD's own contract +implies but nothing before this task actually built -- momentum's own +components advance through the ordinary `step` path with no pressure +term (the predictor), the result is reassembled and handed to whichever +`PressureCoupling` was configured (the corrector), and the corrected +components replace the predictor's own. `docs/planning/roadmap.md` +TASK-034's own Design decisions record the full numerical investigation, +including the Poisson-matrix caching fix this task found needed while +measuring the Lid-Driven Cavity validation's own real runtime. + --- ## Linear Solver diff --git a/docs/implementation/config-template.yaml b/docs/implementation/config-template.yaml index a554ccd..729a25b 100644 --- a/docs/implementation/config-template.yaml +++ b/docs/implementation/config-template.yaml @@ -108,9 +108,9 @@ field_display: # Live, repeatedly-stepped simulation seeding -- distinct from field_display # above, which renders one static frame. simulation: - # Valid: null (no live simulation runs) or "gaussian_blob", the only - # built-in pattern this field currently accepts. Invalid: any other - # string. + # Valid: null (no live simulation runs), "gaussian_blob", or + # "sinusoidal_mode", the two built-in patterns this field currently + # accepts. Invalid: any other string. scalar_pattern: null # Valid: null or "uniform", the only built-in pattern this field currently # accepts. Invalid: any other string. diff --git a/docs/implementation/golden-demos.md b/docs/implementation/golden-demos.md index 79c3e53..d5747b6 100644 --- a/docs/implementation/golden-demos.md +++ b/docs/implementation/golden-demos.md @@ -257,32 +257,105 @@ momentum equations); this demo's velocity is a fixed, uniform vector from configuration, transporting the scalar the same way a wind field transports smoke without itself being computed from the smoke. -## Initial Golden Demo +## Lid-Driven Cavity -A 2D air-current simulation, corresponding to the MVP -(`docs/implementation/mvp.md`). It must: +TASK-034's own golden demo (`docs/planning/roadmap.md`, Stage 5 +Completion Criterion 8) -- **the Initial Golden Demo described below, +now built.** A square cavity, no-slip on every wall, the top wall +moving tangentially at a constant speed: the classic incompressible +Navier-Stokes benchmark, and the demo the MVP's own Definition of Done +refers to as "golden demo exists." -- construct the domain (structured 2D Cartesian mesh); -- configure the numerical components (via `src/pyflow/configuration/`, - not hardcoded -- see `adr/ADR-003-modular-numerical-strategies.md`); -- execute timesteps; -- produce measurable velocity fields; -- render the result. +"Working" means, concretely: + +- the demo *is* `examples/golden-demos/lid_driven_cavity.yaml` -- a + `numerics.boundary_conditions.north.field_values` entry + (`velocity.0`/`velocity.1`, `VectorField.component_name`) sets the + lid's own tangential-only prescribed velocity, every other wall keeps + the schema's own no-slip default; `simulation.velocity_solved: true` + with no `simulation.scalar_pattern` selects `bootstrap.py`'s own + velocity-only live path (`_add_solved_velocity_rendering`); run via + `uv run python -m pyflow run --config examples/golden-demos/lid_driven_cavity.yaml`; +- every rendered frame is a real `simulation.navier_stokes_step()` call + -- predictor, corrector, corrected state -- not only `simulation.step()` + (`_add_passive_scalar_transport`'s own `velocity_solved` path never + pressure-corrects, a genuine pre-existing gap this task's own + `navier_stokes_step` is what actually closes for a live run); +- the velocity field is rendered as arrows (`build_vector_field_arrows`, + TASK-017, rebuilt every frame the same "remove the old object, build a + new one" way the scalar demo's own mesh is) -- **the first velocity + field PyFlow has ever rendered that was *solved*, not prescribed or + seeded** (`docs/implementation/mvp.md`'s own "visualisation shows the + result", true here for the first time); +- `tests/golden/test_lid_driven_cavity.py` checks the demo is + reproducible via the real CLI, that the rendered velocity has genuinely + moved away from its at-rest initial condition, and that the same + configuration run twice is bit-identical; +- it runs headlessly via `--backend offscreen`, same as every other demo. + +**The quantitative comparison against Ghia, Ghia & Shin (1982) is not +this file's own regression test.** `tests/features/ +navier_stokes_timestep.feature`'s own scenario runs the comparison +directly against the engine at three mesh resolutions to a measured +steady state -- the most computationally expensive check in this +project, deliberately kept separate from a demo's own lightweight +reproducibility smoke test (`docs/planning/roadmap.md` TASK-034's own +Design decision). This demo's own regression test does not assert an +absolute divergence bound either, for a related, measured reason: at +this demo's own coarse mesh and early frame count, `GreenGaussDivergence`'s +naive face-averaged divergence is not the Rhie-Chow-consistent measure +`PISO`'s own corrector loop actually drives to tolerance, and is +additionally distorted near the lid's own two corner singularities (a +genuine, well-documented property of this exact benchmark, not an +artefact) -- see `tests/golden/test_lid_driven_cavity.py`'s own module +docstring for the measured finding. + +## Heat Diffusion + +TASK-034's second golden demo -- Stage 5's own reconciliation of +`docs/implementation/mvp.md`'s Validation section (2026-08-28, +maintainer's call, `docs/planning/roadmap.md` Stage 5 Completion +Criterion 8): heat diffusion *is* the diffusion equation on a +transported scalar, with no named Temperature field needed (that field, +and its buoyancy coupling, is Stage 6's TASK-035 -- a genuinely +different claim, not this one repeated). `docs/planning/ +implementation-plan.md` and `planning/data/demos.yaml` are amended in +the same change as that decision, not left to diverge. + +"Working" means, concretely: -This is the demo the MVP's own Definition of Done refers to as "golden -demo exists." +- the demo *is* `examples/golden-demos/heat_diffusion.yaml` -- a + `simulation.scalar_pattern: sinusoidal_mode` initial condition (TASK-034's + own new pattern, `bootstrap.py`'s `_simulation_scalar_initializer`): a + single spatial Fourier mode, one full wavelength across the mesh's own + x-extent, on a domain periodic on every edge and no prescribed + velocity at all -- pure diffusion, no advection; run via + `uv run python -m pyflow run --config examples/golden-demos/heat_diffusion.yaml`; +- **it validates something quantitative, per `mvp.md`'s own Validation + section rather than its Components one**: a single mode decays + exponentially at a rate set only by `fluid.diffusion_coefficient` and + the mode's own wavenumber -- an exact, closed-form answer, checked + directly (`tests/golden/test_heat_diffusion.py`, the same "bootstrap + at two frame counts, measure across real elapsed time" + shape `test_passive_scalar_transport.py` already established), not + only "heat visibly spread". Distinct from Stage 4's own diffusion + criteria (`central_difference_diffusion.feature`), which measured + spatial convergence order and conservation -- neither of which is a + decay *rate*; +- it runs headlessly via `--backend offscreen`, same as every other demo. ## Future Demos Add an entry here when a new capability is implemented, per -`docs/planning/implementation-plan.md`'s Golden Demos table (Scalar -Transport, Heat Diffusion, Poiseuille Flow, Lid-Driven Cavity, -Rayleigh-Bénard Convection, Taylor-Green Vortex, Kelvin-Helmholtz -Instability, Flow Around Cylinder, Vortex, Dam Break, 3D Cavity -- the -four added 2026-08-20, `docs/planning/backlog.md` "physical correctness -validation"). Do not add a demo entry for a capability that doesn't -exist yet -- these get written when the corresponding capability level -is reached, not speculatively ahead of it. +`docs/planning/implementation-plan.md`'s Golden Demos table (Poiseuille +Flow, Rayleigh-Bénard Convection, Taylor-Green Vortex, Kelvin-Helmholtz +Instability, Flow Around Cylinder, Vortex, Dam Break, 3D Cavity -- +**Scalar Transport, Heat Diffusion and Lid-Driven Cavity are built and +have their own sections above, not listed here any more.** The four +added 2026-08-20, `docs/planning/backlog.md` "physical correctness +validation". Do not add a demo entry for a capability that doesn't exist +yet -- these get written when the corresponding capability level is +reached, not speculatively ahead of it. **The four added 2026-08-20 exist specifically to validate physical correctness, not just demonstrate a capability** -- unlike every demo @@ -298,6 +371,7 @@ quantitative pass/fail check against that known answer -- see This file defines *what* each golden demo must do and how it's verified. `examples/golden-demos/` holds each demo's actual configuration file -- not code, per the public-API rule above. Empty Window's is -`empty_window.yaml`; the Initial Golden Demo's doesn't exist yet, since -the demo it configures doesn't either -- this is the specification, not -the implementation. +`empty_window.yaml`; the Lid-Driven Cavity's (the demo this file used to +call "Initial Golden Demo" before TASK-034 built it) is +`lid_driven_cavity.yaml` -- this is the specification, not the +implementation. diff --git a/docs/planning/roadmap.md b/docs/planning/roadmap.md index c85e25b..8134585 100644 --- a/docs/planning/roadmap.md +++ b/docs/planning/roadmap.md @@ -182,14 +182,22 @@ TASK-000..010 rows below reading **Done**. | TASK-006 Logging Framework | **Done** 2026-08-16 -- stdlib `logging`, centralised on the `pyflow` logger; every subsystem gets its logger via `get_logger(__name__)` and inherits level/formatting through the hierarchy | | TASK-007 Rendering Framework | **Done** 2026-08-16 -- wgpu/pygfx (`adr/ADR-005`) window creation, render loop, clean shutdown; canvas backend (glfw interactive / offscreen headless) selected via configuration, both behind one interface (`src/pyflow/rendering/canvas.py`) | | TASK-008 Repository Documentation | **Done** -- this row previously read "Partial -- core documents drafted; the Handbook is largely empty", stale since 2026-08-17 when all sixteen Handbook entries (E3/E4) were written; corrected 2026-08-19. All nine artifacts TASK-008 names (README, Handbook, ADRs, Capability Map, Implementation Plan, Engineering Principles, Documentation Guidelines, Practices, Dreams) exist with real content, verified directly by line count, not assumed | -| TASK-009 CLAUDE.md Hierarchy | **Done** 2026-08-19, count kept current since -- 45 files exist as of 2026-08-23 (up from 42 as of 2026-08-22: `tests/features/CLAUDE.md` joined the same day as ADR-007, missed by that day's own consistency sweep; `src/pyflow/engine/numerics/CLAUDE.md` and `tests/unit/numerics/CLAUDE.md` joined with TASK-018, 2026-08-23 -- all three real content, found and fixed while drafting TASK-018, the same "count restated in three places, one file added, count not touched" failure this row exists to warn about. 42 itself up from 40: F2 found `.claude/` and `.claude/hooks/` had no `CLAUDE.md` at all and were untracked by both inventories, fixed with real content, not placeholders; 40 itself down from 43: `assets/icons/`, `assets/shaders/`, `assets/textures/` retired 2026-08-19, E9, no document anywhere having ever stated what they were for, the same test that retired `tools/planner/`/`tools/scripts/`, 2026-08-17, E10; 43 itself down from 45 for that earlier retirement); **4** are still generic placeholders, 41 carry real content. E9's *Done when* was revised the same day it closed: no placeholder may remain in a directory that has content, not no placeholder anywhere -- all 4 remaining (`docs/tutorials/`, `examples/experiments/`, `examples/tutorials/`, `tests/performance/`) sit in directories with no real content yet, verified directly. `docs/planning/backlog.md` E9/F2 hold the file-by-file breakdown and are the authoritative count | +| TASK-009 CLAUDE.md Hierarchy | **Done** 2026-08-19, count kept current since -- 46 files exist as of 2026-08-29 (up from 45 as of 2026-08-23: `tests/fixtures/CLAUDE.md` joined with TASK-034, real content from the day it was created, the same "one file added, count updated in the same change" discipline this row exists to model; 45 itself up from 42 as of 2026-08-22: `tests/features/CLAUDE.md` joined the same day as ADR-007, missed by that day's own consistency sweep; `src/pyflow/engine/numerics/CLAUDE.md` and `tests/unit/numerics/CLAUDE.md` joined with TASK-018, 2026-08-23 -- all three real content, found and fixed while drafting TASK-018, the same "count restated in three places, one file added, count not touched" failure this row exists to warn about. 42 itself up from 40: F2 found `.claude/` and `.claude/hooks/` had no `CLAUDE.md` at all and were untracked by both inventories, fixed with real content, not placeholders; 40 itself down from 43: `assets/icons/`, `assets/shaders/`, `assets/textures/` retired 2026-08-19, E9, no document anywhere having ever stated what they were for, the same test that retired `tools/planner/`/`tools/scripts/`, 2026-08-17, E10; 43 itself down from 45 for that earlier retirement); **4** are still generic placeholders, 41 carry real content. E9's *Done when* was revised the same day it closed: no placeholder may remain in a directory that has content, not no placeholder anywhere -- all 4 remaining (`docs/tutorials/`, `examples/experiments/`, `examples/tutorials/`, `tests/performance/`) sit in directories with no real content yet, verified directly. `docs/planning/backlog.md` E9/F2 hold the file-by-file breakdown and are the authoritative count | | TASK-010 Engine Bootstrap | **Done** 2026-08-16 -- `pyflow run` loads configuration, initialises logging, opens the render window, runs the loop, exits cleanly; verified with both the offscreen backend (automated, `tests/integration/test_bootstrap.py`) and the real interactive glfw backend (manual run, a real window opened and closed cleanly). `make ci`'s pass is what TASK-010 means by "the CI pipeline passes" here, per the C2 scope decision above -- not a claim that GitHub Actions itself has run it | This paragraph previously said `make install` and `make test` were still expected to fail, pending `uv.lock` and a test suite (B2/C1) -- stale since 2026-08-16 and corrected 2026-08-19. Both now succeed: `uv.lock` is committed (B2) and `make test` runs the suite with coverage -(C1a/C1b): **653 tests at 99% as of 2026-08-29**, having been 64 when +(C1a/C1b): **672 tests at 99% as of 2026-08-29 (Stage 5 complete)**, up +from 653 after TASK-033 -- TASK-034's own ten new Gherkin scenarios +in `tests/unit/test_navier_stokes_timestep.py`, two new plain (non-BDD) +unit tests in `tests/unit/test_piso_pressure_coupling.py` proving the +new `_poisson_matrix` cache, one new periodic-aware test each in +`test_gradient_contract.py`/`test_divergence_contract.py`, and five new +Gherkin scenarios across the two new golden-demo modules +(`test_heat_diffusion.py`, `test_lid_driven_cavity.py`). Having been 64 +when this paragraph was rewritten on 2026-08-19, 202 earlier the same day, 212 after TASK-014, 226 after TASK-015, 250 after TASK-016, 287 after TASK-017, 297 after TASK-039, 315 after the Stage 2 exit audit and 337 @@ -498,11 +506,21 @@ corrector loop converging with a non-increasing recorded divergence sequence, a deliberately halving solver forcing and proving multiple genuine (strictly decreasing) passes, and exhausting the iteration limit raising `DivergenceDidNotConvergeError` rather than returning a -best-effort result. **653 tests overall** -- the feature file's own -three scenarios plus three new `NumericsConfig.pressure_correction_ -tolerance`/`pressure_correction_max_iterations` load/reject tests in -`test_configuration.py`, the same config-section-addition shape every -prior task in this run used. +best-effort result; and to 94 with TASK-034's own three feature files, +Stage 5's fifth and last task: ten scenarios in +`navier_stokes_timestep.feature` (the predictor/corrector/corrected +sequence, both null tests, determinism, the ADR-003 substitution check, +Couette flow, the Ghia cavity comparison, the Taylor-Green matched/ +mismatched pair, and kinetic-energy conservation), three in +`lid_driven_cavity.feature`, and two in `heat_diffusion.feature` -- the +MVP's own two golden demos. **672 tests overall**: 653 after TASK-033 +(recorded above), plus TASK-034's own fifteen new Gherkin scenarios +across those three feature files, plus four more plain pytest +functions -- two new (non-BDD) unit tests in `test_piso_pressure_ +coupling.py` proving the new `PISO._poisson_matrix` cache is reused +across calls and safely recomputed for a different mesh, and one new +periodic-aware test each in `test_gradient_contract.py`/ +`test_divergence_contract.py`. **All** `make ci` targets pass, verified via the Makefile itself, not only via `uv tool run` in isolation -- that is `lint`, `typecheck`, `test`, `check-docs`, @@ -7508,6 +7526,121 @@ the reading that replaces it. Navier-Stokes Timestep +**Status: Done, 2026-08-29, Stage 5's fifth and last task in build +order -- this defines the MVP.** `pyflow.engine.simulation. +navier_stokes_step` assembles TASK-031/032/033 into one real +predictor/corrector/corrected-state timestep, reached only through +whichever `PressureCoupling`/`LinearSolver` `assemble_numerics` +resolved (Criterion 13's own substitution check passes: a registered +test double is demonstrably what the timestep calls, not a hardcoded +`PISO`). Eleven scenarios in `tests/features/navier_stokes_timestep.feature`, +bound by `tests/unit/test_navier_stokes_timestep.py`, all real-engine +(no config file, no CLI run -- this is the mechanism, not a demo): + +- The predictor/corrector/corrected sequence, each part observable via + `NavierStokesStepResult`'s own three fields. +- Both null tests: a uniform, non-axis-aligned velocity field on a + fully periodic, zero-viscosity domain stays exactly the same value and + at solver-tolerance divergence over 20 real steps; fluid initially at + rest in a closed no-slip domain stays at rest to floating-point + tolerance over the same 20 steps. +- Determinism: two independent runs from the same initial state produce + bit-identical corrected velocity and pressure. +- The ADR-003 substitution check: a `PressureCoupling` test double + registered under its own name and selected by configuration is + demonstrably what `navier_stokes_step` calls. +- Couette flow: an impulsively-started channel (periodic in the flow + direction, no-slip walls, one stationary and one moving tangentially) + reaches a measured steady state (residual-based, not a step count) + whose streamwise profile matches the exact linear solution to + `abs=1e-6`, with the wall-normal component staying zero throughout. +- **The Ghia cavity comparison -- this project's most computationally + expensive test, deliberately.** Three real runs (resolutions 9, 13, + 17, chosen odd so the centreline always lands exactly on a cell-centre + column/row) to a measured steady state at Re = 100: the RMS error + against Ghia, Ghia & Shin (1982)'s own Table I (`tests/fixtures/ + ghia_1982_re100.py`, cross-checked against two independent public + reproductions of the paper's own table, not transcribed from memory) + decreases strictly across all three resolutions; the finest + resolution's own primary vortex (found by minimum velocity magnitude + in a central sub-region, avoiding the near-wall low-velocity artefact + a naive whole-domain search picks up) lands within 0.1 of Ghia's own + reference point in unit-cavity coordinates; both downstream secondary + corner vortices are detected via opposite-sign discrete vorticity in + each bottom corner's own sub-region. **Confirmed passing on a real run + before this status was written**: 1 passed in 683.82s (11m24s). +- The emergent-phenomenon pair, **Taylor-Green vortex decay, chosen over + Kelvin-Helmholtz roll-up by measurement**: at a viscosity where + physical diffusion dominates, the measured decay rate matched the + exact `2 * wavenumber**2 * viscosity` rate to within ~0.3%; at a 100x + smaller viscosity (mesh and advection scheme held fixed), the measured + rate was off by a factor of roughly 3.8 -- upwind's own numerical + diffusion dominating once the physical rate no longer does, exactly + the failure this bullet exists to catch. +- Kinetic-energy conservation: zero increase, to floating-point + precision, over 20 real steps on a divergent-initial-condition, + near-zero-viscosity, closed no-slip fixture -- checked step by step, + not net over the run. + +Two golden demos, each its own config file, feature file and +CLI-subprocess-tested regression suite: **Lid-Driven Cavity** +(`examples/golden-demos/lid_driven_cavity.yaml`, `tests/golden/ +test_lid_driven_cavity.py`) -- the MVP's own golden demo, the first +velocity field PyFlow has ever rendered that was *solved*, live, via +`bootstrap.py`'s new `_add_solved_velocity_rendering` and a real +`navier_stokes_step` every frame; and **Heat Diffusion** +(`examples/golden-demos/heat_diffusion.yaml`, `tests/golden/ +test_heat_diffusion.py`) -- a single sinusoidal mode +(`SimulationConfig.scalar_pattern`'s new `"sinusoidal_mode"` value) on a +fully periodic domain, decaying at the exact analytic rate to within +~0.6% on a real measured run. + +**Three real findings, made and closed within this task rather than +discovered afterward:** + +1. **`GreenGaussGradient`/`GreenGaussDivergence`/`PISO` had no periodic- + boundary support at all** -- found while building the periodic null + test, which cannot reach `PISO` without it (even *measuring* an + already divergence-free field's divergence raised + `UnconfiguredBoundaryFaceError` unconditionally for any periodic + face). Both operators gained a `periodic_pairs` constructor + parameter, the same shape `CentralDifferenceDiffusion` already had + since TASK-030; `PISO` gained a fifth, defaulted parameter threading + it to `_diffusion` (which had been silently passed a hardcoded `{}` + regardless of what `PISO` itself was told), `_gradient`, + `_divergence`, and its own `_rhie_chow_divergence` correction loop. + Verified directly: a uniform field measures exactly `0.0` divergence + through this path; a hand-assigned non-uniform field gives a real + nonzero value matching a hand-derivation built directly from + `mesh.wrapped_neighbour_cell`. No ADR -- the same registry-level + widening, no-interface-change shape TASK-030's own periodic addition + used. +2. **`PISO._poisson_matrix` was rebuilt from scratch every single + timestep**, even though it depends only on the fixed mesh and + pressure boundary treatment -- found while timing the cavity + validation's own first real run (n=8: 285ms/step; n=16: 3607ms/step, + dominated 70-92% by matrix construction). Caching it per `PISO` + instance (by mesh identity) cut those to 81ms and 469ms respectively + -- a 3.5x and 7.7x reduction -- and is what made the three-resolution + comparison fit inside an 11-minute test rather than an estimated + multi-hour one. `tests/unit/test_piso_pressure_coupling.py` gained + two new plain (non-BDD) unit tests proving the cache is reused across + calls on the same mesh and safely recomputed for a different one. +3. **`velocity_tangential` (Stage 5's own design question two, + resolved 2026-08-28) was never built -- `BoundaryFaceConfig. + field_values`/`field_gradients` (TASK-031c, landed the very next day) + already supply the exact mechanism that question needed.** A + per-field-name override at one wall (`velocity.0 = U, velocity.1 = 0` + at a moving lid) is already fully general, needs no new config field, + and no per-wall tangential-axis wiring inside `assembly.py` (which + stays field-name-agnostic by design). Recorded here explicitly, per + root `CLAUDE.md`'s Validation section, rather than silently building + something different from what was decided without saying so. + +`make ci` green (see this session's own run for the final test/coverage +figures, folded into the count paragraph above Stage 0); `mypy --strict` +clean; `ruff` clean. + **Pause/rewind/replay, noted here as future scope (2026-08-20, raised by the maintainer while scoping TASK-013's live zoom/pan):** not an acceptance criterion of this task, but the natural place to build it once @@ -7558,33 +7691,75 @@ replaces. ### Design questions: the shape is settled, two pieces land here -**Two above is answered** (`velocity_tangential`, 2026-08-28) -- but -whether this task or TASK-031 *builds* it depends on which first needs a -no-slip wall. Both of this task's validation cases do, so if TASK-031 -did not need one, it lands here. - -**Four's timestep half**, which lands squarely on this task even though -the rest of the configuration surface was settled at TASK-031: -Criterion 5 -runs the cavity at three resolutions, and explicit RK4's stability limit -tightens with the mesh (with `dx` for advection, with `dx` squared for -diffusion). A single fixed `numerics.timestep` cannot serve all three, -so either it becomes derivable or each resolution carries its own -- and -the derivation is stated in the scenario either way, since a silently -hand-tuned timestep per resolution is a convergence study measuring the -tuning. +**Two above is answered, but resolved differently than drafted -- a +finding, recorded rather than silently substituted.** `velocity_tangential` +(2026-08-28) was Stage 5's own answer to "a lid-driven cavity's lid is +tangential", but `BoundaryFaceConfig.field_values`/`field_gradients` +(TASK-031c) landed the very next day, after that decision, and already +supply the exact general mechanism the question needed: a per-field-name +override at one wall (`velocity.0 = U, velocity.1 = 0` at the lid, +`VectorField.component_name`), with no new config field, no per-wall +tangential-axis wiring inside `assembly.py` (which stays field-name- +agnostic by design), and no new concept to document. A no-slip +*stationary* wall needs no override at all -- `scalar_value = 0.0` +already zeroes both components identically, since normal and tangential +are both zero there. `velocity_tangential` itself was never built. +`tests/unit/test_navier_stokes_timestep.py`'s own module docstring +records the same finding next to the code it governs. + +**Four's timestep half, resolved**: `pyflow.engine.simulation. +stable_timestep(mesh, viscosity, velocity_scale, safety_factor=0.25)` -- +`min(dx/velocity_scale, dx**2/viscosity) * safety_factor`, the tighter of +the CFL and diffusive stability limits this scheme combination (explicit +RK4, first-order upwind, central-difference diffusion) actually needs. +**`0.25` is measured, not derived**: a disposable prototype swept safety +factors across a mixed advection/diffusion regime, a diffusion-dominated +one, and an advection-dominated one; `0.3` was the largest factor that +stayed stable for 500 steps in every regime tried, `0.35` already blew +up in the mixed one, and `0.25` keeps real margin below that measured +edge. Used directly by every resolution's own dt in the Couette, +Taylor-Green and Ghia cavity scenarios below -- no hand-tuned per- +resolution timestep anywhere, so the cavity's own three-resolution +comparison is a genuine convergence study, not one measuring how well +each resolution's own dt was picked. ### Artifacts Produced - `tests/features/navier_stokes_timestep.feature` -- this task's own - Acceptance Criteria. -- One `.feature` file per demo Criterion 8 requires, named for its demo, - the same pairing every golden demo already has - (`tests/features/passive_scalar_transport.feature` being the most - recent). -- A config file per demo under `examples/golden-demos/`, and a - regression test per demo under `tests/golden/` invoking it through the - real CLI as a subprocess. + Acceptance Criteria, eleven scenarios; bound by `tests/unit/ + test_navier_stokes_timestep.py`. +- `tests/features/lid_driven_cavity.feature` and `tests/features/ + heat_diffusion.feature` -- one `.feature` file per demo Criterion 8 + requires, the same pairing every golden demo already has; bound by + `tests/golden/test_lid_driven_cavity.py`/`test_heat_diffusion.py`. +- `examples/golden-demos/lid_driven_cavity.yaml` and `examples/golden-demos/ + heat_diffusion.yaml`. +- `pyflow.engine.simulation.navier_stokes_step`/`NavierStokesStepResult`/ + `stable_timestep` -- the predictor/corrector/corrected-state assembly + and its own timestep-derivation helper. +- `pyflow.bootstrap._add_solved_velocity_rendering` -- the velocity-only + live rendering path `_add_passive_scalar_transport`'s own docstring + named as this task's likely first consumer, now built. +- `SimulationConfig.scalar_pattern`'s new `"sinusoidal_mode"` value + (`ScalarTransportPattern`), the Heat Diffusion demo's own initial + condition. +- **`tests/fixtures/`, a new top-level test-data convention** (this + repository had none) -- `ghia_1982_re100.py`, the committed Ghia, + Ghia & Shin (1982) Table I reference data (u along the vertical + centreline, v along the horizontal one, the primary vortex centre), + cited and cross-checked against two independent public reproductions + of the paper's own table before being trusted (this environment has no + direct access to the original print journal -- the fixture's own + docstring states that limit explicitly rather than silently). Recorded + in `docs/repository-manifest.md` and this directory's own new + `CLAUDE.md`, per the Blast Radius rule. +- **A real periodic-boundary extension to `GreenGaussGradient`/ + `GreenGaussDivergence`/`PISO`**, not anticipated when this task was + drafted -- found necessary to make the periodic null test reachable at + all; see this task's own Design decisions below. +- **A `PISO._poisson_matrix` caching fix**, found necessary while + measuring the Ghia cavity validation's own real runtime; see this + task's own Design decisions below. - Entries in `docs/implementation/golden-demos.md` for each demo built, written when it exists rather than ahead of it, per that document's own rule. @@ -7617,18 +7792,41 @@ feature file, are the criteria. Written to cover, at minimum: - The emergent-phenomenon pair: the instability under a configuration that should produce it, and its absence under one that should not -- with the candidate chosen by measurement, per Criterion 5's note on - what the MVP's numerical diffusion may suppress. + what the MVP's numerical diffusion may suppress. **Taylor-Green vortex + decay, not Kelvin-Helmholtz roll-up -- chosen by measurement, per + Criterion 5's own instruction, not by preference.** Reuses this task's + own periodic-domain infrastructure directly and has a closed-form decay + rate (`2 * wavenumber**2 * viscosity`) to measure against rather than + needing a roll-up detector. Measured directly before being trusted: at + a viscosity where physical diffusion dominates, the measured rate + agreed with the exact rate to within ~0.3%; at a 100x smaller + viscosity (mesh and advection scheme held fixed), the measured rate + was off by a factor of roughly 3.8 -- upwind's own numerical diffusion + dominating, exactly the failure mode this bullet exists to catch. - No single step increases total kinetic energy for an inviscid, unforced, closed-domain flow -- step by step, not net over the run, since upwind's own dissipation makes the net check pass regardless. + Measured directly (not merely a priori true even for a correct + implementation): zero increase, to floating-point precision, over 20 + real steps on a divergent-initial-condition fixture at near-zero + viscosity. ### Discharges Criteria 4, 8, 9, 10, 11 and 13, entirely. Criterion 5, all bullets, including the Couette one entirely -- TASK-033 supplies the corrector loop it depends on but does not itself scenario-test it (see that -task's own Discharges). Criterion 12, its tangential-boundary and -run-length share. Criterion 6 and Criterion 7, its own share. +task's own Discharges). Criterion 12, its tangential-boundary share -- +**discharged by the `field_values` finding above, not by building +`velocity_tangential`** -- and its run-length/steadiness share: +**deliberately not a new config field.** The Ghia cavity scenario's own +residual-based steadiness detection runs directly against the engine +(`AssembledNumerics` constructed by hand, not through `PyFlowConfig` at +all), so there is no live-run config surface for it to occupy; the two +golden demos themselves never claim to reach steady state (a live +`pyflow run` is stopped by a user, or by `--max-frames`, an existing CLI +flag, not a new config field), so neither needs one either. Criterion 6 +and Criterion 7, its own share. Golden Demo @@ -7636,6 +7834,38 @@ Lid-driven cavity. This defines the MVP of PyFlow. +### Stage 5 Completion Criteria — Exit Audit + +Written 2026-08-29, TASK-034 done, the same "read every criterion back +against what actually landed" discipline Stage 3/4's own exit audits +used. Each row points at the task's own Discharges section (or, for +Criterion 5, this task's own Acceptance Criteria bullets above) for the +full record rather than re-narrating it. + +| Criterion | Verdict | +|-----------|---------| +| 1. Velocity transported by the same mechanism as every other field | **Met.** TASK-031's own subtasks a/c/d; `test_navier_stokes_timestep.py`'s own no-special-casing check still passes with `navier_stokes_step` added. | +| 2. Pressure solved from the constraint, not transported | **Met.** TASK-032. | +| 3. Divergence decreases monotonically to the configured tolerance | **Met.** TASK-033, `PISO` genuinely multi-pass. | +| 4. One timestep solves momentum and continuity together | **Met.** `navier_stokes_step`; predictor/corrector/corrected sequence, both null tests, determinism -- all real-engine scenarios, all passing. | +| 5. Physical correctness against a known answer, per case | **Met.** Couette at solver tolerance; Ghia cavity, monotonic convergence across three real resolutions plus the finest's own vortex structure (confirmed on a real 11m24s run); Taylor-Green matched/mismatched pair; kinetic energy never increasing. | +| 6. Rejection paths exercised against real bad input | **Met.** `UnconfiguredBoundaryFaceError`'s own periodic case now genuinely reachable and tested (this task's own periodic-support fix); every other rejection path already covered by TASK-041/031-033. | +| 7. Executable Gherkin criteria, `make check-scenarios` gates | **Met.** 94 scenarios across 21 feature files, `make check-scenarios` passing. | +| 8. Demonstrations: Lid-Driven Cavity and Heat Diffusion | **Met.** Both built, both with a CLI-subprocess regression test and a quantitative physical check. | +| 9. `make ci` green on a real runner | **Pending this branch's own CI run** -- green locally (`make ci`, this session), not yet checked against a real `ubuntu-latest`/`windows-latest` run, per this project's own standard of evidence (a merged PR's own `gh run` output, not a local pass alone). | +| 10. Documentation matches the tree, capability map included | **Met.** `icds.md`/`engine.md` (Pressure-Velocity Coupling), `golden-demos.md` (two new sections, "Initial Golden Demo" retired), every touched `CLAUDE.md`, both inventories, `tests/fixtures/` recorded as a new convention. Capability map: `planning/data/demos.yaml`/`capabilities.yaml` already named `demo-heat-diffusion`/`demo-lid-driven-cavity` with `validates -> capability-level-2` edges before this task landed either -- nothing to add, verified directly rather than assumed. | +| 11. `mvp.md`'s Definition of Done discharged item by item | **Met**, reading its own table back: simulation runs end-to-end (Criteria 4, 8); physical fields evolve (1, 2); boundary conditions operate (5's Couette/cavity bullets); pressure/velocity coupling works (3, in the strong sense); numerical solution is measurable (5); visualisation shows the result (8's cavity bullet, a solved field rendered live for the first time); golden demo exists (8); documentation describes the implemented functionality (10); tests verify the core behaviour (7); capability map is updated (10's own share, already current). | +| 12. Everything this stage adds is configuration-driven, validated, documented | **Met**, with one deliberate exception recorded rather than silently narrowed: run-length/steadiness stayed a validation-scenario constant, not a config field (this task's own Discharges above explain why neither the demos nor the direct-engine Ghia scenario need one). `fluid:`, the corrector-loop tunables, solved-vs-prescribed velocity, and per-field wall values (superseding `velocity_tangential`) are all real, validated, documented config surface. | +| 13. The solver runs through ADR-003's seams, checked by substitution | **Met.** `navier_stokes_step`'s own substitution scenario: a `PressureCoupling` test double registered under its own name and selected by configuration is demonstrably what gets called. | + +**Criterion 9's own caveat is the one honest gap this audit found**: +this session's own `make ci` is green, but per this project's own +Merge Gate (root `CLAUDE.md`), "mechanically green" means a real run on +both platforms, checked from the actual `gh run` output once this +branch's PR exists -- not inferred from a local pass. Recorded here +rather than silently assumed, per the Merge Gate's own fourth +requirement ("said honestly"). + --- # Stage 6 — Additional Physical Fields diff --git a/docs/planning/status.md b/docs/planning/status.md index cd5fd2f..2043516 100644 --- a/docs/planning/status.md +++ b/docs/planning/status.md @@ -15,14 +15,14 @@ demand, not part of this file. ## Progress -**37/42 tasks complete (88%)** across 14 planned stages. For the full plan, including +**38/42 tasks complete (90%)** across 14 planned stages. For the full plan, including stages below not yet broken into tasks: [roadmap.md](roadmap.md). ```mermaid pie showData title "Tasks across the roadmap" - "Done" : 37 - "Not started" : 5 + "Done" : 38 + "Not started" : 4 ``` ### Milestones @@ -35,13 +35,13 @@ pie showData ### Up next -**Stage 5 -- First Fluid Solver** is next, starting with TASK-034 (Navier-Stokes Timestep). +**Stage 5 -- First Fluid Solver** has no pending tasks recorded, but isn't marked complete -- likely awaiting its exit audit. ## Live repository facts -- **45** `CLAUDE.md` files -- **653** tests collected -- **79** Gherkin scenarios (`tests/features/*.feature`) +- **46** `CLAUDE.md` files +- **672** tests collected +- **94** Gherkin scenarios (`tests/features/*.feature`) ## Stages @@ -115,7 +115,7 @@ pie showData ### Stage 5 -- First Fluid Solver -**no status recorded** -- `████████░░` 4/5 tasks; 13 criteria defined, no status line yet +**no status recorded** -- `██████████` 5/5 tasks; 13 criteria defined, no status line yet | Task | Status | Date | Artifact | |------|--------|------|----------| @@ -123,7 +123,7 @@ pie showData | TASK-031 | Done | 2026-08-29 | `advection.py` | | TASK-032 | Done | 2026-08-29 | `src/pyflow/engine/scalar_field.py` | | TASK-033 | Done | 2026-08-29 | `PressureCoupling.correct` | -| TASK-034 -- Navier-Stokes Timestep | Not started | | | +| TASK-034 | Done | 2026-08-29 | `tests/features/navier_stokes_timestep.feature` | ### Stage 6 -- Additional Physical Fields diff --git a/docs/repository-inventory.md b/docs/repository-inventory.md index b28c3a0..7c91e3a 100644 --- a/docs/repository-inventory.md +++ b/docs/repository-inventory.md @@ -16,7 +16,7 @@ reading job and lives in the manifest. Test counts and coverage are not here either -- those come from running the suite, not from listing files. -**278 tracked files** across 45 directories; +**289 tracked files** across 46 directories; 4 are empty. ## (root) @@ -176,6 +176,8 @@ listing files. - `empty_mesh.yaml` - `empty_window.yaml` - `field_display.yaml` +- `heat_diffusion.yaml` +- `lid_driven_cavity.yaml` - `numerics_assembly.yaml` - `passive_scalar_transport.yaml` @@ -309,6 +311,9 @@ listing files. - `field_display.feature` - `first_order_upwind_advection.feature` - `fluid_configuration.feature` +- `heat_diffusion.feature` +- `lid_driven_cavity.feature` +- `navier_stokes_timestep.feature` - `neumann_boundary.feature` - `numerics_assembly.feature` - `passive_scalar_transport.feature` @@ -320,6 +325,12 @@ listing files. - `simulation_orchestrator.feature` - `velocity_field_support.feature` +## tests/fixtures + +- `CLAUDE.md` +- `__init__.py` +- `ghia_1982_re100.py` + ## tests/golden - `CLAUDE.md` @@ -329,6 +340,8 @@ listing files. - `test_empty_mesh.py` - `test_empty_window.py` - `test_field_display.py` +- `test_heat_diffusion.py` +- `test_lid_driven_cavity.py` - `test_numerics_assembly.py` - `test_passive_scalar_transport.py` @@ -380,6 +393,7 @@ listing files. - `test_main.py` - `test_mesh_contract.py` - `test_mesh_visualization.py` +- `test_navier_stokes_timestep.py` - `test_neumann_boundary.py` - `test_periodic_boundary.py` - `test_piso_pressure_coupling.py` diff --git a/docs/repository-manifest.md b/docs/repository-manifest.md index a2eb6ad..670c921 100644 --- a/docs/repository-manifest.md +++ b/docs/repository-manifest.md @@ -595,7 +595,7 @@ stage boundary, not only when something here is being edited. # tests/ -`tests/` with `unit/`, `integration/`, `golden/`, `performance/`. +`tests/` with `unit/`, `integration/`, `golden/`, `performance/`, `fixtures/`. 🟨 — 56 test modules, **605 tests, 99% coverage** (2026-08-28; two of those from a same-day CLI help-message accuracy fix in existing modules @@ -821,6 +821,25 @@ centroid moves downstream at approximately the prescribed velocity over real elapsed time (measured over two independent real `bootstrap()` runs, tolerance derived from an actual run rather than guessed), not only that rendered pixels changed. +`tests/golden/test_heat_diffusion.py` (TASK-034) binds +`tests/features/heat_diffusion.feature` -- Stage 5's own reconciliation +of `mvp.md`'s Validation section (heat diffusion as the diffusion +equation on a transported scalar): the required CLI-subprocess scenario, +plus a quantitative scenario checking a single sinusoidal mode's own RMS +amplitude decays at the exact analytic rate `Gamma * wavenumber**2` +predicts (measured over two independent real `bootstrap()` runs, +tolerance measured from an actual run: ~0.6% agreement). +`tests/golden/test_lid_driven_cavity.py` (TASK-034) binds +`tests/features/lid_driven_cavity.feature` -- the MVP's own golden demo, +the first velocity field PyFlow has ever rendered that was *solved*: the +required CLI-subprocess scenario, plus scenarios checking the rendered +velocity has genuine nonzero motion away from the lid and that the same +configuration run twice is bit-identical. Deliberately does not assert +an absolute divergence bound (found to fail even away from the lid's own +two corner singularities, since `GreenGaussDivergence`'s naive +divergence is not the Rhie-Chow-consistent measure `PISO` itself uses) +-- the real, tolerance-gated divergence claim is +`tests/features/navier_stokes_timestep.feature`'s own scenarios. `unit/` otherwise holds config/logging/rendering (D1/D2/D3), the tooling tests @@ -907,6 +926,19 @@ the configured `numerics` section, and that adding one to `field_display.yaml` renders pixel-identical output. `performance/` still empty (nothing to benchmark yet). +`fixtures/` (TASK-034, 2026-08-29) is a new top-level test-data +convention -- this repository had none until the Lid-Driven Cavity +validation needed committed reference data external to the project. +`fixtures/ghia_1982_re100.py` is the first occupant: U. Ghia, K. N. +Ghia and C. T. Shin (1982)'s own Table I at Reynolds number 100 (u +along the vertical centreline, v along the horizontal one, the primary +vortex centre), cited and cross-checked against two independent public +reproductions of the paper's own table (this environment has no direct +access to the original print journal, stated explicitly in the module's +own docstring rather than left implicit). Distinct from `unit/ +_numerics.py`/`golden/_demo.py`, which stay local machinery this project +derives itself -- see `fixtures/CLAUDE.md` for the full distinction. + The count above read "42 tests, 87% coverage" from 2026-08-16 until 2026-08-21 -- see the `src/` section above for why neither `make ci` nor `make check-claims` was ever going to notice. See `pyproject.toml`'s @@ -926,9 +958,11 @@ and collided. Roadmap TASK-003, done. 🟨 — `golden-demos/empty_window.yaml` (D5, 2026-08-16), `golden-demos/empty_mesh.yaml` (TASK-013, 2026-08-20), `golden-demos/field_display.yaml` (TASK-017, 2026-08-21), -`golden-demos/numerics_assembly.yaml` (TASK-021, 2026-08-23), and -`golden-demos/passive_scalar_transport.yaml` (TASK-030, 2026-08-28) are -the five demos so far: plain configuration files, no Python -- golden demos run +`golden-demos/numerics_assembly.yaml` (TASK-021, 2026-08-23), +`golden-demos/passive_scalar_transport.yaml` (TASK-030, 2026-08-28), +`golden-demos/heat_diffusion.yaml` and `golden-demos/lid_driven_cavity.yaml` +(both TASK-034, 2026-08-29 -- the latter the MVP's own golden demo) are +the seven demos so far: plain configuration files, no Python -- golden demos run through the public `pyflow run --config ` CLI, per `docs/implementation/golden-demos.md`'s public-API rule, so there is no demo-specific script here (an earlier `empty_window.py` was replaced by diff --git a/examples/golden-demos/CLAUDE.md b/examples/golden-demos/CLAUDE.md index 547f65e..b396ab6 100644 --- a/examples/golden-demos/CLAUDE.md +++ b/examples/golden-demos/CLAUDE.md @@ -17,9 +17,9 @@ was written, once the one thing that made it "Empty Window" (a solid background colour) became `RenderingConfig.background_color`, a real configuration option instead of code. -Five demos live here as of 2026-08-28, one per stage that has produced a -visible capability, plus one that deliberately has nothing new to -render: +Seven demos live here as of 2026-08-29 (TASK-034, Stage 5), one per +stage that has produced a visible capability, plus one that deliberately +has nothing new to render: - `empty_window.yaml` (D5, 2026-08-16), Capability Level 0's: sets `rendering.background_color`, nothing else -- everything about running @@ -46,10 +46,20 @@ render: section whose east/west edges are `periodic`. See `docs/implementation/golden-demos.md`'s own section for what "working" means concretely. - -The 2D air-current simulation is still ahead of all four, waiting on -the MVP to exist -- `docs/implementation/golden-demos.md`'s "Initial -Golden Demo" section. +- `heat_diffusion.yaml` (TASK-034, 2026-08-29), Stage 5's own + reconciliation of `mvp.md`'s Validation section: a single sinusoidal + mode (`simulation.scalar_pattern: sinusoidal_mode`, this task's own new + pattern) decaying on a fully periodic domain with no velocity at all -- + pure diffusion, checked against its own exact closed-form decay rate. +- `lid_driven_cavity.yaml` (TASK-034, 2026-08-29), Stage 5's own -- the + MVP's own golden demo (this document used to call it "the Initial + Golden Demo" before this task built it). `simulation.velocity_solved: + true` with no `scalar_pattern` selects `bootstrap.py`'s own new + velocity-only live path; the moving lid is a `numerics. + boundary_conditions.north.field_values` entry + (`velocity.0`/`velocity.1`), not a new config field -- see this task's + own roadmap entry for why `velocity_tangential` (Stage 5's own design + question two) was never built. Every demo here should follow the same shape: diff --git a/examples/golden-demos/heat_diffusion.yaml b/examples/golden-demos/heat_diffusion.yaml new file mode 100644 index 0000000..3beab4f --- /dev/null +++ b/examples/golden-demos/heat_diffusion.yaml @@ -0,0 +1,57 @@ +# Heat Diffusion golden demo (Stage 5, TASK-034) -- Stage 5's own +# reconciliation of `docs/implementation/mvp.md`'s Validation section +# (2026-08-28, maintainer's call): heat diffusion *is* the diffusion +# equation on a transported scalar, with no named Temperature field +# needed (that field, and its buoyancy coupling, is Stage 6's TASK-035, +# a genuinely different claim). +# +# A single sinusoidal mode, one full wavelength across the mesh's own +# x-extent, decays exponentially at a rate set only by the diffusion +# coefficient and the mode's own wavenumber -- an exact, checkable +# answer (`tests/features/heat_diffusion.feature`), not just "heat +# visibly spreads". No velocity at all (`SimulationConfig.velocity_ +# pattern` unset): this demo is pure diffusion, the reading +# `docs/handbook/numerical-methods/diffusion.md` describes for a single +# mode on a periodic domain. +# +# Run it exactly the way any user would: +# +# uv run python -m pyflow run --config examples/golden-demos/heat_diffusion.yaml +# +# Every boundary is periodic: the mode has no y-dependence, so north/ +# south's own periodicity is never exercised by the flow, but a +# transported field still needs *some* configured condition on every +# edge regardless, and periodic is what keeps the domain's own physical +# picture (an unbounded medium carrying one Fourier mode) honest rather +# than introducing an artificial wall the closed-form solution assumes +# away. + +mesh: + extent: [24, 8] + spacing: [0.125, 0.125] + +numerics: + timestep: 0.01 + boundary_conditions: + north: + type: periodic + south: + type: periodic + east: + type: periodic + west: + type: periodic + +simulation: + scalar_pattern: sinusoidal_mode + +fluid: + diffusion_coefficient: 0.05 + +field_display: + low_color: "#0a0a2a" + high_color: "#ff4400" + value_range: [-1.0, 1.0] + +rendering: + background_color: "#1a1a2e" diff --git a/examples/golden-demos/lid_driven_cavity.yaml b/examples/golden-demos/lid_driven_cavity.yaml new file mode 100644 index 0000000..c4444ed --- /dev/null +++ b/examples/golden-demos/lid_driven_cavity.yaml @@ -0,0 +1,61 @@ +# Lid-Driven Cavity golden demo (Stage 5, TASK-034) -- the MVP's own +# golden demo (`docs/implementation/mvp.md`, `docs/implementation/ +# golden-demos.md`'s "Initial Golden Demo"). A square cavity, no-slip on +# every wall, the top ("north") wall moving tangentially at a constant +# speed -- the classic incompressible Navier-Stokes benchmark, rendered +# live as it solves. +# +# Run it exactly the way any user would: +# +# uv run python -m pyflow run --config examples/golden-demos/lid_driven_cavity.yaml +# +# `simulation.velocity_solved: true` with no `simulation.scalar_pattern` +# selects `bootstrap.py`'s own velocity-only live path +# (`_add_solved_velocity_rendering`, added by this task): every frame +# advances one real `navier_stokes_step` -- predictor, corrector, +# corrected state -- and redraws the corrected velocity as arrows. The +# quantitative comparison against Ghia, Ghia & Shin (1982) is +# `tests/features/navier_stokes_timestep.feature`'s own scenario, run +# directly against the engine at several resolutions to reach a measured +# steady state; this config is the reproducible, visible demonstration +# `docs/implementation/golden-demos.md`'s public-API rule requires, not +# a second copy of that validation. +# +# The moving lid is expressed as `numerics.boundary_conditions.north. +# field_values` -- `velocity.0`/`velocity.1` name momentum's own two +# decomposed components (`VectorField.component_name`), so `u = 1.0, +# v = 0.0` at the lid is a tangential-only prescribed velocity, no +# penetration. Every other wall is the schema's own no-slip default +# (`scalar_value: 0.0` on both components alike). + +mesh: + extent: [16, 16] + spacing: [0.0625, 0.0625] + +numerics: + timestep: 0.008 + boundary_conditions: + north: + type: dirichlet + field_values: + velocity.0: 1.0 + velocity.1: 0.0 + south: + type: dirichlet + east: + type: dirichlet + west: + type: dirichlet + +simulation: + velocity_solved: true + +fluid: + viscosity: 0.01 + +field_display: + arrow_color: "#00e0ff" + arrow_scale: 0.05 + +rendering: + background_color: "#1a1a2e" diff --git a/src/pyflow/bootstrap.py b/src/pyflow/bootstrap.py index fed1bae..c083dcb 100644 --- a/src/pyflow/bootstrap.py +++ b/src/pyflow/bootstrap.py @@ -7,7 +7,12 @@ TASK-030 (Stage 4, 2026-08-28):** `_add_passive_scalar_transport` wires a real `simulation.step()` into the render loop, one timestep per rendered frame, whenever `config.simulation.scalar_pattern` is set. -Every other configuration still renders without stepping anything. +**A second live path arrived in TASK-034 (Stage 5, 2026-08-29):** +`_add_solved_velocity_rendering` wires a real, genuinely pressure- +corrected `simulation.navier_stokes_step()` into the render loop instead, +whenever `config.simulation.velocity_solved` is set with no +`scalar_pattern` -- the Lid-Driven Cavity demo's own shape. Every other +configuration still renders without stepping anything. This docstring read "No simulation functionality -- Stage 0's job..." until the 2026-08-28 Stage 4 exit audit, in a module that by then @@ -47,6 +52,7 @@ from pyflow.engine.mesh import Mesh, StructuredCartesianMesh from pyflow.engine.numerics.assembly import assemble_numerics from pyflow.engine.scalar_field import ScalarField +from pyflow.engine.simulation import navier_stokes_step from pyflow.engine.simulation import step as simulation_step from pyflow.engine.vector_field import VectorField from pyflow.rendering import RenderWindow @@ -109,6 +115,18 @@ def _simulation_scalar_initializer( scalar_pattern` (TASK-030) -- the live-simulation counterpart to `_scalar_display_initializer` above, sharing its "derive shape from mesh bounds, don't add a config field for it" reasoning. + + **`"sinusoidal_mode"` (TASK-034, Stage 5) is the Heat Diffusion + golden demo's own initial condition** -- a single spatial Fourier + mode, one full wavelength across the mesh's own x-extent + (`wavenumber = 2*pi / domain_width`, the same "derived from mesh + bounds" precedent `"gaussian_blob"`'s own `sigma` already sets), with + no y-dependence. This is the one initial condition PyFlow's diffusion + equation has a closed-form solution for at all: a single mode decays + exponentially at a rate `Gamma * wavenumber**2`, set by the diffusion + coefficient and the mode's own wavenumber alone -- `tests/features/ + heat_diffusion.feature`'s own criterion measures exactly that rate + against this closed form. """ if pattern == "gaussian_blob": min_x, min_y, max_x, max_y = bounds @@ -117,6 +135,11 @@ def _simulation_scalar_initializer( center_y = (min_y + max_y) / 2 sigma = 0.08 * domain_width return lambda x, y: math.exp(-((x - center_x) ** 2 + (y - center_y) ** 2) / (2 * sigma**2)) + if pattern == "sinusoidal_mode": + min_x, _min_y, max_x, _max_y = bounds + domain_width = max_x - min_x + wavenumber = 2 * math.pi / domain_width + return lambda x, y: math.sin(wavenumber * (x - min_x)) raise ValueError(f"unknown simulation scalar pattern: {pattern!r}") # pragma: no cover @@ -231,6 +254,71 @@ def _advance() -> None: return _advance +def _add_solved_velocity_rendering( + window: RenderWindow, mesh: Mesh, config: PyFlowConfig +) -> Callable[[], None]: + """Wires a real `simulation.navier_stokes_step()` into a live `pyflow + run` (TASK-034, Stage 5) -- the mechanism the Lid-Driven Cavity + golden demo needs and no demo before it does: a *solved* velocity + field, rendered live, with no scalar alongside it at all. + + `_add_passive_scalar_transport`'s own docstring named this gap in + advance (TASK-031, 2026-08-29): "a velocity-only live run has + nothing this function knows how to render yet (no vector-arrow-per- + frame path exists)... revisit when a demo genuinely needs + velocity-only live rendering (TASK-034's own Lid-Driven Cavity is the + likely first)". This is that revisit -- `build_vector_field_arrows` + (TASK-017) already existed for a *static* vector display; this + function is what rebuilds it every frame, the same "remove the old + `gfx.Line`, build a new one" shape `_add_passive_scalar_transport` + already uses for its own scalar mesh. + + **Uses `navier_stokes_step`, not plain `step`** -- the real + difference from `_add_passive_scalar_transport`'s own `velocity_ + solved` path, which only ever transports velocity's components like + an ordinary scalar and never pressure-corrects them (a genuine, + pre-existing gap in that path, out of this task's own scope to + close: nothing before TASK-034 had a corrector loop to call). A + demo using this function is genuinely incompressible, frame by + frame, not merely self-advected. + """ + assert window.assembled_numerics is not None + numerics = window.assembled_numerics + + velocity_initializer = _simulation_velocity_initializer( + config.simulation.velocity_pattern, config.simulation.velocity + ) + velocity_field = VectorField( + mesh, "velocity", num_components=2, initial_value=velocity_initializer + ) + state: dict[str, Field] = {c.name: c for c in velocity_field.decompose()} + window.simulation_fields = state + + rendered_object = build_vector_field_arrows( + velocity_field, config.field_display.arrow_color, config.field_display.arrow_scale + ) + if rendered_object is not None: + rendered_object.local.position = (0.0, 0.0, _ARROWS_Z) + window.scene.add(rendered_object) + + def _advance() -> None: + nonlocal state, rendered_object, velocity_field + result = navier_stokes_step(state, "velocity", numerics, config.numerics.timestep) + state = result.fields + window.simulation_fields = state + velocity_field = result.corrected_velocity + if rendered_object is not None: + window.scene.remove(rendered_object) + rendered_object = build_vector_field_arrows( + velocity_field, config.field_display.arrow_color, config.field_display.arrow_scale + ) + if rendered_object is not None: + rendered_object.local.position = (0.0, 0.0, _ARROWS_Z) + window.scene.add(rendered_object) + + return _advance + + def _add_field_display( window: RenderWindow, mesh: Mesh, field_display: FieldDisplayConfig ) -> _Bounds: @@ -354,7 +442,18 @@ def bootstrap( config.field_display.scalar_pattern is not None or config.field_display.vector_pattern is not None ) - run_simulation = config.simulation.scalar_pattern is not None + run_scalar_simulation = config.simulation.scalar_pattern is not None + # TASK-034 (Stage 5): a velocity-only live run -- solved, rendered as + # arrows, no scalar alongside it (`_add_solved_velocity_rendering`'s + # own docstring). Mutually exclusive with `run_scalar_simulation`, + # the same "one live-simulation path per run" shape TASK-030 already + # established -- `_add_passive_scalar_transport`'s own `velocity_ + # solved` still covers a solved velocity carrying a scalar alongside + # it, unaffected by this addition. + run_velocity_only_simulation = ( + config.simulation.velocity_solved and config.simulation.scalar_pattern is None + ) + run_simulation = run_scalar_simulation or run_velocity_only_simulation on_frame: Callable[[], None] | None = None if config.rendering.show_mesh or show_fields or run_simulation: # TASK-013/017: visualise the configured mesh's grid and/or its @@ -374,12 +473,14 @@ def bootstrap( window.scene.add(build_mesh_grid_line(mesh, config.rendering.grid_color)) if show_fields: bounds = _add_field_display(window, mesh, config.field_display) - if run_simulation: + if run_scalar_simulation: # TASK-030: the first config that wires a real `simulation. # step()` into this run's own render loop, one timestep per # rendered frame -- every capability before it only ever # rendered one static frame. on_frame = _add_passive_scalar_transport(window, mesh, config) + elif run_velocity_only_simulation: + on_frame = _add_solved_velocity_rendering(window, mesh, config) fit_camera_to_bounds(window.camera, bounds) window.apply_camera_config() diff --git a/src/pyflow/configuration/CLAUDE.md b/src/pyflow/configuration/CLAUDE.md index e0a142c..de263cd 100644 --- a/src/pyflow/configuration/CLAUDE.md +++ b/src/pyflow/configuration/CLAUDE.md @@ -153,6 +153,26 @@ alongside any transported scalar and are reassembled (`VectorField.assemble`) after every frame -- still requires a `scalar_pattern` too, since there is no velocity-only live rendering path yet (`docs/planning/roadmap.md` TASK-031's own Status note). +**That gap is closed by TASK-034 (Stage 5, 2026-08-29)**: +`bootstrap.py`'s new `_add_solved_velocity_rendering` is the +velocity-only path this note anticipated, selected when `velocity_solved` +is `True` and `scalar_pattern` is `None` -- the Lid-Driven Cavity demo's +own shape. Uses `navier_stokes_step`, not plain `step`, so this path is +genuinely pressure-corrected every frame; `_add_passive_scalar_transport`'s +own `velocity_solved` path (a scalar *and* a solved velocity together) +is unaffected and still uses plain `step`, a real, pre-existing gap this +task did not close, since nothing before TASK-034 had a corrector loop +to call there either. + +**`ScalarTransportPattern` gained a second value, `"sinusoidal_mode"` +(TASK-034, added 2026-08-29)** -- the Heat Diffusion demo's own initial +condition: a single spatial Fourier mode, one full wavelength across the +mesh's own x-extent (`wavenumber = 2*pi / domain_width`, the same +"derived from mesh bounds" precedent every other pattern here follows), +with no y-dependence. The one initial condition PyFlow's diffusion +equation has a closed-form decay-rate answer for at all -- P-016 permits +this addition directly, per Criterion 12's own "a member is added [to a +pattern set] precisely because a demo needs it". **`generator.py`'s `generate_config_yaml` (TASK-039, added 2026-08-21) is `loader.py` run in reverse**: `load_config` turns YAML into a diff --git a/src/pyflow/configuration/schema.py b/src/pyflow/configuration/schema.py index 7ea6ddd..1090564 100644 --- a/src/pyflow/configuration/schema.py +++ b/src/pyflow/configuration/schema.py @@ -327,7 +327,7 @@ def validate(self) -> None: ) -ScalarTransportPattern = Literal["gaussian_blob"] +ScalarTransportPattern = Literal["gaussian_blob", "sinusoidal_mode"] VelocityPrescriptionPattern = Literal["uniform"] _VALID_SCALAR_TRANSPORT_PATTERNS = frozenset(get_args(ScalarTransportPattern)) diff --git a/src/pyflow/engine/CLAUDE.md b/src/pyflow/engine/CLAUDE.md index 0003815..211f356 100644 --- a/src/pyflow/engine/CLAUDE.md +++ b/src/pyflow/engine/CLAUDE.md @@ -1059,6 +1059,63 @@ than a best-effort result) are `tests/features/ pressure_correction_loop.feature`, bound by `tests/unit/ test_pressure_correction_loop.py`. +**`PISO` gained periodic-boundary support and a real performance fix, +both TASK-034 (Stage 5, 2026-08-29).** `GreenGaussGradient`/ +`GreenGaussDivergence` had no periodic case at all before this -- +`UnconfiguredBoundaryFaceError` unconditionally for any periodic +boundary face, found while building TASK-034's own mandated "uniform +flow on a fully periodic domain" null test: that scenario cannot reach +`PISO` at all without this, since even *measuring* an already +divergence-free field's divergence goes through `GreenGaussDivergence`. +Both gained a `periodic_pairs` constructor parameter -- the same shape +`CentralDifferenceDiffusion` already had since TASK-030 -- substituting +`mesh.wrapped_neighbour_cell` for a periodic face's own missing +mesh-level neighbour before falling through the ordinary interior-face +formula; `GreenGaussGradient`'s own periodic branch needs no distance +term to double, unlike diffusion's central-difference formula, since +Green-Gauss face averaging never divides by distance. **Verified +directly before being trusted**: a uniform, non-axis-aligned velocity +field on a genuinely periodic mesh measures exactly `0.0` divergence +through this path (float-exact); a hand-assigned non-uniform field gives +a real nonzero divergence matching a value hand-derived directly from +`mesh.wrapped_neighbour_cell`, proving the periodic branch reads the +real wrapped neighbour rather than silently skipping the face (which +would also avoid raising). `PISO` itself gained a fifth, defaulted +constructor parameter (`periodic_pairs`, empty by default -- every +existing call site passing only `(linear_solver, boundary_conditions)` +keeps working unchanged), threaded to `_diffusion` (the Poisson matrix, +which had been silently passed a hardcoded `{}` regardless of what +`PISO` itself was told -- only diffusion's own periodic support was ever +reachable before this), `_gradient`, `_divergence`, and +`_rhie_chow_divergence`'s own per-face loop (a periodic face now gets +the same Rhie-Chow correction an interior face does). +`register_pressure_coupling`'s own factory widened from four arguments +to five (`_resolve_with_five_arguments`, a new generic helper alongside +`_resolve_with_four_arguments` in `assembly.py`) to thread the same +`periodic_pairs` mapping advection/diffusion already receive. + +**`_poisson_matrix` is now cached per `PISO` instance, not rebuilt on +every `correct` call -- found while measuring the Lid-Driven Cavity +validation's own real runtime, not a speculative optimisation.** The +construction is `O(num_cells * num_faces)` (one full `self._diffusion. +flux` call per matrix column); timed directly with a disposable +prototype, it dominated 70-90% of measured per-timestep cost at MVP +cavity mesh sizes (8x8 through 16x16) even though the matrix depends +only on `mesh` and this instance's own fixed pressure boundary +treatment -- never on the current velocity, pressure, or `dt` -- so +rebuilding it every timestep was never buying correctness. Cached by +mesh *identity*: a real run always hands `correct` the same mesh object +every timestep (the common case, and a cache hit), while a genuinely +different mesh object safely recomputes rather than serving a stale +matrix. Cut measured per-timestep cost by roughly 3.5x at 8x8 and 7.7x +at 16x16 in the same prototype, which is what made the cavity +validation's own three-resolution comparison fit inside a real (if +still substantial) test budget at all -- +`tests/unit/test_piso_pressure_coupling.py`'s own two new plain (non-BDD) +unit tests prove the reuse and the safe-recompute-on-a-different-mesh +case directly, since caching correctness is an implementation detail, +not a new physical-correctness claim needing its own Gherkin scenario. + **`gradient.py`/`divergence.py`** (TASK-018, Stage 3, interface-only until TASK-027) hold `GradientScheme`/`DivergenceScheme` -- two of the three operators (with `source.py`) that jointly compute the Flux layer @@ -1405,3 +1462,57 @@ is the mechanism a future demo (TASK-030) is built on top of, not a demo itself. `tests/unit/test_simulation.py` binds it directly, per `tests/unit/CLAUDE.md`'s own scope, rather than living under `tests/golden/`. + +**`simulation.py` gained `navier_stokes_step`/`NavierStokesStepResult` +and `stable_timestep` (TASK-034, Stage 5, 2026-08-29) -- the assembly +this task's own Purpose names: "assemble the first three tasks into one +incompressible Navier-Stokes timestep".** `navier_stokes_step(fields, +velocity_field_name, numerics, dt)` is the fractional-step sequence +Stage 5's own design question five settled: momentum's own two +components advance through the *existing* `step` path above (self- +advected by the current velocity, no pressure term at all -- the +predictor), the advanced components are reassembled into a provisional +`VectorField` and handed to `numerics.pressure_coupling.correct` (the +corrector), and the corrected components replace the predictor's own +inside the returned state (the corrected state). Returns a +`NavierStokesStepResult` exposing all three parts separately +(`fields`, `provisional_velocity`, `corrected_velocity`, `pressure`) -- +Stage 5 Completion Criterion 4's own "each part observable, not only the +end state". + +**`velocity_field_name` is an explicit parameter, not a fixed +convention -- the same structural guarantee `step` itself already +carries (Stage 5 Completion Criterion 1) stays true with this function +added.** `tests/features/velocity_field_support.feature`'s own +"no `"velocity"` string literal, no `VectorField` isinstance check, no +hardcoded component-name pair" check runs against this whole module's +source unchanged; `navier_stokes_step` computes both component names +via `VectorField.component_name(velocity_field_name, i)` rather than +ever writing either literally. `bootstrap.py` is still the one place +that legitimately knows a live run's velocity field is conventionally +named `"velocity"` (`src/pyflow/CLAUDE.md`'s own rule) -- it is the +caller that supplies the literal, not this module. + +**`stable_timestep(mesh, viscosity, velocity_scale, safety_factor=0.25)` +resolves Stage 5's own design question four's timestep half**: the +tighter of the CFL limit (`dx / velocity_scale`) and the diffusive limit +(`dx**2 / viscosity`), scaled by `safety_factor`. **`0.25` is a measured +constant, not a derived one** -- a disposable prototype swept safety +factors across three regimes (mixed advection/diffusion, diffusion- +dominated, advection-dominated) before settling on it: `0.3` was the +largest factor that stayed stable for 500 steps in every regime tried, +`0.35` already blew up in the mixed one. Used by every real-engine +scenario in `tests/unit/test_navier_stokes_timestep.py` that needs a +resolution-appropriate timestep (Couette, Taylor-Green, the Ghia cavity +comparison) -- no hand-tuned per-resolution `dt` anywhere, which is what +makes the cavity's own three-resolution comparison a genuine convergence +study rather than one measuring how well each resolution's timestep was +picked. + +**Found while building TASK-034's own mandated periodic null test, not +anticipated in advance: `GreenGaussGradient`/`GreenGaussDivergence`/ +`PISO` had no periodic-boundary support at all.** See +`engine/numerics/CLAUDE.md`'s own `PISO` entry, below, for the full +finding and fix -- summarised here only because it is what makes the +"uniform flow on a fully periodic domain" scenario in `tests/features/ +navier_stokes_timestep.feature` reachable at all. diff --git a/src/pyflow/engine/numerics/CLAUDE.md b/src/pyflow/engine/numerics/CLAUDE.md index 49071da..81b4bb5 100644 --- a/src/pyflow/engine/numerics/CLAUDE.md +++ b/src/pyflow/engine/numerics/CLAUDE.md @@ -185,6 +185,22 @@ composed `Gradient`/`Divergence` pair TASK-027 tried and measured failing), and `DivergenceDidNotConvergeError` as the outer loop's own honest-exhaustion counterpart to `PressureSolveDidNotConvergeError`. +**`pressure_coupling.py`/`gradient.py`/`divergence.py` all changed again +the same day (TASK-034): periodic-boundary support, and a real +performance fix.** `register_pressure_coupling`'s own factory widens +once more, from four arguments to five (`_resolve_with_five_arguments`, +a new generic helper alongside `_resolve_with_four_arguments` rather +than widening it, since diffusion still needs only four) -- +`periodic_pairs`, the same mapping advection/diffusion already receive, +needed because `GreenGaussGradient`/`GreenGaussDivergence` had no +periodic case at all before this task and raised unconditionally for any +periodic boundary face. `PISO` also gained a real performance fix found +while measuring the Ghia cavity validation's own runtime: `_poisson_matrix` +is now cached per instance rather than rebuilt every `correct` call, +since it depends only on the fixed mesh and this instance's own pressure +boundary treatment. See `src/pyflow/engine/CLAUDE.md`'s own `PISO` entry +for the full finding on both. + Full design rationale -- why a subpackage, why every operator takes `Field` rather than a concrete subclass, why the return shapes split face-valued (Advection/Diffusion) from cell-valued diff --git a/src/pyflow/engine/numerics/assembly.py b/src/pyflow/engine/numerics/assembly.py index 7a22d7b..0510dc2 100644 --- a/src/pyflow/engine/numerics/assembly.py +++ b/src/pyflow/engine/numerics/assembly.py @@ -66,9 +66,13 @@ geometry (`StructuredCartesianMesh.wrapped_neighbour_cell`), not a prescribed value, so it never becomes a `BoundaryCondition` instance at all, and a concrete scheme consults `periodic_pairs` itself rather than -`assemble_numerics` (or the orchestrator) special-casing it. Not exposed -on `AssembledNumerics` itself -- only advection/diffusion need it, the -same "no field nothing reads" discipline the rest of that dataclass +`assemble_numerics` (or the orchestrator) special-casing it. **Also +threaded into pressure_coupling's own factory since TASK-034 +(2026-08-29)** -- `PISO`'s pressure treatment had no periodic case at +all before that task (`pressure_coupling.py`'s own entry). Not exposed +on `AssembledNumerics` itself -- advection, diffusion and +pressure_coupling each receive it at construction instead, the same "no +field nothing reads directly" discipline the rest of that dataclass already follows. """ @@ -154,7 +158,11 @@ class AssembledNumerics: _time_integrator_registry: dict[str, Callable[[], TimeIntegrator]] = {} _linear_solver_registry: dict[str, Callable[[float, int], LinearSolver]] = {} _pressure_coupling_registry: dict[ - str, Callable[[LinearSolver, Mapping[str, BoundaryCondition], float, int], PressureCoupling] + str, + Callable[ + [LinearSolver, Mapping[str, BoundaryCondition], float, int, Mapping[str, str]], + PressureCoupling, + ], ] = {} _boundary_condition_registry: dict[str, Callable[[BoundaryFaceConfig], BoundaryCondition]] = {} @@ -235,22 +243,29 @@ def register_linear_solver(name: str, factory: Callable[[float, int], LinearSolv def register_pressure_coupling( name: str, factory: Callable[ - [LinearSolver, Mapping[str, BoundaryCondition], float, int], PressureCoupling + [LinearSolver, Mapping[str, BoundaryCondition], float, int, Mapping[str, str]], + PressureCoupling, ], ) -> None: """Make `name` resolve to `factory(linear_solver, boundary_conditions, - pressure_correction_tolerance, pressure_correction_max_iterations)` - in future `assemble_numerics` calls -- `boundary_conditions` is the - same face-name-keyed mapping `AssembledNumerics.boundary_conditions` - carries (velocity's own boundary conditions), added in TASK-027 so a - concrete strategy's own `DivergenceScheme` can be boundary-aware at - construction, the same "constructed with it, not handed it after the - fact" reasoning `boundary_conditions`/`diffusion_coefficient` already - established for advection/diffusion. The tolerance/iterations pair - (TASK-033, added 2026-08-29) is `NumericsConfig. - pressure_correction_tolerance`/`pressure_correction_max_iterations` -- - a corrector *loop*'s own outer convergence tunables, "outer-loop state - the strategy owns" (Stage 5's own design question three). + pressure_correction_tolerance, pressure_correction_max_iterations, + periodic_pairs)` in future `assemble_numerics` calls -- + `boundary_conditions` is the same face-name-keyed mapping + `AssembledNumerics.boundary_conditions` carries (velocity's own + boundary conditions), added in TASK-027 so a concrete strategy's own + `DivergenceScheme` can be boundary-aware at construction, the same + "constructed with it, not handed it after the fact" reasoning + `boundary_conditions`/`diffusion_coefficient` already established for + advection/diffusion. The tolerance/iterations pair (TASK-033, added + 2026-08-29) is `NumericsConfig.pressure_correction_tolerance`/ + `pressure_correction_max_iterations` -- a corrector *loop*'s own outer + convergence tunables, "outer-loop state the strategy owns" (Stage 5's + own design question three). **`periodic_pairs` (TASK-034, added + 2026-08-29) is the same mapping advection/diffusion already receive** + -- `PISO`'s own pressure treatment had no periodic case at all before + this task (`pressure_coupling.py`'s own entry), needed to make Stage 5 + Completion Criterion 4's "uniform flow on a fully periodic domain" + null test reachable at all. """ _register(_pressure_coupling_registry, name, factory, "pressure_coupling") @@ -303,22 +318,22 @@ def _resolve_with_four_arguments[T, A, B, C, D]( argument_d: D, component: str, ) -> T: - """Same as `_resolve_with_two_arguments`, for the two components whose - factory needs four constructor arguments rather than two: diffusion + """Same as `_resolve_with_two_arguments`, for the one component whose + factory needs exactly four constructor arguments: diffusion (`boundary_conditions`, `periodic_pairs`, `diffusion_coefficient`, - `coefficient_overrides` since TASK-031b, 2026-08-29) and - pressure_coupling (`linear_solver`, `boundary_conditions`, - `pressure_correction_tolerance`, `pressure_correction_max_iterations` - since TASK-033, 2026-08-29 -- the outer corrector loop's own tunables, - "outer-loop state the strategy owns"). Kept as its own generic helper - instead of widening `_resolve_with_two_arguments` itself, since - advection/linear_solver still only need two arguments each and a - shared four-argument signature would force both to pass unused ones. - **Replaces `_resolve_with_three_arguments`, TASK-030's own three-argument - helper** -- diffusion was its only caller, and TASK-031b's own fourth - argument left it with none; deleted in the same change as genuinely - dead code, the same "no remaining caller" reasoning this docstring's - predecessor already applied to `_resolve_with_argument`. + `coefficient_overrides` since TASK-031b, 2026-08-29). Kept as its own + generic helper instead of widening `_resolve_with_two_arguments` + itself, since advection/linear_solver still only need two arguments + each and a shared four-argument signature would force both to pass + unused ones. **Replaces `_resolve_with_three_arguments`, TASK-030's + own three-argument helper** -- diffusion was its only caller, and + TASK-031b's own fourth argument left it with none; deleted in the + same change as genuinely dead code, the same "no remaining caller" + reasoning this docstring's predecessor already applied to + `_resolve_with_argument`. **pressure_coupling used to be a second + user, from TASK-033 (2026-08-29) until TASK-034 (2026-08-29, the same + day) widened its own factory to five arguments** -- see + `_resolve_with_five_arguments` below. """ factory = registry.get(name) if factory is None: @@ -326,6 +341,34 @@ def _resolve_with_four_arguments[T, A, B, C, D]( return factory(argument_a, argument_b, argument_c, argument_d) +def _resolve_with_five_arguments[T, A, B, C, D, E]( + registry: Mapping[str, Callable[[A, B, C, D, E], T]], + name: str, + argument_a: A, + argument_b: B, + argument_c: C, + argument_d: D, + argument_e: E, + component: str, +) -> T: + """Same as `_resolve_with_four_arguments`, for the one component whose + factory needs five constructor arguments: pressure_coupling + (`linear_solver`, `boundary_conditions`, `pressure_correction_tolerance`, + `pressure_correction_max_iterations`, `periodic_pairs` since TASK-034, + 2026-08-29 -- the same mapping advection/diffusion already receive, + needed so `PISO`'s own pressure treatment can reach a periodic + boundary face without raising; `pressure_coupling.py`'s own entry has + the full reasoning). Kept as its own generic helper rather than + widening `_resolve_with_four_arguments` itself, since diffusion still + needs only four arguments and a shared five-argument signature would + force it to pass an unused one. + """ + factory = registry.get(name) + if factory is None: + raise UnknownSchemeError(f"no {component} implementation registered under {name!r}") + return factory(argument_a, argument_b, argument_c, argument_d, argument_e) + + def assemble_numerics( config: NumericsConfig, diffusion_coefficient: float = 1.0, @@ -410,13 +453,14 @@ def assemble_numerics( config.linear_solver_max_iterations, "linear_solver", ) - pressure_coupling = _resolve_with_four_arguments( + pressure_coupling = _resolve_with_five_arguments( _pressure_coupling_registry, config.pressure_coupling, linear_solver, boundary_conditions, config.pressure_correction_tolerance, config.pressure_correction_max_iterations, + periodic_pairs, "pressure_coupling", ) diff --git a/src/pyflow/engine/numerics/divergence.py b/src/pyflow/engine/numerics/divergence.py index 5849d6d..9fcaa6b 100644 --- a/src/pyflow/engine/numerics/divergence.py +++ b/src/pyflow/engine/numerics/divergence.py @@ -79,10 +79,35 @@ class GreenGaussDivergence(DivergenceScheme): positive outward); a Neumann (`"gradient"`) condition extrapolates zero-order from the owner's own normal-component velocity, the same convention `FirstOrderUpwindAdvection`'s own Neumann handling uses. + + **Periodic-aware the same way `CentralDifferenceDiffusion` is + (TASK-030), added by TASK-034 (Stage 5) once `PISO` needed it.** A + fully periodic domain's own null test (`docs/planning/roadmap.md` + Stage 5 Completion Criterion 4) routes a real, already divergence-free + velocity field through `PISO`'s own pressure solve, which measures + divergence through this class at every boundary face -- until this + addition, unconditionally `UnconfiguredBoundaryFaceError`, since a + periodic face is never given a `BoundaryCondition` object at all + (`assemble_numerics`'s own `periodic_pairs`/`boundary_conditions` + split). **Verified directly before being trusted, not assumed**: a + uniform, non-axis-aligned velocity field on a genuinely periodic mesh + measures exactly `0.0` divergence at every cell through this path + (float-exact, not merely small), while a non-uniform field on the same + mesh still measures a real nonzero divergence -- confirming the wrap + reads real neighbour values rather than silently zeroing every + boundary face's own contribution. At a face named in `periodic_pairs`, + `divergence` substitutes `mesh.wrapped_neighbour_cell` for `neighbour` + before falling through to the ordinary interior-face averaging; + `boundary_conditions` is never consulted for a periodic face. """ - def __init__(self, boundary_conditions: Mapping[str, BoundaryCondition]) -> None: + def __init__( + self, + boundary_conditions: Mapping[str, BoundaryCondition], + periodic_pairs: Mapping[str, str], + ) -> None: self._boundary_conditions = boundary_conditions + self._periodic_pairs = periodic_pairs def _check_field(self, field: CollocatedField[Any]) -> None: if field.component_shape != (_SPATIAL_DIMENSIONS,): @@ -102,6 +127,10 @@ def divergence(self, field: Field) -> torch.Tensor: owner, neighbour = mesh.face_neighbours(face) normal_x, normal_y = mesh.face_normal(face) owner_x, owner_y = field.value_at(owner) + if neighbour is None: + boundary_name = mesh.boundary_face_name(face) + if boundary_name in self._periodic_pairs: + neighbour = mesh.wrapped_neighbour_cell(face) if neighbour is not None: neighbour_x, neighbour_y = field.value_at(neighbour) value_x, value_y = (owner_x + neighbour_x) / 2, (owner_y + neighbour_y) / 2 diff --git a/src/pyflow/engine/numerics/gradient.py b/src/pyflow/engine/numerics/gradient.py index 7a8b233..46ec6fb 100644 --- a/src/pyflow/engine/numerics/gradient.py +++ b/src/pyflow/engine/numerics/gradient.py @@ -77,10 +77,28 @@ class GreenGaussGradient(GradientScheme): (`owner_value + gradient * distance`) -- exact for a linear field, and reduces to zero-order extrapolation for a zero-gradient condition (the impermeable-wall assumption `PISO` uses for pressure). + + **Periodic-aware the same way `CentralDifferenceDiffusion` is + (TASK-030), added by TASK-034 (Stage 5) once `PISO` needed it: a + fully periodic domain's own pressure Poisson solve reaches every + boundary face through this class, and until this addition it had no + periodic case at all -- `UnconfiguredBoundaryFaceError` unconditionally, + the same gap `divergence.py`'s own entry describes.** At a face named + in `periodic_pairs`, `gradient` substitutes `mesh.wrapped_neighbour_cell` + for `neighbour` before falling through to the ordinary interior-face + averaging -- no distance term to double here, unlike diffusion's own + central-difference formula, since Green-Gauss face averaging never + divides by distance; `boundary_conditions` is never consulted for a + periodic face. """ - def __init__(self, boundary_conditions: Mapping[str, BoundaryCondition]) -> None: + def __init__( + self, + boundary_conditions: Mapping[str, BoundaryCondition], + periodic_pairs: Mapping[str, str], + ) -> None: self._boundary_conditions = boundary_conditions + self._periodic_pairs = periodic_pairs def gradient(self, field: Field) -> torch.Tensor: assert isinstance(field, CollocatedField) @@ -94,6 +112,10 @@ def gradient(self, field: Field) -> torch.Tensor: owner, neighbour = mesh.face_neighbours(face) normal_x[face], normal_y[face] = mesh.face_normal(face) owner_value = float(field.value_at(owner)) + if neighbour is None: + boundary_name = mesh.boundary_face_name(face) + if boundary_name in self._periodic_pairs: + neighbour = mesh.wrapped_neighbour_cell(face) if neighbour is not None: neighbour_value = float(field.value_at(neighbour)) face_values[face] = (owner_value + neighbour_value) / 2 diff --git a/src/pyflow/engine/numerics/pressure_coupling.py b/src/pyflow/engine/numerics/pressure_coupling.py index 34741c7..875a1cb 100644 --- a/src/pyflow/engine/numerics/pressure_coupling.py +++ b/src/pyflow/engine/numerics/pressure_coupling.py @@ -191,6 +191,34 @@ class PISO(PressureCoupling): boundary-aware logic); `GreenGaussGradient` still computes both the cell-centred pressure gradient the velocity correction is built from and the per-cell gradients the Rhie-Chow correction term needs. + + **`periodic_pairs` (added TASK-034, Stage 5) is a new fifth, + defaulted constructor parameter -- `PISO`'s own pressure treatment + had no periodic case at all before this, unconditionally + `UnconfiguredBoundaryFaceError` for any periodic boundary face (see + `gradient.py`/`divergence.py`'s own entries).** Found while building + TASK-034's mandated "uniform flow on a fully periodic domain" null + test (Stage 5 Completion Criterion 4): that scenario cannot reach + `PISO` at all without this, since even *measuring* an already + divergence-free field's divergence goes through `GreenGaussDivergence`, + which raised for every periodic face regardless of the field's actual + values. Threaded to `_diffusion` (the Poisson matrix, which already + knew how to be periodic via `CentralDifferenceDiffusion`'s own + TASK-030 support -- only `PISO` was passing it a hardcoded `{}`), + `_gradient`, `_divergence`, and `_rhie_chow_divergence`'s own + per-face loop (a periodic face gets the same Rhie-Chow correction an + interior face does, via `mesh.wrapped_neighbour_cell` and the same + doubled-distance convention `CentralDifferenceDiffusion` already + established). Defaults to an empty mapping, so every existing call + site that only passes `(linear_solver, boundary_conditions)` keeps + working unchanged, the same courtesy `tolerance`/`max_iterations` + already extend. **Verified directly, not assumed**: a uniform, + non-axis-aligned velocity field on a genuinely periodic mesh measures + exactly `0.0` divergence through this path, so `correct`'s very first + iteration returns without ever calling the linear solver -- the + physically correct answer for a flow that is already steady, and the + reason this task's own periodic null test needed no further PISO + correctness work beyond making the measurement possible at all. """ def __init__( @@ -199,21 +227,25 @@ def __init__( boundary_conditions: Mapping[str, BoundaryCondition], tolerance: float = 1e-6, max_iterations: int = 50, + periodic_pairs: Mapping[str, str] = MappingProxyType({}), ) -> None: super().__init__(linear_solver) self._tolerance = tolerance self._max_iterations = max_iterations + self._periodic_pairs = periodic_pairs pressure_boundary_conditions = MappingProxyType( {name: _ZeroGradientPressureCondition() for name in _PRESSURE_BOUNDARY_FACE_NAMES} ) self._diffusion = CentralDifferenceDiffusion( pressure_boundary_conditions, - {}, # no periodic pressure boundaries in this task's own scope (TASK-030) + periodic_pairs, diffusion_coefficient=1.0, ) - self._gradient = GreenGaussGradient(pressure_boundary_conditions) - self._divergence = GreenGaussDivergence(boundary_conditions) + self._gradient = GreenGaussGradient(pressure_boundary_conditions, periodic_pairs) + self._divergence = GreenGaussDivergence(boundary_conditions, periodic_pairs) self.last_divergence_history: tuple[float, ...] = () + self._cached_poisson_mesh: StructuredCartesianMesh | None = None + self._cached_poisson_matrix: torch.Tensor | None = None def correct( self, provisional_velocity: VectorField, dt: float @@ -261,12 +293,36 @@ def correct( ) def _poisson_matrix(self, mesh: StructuredCartesianMesh) -> torch.Tensor: + """Built once per distinct `mesh` and cached for the rest of this + `PISO` instance's own lifetime (TASK-034, Stage 5), not rebuilt on + every `correct` call as before -- found while measuring the Lid- + Driven Cavity validation's own real runtime (Stage 5 Completion + Criterion 5's "the runtime this implies is part of the criterion, + not a surprise to discover in CI"): this construction is + `O(num_cells * num_faces)` (one full `self._diffusion.flux` call + per column), dominating measured per-timestep cost by roughly + 70-90% at MVP cavity mesh sizes (8x8 through 16x16, timed + directly with a disposable prototype before this change, not + assumed), even though the matrix depends only on `mesh` and this + instance's own fixed pressure boundary treatment -- never on the + current velocity, pressure, or `dt` -- so nothing about repeating + it across timesteps was ever buying correctness. Cached by mesh + *identity*, not equality: a real run always hands `correct` the + same mesh object every timestep, so the common case is a cache + hit; a genuinely different mesh object (unusual -- no code path + in this repository reuses one `PISO` instance across meshes today) + safely recomputes rather than serving a stale matrix. + """ + if self._cached_poisson_mesh is mesh and self._cached_poisson_matrix is not None: + return self._cached_poisson_matrix num_cells = mesh.num_cells matrix = torch.zeros((num_cells, num_cells), dtype=torch.float64) for column in range(num_cells): basis = ScalarField(mesh, "e") basis.values[column] = 1.0 matrix[:, column] = -accumulate_flux_to_cells(mesh, self._diffusion.flux(basis)) + self._cached_poisson_mesh = mesh + self._cached_poisson_matrix = matrix return matrix def _rhie_chow_divergence( @@ -296,10 +352,15 @@ def _rhie_chow_divergence( correction_face = torch.zeros(mesh.num_faces, dtype=torch.float64) for face in range(mesh.num_faces): owner, neighbour = mesh.face_neighbours(face) + distance = mesh.face_centroid_distance(face) if neighbour is None: - continue + boundary_name = mesh.boundary_face_name(face) + if boundary_name in self._periodic_pairs: + neighbour = mesh.wrapped_neighbour_cell(face) + distance = 2 * distance + else: + continue normal_x, normal_y = mesh.face_normal(face) - distance = mesh.face_centroid_distance(face) direct = (pressure.value_at(neighbour) - pressure.value_at(owner)) / distance gx_o, gy_o = float(pressure_gradient[owner, 0]), float(pressure_gradient[owner, 1]) gx_n, gy_n = ( diff --git a/src/pyflow/engine/simulation.py b/src/pyflow/engine/simulation.py index 3dc2a3a..3fbe307 100644 --- a/src/pyflow/engine/simulation.py +++ b/src/pyflow/engine/simulation.py @@ -15,13 +15,14 @@ from __future__ import annotations from collections.abc import Mapping +from dataclasses import dataclass from typing import TYPE_CHECKING import torch from pyflow.engine.field import Field -from pyflow.engine.mesh import Mesh -from pyflow.engine.scalar_field import PressureField +from pyflow.engine.mesh import Mesh, StructuredCartesianMesh +from pyflow.engine.scalar_field import PressureField, ScalarField from pyflow.engine.vector_field import VectorField if TYPE_CHECKING: @@ -159,3 +160,158 @@ def derivative(state: Mapping[str, Field]) -> dict[str, torch.Tensor]: return result return numerics.time_integration.advance(fields, derivative, dt) + + +@dataclass(frozen=True) +class NavierStokesStepResult: + """The three parts of one `navier_stokes_step` call, each on its own + field rather than folded into a single returned mapping -- Stage 5 + Completion Criterion 4's own "each part observable, not only the end + state" (`docs/planning/roadmap.md` TASK-034). + + `fields` is the fully advanced state: every entry `step` advanced, + with momentum's own two components replaced by their corrected + values. `provisional_velocity` is the predictor's own output, before + correction -- divergent in general. `corrected_velocity`/`pressure` + are exactly `numerics.pressure_coupling.correct`'s own return values. + """ + + fields: dict[str, Field] + provisional_velocity: VectorField + corrected_velocity: VectorField + pressure: ScalarField + + +def navier_stokes_step( + fields: Mapping[str, Field], + velocity_field_name: str, + numerics: AssembledNumerics, + dt: float, +) -> NavierStokesStepResult: + """One incompressible Navier-Stokes timestep (TASK-034, Stage 5): + predictor, corrector, corrected state -- the fractional-step sequence + `docs/handbook/numerical-methods/pressure-velocity-coupling.md` + describes. **The projection sits once per timestep, outside + `TimeIntegrator` entirely**, not inside any of RK4's own four stage + evaluations -- Stage 5's own design question five, resolved by the + maintainer in favour of the classical arrangement: the momentum + predictor is fully explicit, with no pressure term at all, and the + corrector loop (`PressureCoupling.correct`) projects the *result* of + that predictor once, not each of RK4's own intermediate states. + + `velocity_field_name` names which two entries of `fields` + (`VectorField.component_name(velocity_field_name, 0)`/`(..., 1)`) are + momentum's own components -- an explicit parameter, not a fixed + convention baked into this module, so this function keeps Stage 5 + Completion Criterion 1's own structural guarantee intact: no + hardcoded component-name pair anywhere in this file's own source + (`tests/features/velocity_field_support.feature`'s own check, + `inspect.getsource(simulation)`, still passes unchanged with this + function added). + + **Predictor:** every entry of `fields` -- momentum's own two + components and any other transported field alongside them (a scalar + transported by the same velocity, say) -- advances through the + ordinary `step` path above, self-advected by the *current* (not yet + corrected) velocity, with no pressure term. **Corrector:** the two + advanced momentum components are reassembled into a provisional + `VectorField` and handed to `numerics.pressure_coupling.correct`, + which projects it onto a divergence-free field and reports the + pressure consistent with that projection -- reached only through the + configured `PressureCoupling`/`LinearSolver`, never a hardcoded + concrete class (Stage 5 Completion Criterion 13's own substitution + check). **Corrected state:** the provisional momentum components + inside the predictor's own result are overwritten with the corrected + ones; every other field is left exactly as the predictor advanced it, + since nothing pressure-corrects a scalar. + + Raises whatever `step`/`numerics.pressure_coupling.correct` raise -- + `MismatchedMeshError`/`PressureFieldTransportError` from the former, + `PressureSolveDidNotConvergeError`/`DivergenceDidNotConvergeError` + from the latter (`pressure_coupling.py`) -- rather than catching and + re-wrapping either. + """ + u_name = VectorField.component_name(velocity_field_name, 0) + v_name = VectorField.component_name(velocity_field_name, 1) + u_field = fields[u_name] + v_field = fields[v_name] + assert isinstance(u_field, ScalarField) + assert isinstance(v_field, ScalarField) + current_velocity = VectorField.assemble([u_field, v_field], velocity_field_name) + + predicted = step(fields, current_velocity, numerics, dt) + + predicted_u = predicted[u_name] + predicted_v = predicted[v_name] + assert isinstance(predicted_u, ScalarField) + assert isinstance(predicted_v, ScalarField) + provisional_velocity = VectorField.assemble([predicted_u, predicted_v], velocity_field_name) + + corrected_velocity, pressure = numerics.pressure_coupling.correct(provisional_velocity, dt) + + new_fields = dict(predicted) + for component in corrected_velocity.decompose(): + new_fields[component.name] = component + + return NavierStokesStepResult( + fields=new_fields, + provisional_velocity=provisional_velocity, + corrected_velocity=corrected_velocity, + pressure=pressure, + ) + + +_STABILITY_SAFETY_FACTOR = 0.25 +"""Multiplies the tighter of the CFL/diffusive stability limits below to +get a timestep this project's own explicit RK4 predictor -- first-order +upwind advection plus central-difference diffusion, `navier_stokes_step`'s +own scheme combination -- actually stays stable at. **Measured directly +with a disposable prototype script before being trusted, not derived +analytically** (`docs/planning/roadmap.md` TASK-034's own design +question four, "the derivation is stated in the scenario either way"): +swept across a mixed advection/diffusion regime, a diffusion-dominated +regime, and an advection-dominated regime alike, `0.3` was the largest +factor that stayed stable for 500 steps in every regime tried and `0.35` +already blew up in the mixed one; `0.25` keeps real margin below that +measured edge rather than shipping a value discovered exactly at the +boundary of blowing up. +""" + + +def stable_timestep( + mesh: StructuredCartesianMesh, + viscosity: float, + velocity_scale: float, + safety_factor: float = _STABILITY_SAFETY_FACTOR, +) -> float: + """A timestep within this scheme combination's own explicit stability + limit on `mesh`, for a flow whose characteristic speed is + `velocity_scale` and whose viscosity is `viscosity` (TASK-034, Stage + 5's own design question four -- "the timestep becomes derivable from + the mesh and the configured viscosity", chosen over a separately + hand-tuned timestep per mesh resolution, which `docs/planning/ + roadmap.md` TASK-034 itself names as "a convergence study measuring + the tuning" rather than a real one). + + Explicit RK4 is stability-limited two ways at once, and the tighter + one governs: the CFL condition (`dx / velocity_scale`, advection) and + the diffusive limit (`dx**2 / viscosity`, central-difference + diffusion) -- both scale differently under mesh refinement (the first + linearly in `dx`, the second quadratically), so a single fixed + timestep that is stable at a coarse resolution is guaranteed unstable + at a finer one. `dx` is this mesh's own smallest cell spacing (`min` + of its two axes, so a non-square mesh is governed by its tighter + axis); `velocity_scale <= 0` (the periodic null test's own zero- + forcing, non-advecting fixtures never call this at all, but a caller + that does pass one gets a well-defined answer rather than a division + by zero) is treated as "no advective limit", leaving the diffusive + one alone to govern. + """ + north_face = next(f for f in range(mesh.num_faces) if mesh.boundary_face_name(f) == "north") + west_face = next(f for f in range(mesh.num_faces) if mesh.boundary_face_name(f) == "west") + dx = min(float(mesh.face_area(north_face)), float(mesh.face_area(west_face))) + diffusive_limit = dx**2 / viscosity + if velocity_scale <= 0: + return safety_factor * diffusive_limit + convective_limit = dx / velocity_scale + return safety_factor * min(convective_limit, diffusive_limit) diff --git a/tests/CLAUDE.md b/tests/CLAUDE.md index f168103..3e68fdb 100644 --- a/tests/CLAUDE.md +++ b/tests/CLAUDE.md @@ -14,6 +14,11 @@ Four kinds of test, split by what they exercise: here still without a real test; its conventions get written when the first benchmark sets a precedent, not ahead of it (`docs/planning/backlog.md` E9). +- `fixtures/` -- committed reference data external to this repository + (a published paper's own tabulated numbers), not machinery this + project derives itself. New as of TASK-034 (Stage 5, 2026-08-29); see + its own `CLAUDE.md` for the distinction from `unit/_numerics.py`/ + `golden/_demo.py`, which stay local machinery. This split was undocumented until 2026-08-15, when the first real test (`integration/test_cli.py`) gave it a concrete precedent to write down diff --git a/tests/features/heat_diffusion.feature b/tests/features/heat_diffusion.feature new file mode 100644 index 0000000..12bcd8e --- /dev/null +++ b/tests/features/heat_diffusion.feature @@ -0,0 +1,28 @@ +# The acceptance criteria for the Heat Diffusion golden demo (TASK-034, +# Stage 5 Completion Criterion 8, and `docs/implementation/mvp.md`'s own +# Validation section, reconciled 2026-08-28 -- heat diffusion is the +# diffusion equation on a transported scalar, no named Temperature field +# needed). `tests/golden/test_heat_diffusion.py` binds these scenarios. + +Feature: Heat Diffusion + A single sinusoidal mode, one full wavelength across the mesh's own + x-extent, on an otherwise unforced, unadvected, fully periodic domain + -- the one initial condition PyFlow's diffusion equation has a + closed-form answer for: exponential decay at a rate set only by the + diffusion coefficient and the mode's own wavenumber. + + Background: + Given the golden demo "heat_diffusion" + + Scenario: A user can run it with the documented command + When it is run through the public CLI, headless + Then the command exits cleanly + + # "Numerical solution is measurable" (`docs/implementation/mvp.md`'s + # own Definition of Done), checked against an exact analytic answer + # rather than only "heat visibly spread" -- distinct from Stage 4's + # own diffusion criteria, which measured spatial convergence order and + # conservation, neither of which is a decay *rate*. + Scenario: The mode's own amplitude decays at the exact analytic rate set by the diffusion coefficient and its wavenumber + When it is bootstrapped once after a few real timesteps and again after many more + Then the measured decay rate matches the analytic rate for the configured diffusion coefficient and wavenumber diff --git a/tests/features/lid_driven_cavity.feature b/tests/features/lid_driven_cavity.feature new file mode 100644 index 0000000..8084955 --- /dev/null +++ b/tests/features/lid_driven_cavity.feature @@ -0,0 +1,30 @@ +# The acceptance criteria for the Lid-Driven Cavity golden demo +# (TASK-034, Stage 5 Completion Criterion 8 -- the MVP's own golden +# demo). `tests/golden/test_lid_driven_cavity.py` binds these scenarios. +# The quantitative comparison against Ghia, Ghia & Shin (1982) is +# `tests/features/navier_stokes_timestep.feature`'s own scenario, run +# directly against the engine, not repeated here -- this file is the +# reproducible, visible demonstration the public-API rule requires. + +Feature: Lid-Driven Cavity + A square cavity, no-slip on every wall, the top wall moving + tangentially at a constant speed -- rendered live as a real, + incompressible `navier_stokes_step` solves it, one timestep per frame. + The first velocity field PyFlow has ever rendered that was *solved*, + not prescribed or seeded (`docs/implementation/mvp.md`'s own + "visualisation shows the result"). + + Background: + Given the golden demo "lid_driven_cavity" + + Scenario: A user can run it with the documented command + When it is run through the public CLI, headless + Then the command exits cleanly + + Scenario: The rendered velocity field is genuinely solved, not the zero it started from + When it is bootstrapped for a few real timesteps + Then the velocity field has real nonzero motion away from the lid + + Scenario: The same configuration run twice produces identical state + When it is bootstrapped for a few real timesteps twice + Then both runs produce identical velocity fields diff --git a/tests/features/navier_stokes_timestep.feature b/tests/features/navier_stokes_timestep.feature new file mode 100644 index 0000000..492c403 --- /dev/null +++ b/tests/features/navier_stokes_timestep.feature @@ -0,0 +1,94 @@ +# The acceptance criteria for Navier-Stokes Timestep (TASK-034, Stage +# 5's fifth and last task in build order). Assembles TASK-031/032/033 +# into one incompressible Navier-Stokes timestep, then validates it. +# Not a golden demo -- no config file under `examples/golden-demos/`, no +# CLI subprocess run, since every claim here is checked against the +# engine mechanism directly, the same `tests/unit/` shape every prior +# Stage 4/5 numerical-scheme feature file already established. Each of +# this stage's two golden demos (Lid-Driven Cavity, Heat Diffusion) has +# its own separate feature file, per `docs/planning/roadmap.md` +# TASK-034's own Artifacts Produced bullet. `tests/unit/ +# test_navier_stokes_timestep.py` binds these scenarios. + +Feature: Navier-Stokes Timestep + + # -- Criterion 4: predictor, corrector, corrected state, each observable + + Scenario: One timestep produces an observable predictor, corrector, and corrected state + Given a closed, no-slip domain with a divergent initial velocity field + When one Navier-Stokes timestep is taken + Then the provisional velocity, the corrected velocity, and the pressure field are all present + And the corrected velocity differs from the provisional velocity + And the corrected velocity's own divergence is smaller than the provisional velocity's + + # -- Criterion 4: the two null tests neither of the others can substitute for + + Scenario: Uniform flow on a fully periodic, inviscid, unforced domain stays divergence-free and unchanged over many steps + Given a fully periodic domain with a uniform, non-axis-aligned velocity field and zero viscosity + When many Navier-Stokes timesteps are taken + Then the velocity field is exactly the same uniform value at every step + And the velocity field's own divergence never leaves solver tolerance at any step + + Scenario: Fluid at rest in a closed, no-slip domain stays at rest over many steps + Given a closed, no-slip domain with the fluid initially at rest + When many Navier-Stokes timesteps are taken + Then the velocity field stays at rest to floating-point tolerance at every step + + # -- Criterion 4: determinism + + Scenario: The same configuration run twice produces identical state + Given a closed, no-slip domain with a divergent initial velocity field + When one Navier-Stokes timestep is taken twice from the same initial state + Then both runs produce identical corrected velocity and pressure fields + + # -- Criterion 13: the solver runs through ADR-003's seams, not around them + + Scenario: The timestep calls the configured PressureCoupling strategy, not a hardcoded one + Given a PressureCoupling test double registered under its own name and selected by configuration + When one Navier-Stokes timestep is taken + Then the test double's own distinctive pressure value appears in the result, not a real solve's + + # -- Criterion 5: Couette flow, at solver tolerance rather than a loose one + + Scenario: Couette flow reaches the exact linear steady velocity profile + Given a channel periodic in the flow direction, no-slip walls, one stationary and one moving tangentially + When Navier-Stokes timesteps are taken until the flow reaches steady state + Then the steady streamwise velocity profile matches the exact linear Couette solution at solver tolerance + And the wall-normal velocity component stays zero everywhere + + # -- Criterion 5: Lid-driven cavity against Ghia, Ghia & Shin (1982) -- + # the criterion is convergence across resolutions, not a fixed + # percentage at one; this is this project's most computationally + # expensive scenario (three real runs to a measured steady state), + # deliberately, per this task's own Design decision. + + Scenario: The Ghia comparison error decreases monotonically across three mesh resolutions, and the finest shows the right vortex structure + Given three lid-driven cavity meshes at increasing resolution, at Reynolds number 100 + When each is run to a measured steady state + Then the error against Ghia's centreline profiles decreases monotonically across the three resolutions + And the finest resolution's primary vortex centre is within a stated distance of Ghia's own + And the finest resolution shows both downstream secondary corner vortices, rotating opposite the primary + + # -- Criterion 5: the emergent-phenomenon pair, and its negative control. + # Taylor-Green vortex decay, chosen over Kelvin-Helmholtz by measurement + # (`docs/planning/roadmap.md` TASK-034's own Design decision): it + # reuses this task's own periodic-domain infrastructure directly and + # has a closed-form decay rate to measure against, rather than needing + # a roll-up detector. + + Scenario: Taylor-Green vortex decay matches the exact rate when physical viscosity dominates numerical diffusion + Given a Taylor-Green vortex on a periodic domain at a viscosity where physical diffusion dominates + When the vortex is advanced and its own decay rate is measured + Then the measured decay rate matches the exact closed-form rate closely + + Scenario: Taylor-Green vortex decay does not match the exact rate when the advection scheme's own numerical diffusion dominates + Given a Taylor-Green vortex on a periodic domain at a viscosity where numerical diffusion dominates + When the vortex is advanced and its own decay rate is measured + Then the measured decay rate does not match the exact closed-form rate + + # -- Criterion 5: conservation, a claim none of the scenarios above make + + Scenario: No single step increases total kinetic energy for an inviscid, unforced, closed-domain flow + Given a closed, no-slip domain with a divergent initial velocity field and negligible viscosity + When many Navier-Stokes timesteps are taken + Then total kinetic energy never increases from one step to the next diff --git a/tests/fixtures/CLAUDE.md b/tests/fixtures/CLAUDE.md new file mode 100644 index 0000000..98ad9c2 --- /dev/null +++ b/tests/fixtures/CLAUDE.md @@ -0,0 +1,37 @@ +# CLAUDE + +Committed reference data external to this repository -- published +results, not code. **New as of TASK-034 (Stage 5, 2026-08-29)**: this +repository had no test-data directory before this task needed one, per +`docs/planning/roadmap.md` Stage 5 Completion Criterion 5's own +instruction ("wherever it lands is a new convention... updated in the +same change, per the Blast Radius rule"). + +Distinct from `tests/unit/_numerics.py`/`tests/golden/_demo.py`: those +are machinery a step definition is *built from* (fixture constants, test +doubles, geometry helpers) and live next to the tests that use them, +each with its own module-level "local by default" scope. What lives here +instead is data with an external citation -- a published paper's own +tabulated numbers, not anything this project derived -- and is imported +by whichever test module needs it (`from fixtures. import ...`, +reachable because `tests/` itself has no `__init__.py` and pytest's own +rootdir insertion puts it on `sys.path`), not copied into each one. + +**`ghia_1982_re100.py`** is the first occupant: U. Ghia, K. N. Ghia and +C. T. Shin's own Table I (Re = 100), the Lid-Driven Cavity validation's +reference data (`tests/unit/test_navier_stokes_timestep.py`, `docs/ +planning/roadmap.md` TASK-034). Every number carries the paper, table, +and column it came from, per Stage 5 Completion Criterion 5's own "the +reference values are committed data with their citation attached, not +literals typed into an assertion" -- see the module's own docstring for +exactly what was cross-checked against what, and the honest limit +recorded there: this environment has no direct access to the original +1982 print journal, so the tables were cross-checked against two +independent public reproductions of the paper's own Table I rather than +transcribed from it directly. + +**A file here should be added only when a real committed reference +value needs a home, the same "real content first" discipline +`tools/CLAUDE.md` states for `generators/`/`validators/`** -- not +speculatively, and not for anything this project derives itself (that +stays local machinery, in `_numerics.py`/`_demo.py`). diff --git a/tests/fixtures/__init__.py b/tests/fixtures/__init__.py new file mode 100644 index 0000000..bb9861c --- /dev/null +++ b/tests/fixtures/__init__.py @@ -0,0 +1,10 @@ +"""Makes `fixtures` a package rather than a bare directory. + +Needed so `mypy` can resolve `from fixtures.ghia_1982_re100 import ...` +(`tests/unit/test_navier_stokes_timestep.py`) without ambiguity -- +without this, `mypy` finds `ghia_1982_re100.py` under two different +module names (`ghia_1982_re100` and `fixtures.ghia_1982_re100`) and +refuses to proceed, the same "no `__init__.py`" collision +`tests/unit/__init__.py`'s own docstring records for pytest/mypy module +identification generally. +""" diff --git a/tests/fixtures/ghia_1982_re100.py b/tests/fixtures/ghia_1982_re100.py new file mode 100644 index 0000000..2d178af --- /dev/null +++ b/tests/fixtures/ghia_1982_re100.py @@ -0,0 +1,94 @@ +"""Reference data for the Lid-Driven Cavity validation (TASK-034, Stage +5 Completion Criterion 5): U. Ghia, K. N. Ghia and C. T. Shin, "High-Re +Solutions for Incompressible Flow Using the Navier-Stokes Equations and +a Multigrid Method", Journal of Computational Physics 48, 387-411 +(1982), Table I, Reynolds number 100. + +`docs/planning/roadmap.md` TASK-034's own Criterion 5 bullet requires +the reference values to be "committed data with their citation attached, +not literals typed into an assertion" and "read off the paper itself, +not from memory or a secondary source". **The honest limit on that +second half, stated here rather than left implicit**: this environment +has no direct access to the original 1982 print journal. The two tables +below were cross-checked against two independent public reproductions of +the paper's own Table I (not each other) before being committed -- + (u +along the vertical centreline) and + (v +along the horizontal centreline), both explicitly labelled as +transcriptions of Ghia, Ghia & Shin (1982)'s own published table -- and +agree with each other to every digit. `PRIMARY_VORTEX_CENTER` is the +figure this paper's own results are most commonly cited by (cross- +checked the same way, via web search rather than a single source) and +matches this module's own docstring note below on what that figure +means. Nothing here is typed from memory alone. + +**`PRIMARY_VORTEX_CENTER` is the widely-quoted approximate figure for +this paper's own Re = 100 result, used as a sanity check on this +fixture's own scale (`docs/planning/roadmap.md` TASK-034's own +Criterion 5 bullet says so explicitly), not a literal digit-for-digit +transcription of a table cell the way the two velocity profiles above +are** -- the original paper reports it via a streamfunction/vorticity +contour plot and an accompanying table of local extrema, not as a single +coordinate pair printed in text. + +**Secondary corner vortices are checked for presence, not against +Ghia's own tabulated coordinates** -- `docs/planning/roadmap.md` +TASK-034's own Criterion 5 bullet asks only that "both downstream +secondary corner vortices [are] present" at the finest resolution, which +this fixture does not need a number for: the validating scenario detects +a genuine local recirculation near each bottom corner directly from +PyFlow's own computed velocity field. +""" + +from __future__ import annotations + +# u-velocity along the vertical centreline (x = 0.5 in unit-cavity +# coordinates), y from the moving lid (y=1) down to the stationary floor +# (y=0). 17 points, Table I. +U_VELOCITY_ALONG_VERTICAL_CENTERLINE: tuple[tuple[float, float], ...] = ( + (1.0000, 1.00000), + (0.9766, 0.84123), + (0.9688, 0.78871), + (0.9609, 0.73722), + (0.9531, 0.68717), + (0.8516, 0.23151), + (0.7344, 0.00332), + (0.6172, -0.13641), + (0.5000, -0.20581), + (0.4531, -0.21090), + (0.2813, -0.15662), + (0.1719, -0.10150), + (0.1016, -0.06434), + (0.0703, -0.04775), + (0.0625, -0.04192), + (0.0547, -0.03717), + (0.0000, 0.00000), +) + +# v-velocity along the horizontal centreline (y = 0.5 in unit-cavity +# coordinates), x from the right wall (x=1) to the left wall (x=0). 17 +# points, Table I. +V_VELOCITY_ALONG_HORIZONTAL_CENTERLINE: tuple[tuple[float, float], ...] = ( + (1.0000, 0.00000), + (0.9688, -0.05906), + (0.9609, -0.07391), + (0.9531, -0.08864), + (0.9453, -0.10313), + (0.9063, -0.16914), + (0.8594, -0.22445), + (0.8047, -0.24533), + (0.5000, 0.05454), + (0.2344, 0.17527), + (0.2266, 0.17507), + (0.1563, 0.16077), + (0.0938, 0.12317), + (0.0781, 0.10890), + (0.0703, 0.10091), + (0.0625, 0.09233), + (0.0000, 0.00000), +) + +# (x, y) in unit-cavity coordinates -- see this module's own docstring +# for what this figure is and is not. +PRIMARY_VORTEX_CENTER: tuple[float, float] = (0.6172, 0.7344) diff --git a/tests/golden/CLAUDE.md b/tests/golden/CLAUDE.md index 2fbd9ce..c0ad253 100644 --- a/tests/golden/CLAUDE.md +++ b/tests/golden/CLAUDE.md @@ -81,3 +81,36 @@ from a real run, not guessed** -- a real run agrees with the closed-form prediction to within ~4%; confirmed to actually fail under a mutation that froze the simulation state every frame (never calling `simulation.step`) before being trusted. + +**`test_heat_diffusion.py` (TASK-034, added 2026-08-29) is the sixth +demo module** -- Stage 5's own reconciliation of `mvp.md`'s Validation +section: heat diffusion as the diffusion equation on a transported +scalar, no named Temperature field needed. Same shape as +`test_passive_scalar_transport.py`'s own join: the required +CLI-subprocess scenario, plus one demo-specific step (bootstraps twice, +at two frame counts, measures the transported field's own RMS amplitude +at each) proving a genuine physical claim -- a single sinusoidal mode's +own amplitude decays at the exact rate `Gamma * wavenumber**2` predicts, +not only that the field changed. **The tolerance (`rel=0.1`) was +measured from a real run** -- a real run agrees with the closed-form +rate to within ~0.6%. + +**`test_lid_driven_cavity.py` (TASK-034, added 2026-08-29) is the +seventh demo module** -- the MVP's own golden demo +(`docs/implementation/mvp.md`, `docs/implementation/golden-demos.md`'s +"Lid-Driven Cavity" section), and the first velocity field PyFlow has +ever rendered that was *solved*, not prescribed or seeded. Reads +`RenderWindow.simulation_fields` back the same way +`test_passive_scalar_transport.py` does, reassembling velocity's own two +decomposed components via `VectorField.assemble` since this demo has no +scalar to read at all. **Deliberately does not assert an absolute +divergence bound** -- tried first, and found to fail even at cells well +away from the lid's own two corner singularities, because +`GreenGaussDivergence`'s own naive face-averaged divergence is not the +Rhie-Chow-consistent measure `PISO`'s own corrector loop actually drives +to tolerance (`src/pyflow/engine/numerics/pressure_coupling.py`'s own +`_rhie_chow_divergence` docstring). The real, tolerance-gated divergence +claim is `tests/features/navier_stokes_timestep.feature`'s own scenarios, +measured the way `PISO` itself measures; this module's own two physical +checks are lighter (genuine nonzero motion away from the lid, +determinism) -- see this module's own docstring for the full finding. diff --git a/tests/golden/test_heat_diffusion.py b/tests/golden/test_heat_diffusion.py new file mode 100644 index 0000000..02e0e31 --- /dev/null +++ b/tests/golden/test_heat_diffusion.py @@ -0,0 +1,102 @@ +"""Heat Diffusion golden demo (TASK-034, Stage 5). + +The acceptance criteria are `tests/features/heat_diffusion.feature` +(`adr/ADR-007-executable-acceptance-criteria.md`). This module binds +them and supplies the one step only this demo needs -- measuring the +transported field's own decay rate between two frame counts, the same +"two independent bootstraps, not one run read back twice" shape +`test_passive_scalar_transport.py` already established, reused here +rather than re-derived. +""" + +from __future__ import annotations + +import math + +import pytest +from pytest_bdd import scenarios, then, when + +from pyflow.bootstrap import bootstrap +from pyflow.engine.scalar_field import ScalarField + +from ._demo import DemoRun + +scenarios("heat_diffusion.feature") + +_EARLY_FRAMES = 1 +_LATE_FRAMES = 201 +"""200 real RK4 timesteps apart -- at this demo's own configured +`numerics.timestep` (0.01), that is 2.0 time units of decay, enough for +the mode's own amplitude to drop by roughly 35% at the configured +diffusion coefficient -- measured directly (`docs/planning/roadmap.md` +TASK-034's own Design decision) before choosing this frame count, not +guessed: too few steps leaves too little decay to measure the rate from +accurately, too many costs test runtime for no added precision. +""" + + +def _rms_amplitude(field: ScalarField) -> float: + """The field's own root-mean-square value over every cell -- a + single number describing "how much mode is left" that does not + depend on knowing which cell the mode's own peak currently sits in, + unlike sampling one specific cell. + """ + mesh = field.mesh + total = 0.0 + for cell in range(mesh.num_cells): + total += float(field.value_at(cell)) ** 2 + return math.sqrt(total / mesh.num_cells) + + +@when( + "it is bootstrapped once after a few real timesteps and again after many more", + target_fixture="amplitudes", +) +def _when_bootstrapped_at_two_frame_counts(demo: DemoRun) -> tuple[float, float]: + early_window = bootstrap(demo.config_path, backend="offscreen", max_frames=_EARLY_FRAMES) + assert early_window.simulation_fields is not None + early_tracer = early_window.simulation_fields["tracer"] + assert isinstance(early_tracer, ScalarField) + + late_window = bootstrap(demo.config_path, backend="offscreen", max_frames=_LATE_FRAMES) + assert late_window.simulation_fields is not None + late_tracer = late_window.simulation_fields["tracer"] + assert isinstance(late_tracer, ScalarField) + + return (_rms_amplitude(early_tracer), _rms_amplitude(late_tracer)) + + +@then( + "the measured decay rate matches the analytic rate for the configured diffusion " + "coefficient and wavenumber" +) +def _then_decay_rate_matches_analytic(demo: DemoRun, amplitudes: tuple[float, float]) -> None: + early_amplitude, late_amplitude = amplitudes + dt = demo.config.numerics.timestep + elapsed_steps = _LATE_FRAMES - _EARLY_FRAMES + elapsed_time = dt * elapsed_steps + measured_rate = -math.log(late_amplitude / early_amplitude) / elapsed_time + + min_x, _min_y, max_x, _max_y = _mesh_x_extent(demo) + domain_width = max_x - min_x + wavenumber = 2 * math.pi / domain_width + analytic_rate = demo.config.fluid.diffusion_coefficient * wavenumber**2 + + # Measured directly before choosing this bound, not guessed + # (`docs/planning/roadmap.md` TASK-034's own Design decision): a real + # run agrees with the closed-form rate to within ~0.6% at this + # demo's own resolution/timestep -- rel=0.1 stays comfortably above + # that margin without being so loose a genuinely broken diffusion + # scheme (wrong coefficient, wrong sign, no diffusion at all) could + # still pass. + assert measured_rate == pytest.approx(analytic_rate, rel=0.1), ( + f"expected the mode to decay at roughly {analytic_rate} per unit time, " + f"measured {measured_rate} instead" + ) + + +def _mesh_x_extent(demo: DemoRun) -> tuple[float, float, float, float]: + nx, ny = demo.config.mesh.extent + dx, dy = demo.config.mesh.spacing + origin_x, origin_y = demo.config.mesh.origin + return (origin_x, origin_y, origin_x + nx * dx, origin_y + ny * dy) diff --git a/tests/golden/test_lid_driven_cavity.py b/tests/golden/test_lid_driven_cavity.py new file mode 100644 index 0000000..e1f6ac3 --- /dev/null +++ b/tests/golden/test_lid_driven_cavity.py @@ -0,0 +1,89 @@ +"""Lid-Driven Cavity golden demo (TASK-034, Stage 5) -- the MVP's own +golden demo. + +The acceptance criteria are `tests/features/lid_driven_cavity.feature` +(`adr/ADR-007-executable-acceptance-criteria.md`). This module binds +them and supplies the steps only this demo needs: reading a *solved* +velocity field back from `RenderWindow.simulation_fields` +(`_add_solved_velocity_rendering`, `src/pyflow/bootstrap.py`), which +`conftest.py`'s shared vocabulary has no step for since no earlier demo +rendered a solved vector field at all. + +**Deliberately does not assert an absolute divergence bound at this +demo's own coarse (16x16), early-frame (10 steps) resolution.** Tried +first, and found to fail even well away from the lid's own two corner +singularities (`u` jumps from 1 to 0 discontinuously where the moving +lid meets a stationary wall, a genuine, well-documented property of this +exact benchmark, not an artefact): the reason is that `GreenGaussDivergence`'s +own naive face-averaged divergence is, by construction, not the +Rhie-Chow-consistent measure `PISO`'s own corrector loop actually drives +to tolerance (`src/pyflow/engine/numerics/pressure_coupling.py`'s own +`_rhie_chow_divergence` docstring: "precisely the measure that does +*not* see this"). The real, tolerance-gated divergence claim is +`tests/features/navier_stokes_timestep.feature`'s own predictor/ +corrector scenario and the Ghia comparison, both measured the way `PISO` +itself measures; this module stays a lighter reproducibility smoke test. +""" + +from __future__ import annotations + +from dataclasses import dataclass + +import torch +from pytest_bdd import scenarios, then, when + +from pyflow.bootstrap import bootstrap +from pyflow.engine.numerics.boundary_condition import BoundaryCondition +from pyflow.engine.scalar_field import ScalarField +from pyflow.engine.vector_field import VectorField + +from ._demo import DemoRun + +scenarios("lid_driven_cavity.feature") + +_FEW_FRAMES = 10 + + +@dataclass +class _SolvedRun: + velocity: VectorField + boundary_conditions: dict[str, BoundaryCondition] + + +def _run(demo: DemoRun) -> _SolvedRun: + window = bootstrap(demo.config_path, backend="offscreen", max_frames=_FEW_FRAMES) + assert window.simulation_fields is not None + assert window.assembled_numerics is not None + u = window.simulation_fields[VectorField.component_name("velocity", 0)] + v = window.simulation_fields[VectorField.component_name("velocity", 1)] + assert isinstance(u, ScalarField) + assert isinstance(v, ScalarField) + velocity = VectorField.assemble([u, v], "velocity") + return _SolvedRun( + velocity=velocity, boundary_conditions=dict(window.assembled_numerics.boundary_conditions) + ) + + +@when("it is bootstrapped for a few real timesteps", target_fixture="run") +def _when_bootstrapped(demo: DemoRun) -> _SolvedRun: + return _run(demo) + + +@when("it is bootstrapped for a few real timesteps twice", target_fixture="runs") +def _when_bootstrapped_twice(demo: DemoRun) -> tuple[_SolvedRun, _SolvedRun]: + return (_run(demo), _run(demo)) + + +@then("the velocity field has real nonzero motion away from the lid") +def _then_nonzero_motion(run: _SolvedRun) -> None: + # The lid itself is a boundary, not a cell -- interior cells starting + # at rest must have picked up real motion from the moving wall by + # now, or this demo is rendering the initial condition and nothing + # else. + assert float(run.velocity.magnitude().max()) > 1e-6 + + +@then("both runs produce identical velocity fields") +def _then_identical(runs: tuple[_SolvedRun, _SolvedRun]) -> None: + first, second = runs + torch.testing.assert_close(first.velocity.values, second.velocity.values, rtol=0, atol=0) diff --git a/tests/unit/CLAUDE.md b/tests/unit/CLAUDE.md index 1865daf..669d941 100644 --- a/tests/unit/CLAUDE.md +++ b/tests/unit/CLAUDE.md @@ -282,6 +282,50 @@ distinct from `test_piso_pressure_coupling.py`'s own `_NeverConvergesSolver` (`converged=False`), since this scenario is about the *outer* loop giving up, not the *inner* solve failing. +**`test_navier_stokes_timestep.py` (TASK-034, added 2026-08-29) is the +thirteenth, and Stage 5's fourth and last module in this lineage** -- +Stage 5's fifth task, binding `tests/features/ +navier_stokes_timestep.feature`'s eleven scenarios: `simulation. +navier_stokes_step`'s own predictor/corrector/corrected-state sequence, +both null tests, determinism, the ADR-003 substitution check, Couette +flow, the Ghia cavity comparison, the Taylor-Green emergent-phenomenon +pair, and kinetic-energy conservation. Same shape as every module before +it: its own `_Context` dataclass, its own local doubles +(`_MarkerPressureCoupling`, the substitution check's own test double), +no golden-demo config file or CLI run for this file -- the two golden +demos this task also builds (Lid-Driven Cavity, Heat Diffusion) are +bound separately, in `tests/golden/`, per that directory's own +convention. **Imports `tests/fixtures/ghia_1982_re100.py`** (this +repository's first use of the new top-level `tests/fixtures/` +convention, see that directory's own `CLAUDE.md`) for the Ghia, Ghia & +Shin (1982) reference data, rather than the local-doubles pattern every +earlier module in this list uses for its own fixture data -- committed, +cited reference data is not a test double, and belongs where any other +module needing the same table could import it too. + +**The Ghia cavity scenario is this project's most computationally +expensive test, deliberately** -- three real runs (resolutions 9, 13, +17) to a measured steady state, not a fixed step count. Chosen odd so +the vertical/horizontal centreline always lands exactly on a column/row +of cell centres, no interpolation needed against Ghia's own tabulated +points. **A real, measured performance fix was needed to keep this +runtime tractable at all**: `PISO._poisson_matrix` used to rebuild an +`O(num_cells * num_faces)` matrix every single timestep even though +nothing about it changes between timesteps on a fixed mesh -- caching it +per `PISO` instance (`src/pyflow/engine/CLAUDE.md`'s own `PISO` entry) +cut measured per-timestep cost by roughly 3.5-7.7x at MVP cavity mesh +sizes, found and applied before this scenario was written, not +discovered afterward as a slow-test complaint. Vortex detection (the +finest resolution's own primary-vortex-centre and secondary-corner- +vortex checks) is built directly from the computed velocity field -- +minimum velocity magnitude in a central sub-region for the primary +vortex, opposite-sign discrete vorticity in each bottom corner's own +sub-region for the secondary ones -- both thresholds measured against a +real (disposable-prototype) converged run before being written into the +test, not guessed: the detected primary vortex landed 0.019 away from +Ghia's own reference point, and clear opposite-sign vorticity (magnitude +0.13-0.43) was found throughout each bottom corner. + **The convention is "local by default, shared where genuinely identical" -- amended 2026-08-28 by the Stage 4 exit audit, which found the older blanket form ("each binding test supplies its own local diff --git a/tests/unit/numerics/test_assembly.py b/tests/unit/numerics/test_assembly.py index 4b0b5dd..0822c92 100644 --- a/tests/unit/numerics/test_assembly.py +++ b/tests/unit/numerics/test_assembly.py @@ -103,8 +103,9 @@ class _CapturingPressureCoupling(PressureCoupling): """Records the exact `boundary_conditions` mapping it was constructed with -- the pressure-coupling analogue of `_CapturingAdvection`/ `_CapturingDiffusion` above (TASK-027). Accepts (and ignores) - `tolerance`/`max_iterations` (TASK-033, 2026-08-29), matching - `register_pressure_coupling`'s new four-argument factory shape. + `tolerance`/`max_iterations` (TASK-033, 2026-08-29) and + `periodic_pairs` (TASK-034, 2026-08-29), matching + `register_pressure_coupling`'s new five-argument factory shape. """ def __init__( @@ -113,8 +114,9 @@ def __init__( boundary_conditions: Mapping[str, BoundaryCondition], tolerance: float, max_iterations: int, + periodic_pairs: Mapping[str, str], ) -> None: - del tolerance, max_iterations + del tolerance, max_iterations, periodic_pairs super().__init__(linear_solver) self.received_boundary_conditions = boundary_conditions diff --git a/tests/unit/numerics/test_divergence_contract.py b/tests/unit/numerics/test_divergence_contract.py index 4da1a3f..f08c238 100644 --- a/tests/unit/numerics/test_divergence_contract.py +++ b/tests/unit/numerics/test_divergence_contract.py @@ -54,7 +54,7 @@ def evaluate(self, field: Field, face: int) -> float: def _green_gauss_divergence() -> GreenGaussDivergence: condition = _ZeroGradientCondition() return GreenGaussDivergence( - {"north": condition, "south": condition, "east": condition, "west": condition} + {"north": condition, "south": condition, "east": condition, "west": condition}, {} ) @@ -176,7 +176,7 @@ def evaluate(self, field: Field, face: int) -> float: condition = _LinearDirichlet() scheme = GreenGaussDivergence( - {"north": condition, "south": condition, "east": condition, "west": condition} + {"north": condition, "south": condition, "east": condition, "west": condition}, {} ) field = VectorField( mesh, @@ -195,7 +195,7 @@ def evaluate(self, field: Field, face: int) -> float: def test_green_gauss_divergence_raises_for_an_unconfigured_boundary_face() -> None: mesh = _mesh() - scheme = GreenGaussDivergence({}) + scheme = GreenGaussDivergence({}, {}) field = VectorField(mesh, "velocity", num_components=2, initial_value=(1.0, 0.0)) with pytest.raises(UnconfiguredBoundaryFaceError): @@ -206,9 +206,60 @@ def test_green_gauss_divergence_rejects_a_field_with_the_wrong_component_shape() mesh = _mesh() condition = _ZeroGradientCondition() scheme = GreenGaussDivergence( - {"north": condition, "south": condition, "east": condition, "west": condition} + {"north": condition, "south": condition, "east": condition, "west": condition}, {} ) field = ScalarField(mesh, "temperature", initial_value=1.0) with pytest.raises(IncompatibleVectorFieldError): scheme.divergence(field) + + +def test_green_gauss_divergence_is_periodic_aware() -> None: + # TASK-034 (Stage 5): a face named in `periodic_pairs` must not raise + # `UnconfiguredBoundaryFaceError`. A uniform velocity field is + # trivially periodic and must give exactly zero divergence everywhere + # (verified numerically before being written here) -- this is also + # the mechanism Stage 5 Completion Criterion 4's own "uniform flow on + # a fully periodic domain stays divergence-free" null test depends on + # being true at the `PISO` level. A non-uniform, hand-assigned field + # (not from a smooth function, so nothing is accidentally continuous + # across the wrap) must give a real nonzero divergence matching a + # value hand-derived directly from `mesh.wrapped_neighbour_cell` -- + # proving the periodic branch reads the real wrapped neighbour rather + # than silently skipping the face. + mesh = _mesh() + assert isinstance(mesh, StructuredCartesianMesh) + all_periodic = {"north": "south", "south": "north", "east": "west", "west": "east"} + scheme = GreenGaussDivergence({}, all_periodic) + + uniform = VectorField(mesh, "velocity", num_components=2, initial_value=(1.3, -0.7)) + uniform_result = scheme.divergence(uniform) + assert torch.allclose( + uniform_result, torch.zeros(mesh.num_cells, dtype=torch.float64), atol=1e-9 + ) + + nonuniform = VectorField(mesh, "velocity", num_components=2, initial_value=(0.0, 0.0)) + for cell in range(mesh.num_cells): + nonuniform.set_value_at(cell, (float((cell * 7) % 5) - 2.0, float((cell * 3) % 4) - 1.5)) + nonuniform_result = scheme.divergence(nonuniform) + assert not torch.allclose( + nonuniform_result, torch.zeros(mesh.num_cells, dtype=torch.float64), atol=1e-9 + ) + + cell = 0 + total = 0.0 + for face in range(mesh.num_faces): + owner, neighbour = mesh.face_neighbours(face) + if owner != cell and neighbour != cell: + continue + sign = 1.0 if owner == cell else -1.0 + normal_x, normal_y = mesh.face_normal(face) + if neighbour is None: + neighbour = mesh.wrapped_neighbour_cell(face) + owner_x, owner_y = nonuniform.value_at(owner) + neighbour_x, neighbour_y = nonuniform.value_at(neighbour) + face_x, face_y = (owner_x + neighbour_x) / 2, (owner_y + neighbour_y) / 2 + face_normal_velocity = face_x * normal_x + face_y * normal_y + total += sign * face_normal_velocity * mesh.face_area(face) + total /= mesh.cell_volume(cell) + assert float(nonuniform_result[cell]) == pytest.approx(total, abs=1e-9) diff --git a/tests/unit/numerics/test_gradient_contract.py b/tests/unit/numerics/test_gradient_contract.py index 25317ba..0f9439c 100644 --- a/tests/unit/numerics/test_gradient_contract.py +++ b/tests/unit/numerics/test_gradient_contract.py @@ -57,7 +57,7 @@ def evaluate(self, field: Field, face: int) -> float: def _green_gauss_gradient() -> GreenGaussGradient: condition = _ZeroGradientCondition() return GreenGaussGradient( - {"north": condition, "south": condition, "east": condition, "west": condition} + {"north": condition, "south": condition, "east": condition, "west": condition}, {} ) @@ -175,7 +175,7 @@ def evaluate(self, field: Field, face: int) -> float: condition = _LinearDirichlet() scheme = GreenGaussGradient( - {"north": condition, "south": condition, "east": condition, "west": condition} + {"north": condition, "south": condition, "east": condition, "west": condition}, {} ) field = ScalarField(mesh, "phi", initial_value=lambda x, y: 2.0 * x - 3.0 * y) @@ -187,8 +187,65 @@ def evaluate(self, field: Field, face: int) -> float: def test_green_gauss_gradient_raises_for_an_unconfigured_boundary_face() -> None: mesh = _mesh() - scheme = GreenGaussGradient({}) + scheme = GreenGaussGradient({}, {}) field = ScalarField(mesh, "phi", initial_value=1.0) with pytest.raises(UnconfiguredBoundaryFaceError): scheme.gradient(field) + + +def test_green_gauss_gradient_is_periodic_aware() -> None: + # TASK-034 (Stage 5): a face named in `periodic_pairs` must not raise + # `UnconfiguredBoundaryFaceError` even with no `BoundaryCondition` at + # all configured. A genuinely linear field is not itself periodic (it + # does not agree with itself across the wrap), so "exact for a linear + # field" is not the claim to check here -- a uniform field is + # trivially periodic and must give an exactly zero gradient + # everywhere, checked first; a non-uniform field with hand-assigned + # per-cell values (not from a smooth function, so nothing about it is + # accidentally continuous across the wrap) must give a real nonzero + # gradient that matches a hand-derived value built directly from + # `mesh.wrapped_neighbour_cell` -- proving the periodic branch reads + # the real wrapped neighbour rather than silently skipping the face + # (which would also avoid raising, but would leave every boundary + # cell's own gradient contribution wrong). Verified numerically before + # being written into this test, not assumed. + mesh = _mesh() + assert isinstance(mesh, StructuredCartesianMesh) + all_periodic = {"north": "south", "south": "north", "east": "west", "west": "east"} + scheme = GreenGaussGradient({}, all_periodic) + + uniform = ScalarField(mesh, "phi", initial_value=3.7) + uniform_result = scheme.gradient(uniform) + assert torch.allclose( + uniform_result, torch.zeros((mesh.num_cells, 2), dtype=torch.float64), atol=1e-9 + ) + + nonuniform = ScalarField(mesh, "phi", initial_value=0.0) + for cell in range(mesh.num_cells): + nonuniform.set_value_at(cell, float((cell * 7) % 5) - 2.0) + nonuniform_result = scheme.gradient(nonuniform) + assert not torch.allclose( + nonuniform_result, torch.zeros((mesh.num_cells, 2), dtype=torch.float64), atol=1e-9 + ) + + # Hand-derived directly from the mesh's own wrapped-neighbour lookup, + # not from this scheme -- independent of the implementation under + # test, the same discipline `_numerics.py`'s own `face_normal_velocity` + # helper uses. + cell = 0 + total = torch.zeros(2, dtype=torch.float64) + for face in range(mesh.num_faces): + owner, neighbour = mesh.face_neighbours(face) + if owner != cell and neighbour != cell: + continue + sign = 1.0 if owner == cell else -1.0 + normal_x, normal_y = mesh.face_normal(face) + if neighbour is None: + neighbour = mesh.wrapped_neighbour_cell(face) + face_value = (nonuniform.value_at(owner) + nonuniform.value_at(neighbour)) / 2 + area = mesh.face_area(face) + total[0] += sign * face_value * normal_x * area + total[1] += sign * face_value * normal_y * area + total = total / mesh.cell_volume(cell) + assert torch.allclose(nonuniform_result[cell], total, atol=1e-9) diff --git a/tests/unit/test_navier_stokes_timestep.py b/tests/unit/test_navier_stokes_timestep.py new file mode 100644 index 0000000..2706363 --- /dev/null +++ b/tests/unit/test_navier_stokes_timestep.py @@ -0,0 +1,804 @@ +"""Binds `tests/features/navier_stokes_timestep.feature` (TASK-034, +Stage 5's fifth and last task in build order) -- assembles TASK-031/032/ +033 into one incompressible Navier-Stokes timestep +(`pyflow.engine.simulation.navier_stokes_step`), then validates it. Not +a golden demo -- no config file under `examples/golden-demos/`, no CLI +subprocess run, the same `tests/unit/` shape every prior Stage 4/5 +numerical-scheme feature file already established. Reuses +`tests/unit/_numerics.py`'s shared building blocks where they fit; +supplies its own local doubles and `_Context` otherwise, per this +directory's own "local by default, shared where genuinely identical" +convention. + +**Per-component wall velocities use `BoundaryFaceConfig.field_values` +directly (`velocity.0`/`velocity.1`, `VectorField.component_name`), not +a dedicated `velocity_tangential` field -- a design finding, not an +oversight.** Stage 5's own design question two (resolved 2026-08-28, +`docs/planning/roadmap.md`) named `velocity_tangential` as its answer, +but `field_values`/`field_gradients` (TASK-031c) landed the very next +day and already supply the exact general mechanism that question needed +-- a per-field-name override at one wall -- with no new config field, no +new per-wall tangential-axis wiring inside `assembly.py` (which stays +field-name-agnostic, per its own established discipline), and no new +concept to document. A no-slip *stationary* wall needs no override at +all (`scalar_value=0.0` already zeroes both components identically, +since normal and tangential are both zero); a *moving* wall (Couette's +plate, the cavity's lid) sets `field_values` for both component names +directly. See this module's own commit message and +`docs/planning/roadmap.md` TASK-034's own Design decision for the full +reasoning -- recorded explicitly, per root `CLAUDE.md`'s Validation +section, rather than silently building something different from what +was decided without saying so. +""" + +from __future__ import annotations + +import math +from dataclasses import dataclass, field + +import pytest +import torch +from fixtures.ghia_1982_re100 import ( + PRIMARY_VORTEX_CENTER, + U_VELOCITY_ALONG_VERTICAL_CENTERLINE, + V_VELOCITY_ALONG_HORIZONTAL_CENTERLINE, +) +from pytest_bdd import given, scenarios, then, when + +from pyflow.configuration.schema import ( + BoundaryConditionsConfig, + BoundaryFaceConfig, + NumericsConfig, +) +from pyflow.engine import simulation +from pyflow.engine.field import Field +from pyflow.engine.mesh import StructuredCartesianMesh +from pyflow.engine.numerics.advection import FirstOrderUpwindAdvection +from pyflow.engine.numerics.assembly import ( + AssembledNumerics, + assemble_numerics, + register_pressure_coupling, +) +from pyflow.engine.numerics.boundary_condition import BoundaryCondition, DirichletBoundaryCondition +from pyflow.engine.numerics.diffusion import CentralDifferenceDiffusion +from pyflow.engine.numerics.divergence import GreenGaussDivergence +from pyflow.engine.numerics.linear_solver import ConjugateGradientSolver +from pyflow.engine.numerics.pressure_coupling import PISO, PressureCoupling +from pyflow.engine.numerics.time_integrator import RK4Integrator +from pyflow.engine.scalar_field import PressureField, ScalarField +from pyflow.engine.simulation import NavierStokesStepResult +from pyflow.engine.vector_field import VectorField + +scenarios("navier_stokes_timestep.feature") + +_VELOCITY_NAME = "velocity" +_U_NAME = VectorField.component_name(_VELOCITY_NAME, 0) +_V_NAME = VectorField.component_name(_VELOCITY_NAME, 1) +_DT = 0.02 +_TOLERANCE = 1e-6 + + +def _no_slip_mesh(extent: tuple[int, int] = (4, 3)) -> StructuredCartesianMesh: + # Non-square, non-trivially-origined, non-unit spacing -- the same + # "distinct factors" discipline every other Stage 4/5 fixture in + # this repository follows. + return StructuredCartesianMesh(origin=(0.4, -0.3), spacing=(0.25, 0.2), extent=extent) + + +def _no_slip_boundary_conditions() -> dict[str, BoundaryCondition]: + condition = DirichletBoundaryCondition(0.0) + return {"north": condition, "south": condition, "east": condition, "west": condition} + + +def _divergent_velocity(mesh: StructuredCartesianMesh) -> VectorField: + # Non-axis-aligned, non-uniform -- the same fixture shape + # `test_piso_pressure_coupling.py`/`test_pressure_correction_loop.py` + # already use, so a real divergent flow needs genuine correction. + cx, cy = 1.5, -0.9 + + def value(x: float, y: float) -> tuple[float, float]: + return (0.6 * (x - cx) - 0.2 * (y - cy), 0.3 * (x - cx) + 0.9 * (y - cy)) + + return VectorField(mesh, _VELOCITY_NAME, num_components=2, initial_value=value) + + +def _real_numerics( + boundary_conditions: dict[str, BoundaryCondition], + periodic_pairs: dict[str, str] | None = None, + viscosity: float = 1.0, + pressure_coupling: PressureCoupling | None = None, +) -> AssembledNumerics: + pairs = periodic_pairs or {} + solver = ConjugateGradientSolver(tolerance=1e-10, max_iterations=1000) + return AssembledNumerics( + advection=FirstOrderUpwindAdvection(boundary_conditions, pairs), + diffusion=CentralDifferenceDiffusion(boundary_conditions, pairs, viscosity), + time_integration=RK4Integrator(), + linear_solver=solver, + pressure_coupling=pressure_coupling + or PISO(solver, boundary_conditions, tolerance=_TOLERANCE, periodic_pairs=pairs), + boundary_conditions=boundary_conditions, + names={}, + ) + + +@dataclass +class _Context: + mesh: StructuredCartesianMesh + numerics: AssembledNumerics | None = None + fields: dict[str, Field] = field(default_factory=dict) + result: NavierStokesStepResult | None = None + other_result: NavierStokesStepResult | None = None + history: list[torch.Tensor] = field(default_factory=list) + divergence_history: list[float] = field(default_factory=list) + lid_speed: float = 0.0 + viscosity: float = 0.0 + channel_height: float = 0.0 + channel_bottom: float = 0.0 + steady: bool = False + ke_history: list[float] = field(default_factory=list) + cavity_errors: list[float] = field(default_factory=list) + cavity_finest_fields: dict[str, Field] = field(default_factory=dict) + cavity_finest_mesh: StructuredCartesianMesh | None = None + cavity_finest_extent: tuple[int, int] = (0, 0) + + +# -- Given ------------------------------------------------------------- + + +@given("a closed, no-slip domain with a divergent initial velocity field", target_fixture="ctx") +def _given_closed_divergent() -> _Context: + mesh = _no_slip_mesh() + bcs = _no_slip_boundary_conditions() + ctx = _Context(mesh=mesh, numerics=_real_numerics(bcs)) + velocity = _divergent_velocity(mesh) + ctx.fields = {c.name: c for c in velocity.decompose()} + return ctx + + +@given( + "a fully periodic domain with a uniform, non-axis-aligned velocity field and zero viscosity", + target_fixture="ctx", +) +def _given_periodic_uniform() -> _Context: + mesh = _no_slip_mesh() + pairs = {"north": "south", "south": "north", "east": "west", "west": "east"} + ctx = _Context(mesh=mesh, numerics=_real_numerics({}, pairs, viscosity=0.0)) + # `CentralDifferenceDiffusion` still needs a positive coefficient to + # be constructed meaningfully; zero viscosity is expressed by scaling + # the coefficient itself to exactly 0.0, so the diffusive flux is + # identically zero at every face regardless of the field. + velocity = VectorField(mesh, _VELOCITY_NAME, num_components=2, initial_value=(1.3, -0.7)) + ctx.fields = {c.name: c for c in velocity.decompose()} + return ctx + + +@given("a closed, no-slip domain with the fluid initially at rest", target_fixture="ctx") +def _given_closed_at_rest() -> _Context: + mesh = _no_slip_mesh() + bcs = _no_slip_boundary_conditions() + ctx = _Context(mesh=mesh, numerics=_real_numerics(bcs)) + velocity = VectorField(mesh, _VELOCITY_NAME, num_components=2, initial_value=(0.0, 0.0)) + ctx.fields = {c.name: c for c in velocity.decompose()} + return ctx + + +class _MarkerPressureCoupling(PressureCoupling): + """Returns a distinctive, obviously-not-real-physics pressure field + (a large constant no real solve would produce for this fixture) -- + exists only to prove `navier_stokes_step` calls whichever + `PressureCoupling` `assemble_numerics` resolved, not a hardcoded + `PISO` (Stage 5 Completion Criterion 13's own substitution check). + """ + + _MARKER = 12345.0 + + def correct( + self, provisional_velocity: VectorField, dt: float + ) -> tuple[VectorField, ScalarField]: + del dt + pressure = PressureField(provisional_velocity.mesh, "pressure") + pressure.values[:] = self._MARKER + return provisional_velocity.copy(), pressure + + +@given( + "a PressureCoupling test double registered under its own name and selected by configuration", + target_fixture="ctx", +) +def _given_marker_pressure_coupling() -> _Context: + mesh = _no_slip_mesh() + name = "test_only_marker_pressure_coupling" + register_pressure_coupling( + name, + lambda linear_solver, boundary_conditions, tolerance, max_iterations, periodic_pairs: ( + _MarkerPressureCoupling(linear_solver) + ), + ) + config = NumericsConfig( + pressure_coupling=name, # type: ignore[arg-type] + boundary_conditions=BoundaryConditionsConfig( + north=BoundaryFaceConfig(type="dirichlet", velocity=None, scalar_value=0.0), + south=BoundaryFaceConfig(type="dirichlet", velocity=None, scalar_value=0.0), + east=BoundaryFaceConfig(type="dirichlet", velocity=None, scalar_value=0.0), + west=BoundaryFaceConfig(type="dirichlet", velocity=None, scalar_value=0.0), + ), + ) + numerics = assemble_numerics(config) + ctx = _Context(mesh=mesh, numerics=numerics) + velocity = _divergent_velocity(mesh) + ctx.fields = {c.name: c for c in velocity.decompose()} + return ctx + + +_COUETTE_LID_SPEED = 1.7 +_COUETTE_VISCOSITY = 0.8 +_COUETTE_EXTENT = (3, 8) +_COUETTE_SPACING = (0.2, 0.2) +_COUETTE_ORIGIN = (0.4, -0.3) + + +@given( + "a channel periodic in the flow direction, no-slip walls, one stationary and one moving " + "tangentially", + target_fixture="ctx", +) +def _given_couette_channel() -> _Context: + mesh = StructuredCartesianMesh( + origin=_COUETTE_ORIGIN, spacing=_COUETTE_SPACING, extent=_COUETTE_EXTENT + ) + south = DirichletBoundaryCondition(0.0) + north = DirichletBoundaryCondition(0.0, {_U_NAME: _COUETTE_LID_SPEED, _V_NAME: 0.0}) + bcs: dict[str, BoundaryCondition] = {"south": south, "north": north} + periodic = {"east": "west", "west": "east"} + solver = ConjugateGradientSolver(tolerance=1e-10, max_iterations=1000) + numerics = AssembledNumerics( + advection=FirstOrderUpwindAdvection(bcs, periodic), + diffusion=CentralDifferenceDiffusion(bcs, periodic, _COUETTE_VISCOSITY), + time_integration=RK4Integrator(), + linear_solver=solver, + pressure_coupling=PISO(solver, bcs, tolerance=1e-8, periodic_pairs=periodic), + boundary_conditions=bcs, + names={}, + ) + ctx = _Context( + mesh=mesh, + numerics=numerics, + lid_speed=_COUETTE_LID_SPEED, + viscosity=_COUETTE_VISCOSITY, + channel_height=_COUETTE_EXTENT[1] * _COUETTE_SPACING[1], + channel_bottom=_COUETTE_ORIGIN[1], + ) + velocity = VectorField(mesh, _VELOCITY_NAME, num_components=2, initial_value=(0.0, 0.0)) + ctx.fields = {c.name: c for c in velocity.decompose()} + return ctx + + +_NEGLIGIBLE_VISCOSITY = 1e-6 + + +@given( + "a closed, no-slip domain with a divergent initial velocity field and negligible viscosity", + target_fixture="ctx", +) +def _given_closed_divergent_inviscid() -> _Context: + mesh = _no_slip_mesh() + bcs = _no_slip_boundary_conditions() + ctx = _Context(mesh=mesh, numerics=_real_numerics(bcs, viscosity=_NEGLIGIBLE_VISCOSITY)) + velocity = _divergent_velocity(mesh) + ctx.fields = {c.name: c for c in velocity.decompose()} + return ctx + + +# -- When ---------------------------------------------------------------- + + +@when("one Navier-Stokes timestep is taken") +def _when_one_step(ctx: _Context) -> None: + assert ctx.numerics is not None + ctx.result = simulation.navier_stokes_step(ctx.fields, _VELOCITY_NAME, ctx.numerics, _DT) + + +@when("one Navier-Stokes timestep is taken twice from the same initial state") +def _when_stepped_twice(ctx: _Context) -> None: + bcs = _no_slip_boundary_conditions() + numerics_a = _real_numerics(bcs) + numerics_b = _real_numerics(bcs) + velocity_a = _divergent_velocity(ctx.mesh) + velocity_b = _divergent_velocity(ctx.mesh) + fields_a = {c.name: c for c in velocity_a.decompose()} + fields_b = {c.name: c for c in velocity_b.decompose()} + ctx.result = simulation.navier_stokes_step(fields_a, _VELOCITY_NAME, numerics_a, _DT) + ctx.other_result = simulation.navier_stokes_step(fields_b, _VELOCITY_NAME, numerics_b, _DT) + + +def _velocity_divergence(mesh: StructuredCartesianMesh, velocity: VectorField) -> torch.Tensor: + pairs = {"north": "south", "south": "north", "east": "west", "west": "east"} + return GreenGaussDivergence({}, pairs).divergence(velocity) + + +def _kinetic_energy(mesh: StructuredCartesianMesh, velocity: VectorField) -> float: + total = 0.0 + for cell in range(mesh.num_cells): + vx, vy = velocity.value_at(cell) + total += 0.5 * (vx * vx + vy * vy) * mesh.cell_volume(cell) + return total + + +@when("many Navier-Stokes timesteps are taken") +def _when_many_steps(ctx: _Context) -> None: + assert ctx.numerics is not None + fields = ctx.fields + for _ in range(20): + result = simulation.navier_stokes_step(fields, _VELOCITY_NAME, ctx.numerics, _DT) + fields = result.fields + ctx.history.append(result.corrected_velocity.values.clone()) + max_divergence = float( + _velocity_divergence(ctx.mesh, result.corrected_velocity).abs().max() + ) + ctx.divergence_history.append(max_divergence) + ctx.ke_history.append(_kinetic_energy(ctx.mesh, result.corrected_velocity)) + ctx.result = result + + +_STEADY_RESIDUAL_TOLERANCE = 1e-9 +_STEADY_MAX_STEPS = 5000 + + +@when("Navier-Stokes timesteps are taken until the flow reaches steady state") +def _when_run_to_steady_state(ctx: _Context) -> None: + assert ctx.numerics is not None + dt = simulation.stable_timestep(ctx.mesh, ctx.viscosity, ctx.lid_speed) + fields = ctx.fields + previous_u: torch.Tensor | None = None + # Steadiness is a measured residual, not a step count -- Stage 5 + # Completion Criterion 5's own "fails on not reaching it rather than + # silently comparing an unconverged field". `ctx.steady` stays False + # (its own default) if the cap is exhausted, and the Then step below + # asserts it explicitly rather than only checking the profile. + for _ in range(_STEADY_MAX_STEPS): + result = simulation.navier_stokes_step(fields, _VELOCITY_NAME, ctx.numerics, dt) + fields = result.fields + u_field = fields[_U_NAME] + assert isinstance(u_field, ScalarField) + u_values = u_field.values + if previous_u is not None: + residual = float((u_values - previous_u).abs().max()) + if residual < _STEADY_RESIDUAL_TOLERANCE: + ctx.steady = True + ctx.result = result + return + previous_u = u_values.clone() + ctx.result = result + + +# -- Then ------------------------------------------------------------------ + + +@then("the provisional velocity, the corrected velocity, and the pressure field are all present") +def _then_all_present(ctx: _Context) -> None: + assert ctx.result is not None + assert isinstance(ctx.result.provisional_velocity, VectorField) + assert isinstance(ctx.result.corrected_velocity, VectorField) + assert isinstance(ctx.result.pressure, ScalarField) + + +@then("the corrected velocity differs from the provisional velocity") +def _then_corrected_differs(ctx: _Context) -> None: + assert ctx.result is not None + assert not torch.equal( + ctx.result.corrected_velocity.values, ctx.result.provisional_velocity.values + ) + + +@then("the corrected velocity's own divergence is smaller than the provisional velocity's") +def _then_divergence_reduced(ctx: _Context) -> None: + assert ctx.result is not None + mesh = ctx.mesh + divergence_fn = GreenGaussDivergence(_no_slip_boundary_conditions(), {}) + provisional_divergence = float( + divergence_fn.divergence(ctx.result.provisional_velocity).abs().max() + ) + corrected_divergence = float( + divergence_fn.divergence(ctx.result.corrected_velocity).abs().max() + ) + assert corrected_divergence < provisional_divergence + assert mesh is ctx.result.corrected_velocity.mesh + + +@then("the velocity field is exactly the same uniform value at every step") +def _then_uniform_unchanged(ctx: _Context) -> None: + initial = torch.tensor([1.3, -0.7], dtype=torch.float64) + for values in ctx.history: + assert torch.allclose(values, initial.expand_as(values), atol=1e-9) + + +@then("the velocity field's own divergence never leaves solver tolerance at any step") +def _then_divergence_in_tolerance(ctx: _Context) -> None: + assert ctx.divergence_history + for max_divergence in ctx.divergence_history: + assert max_divergence <= _TOLERANCE + + +@then("the velocity field stays at rest to floating-point tolerance at every step") +def _then_stays_at_rest(ctx: _Context) -> None: + assert ctx.history + for values in ctx.history: + assert torch.allclose(values, torch.zeros_like(values), atol=1e-9) + + +@then("both runs produce identical corrected velocity and pressure fields") +def _then_identical_runs(ctx: _Context) -> None: + assert ctx.result is not None + assert ctx.other_result is not None + torch.testing.assert_close( + ctx.result.corrected_velocity.values, + ctx.other_result.corrected_velocity.values, + rtol=0, + atol=0, + ) + torch.testing.assert_close( + ctx.result.pressure.values, ctx.other_result.pressure.values, rtol=0, atol=0 + ) + + +@then("the test double's own distinctive pressure value appears in the result, not a real solve's") +def _then_marker_pressure_present(ctx: _Context) -> None: + assert ctx.result is not None + assert torch.allclose( + ctx.result.pressure.values, + torch.full_like(ctx.result.pressure.values, _MarkerPressureCoupling._MARKER), + ) + + +@then( + "the steady streamwise velocity profile matches the exact linear Couette solution at " + "solver tolerance" +) +def _then_couette_profile_matches(ctx: _Context) -> None: + assert ctx.steady, ( + f"flow did not reach steady state within {_STEADY_MAX_STEPS} steps " + f"(residual tolerance {_STEADY_RESIDUAL_TOLERANCE})" + ) + assert ctx.result is not None + u_field = ctx.result.fields[_U_NAME] + assert isinstance(u_field, ScalarField) + for cell in range(ctx.mesh.num_cells): + _x, y = ctx.mesh.cell_centroid(cell) + exact = ctx.lid_speed * (y - ctx.channel_bottom) / ctx.channel_height + actual = float(u_field.value_at(cell)) + assert actual == pytest.approx(exact, abs=1e-6), ( + f"cell {cell} at y={y}: expected {exact}, got {actual}" + ) + + +@then("the wall-normal velocity component stays zero everywhere") +def _then_wall_normal_velocity_zero(ctx: _Context) -> None: + assert ctx.result is not None + v_field = ctx.result.fields[_V_NAME] + assert isinstance(v_field, ScalarField) + assert torch.allclose(v_field.values, torch.zeros_like(v_field.values), atol=1e-9) + + +@then("total kinetic energy never increases from one step to the next") +def _then_ke_never_increases(ctx: _Context) -> None: + assert len(ctx.ke_history) >= 2 + for previous, current in zip(ctx.ke_history, ctx.ke_history[1:], strict=False): + assert current <= previous + 1e-12, f"kinetic energy increased: {previous} -> {current}" + + +# -- Taylor-Green vortex decay: the emergent-phenomenon pair --------------- + + +_TG_AMPLITUDE = 0.3 +_TG_EXTENT = (12, 12) +_TG_DOMAIN_LENGTH = 1.0 +_TG_MEASURE_STEPS = (5, 40) +# Measured directly before being trusted (this module's own commit +# message / `docs/planning/roadmap.md` TASK-034's own Design decision): +# at `_TG_VISCOSITY_MATCHED`, the measured decay rate agreed with the +# exact closed form to within ~0.3%; at `_TG_VISCOSITY_MISMATCHED` (100x +# smaller -- the mesh's own advective numerical diffusion stays fixed +# while the physical one shrinks) the measured rate was off by a factor +# of roughly 3.8. Neither bound below was chosen to make a marginal +# result pass -- both keep real margin around the two measured ratios. +_TG_VISCOSITY_MATCHED = 0.05 +_TG_VISCOSITY_MISMATCHED = 0.0005 + + +def _taylor_green_wavenumber() -> float: + + return 2 * math.pi / _TG_DOMAIN_LENGTH + + +@dataclass +class _TaylorGreenContext: + viscosity: float + measured_rate: float = 0.0 + exact_rate: float = 0.0 + + +def _run_taylor_green(viscosity: float) -> _TaylorGreenContext: + + k = _taylor_green_wavenumber() + mesh = StructuredCartesianMesh( + origin=(0.0, 0.0), + spacing=(_TG_DOMAIN_LENGTH / _TG_EXTENT[0], _TG_DOMAIN_LENGTH / _TG_EXTENT[1]), + extent=_TG_EXTENT, + ) + periodic = {"north": "south", "south": "north", "east": "west", "west": "east"} + + def value(x: float, y: float) -> tuple[float, float]: + return ( + _TG_AMPLITUDE * math.cos(k * x) * math.sin(k * y), + -_TG_AMPLITUDE * math.sin(k * x) * math.cos(k * y), + ) + + numerics = _real_numerics({}, periodic, viscosity=viscosity) + dt = simulation.stable_timestep(mesh, viscosity, _TG_AMPLITUDE, safety_factor=0.25) + velocity = VectorField(mesh, _VELOCITY_NAME, num_components=2, initial_value=value) + fields: dict[str, Field] = {c.name: c for c in velocity.decompose()} + probe_cell = mesh.cell_id(_TG_EXTENT[0] // 4, _TG_EXTENT[1] // 4) + + amplitudes: list[float] = [] + times: list[float] = [] + elapsed = 0.0 + max_step = max(_TG_MEASURE_STEPS) + u_field = fields[_U_NAME] + assert isinstance(u_field, ScalarField) + for step in range(max_step + 1): + if step in _TG_MEASURE_STEPS: + amplitudes.append(abs(float(u_field.value_at(probe_cell)))) + times.append(elapsed) + result = simulation.navier_stokes_step(fields, _VELOCITY_NAME, numerics, dt) + fields = result.fields + u_field = fields[_U_NAME] + assert isinstance(u_field, ScalarField) + elapsed += dt + + first_amplitude, last_amplitude = amplitudes[0], amplitudes[-1] + first_time, last_time = times[0], times[-1] + measured_rate = -math.log(last_amplitude / first_amplitude) / (last_time - first_time) + exact_rate = 2 * k * k * viscosity + return _TaylorGreenContext( + viscosity=viscosity, measured_rate=measured_rate, exact_rate=exact_rate + ) + + +@given( + "a Taylor-Green vortex on a periodic domain at a viscosity where physical diffusion dominates", + target_fixture="tg_ctx", +) +def _given_taylor_green_matched() -> float: + return _TG_VISCOSITY_MATCHED + + +@given( + "a Taylor-Green vortex on a periodic domain at a viscosity where numerical diffusion dominates", + target_fixture="tg_ctx", +) +def _given_taylor_green_mismatched() -> float: + return _TG_VISCOSITY_MISMATCHED + + +@when("the vortex is advanced and its own decay rate is measured", target_fixture="tg_result") +def _when_taylor_green_measured(tg_ctx: float) -> _TaylorGreenContext: + return _run_taylor_green(tg_ctx) + + +@then("the measured decay rate matches the exact closed-form rate closely") +def _then_taylor_green_matches(tg_result: _TaylorGreenContext) -> None: + ratio = tg_result.measured_rate / tg_result.exact_rate + assert 0.9 <= ratio <= 1.1, ( + f"measured/exact ratio {ratio} outside the expected close-match band" + ) + + +@then("the measured decay rate does not match the exact closed-form rate") +def _then_taylor_green_mismatches(tg_result: _TaylorGreenContext) -> None: + ratio = tg_result.measured_rate / tg_result.exact_rate + assert ratio > 2.0 or ratio < 0.5, ( + f"measured/exact ratio {ratio} was unexpectedly close to a match" + ) + + +# -- Lid-driven cavity against Ghia, Ghia & Shin (1982) -------------------- +# +# This project's most computationally expensive scenario, deliberately +# (`docs/planning/roadmap.md` TASK-034's own Design decision): three real +# runs to a measured steady state. Resolutions chosen odd (9, 13, 17) so +# the vertical/horizontal centreline (x=0.5, y=0.5 in unit-cavity +# coordinates) always lands exactly on a column/row of cell centres, no +# interpolation needed for the two Ghia comparisons themselves. + +_CAVITY_RESOLUTIONS = (9, 13, 17) +_CAVITY_REYNOLDS_NUMBER = 100 +_CAVITY_LID_SPEED = 1.0 +_CAVITY_VISCOSITY = _CAVITY_LID_SPEED / _CAVITY_REYNOLDS_NUMBER # Re = U*L/nu, L = 1 +_CAVITY_STEADY_RESIDUAL_TOLERANCE = 1e-6 +_CAVITY_MAX_STEPS = 6000 +# Ghia's own primary-vortex distance bound: measured directly on a real +# (coarser, n=14) run before being trusted -- the detected centre landed +# 0.019 away from Ghia's own (0.6172, 0.7344) in unit-cavity coordinates. +# 0.1 keeps an order of magnitude of margin around that measured distance +# without being so loose a genuinely wrong vortex location could still +# pass. +_CAVITY_VORTEX_DISTANCE_TOLERANCE = 0.1 +# The same measured run found clear opposite-sign vorticity (magnitude +# 0.13-0.43) throughout each bottom corner's own sub-region; this +# threshold sits an order of magnitude below that measured range. +_CAVITY_VORTICITY_NOISE_THRESHOLD = 0.02 + + +@dataclass +class _CavityRun: + resolution: int + error: float + u_field: ScalarField | None = None + v_field: ScalarField | None = None + mesh: StructuredCartesianMesh | None = None + steady: bool = False + + +def _run_cavity(n: int) -> _CavityRun: + mesh = StructuredCartesianMesh(origin=(0.0, 0.0), spacing=(1.0 / n, 1.0 / n), extent=(n, n)) + lid = DirichletBoundaryCondition(0.0, {_U_NAME: _CAVITY_LID_SPEED, _V_NAME: 0.0}) + wall = DirichletBoundaryCondition(0.0) + bcs: dict[str, BoundaryCondition] = {"north": lid, "south": wall, "east": wall, "west": wall} + numerics = _real_numerics(bcs, viscosity=_CAVITY_VISCOSITY) + dt = simulation.stable_timestep(mesh, _CAVITY_VISCOSITY, _CAVITY_LID_SPEED, safety_factor=0.25) + + velocity = VectorField(mesh, _VELOCITY_NAME, num_components=2, initial_value=(0.0, 0.0)) + fields: dict[str, Field] = {c.name: c for c in velocity.decompose()} + previous_u: torch.Tensor | None = None + steady = False + for _ in range(_CAVITY_MAX_STEPS): + result = simulation.navier_stokes_step(fields, _VELOCITY_NAME, numerics, dt) + fields = result.fields + u_field = fields[_U_NAME] + assert isinstance(u_field, ScalarField) + u_values = u_field.values + if previous_u is not None: + residual = float((u_values - previous_u).abs().max()) / dt + if residual < _CAVITY_STEADY_RESIDUAL_TOLERANCE: + steady = True + break + previous_u = u_values.clone() + + u_field = fields[_U_NAME] + v_field = fields[_V_NAME] + assert isinstance(u_field, ScalarField) + assert isinstance(v_field, ScalarField) + + center = n // 2 + u_errors = [] + for y_ghia, u_ghia in U_VELOCITY_ALONG_VERTICAL_CENTERLINE: + row = min(int(y_ghia * n), n - 1) + cell = mesh.cell_id(center, row) + u_errors.append((float(u_field.value_at(cell)) - u_ghia) ** 2) + v_errors = [] + for x_ghia, v_ghia in V_VELOCITY_ALONG_HORIZONTAL_CENTERLINE: + column = min(int(x_ghia * n), n - 1) + cell = mesh.cell_id(column, center) + v_errors.append((float(v_field.value_at(cell)) - v_ghia) ** 2) + error = math.sqrt((sum(u_errors) + sum(v_errors)) / (len(u_errors) + len(v_errors))) + + return _CavityRun( + resolution=n, error=error, u_field=u_field, v_field=v_field, mesh=mesh, steady=steady + ) + + +@given( + "three lid-driven cavity meshes at increasing resolution, at Reynolds number 100", + target_fixture="cavity_resolutions", +) +def _given_cavity_resolutions() -> tuple[int, ...]: + return _CAVITY_RESOLUTIONS + + +@when("each is run to a measured steady state", target_fixture="cavity_runs") +def _when_cavity_runs(cavity_resolutions: tuple[int, ...]) -> list[_CavityRun]: + return [_run_cavity(n) for n in cavity_resolutions] + + +@then( + "the error against Ghia's centreline profiles decreases monotonically across the three " + "resolutions" +) +def _then_cavity_error_decreases(cavity_runs: list[_CavityRun]) -> None: + for run in cavity_runs: + assert run.steady, ( + f"resolution {run.resolution} did not reach steady state within " + f"{_CAVITY_MAX_STEPS} steps" + ) + errors = [run.error for run in cavity_runs] + for previous, current in zip(errors, errors[1:], strict=False): + assert current < previous, f"error did not decrease monotonically: {errors}" + + +def _finest_run(cavity_runs: list[_CavityRun]) -> _CavityRun: + return max(cavity_runs, key=lambda run: run.resolution) + + +def _vorticity_at(run: _CavityRun, i: int, j: int) -> float: + assert run.mesh is not None + assert run.u_field is not None + assert run.v_field is not None + dx = 1.0 / run.resolution + dv_dx = ( + float(run.v_field.value_at(run.mesh.cell_id(i + 1, j))) + - float(run.v_field.value_at(run.mesh.cell_id(i - 1, j))) + ) / (2 * dx) + du_dy = ( + float(run.u_field.value_at(run.mesh.cell_id(i, j + 1))) + - float(run.u_field.value_at(run.mesh.cell_id(i, j - 1))) + ) / (2 * dx) + return dv_dx - du_dy + + +@then("the finest resolution's primary vortex centre is within a stated distance of Ghia's own") +def _then_primary_vortex_near_ghia(cavity_runs: list[_CavityRun]) -> None: + run = _finest_run(cavity_runs) + assert run.mesh is not None + assert run.u_field is not None + assert run.v_field is not None + n = run.resolution + margin = n // 4 + best: tuple[float, int, int] | None = None + for i in range(margin, n - margin): + for j in range(margin, n - margin): + cell = run.mesh.cell_id(i, j) + u_val = float(run.u_field.value_at(cell)) + v_val = float(run.v_field.value_at(cell)) + magnitude = math.hypot(u_val, v_val) + if best is None or magnitude < best[0]: + best = (magnitude, i, j) + assert best is not None + _magnitude, i, j = best + x, y = run.mesh.cell_centroid(run.mesh.cell_id(i, j)) + ghia_x, ghia_y = PRIMARY_VORTEX_CENTER + distance = math.hypot(x - ghia_x, y - ghia_y) + assert distance < _CAVITY_VORTEX_DISTANCE_TOLERANCE, ( + f"detected primary vortex at ({x}, {y}), {distance} from Ghia's own {PRIMARY_VORTEX_CENTER}" + ) + + +@then( + "the finest resolution shows both downstream secondary corner vortices, rotating opposite " + "the primary" +) +def _then_secondary_vortices_present(cavity_runs: list[_CavityRun]) -> None: + run = _finest_run(cavity_runs) + assert run.mesh is not None + n = run.resolution + margin = n // 4 + best: tuple[float, int, int] | None = None + for i in range(margin, n - margin): + for j in range(margin, n - margin): + assert run.u_field is not None + assert run.v_field is not None + cell = run.mesh.cell_id(i, j) + magnitude = math.hypot( + float(run.u_field.value_at(cell)), float(run.v_field.value_at(cell)) + ) + if best is None or magnitude < best[0]: + best = (magnitude, i, j) + assert best is not None + _magnitude, pi, pj = best + primary_vorticity = _vorticity_at(run, pi, pj) + primary_positive = primary_vorticity > 0 + + def _has_opposite_sign_vorticity(i_range: range, j_range: range) -> bool: + for i in i_range: + for j in j_range: + vorticity = _vorticity_at(run, i, j) + if (vorticity > 0) != primary_positive and abs( + vorticity + ) > _CAVITY_VORTICITY_NOISE_THRESHOLD: + return True + return False + + bottom_left = _has_opposite_sign_vorticity(range(1, n // 3), range(1, n // 3)) + bottom_right = _has_opposite_sign_vorticity(range(2 * n // 3, n - 1), range(1, n // 3)) + assert bottom_left, "no opposite-sign (secondary) vorticity found near the bottom-left corner" + assert bottom_right, "no opposite-sign (secondary) vorticity found near the bottom-right corner" diff --git a/tests/unit/test_piso_pressure_coupling.py b/tests/unit/test_piso_pressure_coupling.py index e5fd3c6..d5567da 100644 --- a/tests/unit/test_piso_pressure_coupling.py +++ b/tests/unit/test_piso_pressure_coupling.py @@ -157,7 +157,7 @@ def _when_corrected(ctx: _Context) -> None: def _divergence(ctx: _Context, velocity: VectorField) -> torch.Tensor: from pyflow.engine.numerics.divergence import GreenGaussDivergence - return GreenGaussDivergence(ctx.boundary_conditions).divergence(velocity) + return GreenGaussDivergence(ctx.boundary_conditions, {}).divergence(velocity) @then( @@ -187,3 +187,71 @@ def _then_max_smaller(ctx: _Context) -> None: @then("a pressure solve non-convergence error is raised") def _then_non_convergence_raised(ctx: _Context) -> None: assert isinstance(ctx.raised, PressureSolveDidNotConvergeError) + + +# -- A plain (non-BDD) unit test, not an acceptance criterion of its own -- +# +# TASK-034 (Stage 5): `_poisson_matrix` is now cached per `PISO` instance +# rather than rebuilt every `correct` call -- an implementation-detail +# performance fix (found while measuring the Lid-Driven Cavity +# validation's own real runtime, `pressure_coupling.py`'s own entry for +# the full reasoning), not a new physical-correctness claim, so a plain +# pytest test rather than a new Gherkin scenario. + + +def test_poisson_matrix_is_cached_across_repeated_correct_calls_on_the_same_mesh() -> None: + from pyflow.engine.numerics.linear_solver import ConjugateGradientSolver + + mesh = default_mesh() + condition = _ZeroNormalVelocity() + boundary_conditions: dict[str, BoundaryCondition] = { + "north": condition, + "south": condition, + "east": condition, + "west": condition, + } + solver = ConjugateGradientSolver(tolerance=1e-10, max_iterations=500) + piso = PISO(solver, boundary_conditions, tolerance=1e-8) + + velocity = VectorField( + mesh, "velocity", num_components=2, initial_value=lambda x, y: (0.6 * x, 0.3 * y) + ) + piso.correct(velocity, dt=0.1) + first_matrix = piso._cached_poisson_matrix + assert first_matrix is not None + + piso.correct(velocity, dt=0.1) + second_matrix = piso._cached_poisson_matrix + assert second_matrix is first_matrix, "the matrix was rebuilt on a second call, not reused" + + +def test_poisson_matrix_recomputes_for_a_genuinely_different_mesh() -> None: + from pyflow.engine.numerics.linear_solver import ConjugateGradientSolver + + condition = _ZeroNormalVelocity() + boundary_conditions: dict[str, BoundaryCondition] = { + "north": condition, + "south": condition, + "east": condition, + "west": condition, + } + solver = ConjugateGradientSolver(tolerance=1e-10, max_iterations=500) + piso = PISO(solver, boundary_conditions, tolerance=1e-8) + + mesh_a = default_mesh(extent=(3, 2)) + velocity_a = VectorField( + mesh_a, "velocity", num_components=2, initial_value=lambda x, y: (0.6 * x, 0.3 * y) + ) + piso.correct(velocity_a, dt=0.1) + matrix_a = piso._cached_poisson_matrix + + mesh_b = default_mesh(extent=(4, 3)) + velocity_b = VectorField( + mesh_b, "velocity", num_components=2, initial_value=lambda x, y: (0.6 * x, 0.3 * y) + ) + piso.correct(velocity_b, dt=0.1) + matrix_b = piso._cached_poisson_matrix + + assert matrix_a is not None + assert matrix_b is not None + assert matrix_a.shape != matrix_b.shape diff --git a/tests/unit/test_pressure_field.py b/tests/unit/test_pressure_field.py index d5742ba..16f75a1 100644 --- a/tests/unit/test_pressure_field.py +++ b/tests/unit/test_pressure_field.py @@ -187,7 +187,7 @@ def _when_constant_added(ctx: _Context) -> None: pressure_boundary_conditions = { name: FixedGradientCondition(0.0) for name in ("north", "south", "east", "west") } - gradient_scheme = GreenGaussGradient(pressure_boundary_conditions) + gradient_scheme = GreenGaussGradient(pressure_boundary_conditions, {}) ctx.shifted_correction = ctx.provisional_velocity.values - ctx.dt * gradient_scheme.gradient( shifted ) diff --git a/tools/generators/generate_config_template.py b/tools/generators/generate_config_template.py index 3944903..551fd9e 100644 --- a/tools/generators/generate_config_template.py +++ b/tools/generators/generate_config_template.py @@ -177,9 +177,9 @@ "field_display.arrow_scale": "Valid: a positive number. Invalid: zero or negative.", "field_display.show_legend": "Valid: true or false.", "simulation.scalar_pattern": ( - 'Valid: null (no live simulation runs) or "gaussian_blob", the ' - "only built-in pattern this field currently accepts. Invalid: " - "any other string." + 'Valid: null (no live simulation runs), "gaussian_blob", or ' + '"sinusoidal_mode", the two built-in patterns this field ' + "currently accepts. Invalid: any other string." ), "simulation.velocity_pattern": ( 'Valid: null or "uniform", the only built-in pattern this field ' diff --git a/tools/validators/check_references.py b/tools/validators/check_references.py index 328bbd0..cca4885 100644 --- a/tools/validators/check_references.py +++ b/tools/validators/check_references.py @@ -109,9 +109,7 @@ # No Stage 5 *source* module is listed, deliberately: which modules # change is what that Stage's own design question one decides, and a # guess here would be the speculation P-016 refuses. -PLANNED: dict[str, str] = { - "tests/features/navier_stokes_timestep.feature": "TASK-034", -} +PLANNED: dict[str, str] = {} EXTS = (".md", ".py", ".yaml", ".yml", ".toml", ".cfg", ".txt", ".lock", ".json", ".ini") SPAN = re.compile(r"`([^`\n]+)`")