diff --git a/adr/ADR-003-modular-numerical-strategies.md b/adr/ADR-003-modular-numerical-strategies.md index 55ec10f..66aa1fb 100644 --- a/adr/ADR-003-modular-numerical-strategies.md +++ b/adr/ADR-003-modular-numerical-strategies.md @@ -115,6 +115,19 @@ interface change, no new ADR. `BoundaryFaceConfig.scalar_gradient` is `scalar_value`'s exact mirror, the gap TASK-028's own drafting had already named in advance for this task to resolve. +**Periodic Boundary followed the next day** (TASK-030, Stage 4, +2026-08-28, Stage 4's last task): not a seventh `BoundaryCondition` +implementation -- a periodic face is mesh geometry (`StructuredCartesianMesh. +wrapped_neighbour_cell`), not a prescribed value, so it bypasses this +ADR's own registry mechanism entirely rather than joining it. What *did* +land through this ADR's own pattern: `register_advection_scheme`/ +`register_diffusion_scheme`'s factory types widened to also receive a +`periodic_pairs` mapping (mirrors `diffusion_coefficient`'s own +TASK-024 precedent -- no interface change, no new ADR), and +`FirstOrderUpwindAdvection`/`CentralDifferenceDiffusion` both gained the +logic to consult it. `mvp.md`'s "Periodic (where practical)" bullet is +now real. + **What no longer exists, and is worth stating plainly: every component this ADR names now has a real concrete implementation.** All six -- Advection, Diffusion, Time Integration, Linear Solver, Pressure-Velocity @@ -122,11 +135,12 @@ Coupling, Boundary Condition (both its Dirichlet and Neumann shapes) -- went real across TASK-023 through TASK-029, and `assembly.py`'s own `_Null*` reference implementations, the whole Stage 3 Completion Criterion 1 carve-out this section has been tracking task by task, are -now zero. Only `source.py` and periodic boundary faces remain unbuilt -among this project's numerical machinery, and neither is one of this -ADR's own six -- `source.py` per TASK-018's own P-016 reasoning (no -second implementation identified yet), periodic per `boundary_condition. -py`'s own deliberately-narrower scope (TASK-019). +now zero. Only `source.py` remains unbuilt among this project's +numerical machinery, and it is not one of this ADR's own six -- +TASK-018's own P-016 reasoning (no second implementation identified +yet). Periodic boundary faces, the other item this paragraph used to +name here, are real as of TASK-030 (above) -- deliberately outside this +ADR's own six-component registry, not merely unbuilt. --- diff --git a/docs/architecture/CLAUDE.md b/docs/architecture/CLAUDE.md index 265d188..3f0dfa7 100644 --- a/docs/architecture/CLAUDE.md +++ b/docs/architecture/CLAUDE.md @@ -30,18 +30,24 @@ Grounded directly in `bootstrap.py`, `engine/simulation.py`, `engine/collocated_field.py` -- read those files, not this note, for anything beyond orientation. -**Two of its four sections carry a `Planned` subsection for a mechanism -that doesn't exist yet** (driving `simulation.step()` from a live render -loop, and checkpointing simulation state) -- each anchored to the -specific roadmap task that will build it (TASK-030, TASK-034) rather than -left open-ended, per the maintainer's direction that an unbuilt piece -gets a placeholder and a backlog anchor, not silence or a fabricated -mechanism. Both of those tasks' own `docs/planning/roadmap.md` entries -carry a matching note asking for `sequences.md` to be updated in the same -change that lands them -- check both notes still agree with reality -whenever either task is touched, the same "a diagram makes claims too" -discipline this directory already applies to `overview.md`'s system -diagram (below). +**One of its four sections still carries a `Planned` subsection for a +mechanism that doesn't exist yet** (checkpointing simulation state, +Section 3) -- anchored to the specific roadmap task that will build it +(TASK-034) rather than left open-ended, per the maintainer's direction +that an unbuilt piece gets a placeholder and a backlog anchor, not +silence or a fabricated mechanism. That task's own `docs/planning/ +roadmap.md` entry carries a matching note asking for `sequences.md` to +be updated in the same change that lands it -- check the note still +agrees with reality whenever that task is touched, the same "a diagram +makes claims too" discipline this directory already applies to +`overview.md`'s system diagram (below). + +**Section 2's own `Planned` subsection (driving `simulation.step()` from +a live render loop) was replaced with the real, built sequence 2026-08-28 +(TASK-030)** -- `bootstrap.py`'s `_add_passive_scalar_transport`, the +first config to wire a real timestepping loop into an actual `pyflow +run`. Do not assume it is still a placeholder from an older reading of +this note. `icds.md` (KA-030, Interface Contract Definitions -- the user/configuration-facing interfaces PyFlow's components expose, *not* diff --git a/docs/architecture/engine.md b/docs/architecture/engine.md index 10d7c10..34abf8e 100644 --- a/docs/architecture/engine.md +++ b/docs/architecture/engine.md @@ -305,15 +305,21 @@ condition type. schemes, `src/pyflow/engine/numerics/boundary_condition.py` (`BoundaryCondition`, TASK-019 Boundary Condition Interface, Stage 3; `DirichletBoundaryCondition`, TASK-028; `NeumannBoundaryCondition`, -TASK-029; both Stage 4, 2026-08-28). Periodic -- TASK-030 Periodic -Boundary (Stage 4), not yet built. The interface covers only the +TASK-029; both Stage 4, 2026-08-28). The interface covers only the Dirichlet/Neumann shapes; periodic fits neither and is deliberately not -modelled (see that module's own docstring). +modelled there (see that module's own docstring). `src/pyflow/engine/numerics/assembly.py` registers `DirichletBoundaryCondition`/`NeumannBoundaryCondition` under `"dirichlet"`/`"neumann"` -- the last two of the six `adr/ADR-003` components to go real, retiring the module's final `_Null*` reference -implementation. +implementation. **Periodic -- TASK-030 Periodic Boundary (Stage 4), +built 2026-08-28.** Not a `BoundaryCondition` implementation at all: a +periodic face is mesh geometry, not a prescribed value, so +`StructuredCartesianMesh.wrapped_neighbour_cell(face) -> int` +(`src/pyflow/engine/mesh.py`) is the real mechanism, and +`assemble_numerics` threads a second mapping (`periodic_pairs`) into the +advection/diffusion factories alongside `boundary_conditions` for a +concrete scheme to consult at a periodic face instead. **Upgrade path:** basic edge boundaries → mixed conditions → internal boundaries → arbitrary surfaces/geometries (`upgrade-paths.md` "Boundary diff --git a/docs/architecture/icds.md b/docs/architecture/icds.md index eb8d390..bce5686 100644 --- a/docs/architecture/icds.md +++ b/docs/architecture/icds.md @@ -362,6 +362,21 @@ Completion Criterion 1's carve-out for good**: all six `adr/ADR-003` components now have a real concrete scheme, and zero `_Null*` reference implementations remain in `assembly.py`. +**Done, TASK-030, 2026-08-28, periodic's own half -- Stage 4's last +task.** Not a `BoundaryCondition` implementation, as this document's own +"Expected behaviour" already anticipated (a "wrapped-neighbour reference" +is not a value or a gradient): `StructuredCartesianMesh. +wrapped_neighbour_cell(face) -> int` (`src/pyflow/engine/mesh.py`) is +the real mechanism, additive and off the abstract `Mesh` interface, since +"the opposite edge of the domain" has no meaning for a mesh with no +`(i, j)` structure. `assemble_numerics` still resolves no `Boundary +Condition` instance for a periodic face -- it now also builds a second, +separate mapping (`periodic_pairs`, `{face_name: opposite_face_name}`) +threaded into the advection/diffusion factories alongside +`boundary_conditions`, which is what a concrete scheme consults at a +periodic face instead. `mvp.md`'s "Periodic (where practical)" bullet is +now real, not aspirational. + --- ## Not Yet Addressed: Plugin / Component Discovery diff --git a/docs/architecture/sequences.md b/docs/architecture/sequences.md index 8d15da6..e7ec9c2 100644 --- a/docs/architecture/sequences.md +++ b/docs/architecture/sequences.md @@ -131,26 +131,58 @@ strategies.md` components `assemble_numerics` resolved. See `engine.md` for why that's true and what it buys; this document only shows that it's true, in sequence. -### Planned: driving `step()` from a live run - -**Not built yet.** Nothing in a real `pyflow run` calls `step()` today -- -`bootstrap()` assembles the six numerics components and hands them to -`window.assembled_numerics` for inspection, but never advances a field -over time. `RenderWindow.run(on_frame=...)` (`rendering/window.py`) is -the existing seam a future run will use: it already calls a caller-supplied -callback once per rendered frame, and its own docstring names this -exact purpose ("exactly what a future real-time simulation loop will -need"). What is missing is a caller that closes `on_frame` over a mutable -`fields` mapping and calls `simulation.step(fields, velocity, numerics, -dt)` inside it, then feeds the result back into whatever the render loop -draws. - -This is **TASK-030**'s obligation, not a gap left open-ended: Stage 4's -own Completion Criterion 1 states directly that "TASK-030's golden demo -cannot be assembled at all" without a real simulation-stepping mechanism -running (`docs/planning/roadmap.md`). Update this subsection with the -real sequence once TASK-030 lands -- a note on that task's own roadmap -entry asks for the same thing in the same change. +### Built today: driving `step()` from a live run + +**Built 2026-08-28, TASK-030 -- Stage 4's own Passive Scalar Transport +golden demo.** `bootstrap.py`'s `_add_passive_scalar_transport` is the +caller Section 2's earlier "Planned" note above described in advance: it +builds the transported field and prescribed velocity from +`config.simulation`, renders the first frame, and returns a closure that +`RenderWindow.run(on_frame=...)` (`rendering/window.py`) calls once per +rendered frame thereafter -- the exact seam that subsection's own +docstring anticipated ("exactly what a future real-time simulation loop +will need"). + +```mermaid +sequenceDiagram + participant bootstrap as bootstrap() + participant Add as _add_passive_scalar_transport + participant Window as RenderWindow.run() + participant Advance as on_frame closure + participant Step as simulation.step() + participant Viz as field_visualization + + bootstrap->>Add: _add_passive_scalar_transport(window, mesh, config) + Add->>Viz: build_scalar_field_mesh(scalar_field, colors) + Viz-->>Add: gfx.Mesh (frame 0) + Add->>Window: scene.add(rendered_object) + Add-->>bootstrap: on_frame closure + bootstrap->>Window: window.run(max_frames, on_frame) + loop each rendered frame + Window->>Advance: on_frame() + Advance->>Step: step(state, velocity, numerics, dt) + Step-->>Advance: new state + Advance->>Window: simulation_fields = new state + Advance->>Viz: scalar_field_colors(new tracer, low, high, range) + Viz-->>Advance: per-cell RGBA colors + Advance->>Window: scene.remove(old object); scene.add(new object) + end +``` + +**Each frame's rendered `gfx.Mesh` is rebuilt from scratch, not mutated +in place.** `build_scalar_field_mesh`/`scalar_field_colors` are already +proven correct (TASK-017); an in-place colour-buffer mutation +(`geometry.colors.data[:] = ...`) would be new, unverified pygfx-API +surface for a small win on a small demo mesh -- a deliberate, recorded +deferral (`docs/planning/roadmap.md` TASK-030's own Design decision), not +an oversight. `gfx.Scene.remove` was verified directly (add then remove +leaves `len(scene.children) == 0`) before being relied on. + +`RenderWindow.simulation_fields` is what lets a caller -- the golden +demo's own regression test, most directly -- read back the real field +state a rendered frame came from, the same "`bootstrap()` populates it, +`RenderWindow` itself holds no simulation content" shape +`assembled_numerics` already established (Section 3, below). --- diff --git a/docs/implementation/golden-demos.md b/docs/implementation/golden-demos.md index 9e3a33b..79c3e53 100644 --- a/docs/implementation/golden-demos.md +++ b/docs/implementation/golden-demos.md @@ -198,7 +198,64 @@ implementation (`engine/numerics/assembly.py`'s own docstring explains why one exists under `src/` at all: a real CLI subprocess needs *something* to assemble into). This demo proves the assembly mechanism works, not that PyFlow computes anything yet; the first demo that -computes real physics is Scalar Transport (Stage 4). +computes real physics is Passive Scalar Transport, below (Stage 4). + +## Passive Scalar Transport + +TASK-030's own golden demo (`docs/planning/roadmap.md`, Stage 4 +Completion Criterion 1) -- PyFlow's first demo that computes real +physics, and the first `pyflow run` that steps a real simulation +forward *live*, one timestep per rendered frame, rather than rendering +one static picture. Called "Scalar Transport" in +`docs/planning/implementation-plan.md`'s own Golden Demos table and +`planning/data/demos.yaml` (`demo-scalar-transport`) -- the same demo, +named there before it was built; `mvp.md`'s own validation bullet and +this task's own roadmap entry both say "Passive scalar transport", +which is the name used here and for the real artifact +(`passive_scalar_transport.yaml`). Not reconciled across those two +documents in this change -- noted here rather than left silent, per the +Blast Radius rule's "if something in the radius cannot be updated now, +say so explicitly." + +"Working" means, concretely: + +- the demo *is* `examples/golden-demos/passive_scalar_transport.yaml` -- + a `simulation` section naming a `gaussian_blob` initial condition and + a `uniform` prescribed velocity (`SimulationConfig`, + `src/pyflow/configuration/schema.py`), a `numerics` section whose + east/west boundaries are `periodic` and north/south are `neumann` + (zero gradient -- the prescribed velocity is purely horizontal, so + nothing crosses them, but diffusion still needs some condition there + regardless of flow direction), run via + `uv run python -m pyflow run --config examples/golden-demos/passive_scalar_transport.yaml`; +- a real `simulation.step()` call advances the field once per rendered + frame (`src/pyflow/bootstrap.py`'s `_add_passive_scalar_transport`, + wired through `RenderWindow.run(on_frame=...)`) -- checked directly, + not only by pixel-diffing: the field's own mass-weighted centroid + moves downstream by approximately the prescribed velocity times the + elapsed real time, measured across two independent real runs at + different frame counts (`tests/golden/test_passive_scalar_transport.py`), + within a tolerance derived from an actual measured run (~4% agreement), + not guessed; +- the periodic wrap is exercised by this live run, not only proven in + isolation: rendered offscreen at increasing frame counts, the blob is + seen translating downstream and, once total elapsed travel reaches one + full domain width, reappearing spread across both the east and west + edges -- verified visually during this task's own build, not asserted + by the regression test above (which checks displacement over a + shorter interval, before any wrap occurs, per its own module + docstring); the wrap's own correctness claim belongs to + `periodic_boundary.feature`, checked in isolation as a convergence + property (see that task's own Design decisions, + `docs/planning/roadmap.md` TASK-030); +- it runs headlessly via `--backend offscreen`, same as every other + demo. + +**The velocity field is prescribed, not solved.** Stage 5 is what +eventually solves Navier-Stokes for real (`PressureCoupling`, real +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 diff --git a/docs/planning/backlog.md b/docs/planning/backlog.md index dd9cefd..bcb0731 100644 --- a/docs/planning/backlog.md +++ b/docs/planning/backlog.md @@ -1920,7 +1920,7 @@ here.): - **Advection scheme** -- **done, TASK-023 (Stage 4, 2026-08-27)**: total transported quantity is conserved on a closed domain (every boundary cell's velocity exactly zero, interior cells nonzero -- - no periodic boundary exists yet to test the alternative reading), + no periodic boundary existed yet to test the alternative reading), summing the field over every cell before and after many timesteps to floating-point tolerance -- `tests/features/ first_order_upwind_advection.feature`'s own "Conservation on a @@ -1928,7 +1928,16 @@ here.): own Advection bullet in `docs/planning/roadmap.md`, renumbered 2026-08-26 when TASK-040/Simulation Orchestrator was added as this Stage's own new Criterion 1; discharged when TASK-023 - landed rather than staying a backlog note.) + landed rather than staying a backlog note.) **A periodic boundary + now exists (TASK-030, Stage 4, 2026-08-28) -- the "alternative + reading" this note flagged is still not itself checked**: + `periodic_boundary.feature`'s own round-trip scenario measures + convergence of a field's *shape* toward its starting distribution + under mesh refinement, not total transported quantity summed + before/after under a periodic wrap specifically. A genuinely + distinct, still-open claim, found while closing this note rather + than left implied-done by periodic boundaries existing at all -- + not added to this item's own scope without a decision to do so. - **Diffusion scheme** -- **done, TASK-024 (Stage 4, 2026-08-27)**: same conservation check under zero-flux (Neumann) boundaries as Advection's above -- an insulated domain's field total is diff --git a/docs/planning/roadmap.md b/docs/planning/roadmap.md index d342598..0759b84 100644 --- a/docs/planning/roadmap.md +++ b/docs/planning/roadmap.md @@ -189,7 +189,7 @@ 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): **588 tests at 99% as of 2026-08-28**, having been 64 when +(C1a/C1b): **603 tests at 99% as of 2026-08-28**, 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 @@ -308,12 +308,70 @@ distinct from `scalar_gradient` (a fixture edit, not a new test) -- found necessary by mutation testing: the previous fixture's shared `0.0` for both fields meant a regression reading `velocity` instead of `scalar_gradient` would have passed unnoticed. +603 after TASK-030 (Periodic Boundary, Stage 4's ninth and last task): +four new tests in `test_structured_cartesian_mesh.py` for +`wrapped_neighbour_cell` (correct pairing on all four edges, split into +one rejection test per face orientation -- vertical and horizontal -- +rather than one, found necessary while confirming coverage: a single +`next(...)`-selected interior face always picked the vertical branch +first, leaving the horizontal `raise` line untested); three new Gherkin +scenarios in `tests/unit/test_periodic_boundary.py` (`periodic_boundary. +feature`, this task's own Acceptance Criteria -- advection reading the +wrapped neighbour, diffusion computing the gradient at one full cell +width, and a convergence-based round-trip invariant, below); two new +capture tests in `test_assembly.py` for `periodic_pairs` threading into +the advection/diffusion factories, the periodic analogue of TASK-040's/ +TASK-024's own boundary-conditions capture tests; and four new +`SimulationConfig` load/reject tests in `test_configuration.py` +(reads the section, rejects an unknown `scalar_pattern`/`velocity_ +pattern`, rejects a non-numeric `velocity`), the same shape every prior +config-section addition in this run used. Two new Gherkin scenarios in +`tests/golden/test_passive_scalar_transport.py` (`passive_scalar_ +transport.feature`) for the golden demo this task also builds (Stage 4 +Completion Criterion 1's own "a real simulation-stepping mechanism +running live", not this task's own Acceptance Criteria -- see this +task's own Intent above). +15 overall. + +**The round-trip scenario's own criterion needed a genuine correctness +finding to get right, not just "matches exactly", found by running the +real numbers before writing the assertion (the same discipline TASK-026's +null-space check and TASK-027's Rhie-Chow investigation used).** A plain +"advect a field once around a periodic domain, assert it returns to its +starting distribution" fails even for a *correct* wrap: first-order +upwind's own O(dx) numerical diffusion smooths any field over the +distance it travels, and refining the timestep alone does not shrink +that error (the RK4-integrated semi-discrete system converges to a fixed, +spatially-truncation-dominated limit as `dt -> 0`, verified directly: +identical results at `num_steps` ranging 10-160 on a fixed mesh). What a +*wrong* wrap (mirrored or clamped to the owner's own cell at the +periodic boundary, instead of the opposite edge) cannot reproduce is +refinement actually closing the gap -- measured directly on a genuinely +periodic-compatible fixture (a sine wave in x, period equal to the +domain width, plus a `100*y` term so every row is checked -- a plain +linear ramp was tried first and rejected: wrapping it creates an +artificial discontinuity at the seam that numerical diffusion then +smooths, confounding this specific claim), a real wrap's own round-trip +error drops by roughly 62% over a 4x mesh refinement (16 to 64 cells), +while a mirrored/clamped one (a throwaway mutation, built and run +specifically to check this) drops by only roughly 16% and stays several +times larger throughout. The scenario asserts the fine-resolution error +stays under two thirds of the coarse one -- comfortably separates the two +outcomes measured, confirmed to actually fail under the mirrored/clamped +mutation before being trusted. +**The demo's own centroid-displacement tolerance was measured the same +way, not guessed**: a real run of `passive_scalar_transport.yaml` agrees +with the closed-form `velocity * dt * steps` prediction to within ~4%; +the scenario's own `rel=0.15` bound stays comfortably above that margin +without being so loose a genuinely broken stepping loop (frozen, +backwards, or off by a large factor) could still pass -- confirmed +directly: a mutation that froze `state` every frame (never actually +calling `simulation.step`) fails this scenario. The rest of the climb to 508 that same day is `tests/unit/test_generate_status_report.py` -- the new tool's own test suite -- growing from 23 to 35 tests as that tool itself grew, a live demonstration that this count moves for reasons having nothing to do with the fluid solver and everything to do with why it needs checking -rather than re-reading. **48 of those 588 are Gherkin scenarios +rather than re-reading. **53 of those 603 are Gherkin scenarios rather than pytest functions** (`adr/ADR-007-executable-acceptance-criteria.md`; up from fourteen with `field_display.feature` gaining scenarios and `numerics_assembly.feature` @@ -344,7 +402,18 @@ own Intent; and to 48 with TASK-029's own `neumann_boundary.feature`, two scenarios each building a real interior scheme together with a real `NeumannBoundaryCondition`, both prescribing a nonzero gradient throughout -- diffusion's own reads the gradient's numeric value -directly, advection's own proves the opposite, that it is never read). +directly, advection's own proves the opposite, that it is never read; +and to 51 with TASK-030's own `periodic_boundary.feature`: advection +reading the wrapped neighbour, diffusion computing the gradient at one +full cell width, and a round-trip invariant measured as convergence +under mesh refinement rather than exact equality at one resolution +(this task's own Design decisions record the numerical finding that +makes "matches exactly" the wrong claim to check); and to 53 with +TASK-030's own golden demo, `passive_scalar_transport.feature` -- the +required CLI-subprocess scenario every demo carries, and a quantitative +claim that the transported field's own mass-weighted centroid moves +downstream at approximately the prescribed velocity over real elapsed +time, not only that rendered pixels changed). **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`, @@ -3771,6 +3840,30 @@ were first sketched. Numbered out of sequence rather than renumbered into the 023-030 run, and built first regardless, per the Build order note under the Discharge map above. +### Status as of 2026-08-28: Stage 4 complete, ten of ten criteria met + +| Criterion | Verdict | +|-----------|---------| +| 1. Simulation-stepping mechanism exists, face-flux accumulation uniform across every face | **Met** (TASK-040). `engine/simulation.py`'s `step()`/`accumulate_flux_to_cells` are real, unit-tested (`tests/unit/test_simulation.py`, `simulation_orchestrator.feature`), and never branch on `Mesh.is_boundary_face` -- a concrete Advection/Diffusion scheme handles a boundary face itself. | +| 2. Real implementation replaces reference, under the existing MVP name | **Met.** All seven names TASK-023..029 own resolve to a real scheme (`FirstOrderUpwindAdvection`, `CentralDifferenceDiffusion`, `RK4Integrator`, `ConjugateGradientSolver`, `PISO`, `DirichletBoundaryCondition`, `NeumannBoundaryCondition`), each checked by `isinstance` in `tests/unit/numerics/test_assembly.py`, not just by the name still validating. | +| 3. Contract suite still holds, shown insufficient alone | **Met.** Every real scheme joined its own interface's contract suite (`test_advection_contract.py` etc.) with no edit to any existing test body except where a real interface widening required one (`TimeIntegrator.advance`, `PressureCoupling.correct`, each its own recorded ADR); each scheme's own `.feature` file is what actually proves physical correctness, per this criterion's own "necessary and explicitly not sufficient". | +| 4. Physical correctness, per task | **Met.** Each of TASK-023..030's own Intent lines is discharged by that task's own `.feature` file -- TASK-030's own (the round-trip invariant) is checked as convergence under mesh refinement rather than exact equality at one resolution, a genuine numerical finding recorded in that task's own Design decisions, not a weaker check chosen for convenience. | +| 5. Real implementations' own rejection paths tested | **Met.** Every `UnconfiguredBoundaryFaceError`/`IncompatibleVelocityFieldError`/`IncompatibleVectorFieldError`/`NotABoundaryFaceError` (TASK-030's own, on `wrapped_neighbour_cell`) is exercised directly against real bad input, not only inherited-untested from a shared helper -- `docs/practices.md`'s "rejection criteria stop at the constructor" checked task by task. | +| 6. Executable Gherkin criteria, `make check-scenarios` gates | **Met.** `make check-scenarios`: "All 53 scenario(s) across 14 feature file(s) are bound and run," verified directly, not assumed from the file count. | +| 7. No `_Null*` registration survives under an implemented name | **Met, closed at TASK-029, unaffected by TASK-030** (which retires one more genuinely-dead helper, `_resolve_with_argument`, but no `_Null*` class -- there were none left). `assembly.py`'s own registration calls at the bottom of the file name only real classes. | +| 8. Demonstration: Passive Scalar Transport | **Met** (TASK-030). `examples/golden-demos/passive_scalar_transport.yaml`, run via the real CLI; `tests/golden/test_passive_scalar_transport.py`'s own quantitative scenario (mass-weighted centroid displacement, tolerance measured from a real run); verified visually beyond the regression test -- rendered offscreen at increasing frame counts, the blob is seen translating and, by one full domain width of travel, wrapping around the periodic boundary. | +| 9. `make ci` green on a real runner | **Met.** PR #38 (`feat/task-030-periodic-boundary`), run 33159480722: `ci (ubuntu-latest)` green in 2m57s, `ci (windows-latest)` green in 5m30s -- checked against the actual run via `gh pr checks --watch`, not inferred from the PR merging. | +| 10. Documentation matches the tree | **Met.** `make check-references`/`check-manifest`/`check-inventory`/`check-dependency-tree`/`check-docs`/`check-docs-index`/`check-graph` all pass against the tree as this task leaves it; every stale forward-reference this sweep found (two in `src/pyflow/engine/CLAUDE.md` describing periodic as still raising `UnconfiguredBoundaryFaceError`, one in `docs/planning/backlog.md` saying "no periodic boundary exists yet") was corrected in this same change, not left for a future exit audit to find. | + +**Criterion 9 could not be marked Met from a local checkout alone when +this table was first drafted** -- a local `make ci` pass is not the same +claim as a real CI run, per this repository's own standing distrust of a +green-CI claim nothing has actually watched (`docs/practices.md`, +`CLAUDE.md`'s Merge Gate), the same reasoning Stage 3's own table +applied. Left honestly pending until PR #38's own run actually +completed and was checked directly (above), rather than assumed the +moment the branch was pushed. + ## TASK-040 Simulation Orchestrator @@ -5291,6 +5384,8 @@ contract suite once `NeumannBoundaryCondition` joined it). Periodic Boundary +**Status:** Done, 2026-08-28, Stage 4's ninth and last task. + **Intent:** the claim is that a field advected once around a periodic domain returns to its starting distribution -- a round-trip invariant, which is the only check that distinguishes a genuine wrapped-neighbour @@ -5394,6 +5489,149 @@ from a live run" subsection with the real, built sequence, in the same change -- that document names this task as the anchor for exactly this, and should not still say "not built yet" once it is. +**Done, 2026-08-28, in the same change as the rest of this task.** +`docs/architecture/sequences.md` Section 2 now describes the real, +built sequence -- see that document's own updated section. The demo +itself is `examples/golden-demos/passive_scalar_transport.yaml`: a +prescribed uniform velocity carries a Gaussian scalar blob across a +mesh whose east/west edges are periodic (north/south are `neumann` with +a zero gradient -- an insulated wall, since the prescribed velocity is +purely horizontal and never crosses them). Verified visually, not only +by its own regression test: rendered offscreen at several frame counts, +the blob is seen translating downstream and, by the frame count +corresponding to one full domain width of travel, reappearing spread +across both the east and west edges -- the periodic wrap, genuinely +exercised by a live run, not only by `periodic_boundary.feature`'s own +isolated checks. + +### Artifacts Produced + +- `src/pyflow/engine/mesh.py`: `NotABoundaryFaceError`; + `StructuredCartesianMesh.wrapped_neighbour_cell(face) -> int` -- the + periodic wrap's own geometry (Design decision above). +- `src/pyflow/engine/numerics/advection.py`, + `src/pyflow/engine/numerics/diffusion.py`: both concrete schemes gain a + `periodic_pairs: Mapping[str, str]` constructor parameter and consult + it before falling through their existing interior-neighbour formula -- + no interface change, no new ADR (mirrors `diffusion_coefficient`'s own + TASK-024 precedent). +- `src/pyflow/engine/numerics/assembly.py`: `_PAIRED_BOUNDARY`; + `periodic_pairs` built alongside `boundary_conditions` in + `assemble_numerics`'s existing per-face loop; `register_advection_ + scheme`/`register_diffusion_scheme`'s factory types widen accordingly; + `_resolve_with_three_arguments` (diffusion's own factory now takes + three constructor arguments); `_resolve_with_argument` (the one-argument + helper) deleted as genuinely dead code, its only caller (advection) + having moved to `_resolve_with_two_arguments`. +- `src/pyflow/configuration/schema.py`: `SimulationConfig` + (`scalar_pattern`, `velocity_pattern`, `velocity`) -- the Passive + Scalar Transport demo's own configuration surface, distinct from + `FieldDisplayConfig` (static single-frame display) since this seeds a + real, repeatedly-stepped run. +- `src/pyflow/bootstrap.py`: `_simulation_scalar_initializer`/ + `_simulation_velocity_initializer`; `_add_passive_scalar_transport`, + wiring a real `simulation.step()` call into `RenderWindow.run( + on_frame=...)` -- the first config in this project's history to do so. +- `src/pyflow/rendering/window.py`: `RenderWindow.simulation_fields` -- + the live simulation's own field state, read back by the golden demo's + regression test the same way `assembled_numerics` already lets a + caller read back what got assembled. +- `tests/features/periodic_boundary.feature`, + `tests/features/passive_scalar_transport.feature` -- see Acceptance + Criteria. + +### Implementation + +**The periodic distance is `2 * mesh.face_centroid_distance(face)`, not +a new geometry accessor** -- verified numerically before relying on it: +on a mesh with distinct `dx`/`dy` and a non-trivial origin, doubling a +boundary face's own owner-to-face distance reproduces the true uniform +grid spacing exactly, matching an ordinary interior face's distance on +the same mesh to float precision. The periodic "neighbour" is one full +cell-width away, not the real wrapped cell's actual (far-side-of-the- +domain) centroid, so the plain interior formula `hypot(neighbour.xy - +owner.xy)` would be wildly wrong if applied naively across the wrap. + +`FirstOrderUpwindAdvection.flux`/`CentralDifferenceDiffusion.flux`: at a +boundary face, if its named edge is a key in `periodic_pairs`, +`mesh.wrapped_neighbour_cell(face)` stands in for `neighbour` (and, for +diffusion, the doubled distance stands in for the boundary formula's +own owner-to-face distance) before either scheme's existing +`neighbour is not None` code path runs -- literally the same formula +already used for a real interior neighbour. `boundary_conditions` is +never consulted for a periodic face; `periodic_pairs` and +`boundary_conditions` are deliberately two separate mappings, keeping +"prescribed value/gradient" and "wrapped-neighbour" explicit rather than +inferred from one mapping's absence. + +`bootstrap.py`'s live-stepping branch rebuilds the rendered `gfx.Mesh` +from scratch each frame (`window.scene.remove` the old one, +`build_scalar_field_mesh` a new one) rather than mutating the geometry's +own colour buffer in place -- `build_scalar_field_mesh`/ +`scalar_field_colors` are already proven correct (TASK-017); an in-place +buffer mutation would be new, unverified pygfx-API surface for a small +win on a small demo mesh. `gfx.Scene.remove` was verified to behave as +expected (add then remove leaves `len(scene.children) == 0`) before +being relied on, the same "check implementation details directly" +discipline every prior rendering addition in this codebase has used. + +### Acceptance Criteria + +- **Criterion 4, this task's own claim, per the Intent above:** the + round-trip invariant is checked as *convergence under mesh refinement*, + not exact equality at one resolution -- found necessary by running the + real numbers first: first-order upwind's own O(dx) numerical diffusion + smooths any field over the distance it travels regardless of whether + the wrap itself is correct, and refining the timestep alone does not + shrink that error (verified: near-identical results at `num_steps` + 10-160 on a fixed mesh, since the RK4-integrated semi-discrete system + converges to a fixed, spatially-truncation-dominated limit as + `dt -> 0`). A real wrap's own round-trip error drops by roughly 62% + over a 4x mesh refinement; a mirrored/clamped one (built and run as a + throwaway mutation specifically to check this) drops by only roughly + 16% and stays several times larger throughout -- `periodic_boundary. + feature`'s own scenario asserts the fine-resolution error stays under + two thirds of the coarse one, confirmed to actually fail under that + same mutation before being trusted. +- Advection/diffusion's own wiring is checked directly (not only via the + round-trip): a real periodic pairing reads the wrapped neighbour's own + value/gradient at the correct distance, verified by mutation (reverting + the wrap to the owner's own cell, or the distance to the un-doubled + one, both fail the corresponding scenario). +- **An accessor-level rejection criterion, per the Design decision's own + "still owed" note:** `wrapped_neighbour_cell` raises `NotABoundaryFaceError` + on an interior face, exercised directly for both a vertical and a + horizontal interior face (`test_structured_cartesian_mesh.py`) -- found + necessary while confirming coverage, since a single `next(...)`-selected + interior face always picked the vertical branch first, leaving the + horizontal `raise` line genuinely untested by one test alone. +- Config-surface correctness: `test_configuration.py` gained + `SimulationConfig` load/reject tests, the same shape every prior + config-section addition in this run used. +- **Stage 4 Completion Criterion 1, the golden demo:** `passive_scalar_ + transport.feature`'s own quantitative scenario checks the transported + field's mass-weighted centroid moves downstream by approximately + `velocity * dt * steps`, measured over two independent real runs + (`bootstrap()`, offscreen) rather than only by pixel-diffing two + frames -- "physical fields evolve" (`docs/implementation/mvp.md`'s + Definition of Done) checked directly. The tolerance (`rel=0.15`) was + measured from a real run (~4% actual agreement), not guessed, and + confirmed to fail under a mutation that froze the simulation state + every frame. +- **Stage 3 Completion Criterion 1's carve-out**, already closed at + TASK-029, stays closed -- this task adds no new `_Null*` class, and + deletes one more genuinely dead helper (`_resolve_with_argument`). + +### Discharges + +Stage 4 Completion Criteria 1 (the golden demo, jointly with TASK-040's +own orchestration mechanism), 3, 4 and 5's rejection-path share +(`NotABoundaryFaceError`), Periodic Boundary's own share; Stage 4 +Completion Criterion 6, this task's own two feature files +(`periodic_boundary.feature`, `passive_scalar_transport.feature`), both +bound (`make check-scenarios`); Stage 4 Completion Criterion 7, restated +closed (already true since TASK-029, unaffected by this task). + --- # Stage 5 — First Fluid Solver diff --git a/docs/planning/status.md b/docs/planning/status.md index adb4079..6e8143d 100644 --- a/docs/planning/status.md +++ b/docs/planning/status.md @@ -31,16 +31,17 @@ pie showData - **Stage 1 -- Representing Space** complete (2026-08-21) - **Stage 2 -- Representing Fields** complete (2026-08-22) - **Stage 3 -- Numerical Engine** complete (2026-08-23) +- **Stage 4 -- First Numerical Methods** complete (2026-08-28) ### Up next -**Stage 4 -- First Numerical Methods** is next, starting with TASK-027 (PISO Pressure Coupling), 3 more not yet started in this stage. +**Stage 5 -- First Fluid Solver** is next, starting with TASK-031 (Velocity Field Support), 3 more not yet started in this stage. ## Live repository facts - **45** `CLAUDE.md` files -- **588** tests collected -- **48** Gherkin scenarios (`tests/features/*.feature`) +- **603** tests collected +- **53** Gherkin scenarios (`tests/features/*.feature`) ## Stages @@ -98,7 +99,7 @@ pie showData ### Stage 4 -- First Numerical Methods -**no status recorded** -- `██████░░░░` 5/9 tasks; 10 criteria defined, no status line yet +**complete, as of 2026-08-28** -- `██████░░░░` 5/9 tasks; 10/10 criteria met | Task | Status | Date | Artifact | |------|--------|------|----------| diff --git a/docs/repository-inventory.md b/docs/repository-inventory.md index 30e2e0b..743d8fb 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. -**261 tracked files** across 45 directories; +**266 tracked files** across 45 directories; 4 are empty. ## (root) @@ -176,6 +176,7 @@ listing files. - `empty_window.yaml` - `field_display.yaml` - `numerics_assembly.yaml` +- `passive_scalar_transport.yaml` ## examples/tutorials @@ -308,6 +309,8 @@ listing files. - `first_order_upwind_advection.feature` - `neumann_boundary.feature` - `numerics_assembly.feature` +- `passive_scalar_transport.feature` +- `periodic_boundary.feature` - `piso_pressure_coupling.feature` - `rk4_time_integration.feature` - `simulation_orchestrator.feature` @@ -322,6 +325,7 @@ listing files. - `test_empty_window.py` - `test_field_display.py` - `test_numerics_assembly.py` +- `test_passive_scalar_transport.py` ## tests/integration @@ -369,6 +373,7 @@ listing files. - `test_mesh_contract.py` - `test_mesh_visualization.py` - `test_neumann_boundary.py` +- `test_periodic_boundary.py` - `test_piso_pressure_coupling.py` - `test_rendering.py` - `test_rk4_time_integration.py` diff --git a/docs/repository-manifest.md b/docs/repository-manifest.md index 1e6ec49..736d789 100644 --- a/docs/repository-manifest.md +++ b/docs/repository-manifest.md @@ -527,6 +527,34 @@ Dirichlet-side sibling did. With this field's addition every one of the six `adr/ADR-003` components has a real concrete scheme; `assembly.py`'s own `_null_boundary_value` helper, which only these last two retired classes ever called, is deleted alongside them. +**Periodic Boundary (TASK-030, Stage 4's ninth and last task) closes the +final gap without adding a seventh `BoundaryCondition` shape**: a +periodic face is mesh geometry, not a prescribed value, so +`StructuredCartesianMesh.wrapped_neighbour_cell(face) -> int` +(`mesh.py`, additive, off the abstract `Mesh` interface -- the same +precedent `cell_id`/`cell_index`/`boundary_face_name` set) is the real +mechanism. `FirstOrderUpwindAdvection`/`CentralDifferenceDiffusion` both +gain a `periodic_pairs: Mapping[str, str]` constructor parameter +(registry-level widening, no interface change, no new ADR -- mirrors +`diffusion_coefficient`'s own TASK-024 precedent) and consult it before +falling through their existing interior-neighbour formula; diffusion's +own periodic distance is `2 * mesh.face_centroid_distance(face)`, +verified numerically to reproduce the true uniform grid spacing exactly. +`assemble_numerics` builds `periodic_pairs` in the same per-face loop +that already builds `boundary_conditions`; `register_diffusion_scheme` +widens to a three-argument factory, needing a new +`_resolve_with_three_arguments` helper, and the one-argument +`_resolve_with_argument` helper (advection's own former route) is +deleted as genuinely dead code. A new `SimulationConfig` section +(`configuration/schema.py`) seeds the Passive Scalar Transport golden +demo -- a `gaussian_blob` initial condition and a `uniform` prescribed +velocity, distinct from `FieldDisplayConfig`'s own static-display +patterns since this seeds a real, repeatedly-stepped run. `bootstrap.py` +gains `_add_passive_scalar_transport`, the first code in this project's +history to wire a real `simulation.step()` call into +`RenderWindow.run(on_frame=...)`; `RenderWindow` gains +`simulation_fields`, the same "`bootstrap()` populates it, no simulation +content of its own" shape `assembled_numerics` already established. `rendering/` (`canvas.py`, `window.py` -- `RenderWindow.assembled_numerics`, TASK-021's one addition to this package -- `mesh_visualization.py`, `field_visualization.py`), @@ -557,7 +585,7 @@ stage boundary, not only when something here is being edited. `tests/` with `unit/`, `integration/`, `golden/`, `performance/`. -🟨 — 54 test modules, **588 tests, 99% coverage** (2026-08-28; 35 of +🟨 — 56 test modules, **603 tests, 99% coverage** (2026-08-28; 35 of these being `test_generate_status_report.py` itself; the 47th module, `test_simulation.py` (TASK-040), was the first Gherkin feature file bound outside `tests/golden/` -- not a golden demo, so it lives here per @@ -571,8 +599,11 @@ follows it again for Stage 4's fourth, the 51st, Stage 4's fifth, the 52nd, `test_piso_pressure_coupling.py` (TASK-027), follows it again for Stage 4's sixth, the 53rd, `test_dirichlet_boundary.py` (TASK-028), follows it again for Stage 4's -seventh, and the 54th, `test_neumann_boundary.py` (TASK-029), follows it -again for Stage 4's eighth). +seventh, the 54th, `test_neumann_boundary.py` (TASK-029), follows it +again for Stage 4's eighth, the 55th, `test_periodic_boundary.py` +(TASK-030), follows it again for Stage 4's ninth and last -- and the +56th, `tests/golden/test_passive_scalar_transport.py`, binds that same +task's own golden demo, PyFlow's first that computes real physics). The roadmap's own restatement of this count (`docs/planning/roadmap.md`, just above Stage 1) is cross-checked against `pytest --collect-only` by `make check-status`; this row is not machine-checked and needs the same @@ -618,6 +649,11 @@ feature file, `dirichlet_boundary.feature`/`neumann_boundary.feature`, below. Neumann's own join closes Stage 3 Completion Criterion 1's carve-out for good -- every component this interface (and every other `adr/ADR-003` interface) names now has a real concrete scheme. +**`test_advection_contract.py`/`test_diffusion_contract.py` each gained +a second constructor argument at their one existing real-fixture factory +call (TASK-030)** -- not a new fixture, a signature widening: both +schemes' `periodic_pairs` parameter, passed `{}` since neither suite's +own claims are about periodic behaviour. `test_time_integrator_contract.py` also has no separate inert-teeth-check class, for a different reason: its own acceptance criteria (the zero-derivative case and the nonzero scheme-independence case) already @@ -661,8 +697,12 @@ boundary-conditions-evaluate test's Neumann fixture was rewritten again to assert `scalar_gradient`, this time with a `velocity` deliberately distinct from it, a fix mutation testing found necessary; a new real-scheme-resolution test was added against an explicit -`"neumann"`-typed config, since no default face is Neumann) is not a -contract suite -- it's the +`"neumann"`-typed config, since no default face is Neumann -- and again +TASK-030, not a null-to-real swap this time (none remain), but two new +tests proving `assemble_numerics` threads the resolved `periodic_pairs` +mapping into the advection/diffusion factories, the periodic analogue of +the boundary-conditions capture tests TASK-040/TASK-024 already added) +is not a contract suite -- it's the in-process unit suite for `assemble_numerics`/the six registries, covering Stage 3 Completion Criteria 3 (a newly-registered name resolves with no edit under `src/`) and 4 (mutating a `NumericsConfig` after @@ -745,6 +785,28 @@ into the flux formula; advection's own proves the opposite, that the value is never read -- confirmed by mutation, and this task's own last component closes Stage 3 Completion Criterion 1's carve-out: zero `_Null*` reference implementations remain anywhere in `assembly.py`. +`test_periodic_boundary.py` (TASK-030, Stage 4's ninth and last task) +binds `tests/features/periodic_boundary.feature` -- Criterion 4's +Periodic Boundary bullet. A genuinely different shape from every +boundary-condition binding module above it: there is no condition class +under test, since periodic bypasses `BoundaryCondition` entirely; the +round-trip scenario is checked as convergence under mesh refinement +rather than exact equality at one resolution, a genuine numerical +finding (first-order upwind's own O(dx) numerical diffusion smooths any +field over the distance it travels, regardless of whether the wrap is +correct) verified before being trusted, not assumed -- a real wrap's own +round-trip error drops ~62% over a 4x mesh refinement, a mirrored/ +clamped one (a throwaway mutation built specifically to check this) +drops only ~16% and stays several times larger throughout. +`tests/golden/test_passive_scalar_transport.py` (TASK-030) binds +`tests/features/passive_scalar_transport.feature` -- Stage 4 Completion +Criterion 1's own golden demo, PyFlow's first that computes real +physics: the required CLI-subprocess scenario every demo carries, plus a +quantitative scenario checking the transported field's mass-weighted +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. `unit/` otherwise holds config/logging/rendering (D1/D2/D3), the tooling tests @@ -830,9 +892,10 @@ 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), and -`golden-demos/numerics_assembly.yaml` (TASK-021, 2026-08-23) are the four -demos so far: plain configuration files, no Python -- golden demos run +`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 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 9f65492..547f65e 100644 --- a/examples/golden-demos/CLAUDE.md +++ b/examples/golden-demos/CLAUDE.md @@ -17,7 +17,7 @@ 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. -Four demos live here as of 2026-08-23, one per stage that has produced a +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: @@ -38,8 +38,16 @@ render: output" carve-out means this demo proves configuration assembles into real (if physically trivial) instances, not that anything new appears on screen. +- `passive_scalar_transport.yaml` (TASK-030, 2026-08-28), Stage 4's -- + PyFlow's first demo that computes real physics, and the first + `pyflow run` that steps a real simulation forward live rather than + rendering one static frame: a `simulation` section (`gaussian_blob` + initial condition, `uniform` prescribed velocity) plus a `numerics` + 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 three, waiting on +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. diff --git a/examples/golden-demos/passive_scalar_transport.yaml b/examples/golden-demos/passive_scalar_transport.yaml new file mode 100644 index 0000000..44f9dde --- /dev/null +++ b/examples/golden-demos/passive_scalar_transport.yaml @@ -0,0 +1,51 @@ +# Passive Scalar Transport golden demo (Stage 4, TASK-030) -- the first +# `pyflow run` that steps a real simulation forward, live, rather than +# rendering one static frame. A prescribed (not solved) uniform velocity +# field carries a scalar blob across a domain whose east and west edges +# are periodic: the blob translates and wraps around, continuing on the +# other side. +# +# Run it exactly the way any user would: +# +# uv run python -m pyflow run --config examples/golden-demos/passive_scalar_transport.yaml +# +# `simulation:` names the blob's own initial condition and the +# prescribed velocity that transports it (`SimulationConfig`, +# `src/pyflow/configuration/schema.py`); `field_display:` colours the +# live field the same way it colours a static one. `north`/`south` are +# `neumann` with a zero gradient (an insulated wall) rather than +# periodic -- the prescribed velocity is purely horizontal, so nothing +# ever flows across them, but diffusion still needs *some* condition +# there regardless of flow direction. + +mesh: + extent: [20, 6] + spacing: [0.25, 0.25] + +numerics: + timestep: 0.02 + diffusion_coefficient: 0.02 + boundary_conditions: + east: + type: periodic + west: + type: periodic + north: + type: neumann + scalar_gradient: 0.0 + south: + type: neumann + scalar_gradient: 0.0 + +simulation: + scalar_pattern: gaussian_blob + velocity_pattern: uniform + velocity: [1.0, 0.0] + +field_display: + low_color: "#0a0a2a" + high_color: "#ff8c00" + value_range: [0.0, 1.0] + +rendering: + background_color: "#1a1a2e" diff --git a/src/pyflow/bootstrap.py b/src/pyflow/bootstrap.py index 11d00fd..efd42c9 100644 --- a/src/pyflow/bootstrap.py +++ b/src/pyflow/bootstrap.py @@ -29,11 +29,13 @@ from pyflow import __version__ from pyflow.configuration import load_config -from pyflow.configuration.schema import FieldDisplayConfig, RenderBackend +from pyflow.configuration.schema import FieldDisplayConfig, PyFlowConfig, RenderBackend +from pyflow.engine.field import Field from pyflow.engine.logging_setup import configure_logging, get_logger 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 step as simulation_step from pyflow.engine.vector_field import VectorField from pyflow.rendering import RenderWindow from pyflow.rendering.field_visualization import ( @@ -88,6 +90,105 @@ def _vector_display_initializer( raise ValueError(f"unknown vector display pattern: {pattern!r}") # pragma: no cover +def _simulation_scalar_initializer( + pattern: str, bounds: _Bounds +) -> Callable[[float, float], float]: + """A `Field`-style `(x, y) -> value` callable for `SimulationConfig. + 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. + """ + if pattern == "gaussian_blob": + min_x, min_y, max_x, max_y = bounds + domain_width = max_x - min_x + center_x = min_x + 0.2 * domain_width + 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)) + raise ValueError(f"unknown simulation scalar pattern: {pattern!r}") # pragma: no cover + + +def _simulation_velocity_initializer( + pattern: str | None, velocity: tuple[float, float] +) -> Callable[[float, float], tuple[float, float]]: + """A `Field`-style `(x, y) -> (vx, vy)` callable for `SimulationConfig. + velocity_pattern` -- `None` (no pattern configured) prescribes zero + velocity, independent of whether a scalar pattern is configured, the + same "each of the two names its own thing, `None` its own absence" + shape `FieldDisplayConfig.scalar_pattern`/`vector_pattern` already use. + """ + if pattern is None: + return lambda x, y: (0.0, 0.0) + if pattern == "uniform": + return lambda x, y: velocity + raise ValueError(f"unknown simulation velocity pattern: {pattern!r}") # pragma: no cover + + +def _add_passive_scalar_transport( + window: RenderWindow, mesh: Mesh, config: PyFlowConfig +) -> Callable[[], None]: + """Wires a real `simulation.step()` into a live `pyflow run` + (Stage 4 Completion Criterion 1, TASK-030) -- the mechanism the + Passive Scalar Transport golden demo needs and no demo before it + does. Builds the initial transported scalar field and prescribed + velocity field from `config.simulation`, renders the first frame, + and returns an `on_frame` closure that advances the simulation by one + `config.numerics.timestep` and re-renders after every frame + thereafter. + + Rebuilds the rendered `gfx.Mesh` from scratch each frame (removes the + old one from `window.scene`, `build_scalar_field_mesh`s a new one) + rather than mutating the geometry's own colour buffer in place -- + `build_scalar_field_mesh`/`scalar_field_colors` are already proven + correct (TASK-017); an in-place buffer mutation would be new, + unverified pygfx-API surface for a small win on a small demo mesh + (TASK-030's own Design decision). + """ + assert window.assembled_numerics is not None + numerics = window.assembled_numerics + assert config.simulation.scalar_pattern is not None + + bounds = mesh_bounding_box(mesh) + scalar_initializer = _simulation_scalar_initializer(config.simulation.scalar_pattern, bounds) + velocity_initializer = _simulation_velocity_initializer( + config.simulation.velocity_pattern, config.simulation.velocity + ) + scalar_field = ScalarField(mesh, "tracer", initial_value=scalar_initializer) + velocity_field = VectorField( + mesh, "velocity", num_components=2, initial_value=velocity_initializer + ) + + state: dict[str, Field] = {"tracer": scalar_field} + window.simulation_fields = state + + colors = scalar_field_colors( + scalar_field, + config.field_display.low_color, + config.field_display.high_color, + config.field_display.value_range, + ) + rendered_object = build_scalar_field_mesh(scalar_field, colors) + window.scene.add(rendered_object) + + def _advance() -> None: + nonlocal state, rendered_object + state = simulation_step(state, velocity_field, numerics, config.numerics.timestep) + window.simulation_fields = state + tracer = state["tracer"] + assert isinstance(tracer, ScalarField) + colors = scalar_field_colors( + tracer, + config.field_display.low_color, + config.field_display.high_color, + config.field_display.value_range, + ) + window.scene.remove(rendered_object) + rendered_object = build_scalar_field_mesh(tracer, colors) + window.scene.add(rendered_object) + + return _advance + + def _add_field_display( window: RenderWindow, mesh: Mesh, field_display: FieldDisplayConfig ) -> _Bounds: @@ -193,7 +294,9 @@ def bootstrap( config.field_display.scalar_pattern is not None or config.field_display.vector_pattern is not None ) - if config.rendering.show_mesh or show_fields: + run_simulation = config.simulation.scalar_pattern is not None + 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 # fields -- no bespoke code, per the golden-demo public-API rule. # `show_mesh` is gated separately from `grid_color` being set @@ -211,10 +314,16 @@ 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: + # 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) fit_camera_to_bounds(window.camera, bounds) window.apply_camera_config() - window.run(max_frames=max_frames) + window.run(max_frames=max_frames, on_frame=on_frame) logger.info("pyflow exited cleanly") return window diff --git a/src/pyflow/configuration/CLAUDE.md b/src/pyflow/configuration/CLAUDE.md index e8f3b24..0433c67 100644 --- a/src/pyflow/configuration/CLAUDE.md +++ b/src/pyflow/configuration/CLAUDE.md @@ -111,6 +111,30 @@ The legend's screen position is deliberately *not* a config field -- this schema to the fields a demo author actually needs to vary, not every rendering parameter that happens to exist. +**`SimulationConfig`** (TASK-030, added 2026-08-28): `PyFlowConfig. +simulation`, seeding a real, repeatedly `simulation.step()`-advanced run +-- deliberately a separate section from `FieldDisplayConfig` above, not +a widening of it, since the two answer different questions +(`FieldDisplayConfig` seeds one static rendered frame; this seeds a live +one, driven from `RenderWindow.run(on_frame=...)` via `bootstrap.py`). +`scalar_pattern`/`velocity_pattern` follow `FieldDisplayConfig`'s own +closed-`Literal`-pattern-names precedent (`"gaussian_blob"`/`"uniform"`), +each `None` by default meaning "no live simulation" -- every existing +demo (`field_display`, `numerics_assembly`) is unaffected. Colouring the +live field reuses `field_display.low_color`/`high_color`/`value_range` +as-is; this section has no colour fields of its own, since "how a scalar +field is coloured" is a question `FieldDisplayConfig` already answers +and this section has no reason to answer twice. Shape parameters (the +blob's own center, width) are deliberately not configurable here either, +derived from the mesh's own bounds in `bootstrap.py` instead -- the same +"derived from mesh bounds, not a config field" precedent +`_scalar_display_initializer`'s own `center` already set for +`FieldDisplayConfig`'s "radial_gradient" pattern. `velocity` is a +prescribed (not solved) constant vector, `_number_pair`-normalised the +same way `RenderingConfig.pan`/`MeshConfig.origin` are -- Stage 5 is what +eventually solves Navier-Stokes for real, so a prescribed field is the +only kind of "velocity" any Stage 4 demo can legitimately have. + **`generator.py`'s `generate_config_yaml` (TASK-039, added 2026-08-21) is `loader.py` run in reverse**: `load_config` turns YAML into a validated `PyFlowConfig`; `generate_config_yaml` turns a `PyFlowConfig` diff --git a/src/pyflow/configuration/loader.py b/src/pyflow/configuration/loader.py index 88af91b..accd9a1 100644 --- a/src/pyflow/configuration/loader.py +++ b/src/pyflow/configuration/loader.py @@ -23,6 +23,7 @@ NumericsConfig, PyFlowConfig, RenderingConfig, + SimulationConfig, ) _BOUNDARY_NAMES = ("north", "south", "east", "west") @@ -102,6 +103,7 @@ def load_config(path: str | Path | None = None) -> PyFlowConfig: rendering=RenderingConfig(**raw.get("rendering", {})), mesh=MeshConfig(**raw.get("mesh", {})), field_display=FieldDisplayConfig(**raw.get("field_display", {})), + simulation=SimulationConfig(**raw.get("simulation", {})), numerics=_numerics_config_from_raw(raw.get("numerics", {})), ) config.validate() diff --git a/src/pyflow/configuration/schema.py b/src/pyflow/configuration/schema.py index 3927c6c..3f091b6 100644 --- a/src/pyflow/configuration/schema.py +++ b/src/pyflow/configuration/schema.py @@ -327,6 +327,68 @@ def validate(self) -> None: ) +ScalarTransportPattern = Literal["gaussian_blob"] +VelocityPrescriptionPattern = Literal["uniform"] + +_VALID_SCALAR_TRANSPORT_PATTERNS = frozenset(get_args(ScalarTransportPattern)) +_VALID_VELOCITY_PRESCRIPTION_PATTERNS = frozenset(get_args(VelocityPrescriptionPattern)) + + +@dataclass +class SimulationConfig: + """Live simulation stepping (TASK-030) -- distinct from + `FieldDisplayConfig` above, which seeds one static rendered frame. + `scalar_pattern`/`velocity_pattern` here seed a real, repeatedly + `simulation.step()`-advanced run, driven from `RenderWindow.run( + on_frame=...)` (`src/pyflow/bootstrap.py`) -- Stage 4's own Passive + Scalar Transport golden demo, and the first config section to wire a + live timestepping loop into an actual `pyflow run` at all. + + `None` (the default, for both) means no live simulation -- every + existing demo (`field_display`, `numerics_assembly`) is unaffected. + Colouring the live field reuses `field_display.low_color`/ + `high_color`/`value_range`/`show_legend` as-is, deliberately not + duplicated here: those already answer "how is a scalar field + coloured", a question this section has no reason to answer twice. + + Shape parameters (the blob's own center and width) are deliberately + not configurable here, derived from the mesh's own bounds in + `bootstrap.py` instead -- the same "derived from mesh bounds, not a + config field" precedent `FieldDisplayConfig`'s own + `_scalar_display_initializer`'s `center` already set, keeping this + section small the same way `FieldDisplayConfig` stays small. + `velocity` is a prescribed (not solved) constant vector -- Stage 5 is + what eventually solves Navier-Stokes for real, so a prescribed field + is the only kind of "velocity" any Stage 4 demo can legitimately have. + """ + + scalar_pattern: ScalarTransportPattern | None = None + velocity_pattern: VelocityPrescriptionPattern | None = None + velocity: tuple[float, float] = (1.0, 0.0) + + def __post_init__(self) -> None: + self.velocity = _number_pair(self.velocity, "simulation.velocity") + + def validate(self) -> None: + if ( + self.scalar_pattern is not None + and self.scalar_pattern not in _VALID_SCALAR_TRANSPORT_PATTERNS + ): + raise ValueError( + f"simulation.scalar_pattern must be one of " + f"{sorted(_VALID_SCALAR_TRANSPORT_PATTERNS)} or null, got {self.scalar_pattern!r}" + ) + if ( + self.velocity_pattern is not None + and self.velocity_pattern not in _VALID_VELOCITY_PRESCRIPTION_PATTERNS + ): + raise ValueError( + f"simulation.velocity_pattern must be one of " + f"{sorted(_VALID_VELOCITY_PRESCRIPTION_PATTERNS)} or null, " + f"got {self.velocity_pattern!r}" + ) + + AdvectionSchemeName = Literal["first_order_upwind"] DiffusionSchemeName = Literal["central_difference"] TimeIntegrationSchemeName = Literal["rk4"] @@ -594,6 +656,7 @@ class PyFlowConfig: rendering: RenderingConfig = field(default_factory=RenderingConfig) mesh: MeshConfig = field(default_factory=MeshConfig) field_display: FieldDisplayConfig = field(default_factory=FieldDisplayConfig) + simulation: SimulationConfig = field(default_factory=SimulationConfig) numerics: NumericsConfig = field(default_factory=NumericsConfig) def validate(self) -> None: @@ -601,5 +664,6 @@ def validate(self) -> None: self.rendering.validate() self.mesh.validate() self.field_display.validate() + self.simulation.validate() self.numerics.validate() _validate_boundary_conditions_jointly(self.mesh, self.numerics.boundary_conditions) diff --git a/src/pyflow/engine/CLAUDE.md b/src/pyflow/engine/CLAUDE.md index 4a5d1ed..76a0e6c 100644 --- a/src/pyflow/engine/CLAUDE.md +++ b/src/pyflow/engine/CLAUDE.md @@ -176,6 +176,31 @@ distance); `tests/unit/test_structured_cartesian_mesh.py` checks the exact formula against known grid spacing (full spacing for an interior face, half for a boundary one). +**`wrapped_neighbour_cell`, added 2026-08-28 (TASK-030), additive on +`StructuredCartesianMesh` only -- exactly the shape `boundary_face_name` +above already predicted for it.** A periodic boundary face's own +Design decision (`docs/planning/roadmap.md` TASK-030): a wrapped- +neighbour cell is mesh geometry, not a prescribed value, so it does not +become a third `BoundaryCondition` shape. For a west boundary face at +row `j`, returns the far east cell at the same row (`cell_id(nx - 1, +j)`); the other three edges mirror it on the opposite axis or index. +Raises a new `NotABoundaryFaceError` (`mesh.py`'s own -- cannot reuse +`boundary_condition.py`'s identically-named class without inverting the +dependency direction between the two modules) for an interior face, +exercised directly for both a vertical and a horizontal interior face +(found necessary while confirming coverage: a single `next(...)`- +selected interior face always picked the vertical branch first, leaving +the horizontal `raise` line genuinely untested by one test alone). +**The periodic distance is `2 * mesh.face_centroid_distance(face)`, not +a new geometry accessor** -- verified numerically, not just +algebraically, before relying on it: on a mesh with distinct `dx`/`dy` +and a non-trivial origin, doubling a boundary face's own owner-to-face +distance reproduces the true uniform grid spacing exactly, matching an +ordinary interior face's distance on the same mesh to float precision. +The periodic "neighbour" is one full cell-width away, not the real +wrapped cell's actual (far-side-of-the-domain) centroid, so the plain +interior formula would be wildly wrong applied naively across the wrap. + **`field.py`** (TASK-014, done 2026-08-21) is `Field`, the abstract base every physical quantity the engine transports will share -- Variables, in `docs/architecture/engine.md`'s terms. Deliberately carries no @@ -390,14 +415,16 @@ condition's prescribed value is used directly; a Neumann (`kind == advective term -- its face value is the owner's own, zero-order extrapolation, per `docs/handbook/numerical-methods/ boundary-conditions.md`'s "typically extrapolated from the adjacent -cell-centred value". Inflow at a boundary whose named edge has no -`BoundaryCondition` at all (the periodic case -- `boundary_condition.py` -resolves no object for it) raises `UnconfiguredBoundaryFaceError` rather -than silently extrapolating, which would be a plausible-looking wrong -answer for a periodic boundary specifically (it needs the wrapped -neighbour's actual value, not an extrapolation). Outflow at the same -face never raises it or reads the boundary condition at all -- the -upstream value is simply the owner's. +cell-centred value". Inflow at a boundary whose named edge has no `BoundaryCondition` at all +and is not periodic either raises `UnconfiguredBoundaryFaceError` rather +than silently extrapolating. **A periodic face (`periodic_pairs`, +TASK-030, 2026-08-28) is genuinely handled, not merely not-raised**: +`mesh.wrapped_neighbour_cell(face)` stands in for `neighbour` before +either helper below runs, so the rest of `flux` treats it exactly like a +real interior face -- the wrapped cell's own actual value, not an +extrapolation, and `boundary_conditions` is never consulted for it. +Outflow at the same face never raises it or reads the boundary condition +at all -- the upstream value is simply the owner's. `FirstOrderUpwindAdvection` joins `test_advection_contract.py`'s existing parametrised suite (Stage 4 Completion Criterion 3) with no @@ -436,13 +463,19 @@ extrapolation only). The difference is what each interface's own Neumann shape actually means physically: advection's boundary value is extrapolated because advection has no natural use for a prescribed *gradient*, while diffusion's whole boundary contribution at a Neumann -face *is* the prescribed gradient. No condition configured (the periodic -case) raises `UnconfiguredBoundaryFaceError` -- `diffusion.py`'s own -class, not shared with `advection.py`'s identically-named one (each -numerics interface module owns its own exception vocabulary). Unlike -advection, there is no inflow/outflow carve-out: diffusion has no flow -direction, so every boundary face needs a configured condition -unconditionally. +face *is* the prescribed gradient. No condition configured and not +periodic either raises `UnconfiguredBoundaryFaceError` -- `diffusion.py`'s +own class, not shared with `advection.py`'s identically-named one (each +numerics interface module owns its own exception vocabulary). **A +periodic face (TASK-030, 2026-08-28) substitutes the wrapped neighbour's +own value and the correct one-cell-width distance (`2 * +mesh.face_centroid_distance(face)`, not the plain boundary-face +distance) before falling through the ordinary interior formula** -- see +`src/pyflow/engine/numerics/CLAUDE.md`'s own TASK-030 entry for why +doubling is exactly right on a uniform mesh, verified numerically. +Unlike advection, there is no inflow/outflow carve-out: diffusion has no +flow direction, so every non-periodic boundary face needs a configured +condition unconditionally. `CentralDifferenceDiffusion` joins `test_diffusion_contract.py`'s existing parametrised suite (Stage 4 Completion Criterion 3) with no @@ -958,6 +991,18 @@ Linear Solver, Pressure-Velocity Coupling, Boundary Condition -- now have a real concrete scheme under `src/`; Stage 3 Completion Criterion 1's carve-out (this section's own opening paragraphs) is fully retired. +**TASK-030 (2026-08-28, the next day) widens `assemble_numerics` once +more, with no `_Null*` question left to touch.** A second mapping, +`periodic_pairs` (`{face_name: opposite_face_name}` for every face +actually configured periodic), is built in the same per-face loop that +already builds `boundary_conditions`, and threaded into the advection/ +diffusion factories alongside it -- `register_advection_scheme` widens +to a two-argument factory, `register_diffusion_scheme` to three, needing +a new `_resolve_with_three_arguments` generic helper. The one-argument +helper this module used to route advection through, +`_resolve_with_argument`, is deleted in the same change as genuinely +dead code -- its only caller moved to the two-argument helper. + **Registration refuses to overwrite a different factory** (`DuplicateSchemeError`, added 2026-08-24). The registries are module-level and filled by import side effect, so "last import wins" @@ -977,7 +1022,11 @@ no-op, so a module imported twice does not raise. Dirichlet/Neumann shapes, so `assemble_numerics` reports a periodic face's configured type in `.names` but omits it from `.boundary_conditions` entirely, rather than fabricating an object the -interface has no shape for. +interface has no shape for. **Instead (TASK-030), it becomes a key in the +`periodic_pairs` mapping above** -- `assemble_numerics` never fabricates +a `BoundaryCondition` for it, and a concrete advection/diffusion scheme +consults `periodic_pairs` itself (via `StructuredCartesianMesh. +wrapped_neighbour_cell`) rather than the orchestrator special-casing it. **`assemble_numerics` now resolves `boundary_conditions` before advection/diffusion, and `register_advection_scheme`/ diff --git a/src/pyflow/engine/mesh.py b/src/pyflow/engine/mesh.py index 94a5bc4..8a154f1 100644 --- a/src/pyflow/engine/mesh.py +++ b/src/pyflow/engine/mesh.py @@ -33,6 +33,21 @@ """ +class NotABoundaryFaceError(ValueError): + """Raised when `StructuredCartesianMesh.wrapped_neighbour_cell` is + asked for an interior face's wrapped partner -- meaningless the same + way `BoundaryCondition.evaluate` on an interior face is + (`boundary_condition.py`'s own `NotABoundaryFaceError`), but its own + class here rather than imported from there: `mesh.py` is a + foundational engine layer `engine/numerics/` depends on, not the + reverse, so it cannot reuse a numerics-layer exception without + inverting that dependency. Each numerics interface already owns its + own exception vocabulary (`advection.py`/`diffusion.py`'s own, + identically-reasoned, `UnconfiguredBoundaryFaceError`); this extends + the same pattern to `mesh.py` itself (TASK-030). + """ + + class InvalidMeshEntityError(IndexError): """Raised when a cell or face id does not identify an entity of the mesh it was asked of -- the same exception class for every @@ -359,6 +374,42 @@ def boundary_face_name(self, face: int) -> BoundaryFaceName | None: return "north" return None + def wrapped_neighbour_cell(self, face: int) -> int: + """The cell a periodic boundary face wraps to -- the owner cell + of the *opposite* domain edge, at the same row (west/east) or + column (north/south). + + Pure mesh geometry, not a prescribed value: additive, off the + abstract `Mesh` interface, the same precedent `cell_id`/ + `cell_index`/`boundary_face_name` already set for a concept only + a structured, axis-aligned mesh has (TASK-030's own Design + decision, `docs/planning/roadmap.md`). Whether `face` is + *configured* periodic is a numerics-layer concern + (`engine/numerics/assembly.py`'s own `periodic_pairs`); this + method answers the purely geometric question for any boundary + face regardless of its configured type. + + Raises `NotABoundaryFaceError` if `face` is not a boundary face + -- an interior face already has a real neighbour, so "wrapped + partner" is meaningless for it, the same reasoning + `BoundaryCondition._check_boundary_face` applies to `evaluate`. + """ + self._check_face(face) + if face < self._num_vertical_faces: + p, j = self._decode_vertical_face(face) + if p == 0: + return self.cell_id(self._nx - 1, j) + if p == self._nx: + return self.cell_id(0, j) + raise NotABoundaryFaceError(f"face {face} is not a boundary face") + + i, q = self._decode_horizontal_face(face) + if q == 0: + return self.cell_id(i, self._ny - 1) + if q == self._ny: + return self.cell_id(i, 0) + raise NotABoundaryFaceError(f"face {face} is not a boundary face") + # -- face id encoding/decoding --------------------------------------- # Vertical faces (normal in x): position p in [0, nx] (west/east # boundary at p == 0 / p == nx), row j in [0, ny). Horizontal faces diff --git a/src/pyflow/engine/numerics/CLAUDE.md b/src/pyflow/engine/numerics/CLAUDE.md index 27e1a1a..5624df3 100644 --- a/src/pyflow/engine/numerics/CLAUDE.md +++ b/src/pyflow/engine/numerics/CLAUDE.md @@ -105,6 +105,28 @@ a mutation-caught test-quality finding in `test_assembly.py`'s own fixture (a coincidentally-shared `0.0` had let a wrong-field regression pass unnoticed, fixed per `docs/practices.md`'s "distinct factors" rule). +**`boundary_condition.py`'s own scope stays exactly as Stage 3 left it +(TASK-030, 2026-08-28, Stage 4's ninth and last task) -- periodic never +becomes a third `BoundaryCondition` shape.** A wrapped-neighbour cell is +mesh geometry, not a prescribed value, so it lives on `mesh.py` +(`StructuredCartesianMesh.wrapped_neighbour_cell`) instead -- the same +"structured-only concept, off the abstract `Mesh` interface" precedent +`cell_id`/`cell_index`/`boundary_face_name` already set. `advection.py`/ +`diffusion.py` both gain a `periodic_pairs: Mapping[str, str]` +constructor parameter (mirrors `diffusion_coefficient`'s own TASK-024 +precedent: registry-level widening, no interface change, no new ADR) and +consult it before falling through their existing interior-neighbour +formula -- literally the same formula already used for a real interior +neighbour, once the wrapped cell stands in for `neighbour`. Diffusion's +own distance across a periodic face is `2 * mesh.face_centroid_distance +(face)`, verified numerically (not just algebraically) to reproduce the +true uniform grid spacing exactly, not the wrapped cell's actual +(far-side-of-the-domain) centroid distance. See `src/pyflow/engine/ +CLAUDE.md`'s own entry for the real content, including the numerical +finding behind the round-trip scenario's own convergence-based +criterion (first-order upwind's O(dx) numerical diffusion means "matches +exactly" is the wrong claim to check, even for a correct wrap). + **`assembly.py`** (TASK-021) is different in kind from the interfaces above: not an interface, but the registry (`register_advection_scheme` and five siblings) and `assemble_numerics(NumericsConfig) -> AssembledNumerics` @@ -123,6 +145,17 @@ See the module's own docstring and `src/pyflow/engine/CLAUDE.md`'s `numerics/` entry for the full retirement history and the exception's own now-closed record against Stage 3 Completion Criterion 1. +**TASK-030 (2026-08-28) widens `assemble_numerics` once more, without +touching any `_Null*` question** (there are none left to touch): a +second mapping, `periodic_pairs`, is built in the same per-face loop that +already builds `boundary_conditions`, and threaded into the advection/ +diffusion factories alongside it. `register_advection_scheme` widens to +a two-argument factory; `register_diffusion_scheme` widens to three, +needing a new `_resolve_with_three_arguments` generic helper. The +one-argument helper this file's own history used to route advection +through, `_resolve_with_argument`, is deleted as genuinely dead code -- +its only caller moved to the two-argument helper in the same change. + 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/advection.py b/src/pyflow/engine/numerics/advection.py index f691d2b..658af01 100644 --- a/src/pyflow/engine/numerics/advection.py +++ b/src/pyflow/engine/numerics/advection.py @@ -80,9 +80,9 @@ def flux(self, field: Field, velocity: VectorField) -> torch.Tensor: class UnconfiguredBoundaryFaceError(ValueError): """Raised when inflow occurs at a boundary face whose named edge (`StructuredCartesianMesh.boundary_face_name`) has no - `BoundaryCondition` in this scheme's own mapping -- the periodic - case (`assemble_numerics` resolves no `BoundaryCondition` object for - a periodic-type boundary; TASK-030's own concern, not this scheme's). + `BoundaryCondition` in this scheme's own mapping and is not periodic + either (`periodic_pairs`, TASK-030) -- a face genuinely wired to + neither. Outflow at the same face never raises this: the upstream value is the owner cell's own, and the boundary condition is never consulted @@ -103,10 +103,23 @@ class FirstOrderUpwindAdvection(AdvectionScheme): `docs/planning/roadmap.md`): holds the boundary conditions it needs, keyed by named edge, and consults them itself rather than the orchestrator substituting a value into its output afterward. + + **Periodic-aware the same way (TASK-030).** `periodic_pairs` names + which boundary faces wrap to the opposite edge (e.g. `{"west": + "east", "east": "west"}`) -- absence from it is never read as + "periodic" by omission. At a periodic face, `mesh.wrapped_neighbour_cell` + stands in for `neighbour` before either helper below runs, so the + rest of `flux` treats it exactly like a genuine interior face; a + periodic face never consults `boundary_conditions` at all. """ - 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 flux(self, field: Field, velocity: VectorField) -> torch.Tensor: self._check_velocity(velocity) @@ -117,6 +130,10 @@ def flux(self, field: Field, velocity: VectorField) -> torch.Tensor: result = torch.zeros(mesh.num_faces, dtype=torch.float64) for face in range(mesh.num_faces): owner, neighbour = mesh.face_neighbours(face) + if neighbour is None: + boundary_name = mesh.boundary_face_name(face) + if boundary_name in self._periodic_pairs: + neighbour = mesh.wrapped_neighbour_cell(face) normal_x, normal_y = mesh.face_normal(face) velocity_normal = self._face_normal_velocity( velocity, owner, neighbour, normal_x, normal_y diff --git a/src/pyflow/engine/numerics/assembly.py b/src/pyflow/engine/numerics/assembly.py index a28aa83..118f1f7 100644 --- a/src/pyflow/engine/numerics/assembly.py +++ b/src/pyflow/engine/numerics/assembly.py @@ -56,6 +56,20 @@ `AssembledNumerics.names` but omits it from `.boundary_conditions` entirely, rather than fabricating an object the interface has no shape for. + +**Instead, `assemble_numerics` builds a second, separate mapping +(TASK-030): `periodic_pairs`, `{face_name: opposite_face_name}` for every +face actually configured periodic.** Threaded into the advection/ +diffusion factories alongside `boundary_conditions` (their own second and +third constructor arguments respectively, `register_advection_scheme`/ +`register_diffusion_scheme`'s own docstrings) -- a periodic face is mesh +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 +already follows. """ from __future__ import annotations @@ -77,6 +91,9 @@ from pyflow.engine.numerics.time_integrator import RK4Integrator, TimeIntegrator _BOUNDARY_FACE_NAMES = ("north", "south", "east", "west") +_PAIRED_BOUNDARY = {"north": "south", "south": "north", "east": "west", "west": "east"} +"""Local to this module, deliberately, same reasoning as `_BOUNDARY_FACE_NAMES` +above -- `schema.py`'s own identically-shaped dict is private there.""" class UnknownSchemeError(ValueError): @@ -124,9 +141,11 @@ class AssembledNumerics: names: Mapping[str, str] -_advection_registry: dict[str, Callable[[Mapping[str, BoundaryCondition]], AdvectionScheme]] = {} +_advection_registry: dict[ + str, Callable[[Mapping[str, BoundaryCondition], Mapping[str, str]], AdvectionScheme] +] = {} _diffusion_registry: dict[ - str, Callable[[Mapping[str, BoundaryCondition], float], DiffusionScheme] + str, Callable[[Mapping[str, BoundaryCondition], Mapping[str, str], float], DiffusionScheme] ] = {} _time_integrator_registry: dict[str, Callable[[], TimeIntegrator]] = {} _linear_solver_registry: dict[str, Callable[[float, int], LinearSolver]] = {} @@ -151,31 +170,36 @@ def _register[F](registry: dict[str, F], name: str, factory: F, component: str) def register_advection_scheme( - name: str, factory: Callable[[Mapping[str, BoundaryCondition]], AdvectionScheme] + name: str, + factory: Callable[[Mapping[str, BoundaryCondition], Mapping[str, str]], AdvectionScheme], ) -> None: - """Make `name` resolve to `factory(boundary_conditions)` in future - `assemble_numerics` calls -- `boundary_conditions` is the same - face-name-keyed mapping `AssembledNumerics.boundary_conditions` + """Make `name` resolve to `factory(boundary_conditions, periodic_pairs)` + in future `assemble_numerics` calls -- `boundary_conditions` is the + same face-name-keyed mapping `AssembledNumerics.boundary_conditions` carries, resolved before advection/diffusion so a concrete scheme can receive the boundary conditions it needs at construction (TASK-040's own Design decision, `docs/planning/roadmap.md`), rather than the - orchestrator substituting a value after the fact. + orchestrator substituting a value after the fact. `periodic_pairs` + (TASK-030) is the same shape of addition, one call later: which + boundary faces wrap to the opposite edge, containing only faces + actually configured periodic. """ _register(_advection_registry, name, factory, "advection") def register_diffusion_scheme( - name: str, factory: Callable[[Mapping[str, BoundaryCondition], float], DiffusionScheme] + name: str, + factory: Callable[[Mapping[str, BoundaryCondition], Mapping[str, str], float], DiffusionScheme], ) -> None: - """Make `name` resolve to `factory(boundary_conditions, + """Make `name` resolve to `factory(boundary_conditions, periodic_pairs, diffusion_coefficient)` in future `assemble_numerics` calls -- - `boundary_conditions` the same as `register_advection_scheme`'s own, - `diffusion_coefficient` is `NumericsConfig.diffusion_coefficient` - (TASK-024's own Design decision, `docs/planning/roadmap.md`): a - concrete diffusion scheme is constructed with the physical - coefficient (Gamma) it needs, the same "constructed with it, not - handed it after the fact" reasoning `boundary_conditions` already - established. + `boundary_conditions`/`periodic_pairs` the same as + `register_advection_scheme`'s own, `diffusion_coefficient` is + `NumericsConfig.diffusion_coefficient` (TASK-024's own Design + decision, `docs/planning/roadmap.md`): a concrete diffusion scheme is + constructed with the physical coefficient (Gamma) it needs, the same + "constructed with it, not handed it after the fact" reasoning + `boundary_conditions` already established. """ _register(_diffusion_registry, name, factory, "diffusion") @@ -225,42 +249,52 @@ def _resolve[T](registry: Mapping[str, Callable[[], T]], name: str, component: s return factory() -def _resolve_with_argument[T, A]( - registry: Mapping[str, Callable[[A], T]], name: str, argument: A, component: str +def _resolve_with_two_arguments[T, A, B]( + registry: Mapping[str, Callable[[A, B], T]], + name: str, + argument_a: A, + argument_b: B, + component: str, ) -> T: - """Same as `_resolve`, for the three components whose factory needs - one constructor argument -- advection/diffusion (the boundary- - conditions mapping) and pressure_coupling (the resolved - `LinearSolver`) -- rather than three near-identical inline + """Same as `_resolve`, for the components whose factory needs two + constructor arguments -- advection (`boundary_conditions`, + `periodic_pairs`, TASK-030), linear_solver (`tolerance`, + `max_iterations`) and pressure_coupling (the resolved `LinearSolver`, + `boundary_conditions`) -- rather than several near-identical inline get/raise/call blocks repeating the same lookup (found during - TASK-040's own review cycle). + TASK-040's own review cycle). **Not the one-argument helper this + docstring used to describe advection sharing with diffusion** -- + that helper (`_resolve_with_argument`) was retired the same day + advection itself gained a second constructor argument, leaving it + with no remaining caller. """ factory = registry.get(name) if factory is None: raise UnknownSchemeError(f"no {component} implementation registered under {name!r}") - return factory(argument) + return factory(argument_a, argument_b) -def _resolve_with_two_arguments[T, A, B]( - registry: Mapping[str, Callable[[A, B], T]], +def _resolve_with_three_arguments[T, A, B, C]( + registry: Mapping[str, Callable[[A, B, C], T]], name: str, argument_a: A, argument_b: B, + argument_c: C, component: str, ) -> T: - """Same as `_resolve_with_argument`, for diffusion alone -- the one - component whose factory needs two constructor arguments (the - boundary-conditions mapping *and* `config.diffusion_coefficient`, - TASK-024's own Design decision) rather than one. Kept as its own - generic helper instead of widening `_resolve_with_argument` itself, - since advection/pressure_coupling still only need one argument each - and a shared two-argument signature would force both to pass an - unused second one. + """Same as `_resolve_with_two_arguments`, for diffusion alone -- the + one component whose factory needs three constructor arguments + (`boundary_conditions`, `periodic_pairs`, `diffusion_coefficient`, + TASK-030) rather than two. Kept as its own generic helper instead of + widening `_resolve_with_two_arguments` itself, since advection/ + linear_solver/pressure_coupling still only need two arguments each + and a shared three-argument signature would force all three to pass + an unused third. """ factory = registry.get(name) if factory is None: raise UnknownSchemeError(f"no {component} implementation registered under {name!r}") - return factory(argument_a, argument_b) + return factory(argument_a, argument_b, argument_c) def assemble_numerics(config: NumericsConfig) -> AssembledNumerics: @@ -288,21 +322,27 @@ def assemble_numerics(config: NumericsConfig) -> AssembledNumerics: """ names: dict[str, str] = {} boundary_conditions_by_face: dict[str, BoundaryCondition] = {} + periodic_pairs_by_face: dict[str, str] = {} for face_name in _BOUNDARY_FACE_NAMES: face_config: BoundaryFaceConfig = getattr(config.boundary_conditions, face_name) names[f"boundary_conditions.{face_name}"] = face_config.type + if face_config.type == "periodic": + periodic_pairs_by_face[face_name] = _PAIRED_BOUNDARY[face_name] + continue boundary_factory = _boundary_condition_registry.get(face_config.type) if boundary_factory is not None: boundary_conditions_by_face[face_name] = boundary_factory(face_config) boundary_conditions = MappingProxyType(boundary_conditions_by_face) + periodic_pairs = MappingProxyType(periodic_pairs_by_face) - advection = _resolve_with_argument( - _advection_registry, config.advection, boundary_conditions, "advection" + advection = _resolve_with_two_arguments( + _advection_registry, config.advection, boundary_conditions, periodic_pairs, "advection" ) - diffusion = _resolve_with_two_arguments( + diffusion = _resolve_with_three_arguments( _diffusion_registry, config.diffusion, boundary_conditions, + periodic_pairs, config.diffusion_coefficient, "diffusion", ) diff --git a/src/pyflow/engine/numerics/diffusion.py b/src/pyflow/engine/numerics/diffusion.py index 1a6c61a..c5abe83 100644 --- a/src/pyflow/engine/numerics/diffusion.py +++ b/src/pyflow/engine/numerics/diffusion.py @@ -42,15 +42,16 @@ def flux(self, field: Field) -> torch.Tensor: class UnconfiguredBoundaryFaceError(ValueError): """Raised when a `CentralDifferenceDiffusion` boundary face's named edge (`StructuredCartesianMesh.boundary_face_name`) has no - `BoundaryCondition` in this scheme's own mapping -- the periodic - case, mirroring `advection.py`'s identically-named exception for the - identical underlying reason, but its own class in its own module - (each numerics interface owns its own exception vocabulary). + `BoundaryCondition` in this scheme's own mapping and is not periodic + either (`periodic_pairs`, TASK-030) -- mirroring `advection.py`'s + identically-named exception for the identical underlying reason, but + its own class in its own module (each numerics interface owns its + own exception vocabulary). Unlike advection's own version, there is no inflow/outflow carve-out - here: diffusion has no flow direction, so *every* boundary face needs - a configured condition to compute a diffusive flux at all, not only - the faces where flow happens to be entering. + here: diffusion has no flow direction, so *every* non-periodic + boundary face needs a configured condition to compute a diffusive + flux at all, not only the faces where flow happens to be entering. """ @@ -67,14 +68,26 @@ class CentralDifferenceDiffusion(DiffusionScheme): decision): holds the boundary conditions and the diffusion coefficient (Gamma, `NumericsConfig.diffusion_coefficient`) it needs, rather than the orchestrator substituting either in afterward. + + **Periodic-aware the same way `FirstOrderUpwindAdvection` is + (TASK-030).** At a face named in `periodic_pairs`, `flux` substitutes + `mesh.wrapped_neighbour_cell` for `neighbour` and the correct + one-cell-width distance (`2 * mesh.face_centroid_distance(face)` -- + see this module's own TASK-030 note in `docs/planning/roadmap.md` for + why doubling the owner-to-face distance is exactly right on a uniform + mesh, verified numerically before relying on it) before falling + through the ordinary interior-face formula below; `boundary_conditions` + is never consulted for a periodic face. """ def __init__( self, boundary_conditions: Mapping[str, BoundaryCondition], + periodic_pairs: Mapping[str, str], diffusion_coefficient: float, ) -> None: self._boundary_conditions = boundary_conditions + self._periodic_pairs = periodic_pairs self._gamma = diffusion_coefficient def flux(self, field: Field) -> torch.Tensor: @@ -87,6 +100,11 @@ def flux(self, field: Field) -> torch.Tensor: owner, neighbour = mesh.face_neighbours(face) owner_value = float(field.value_at(owner)) distance = mesh.face_centroid_distance(face) + if neighbour is None: + boundary_name = mesh.boundary_face_name(face) + if boundary_name in self._periodic_pairs: + neighbour = mesh.wrapped_neighbour_cell(face) + distance = 2 * distance if neighbour is not None: neighbour_value = float(field.value_at(neighbour)) gradient = (neighbour_value - owner_value) / distance diff --git a/src/pyflow/engine/numerics/pressure_coupling.py b/src/pyflow/engine/numerics/pressure_coupling.py index c0ccfac..fbc4781 100644 --- a/src/pyflow/engine/numerics/pressure_coupling.py +++ b/src/pyflow/engine/numerics/pressure_coupling.py @@ -146,7 +146,9 @@ def __init__( {name: _ZeroGradientPressureCondition() for name in _PRESSURE_BOUNDARY_FACE_NAMES} ) self._diffusion = CentralDifferenceDiffusion( - pressure_boundary_conditions, diffusion_coefficient=1.0 + pressure_boundary_conditions, + {}, # no periodic pressure boundaries in this task's own scope (TASK-030) + diffusion_coefficient=1.0, ) self._gradient = GreenGaussGradient(pressure_boundary_conditions) self._divergence = GreenGaussDivergence(boundary_conditions) diff --git a/src/pyflow/rendering/CLAUDE.md b/src/pyflow/rendering/CLAUDE.md index ea08504..d9a290e 100644 --- a/src/pyflow/rendering/CLAUDE.md +++ b/src/pyflow/rendering/CLAUDE.md @@ -330,3 +330,26 @@ a `BootstrapResult` wrapper type instead would have changed `bootstrap()`'s return type for every existing caller (tests, `__main__.py`, this file's own module) merely to avoid one import, which is a larger blast radius for a smaller problem. + +## Live Simulation Reporting (TASK-030, done 2026-08-28) + +**`RenderWindow.simulation_fields: Mapping[str, Field] | None`, the same +shape `assembled_numerics` above already established** -- another narrow, +deliberate exception to `RenderWindow`'s own "no simulation content" +scope: `RenderWindow` itself never calls `simulation.step()` and knows +nothing about what the mapping holds beyond its type. +`bootstrap.py`'s `_add_passive_scalar_transport` (the Passive Scalar +Transport golden demo's own mechanism -- the first config to wire a real +`simulation.step()` call into `RenderWindow.run(on_frame=...)`, not a +capability of this package's own) sets it once per frame, inside its own +`on_frame` closure, purely so a caller -- the golden demo's own +regression test, most directly -- has one place to read back the real +field state a rendered frame came from, not only its rendered pixels. +`None` for every run before this task's own (every existing demo) and +for a `RenderWindow` built directly without going through `bootstrap()`. + +Importing `pyflow.engine.field` here is the same kind of new-but-narrow +dependency `assembled_numerics` above already accepted for +`pyflow.engine.numerics.assembly` -- `Field` is only ever used as a type +annotation, and the same "avoid a wrapper type that would change +`bootstrap()`'s return type for every caller" reasoning applies. diff --git a/src/pyflow/rendering/window.py b/src/pyflow/rendering/window.py index ccd6d21..785d6b2 100644 --- a/src/pyflow/rendering/window.py +++ b/src/pyflow/rendering/window.py @@ -9,13 +9,14 @@ from __future__ import annotations -from collections.abc import Callable +from collections.abc import Callable, Mapping from typing import Any import pygfx as gfx from pyflow.configuration.schema import RenderingConfig from pyflow.engine import get_logger +from pyflow.engine.field import Field from pyflow.engine.numerics.assembly import AssembledNumerics from pyflow.rendering.canvas import create_canvas, get_loop @@ -113,6 +114,17 @@ def __init__(self, config: RenderingConfig) -> None: content" scope above -- `bootstrap()` sets this after calling `assemble_numerics()`, purely so a caller has one place to read back what got assembled, per Stage 3 Completion Criterion 8.""" + self.simulation_fields: Mapping[str, Field] | None = None + """The live simulation's own transported fields, updated every + frame (TASK-030), or `None` if this run has no live simulation + (every demo before Passive Scalar Transport) or the caller built + a `RenderWindow` directly. Same shape as `assembled_numerics` + above: `RenderWindow` never advances a simulation itself, true to + its own "no simulation content" scope; `bootstrap()`'s own + `on_frame` closure sets this each step, purely so a caller (a + golden-demo regression test, most directly) has one place to read + back the real field state a rendered frame came from, not only + its rendered pixels.""" self._on_frame: Callable[[], None] | None = None self._pan_drag_start_screen: tuple[float, float] | None = None self._pan_drag_start_position: tuple[float, float, float] | None = None diff --git a/tests/features/passive_scalar_transport.feature b/tests/features/passive_scalar_transport.feature new file mode 100644 index 0000000..ea1ce58 --- /dev/null +++ b/tests/features/passive_scalar_transport.feature @@ -0,0 +1,31 @@ +# The acceptance criteria for the Passive Scalar Transport golden demo, +# in the form that executes. `docs/implementation/golden-demos.md` +# describes what this demo is and why; the scenarios below are what +# "working" means -- see `adr/ADR-007-executable-acceptance-criteria.md`. + +Feature: Passive Scalar Transport + Stage 4's own golden demo (TASK-030, roadmap Stage 4 Completion + Criterion 1): the first `pyflow run` that actually steps a real + simulation forward, live, via `simulation.step()` wired into + `RenderWindow.run(on_frame=...)` -- every earlier demo rendered exactly + one static frame. A prescribed (not solved -- Stage 5 solves + Navier-Stokes for real) uniform velocity field carries a scalar blob + across a periodic domain. + + Background: + Given the golden demo "passive_scalar_transport" + + Scenario: A user can run it with the documented command + When it is run through the public CLI, headless + Then the command exits cleanly + + # "physical fields evolve" (`docs/implementation/mvp.md`'s own + # Definition of Done), measured directly rather than only checked by + # pixel-diffing two frames: the field's own mass-weighted centroid + # (the closest thing a scalar field has to "position") has to move + # downstream at roughly the prescribed velocity over real elapsed + # time -- not merely "the pixels changed somewhere", which a stalled + # or frozen simulation loop could not distinguish from a broken one. + Scenario: The transported field's own centroid moves downstream at the prescribed velocity over real elapsed time + When it is bootstrapped once after a few real timesteps and again after many more + Then the field's mass-weighted centroid has moved downstream by approximately the prescribed velocity times the elapsed time diff --git a/tests/features/periodic_boundary.feature b/tests/features/periodic_boundary.feature new file mode 100644 index 0000000..00425a1 --- /dev/null +++ b/tests/features/periodic_boundary.feature @@ -0,0 +1,51 @@ +# The acceptance criteria for Periodic Boundary (TASK-030, Stage 4's +# ninth and last task). Not a golden demo -- no config file under +# `examples/golden-demos/`, no CLI run; `tests/unit/ +# test_periodic_boundary.py` binds these scenarios directly, per that +# directory's own scope (isolated logic, no process boundary). + +Feature: Periodic Boundary + Periodic bypasses `BoundaryCondition` entirely (`docs/planning/ + roadmap.md` TASK-030's own Design decision): a wrapped-neighbour cell + is mesh geometry, not a prescribed value, so a real Advection/Diffusion + scheme wired with a periodic pairing must consult the opposite edge's + own owner cell directly -- through the same formula it already uses for + a genuine interior neighbour -- rather than any `BoundaryCondition`. + + Background: + Given a small, non-square, non-trivially-origined mesh whose cells each hold a distinct value + + Scenario: A real periodic pairing wired into advection reads the wrapped neighbour's own value, not a boundary condition + Given a west boundary face configured periodic with its east partner + And no boundary condition configured for that face at all + When the advective flux is computed with inflow at that face + Then the inflow boundary face's implied value is the wrapped neighbour's own value, not the owner's + + Scenario: A real periodic pairing wired into diffusion computes the gradient across the wrapped neighbour at one full cell width + Given a west boundary face configured periodic with its east partner + And no boundary condition configured for that face at all + When the diffusive flux is computed at that face + Then that boundary face's flux equals the diffusion coefficient times the difference between the wrapped neighbour's and the owner's own value, divided by one full grid spacing + + Scenario: A field advected once fully around a periodic domain converges toward its starting distribution as the mesh is refined + # Advection alone, not diffusion -- diffusion's own periodic wiring is + # already checked directly above, and diffusion is a genuinely + # irreversible process with no reason to undo itself over one lap, so + # it would falsify this specific claim rather than test it. + # + # "Matches exactly" is not the right claim at any one fixed, cheap + # mesh: first-order upwind's own O(dx) numerical diffusion smooths a + # field over the distance it travels regardless of how correctly the + # wrap itself is implemented, verified numerically (not assumed) + # before writing this scenario. What a *wrong* wrap (mirrored or + # clamped to the owner's own cell at the periodic boundary, instead + # of the opposite edge) cannot reproduce is refinement actually + # closing the gap: measured directly, a real wrap's own round-trip + # error drops by roughly 62% over a 4x mesh refinement, while a + # mirrored/clamped one drops by only roughly 16% and stays several + # times larger throughout -- the discriminating property this + # scenario checks. + Given a mesh whose east and west edges are both periodic + And a uniform velocity field that carries the domain's own scalar field once fully around it in a whole number of real timesteps + When the field is advected for exactly that many real timesteps, at two mesh resolutions four times apart + Then the round-trip error at the finer resolution is well under two thirds of the error at the coarser one diff --git a/tests/golden/CLAUDE.md b/tests/golden/CLAUDE.md index 1f631bf..2fbd9ce 100644 --- a/tests/golden/CLAUDE.md +++ b/tests/golden/CLAUDE.md @@ -62,3 +62,22 @@ There is no demo-specific Python module to load any more (no of this test needed): golden demos are plain config files under `examples/golden-demos/`, not scripts, per the public-API rule in `docs/implementation/golden-demos.md`. + +**`test_passive_scalar_transport.py` (TASK-030, added 2026-08-28) is the +fifth demo module, and PyFlow's first that computes real physics** -- +Stage 4 Completion Criterion 1's own golden demo, the first `pyflow run` +that steps a real simulation forward live rather than rendering one +static frame. Binds `tests/features/passive_scalar_transport.feature`: +the required CLI-subprocess scenario every demo carries, plus one +demo-specific step (`conftest.py`'s shared vocabulary only knows how to +render exactly one or two frames, not step a live simulation forward by +a specific count and read its own field state back) that bootstraps the +demo twice, at two different frame counts, and reads +`RenderWindow.simulation_fields` back from each -- a genuine physical +claim (the transported field's own mass-weighted centroid moves +downstream at approximately the prescribed velocity over real elapsed +time), not only a pixel-diff. **The tolerance (`rel=0.15`) was measured +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. diff --git a/tests/golden/test_passive_scalar_transport.py b/tests/golden/test_passive_scalar_transport.py new file mode 100644 index 0000000..b17e298 --- /dev/null +++ b/tests/golden/test_passive_scalar_transport.py @@ -0,0 +1,96 @@ +"""Passive Scalar Transport golden demo (TASK-030). + +The acceptance criteria are `tests/features/passive_scalar_transport.feature` +(`adr/ADR-007-executable-acceptance-criteria.md`). This module binds them +and supplies the one step only this demo needs -- reading back +`RenderWindow.simulation_fields` at two different frame counts, since the +shared vocabulary in `conftest.py` only knows how to render exactly one or +two *rendered* frames, not step a live simulation forward by a specific +count and read its own field state back. +""" + +from __future__ import annotations + +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("passive_scalar_transport.feature") + +_EARLY_FRAMES = 1 +_LATE_FRAMES = 101 +"""100 real RK4 timesteps apart -- at this demo's own configured +`velocity`/`timestep` (1.0, 0.02), that is 2.0 world units of expected +downstream travel, comfortably less than the mesh's own 5.0-unit domain +width (`mesh.extent[0] * mesh.spacing[0]`), so the blob never wraps +around the periodic boundary during this specific measurement -- the wrap +itself is `periodic_boundary.feature`'s own claim, not this one's. +""" + + +def _centroid_x(field: ScalarField) -> float: + """The field's own mass-weighted centroid x-position -- the closest + thing a scalar field has to "position", valid because every value + here is non-negative (a Gaussian blob never goes negative). + """ + mesh = field.mesh + total_mass = 0.0 + weighted_x = 0.0 + for cell in range(mesh.num_cells): + value = float(field.value_at(cell)) + x, _y = mesh.cell_centroid(cell) + total_mass += value + weighted_x += value * x + assert total_mass > 0, "field has no mass to compute a centroid from" + return weighted_x / total_mass + + +@when( + "it is bootstrapped once after a few real timesteps and again after many more", + target_fixture="centroids", +) +def _when_bootstrapped_at_two_frame_counts(demo: DemoRun) -> tuple[float, float]: + # Two independent real runs (through the same public `bootstrap()` + # entry point every other demo's deeper verification uses), not one + # run read back twice -- there is no pre-step "frame 0" state to + # compare against otherwise, since `RenderWindow.simulation_fields` + # only starts existing once the first `on_frame` call has happened. + 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 (_centroid_x(early_tracer), _centroid_x(late_tracer)) + + +@then( + "the field's mass-weighted centroid has moved downstream by approximately the " + "prescribed velocity times the elapsed time" +) +def _then_centroid_moved_at_the_prescribed_velocity( + demo: DemoRun, centroids: tuple[float, float] +) -> None: + early_x, late_x = centroids + velocity_x, _velocity_y = demo.config.simulation.velocity + dt = demo.config.numerics.timestep + elapsed_steps = _LATE_FRAMES - _EARLY_FRAMES + expected_displacement = velocity_x * dt * elapsed_steps + actual_displacement = late_x - early_x + # Measured directly before choosing this bound (not guessed): a real + # run agrees with the closed-form prediction to within ~4% at this + # demo's own resolution/timestep -- rel=0.15 stays comfortably above + # that margin without being so loose a genuinely broken stepping loop + # (frozen, backwards, or off by a large factor) could still pass. + assert actual_displacement == pytest.approx(expected_displacement, rel=0.15), ( + f"expected the centroid to move roughly {expected_displacement} world units " + f"downstream over {elapsed_steps} real timesteps, moved {actual_displacement} instead" + ) diff --git a/tests/integration/test_cli.py b/tests/integration/test_cli.py index e50c015..882b7ca 100644 --- a/tests/integration/test_cli.py +++ b/tests/integration/test_cli.py @@ -42,7 +42,14 @@ def test_generate_config_prints_valid_yaml_to_stdout() -> None: assert result.returncode == 0, result.stderr parsed = yaml.safe_load(result.stdout) - assert list(parsed.keys()) == ["logging", "rendering", "mesh", "field_display", "numerics"] + assert list(parsed.keys()) == [ + "logging", + "rendering", + "mesh", + "field_display", + "simulation", + "numerics", + ] def test_generate_config_output_writes_file_and_round_trips_through_run( @@ -73,6 +80,7 @@ def test_generate_config_output_writes_file_and_round_trips_through_run( "rendering", "mesh", "field_display", + "simulation", "numerics", ] diff --git a/tests/unit/CLAUDE.md b/tests/unit/CLAUDE.md index c886fb7..9dd4025 100644 --- a/tests/unit/CLAUDE.md +++ b/tests/unit/CLAUDE.md @@ -176,3 +176,31 @@ diffusion scenario and the contract suite's own dedicated test, but leaves the advection scenario passing unchanged, exactly the asymmetry `src/pyflow/engine/CLAUDE.md`'s own `FirstOrderUpwindAdvection`/ `CentralDifferenceDiffusion` entries record. + +**`test_periodic_boundary.py` (TASK-030, added 2026-08-28) is the ninth +and, for this pytest-bdd lineage, last of Stage 4's numerical-scheme +modules**, binding `tests/features/periodic_boundary.feature` -- Stage +4's ninth and last task's own physical-correctness claim. **A genuinely +different shape from every prior binding module in this list**: there is +no condition class under test at all -- periodic bypasses +`BoundaryCondition` entirely (`docs/planning/roadmap.md` TASK-030's own +Design decision), so the only real mechanism is the mesh's own +`wrapped_neighbour_cell` (tested directly in +`test_structured_cartesian_mesh.py`, not here) and the two real interior +schemes' own wiring to it. Own `_Context` dataclass, own local +`_FixedGradientCondition` double for the non-periodic boundary faces the +round-trip scenario's own mesh still needs something configured for +(diffusion has no inflow/outflow carve-out), no golden-demo config file +or CLI run. **The round-trip scenario is checked as convergence under +mesh refinement, not exact equality at one resolution** -- a genuine +numerical finding, not assumed: first-order upwind's own O(dx) numerical +diffusion smooths any field over the distance it travels regardless of +whether the wrap is correct, verified directly (near-identical results +at `num_steps` 10-160 on a fixed mesh, since refining the timestep alone +does not shrink a spatial-truncation-dominated error). A real wrap's own +round-trip error drops by roughly 62% over a 4x mesh refinement; a +mirrored/clamped one (a throwaway mutation, built and run specifically +to check this) drops by only roughly 16% and stays several times larger +throughout -- the scenario's own two-thirds bound was chosen to separate +those two measured outcomes, not guessed, and confirmed to actually fail +under that same mutation before being trusted. diff --git a/tests/unit/numerics/CLAUDE.md b/tests/unit/numerics/CLAUDE.md index 041b332..8b0757d 100644 --- a/tests/unit/numerics/CLAUDE.md +++ b/tests/unit/numerics/CLAUDE.md @@ -274,3 +274,29 @@ tests above it already established -- built against an explicit `"neumann"`-typed config, not the default one, since every default face is `"dirichlet"` (`BoundaryFaceConfig`'s own default has no Neumann face to check by default). + +**`test_advection_contract.py`/`test_diffusion_contract.py` each gained +a second constructor argument at their one existing +`FirstOrderUpwindAdvection`/`CentralDifferenceDiffusion` factory call +(TASK-030, 2026-08-28) -- a real signature-widening join, not a new +fixture.** Both concrete schemes' own `__init__` gained a +`periodic_pairs: Mapping[str, str]` parameter (`docs/planning/ +roadmap.md` TASK-030's own Design decision: a second, separate mapping +threaded alongside `boundary_conditions`); the one existing factory call +in each contract suite passes `{}` (no face configured periodic), since +neither suite's own claims are about periodic behaviour -- that is +`test_periodic_boundary.py`'s (`tests/unit/CLAUDE.md`'s own entry). +**`test_assembly.py` gained two new tests +(`test_advection_and_diffusion_factories_receive_the_resolved_ +periodic_pairs`, `test_a_non_periodic_config_resolves_empty_periodic_ +pairs`) -- the periodic analogue of `_CapturingAdvection`/ +`_CapturingDiffusion`'s own existing boundary-conditions capture +tests.** `AssembledNumerics` has no public `periodic_pairs` field of its +own (only advection/diffusion need it), so this is checked the same +indirect way those two tests already check `boundary_conditions` +threading: through what a real factory actually received. +`_resolve_with_argument` (the one-argument resolve helper advection used +to route through) is deleted in the same change as genuinely dead code, +its only caller having moved to `_resolve_with_two_arguments` -- +`src/pyflow/engine/numerics/CLAUDE.md`'s own entry has the full +reasoning. diff --git a/tests/unit/numerics/test_advection_contract.py b/tests/unit/numerics/test_advection_contract.py index 02dfbc0..edcea0c 100644 --- a/tests/unit/numerics/test_advection_contract.py +++ b/tests/unit/numerics/test_advection_contract.py @@ -56,7 +56,7 @@ def evaluate(self, field: Field, face: int) -> float: def _first_order_upwind_advection() -> FirstOrderUpwindAdvection: condition = _ZeroGradientCondition() return FirstOrderUpwindAdvection( - {"north": condition, "south": condition, "east": condition, "west": condition} + {"north": condition, "south": condition, "east": condition, "west": condition}, {} ) diff --git a/tests/unit/numerics/test_assembly.py b/tests/unit/numerics/test_assembly.py index 2e04b10..27c0266 100644 --- a/tests/unit/numerics/test_assembly.py +++ b/tests/unit/numerics/test_assembly.py @@ -53,8 +53,10 @@ class _TestOnlyAdvection(AdvectionScheme): advection factory now receives (TASK-040's Design decision). """ - def __init__(self, boundary_conditions: Mapping[str, object]) -> None: - del boundary_conditions + def __init__( + self, boundary_conditions: Mapping[str, object], periodic_pairs: Mapping[str, str] + ) -> None: + del boundary_conditions, periodic_pairs def flux(self, field: Field, velocity: VectorField) -> torch.Tensor: self._check_velocity(velocity) @@ -66,8 +68,10 @@ class _OtherTestOnlyAdvection(AdvectionScheme): "a different factory under the same name" is expressible. """ - def __init__(self, boundary_conditions: Mapping[str, object]) -> None: - del boundary_conditions + def __init__( + self, boundary_conditions: Mapping[str, object], periodic_pairs: Mapping[str, str] + ) -> None: + del boundary_conditions, periodic_pairs def flux(self, field: Field, velocity: VectorField) -> torch.Tensor: self._check_velocity(velocity) @@ -75,15 +79,20 @@ def flux(self, field: Field, velocity: VectorField) -> torch.Tensor: class _CapturingAdvection(AdvectionScheme): - """Records the exact `boundary_conditions` mapping it was constructed - with -- every other test-only scheme in this module discards it, so - nothing here would fail if `assemble_numerics` silently passed an - empty mapping (or the wrong one) to the advection factory instead of - the one it just resolved. + """Records the exact `boundary_conditions`/`periodic_pairs` it was + constructed with -- every other test-only scheme in this module + discards them, so nothing here would fail if `assemble_numerics` + silently passed an empty mapping (or the wrong one) to the advection + factory instead of the ones it just resolved. """ - def __init__(self, boundary_conditions: Mapping[str, BoundaryCondition]) -> None: + def __init__( + self, + boundary_conditions: Mapping[str, BoundaryCondition], + periodic_pairs: Mapping[str, str], + ) -> None: self.received_boundary_conditions = boundary_conditions + self.received_periodic_pairs = periodic_pairs def flux(self, field: Field, velocity: VectorField) -> torch.Tensor: self._check_velocity(velocity) @@ -110,17 +119,21 @@ def correct( class _CapturingDiffusion(DiffusionScheme): - """Records the exact `boundary_conditions` mapping and - `diffusion_coefficient` it was constructed with -- the diffusion - analogue of `_CapturingAdvection` above, proving `assemble_numerics` - actually threads both resolved values into the diffusion factory, - not stale or empty ones. + """Records the exact `boundary_conditions` mapping, `periodic_pairs` + mapping, and `diffusion_coefficient` it was constructed with -- the + diffusion analogue of `_CapturingAdvection` above, proving + `assemble_numerics` actually threads all three resolved values into + the diffusion factory, not stale or empty ones. """ def __init__( - self, boundary_conditions: Mapping[str, BoundaryCondition], diffusion_coefficient: float + self, + boundary_conditions: Mapping[str, BoundaryCondition], + periodic_pairs: Mapping[str, str], + diffusion_coefficient: float, ) -> None: self.received_boundary_conditions = boundary_conditions + self.received_periodic_pairs = periodic_pairs self.received_diffusion_coefficient = diffusion_coefficient def flux(self, field: Field) -> torch.Tensor: @@ -314,6 +327,44 @@ def test_diffusion_factory_receives_the_resolved_boundary_conditions_and_coeffic assert assembled.diffusion.received_diffusion_coefficient == 3.5 +def test_advection_and_diffusion_factories_receive_the_resolved_periodic_pairs() -> None: + # `AssembledNumerics` has no public `periodic_pairs` field of its own + # (only advection/diffusion need it, TASK-030's own Design decision), + # so this is checked indirectly through what a real factory actually + # received -- the periodic analogue of the two capture tests above. + advection_name = "test_only_capturing_advection_for_periodic_pairs_test" + diffusion_name = "test_only_capturing_diffusion_for_periodic_pairs_test" + register_advection_scheme(advection_name, _CapturingAdvection) + register_diffusion_scheme(diffusion_name, _CapturingDiffusion) + config = NumericsConfig( + advection=advection_name, # type: ignore[arg-type] + diffusion=diffusion_name, # type: ignore[arg-type] + boundary_conditions=BoundaryConditionsConfig( + north=BoundaryFaceConfig(type="periodic", velocity=None, pressure=None), + south=BoundaryFaceConfig(type="periodic", velocity=None, pressure=None), + ), + ) + + assembled = assemble_numerics(config) + + expected = {"north": "south", "south": "north"} + assert isinstance(assembled.advection, _CapturingAdvection) + assert dict(assembled.advection.received_periodic_pairs) == expected + assert isinstance(assembled.diffusion, _CapturingDiffusion) + assert dict(assembled.diffusion.received_periodic_pairs) == expected + + +def test_a_non_periodic_config_resolves_empty_periodic_pairs() -> None: + advection_name = "test_only_capturing_advection_for_non_periodic_test" + register_advection_scheme(advection_name, _CapturingAdvection) + config = NumericsConfig(advection=advection_name) # type: ignore[arg-type] + + assembled = assemble_numerics(config) + + assert isinstance(assembled.advection, _CapturingAdvection) + assert dict(assembled.advection.received_periodic_pairs) == {} + + # -- The reference ("null") implementations' own behaviour ----------------- # # `assemble_numerics` only proves these construct; the module docstring's diff --git a/tests/unit/numerics/test_diffusion_contract.py b/tests/unit/numerics/test_diffusion_contract.py index 42486b2..8102482 100644 --- a/tests/unit/numerics/test_diffusion_contract.py +++ b/tests/unit/numerics/test_diffusion_contract.py @@ -48,6 +48,7 @@ def _central_difference_diffusion() -> CentralDifferenceDiffusion: condition = _ZeroGradientCondition() return CentralDifferenceDiffusion( {"north": condition, "south": condition, "east": condition, "west": condition}, + {}, 1.0, ) diff --git a/tests/unit/test_central_difference_diffusion.py b/tests/unit/test_central_difference_diffusion.py index 966f7e0..28c4103 100644 --- a/tests/unit/test_central_difference_diffusion.py +++ b/tests/unit/test_central_difference_diffusion.py @@ -106,7 +106,7 @@ def _west_face(mesh: StructuredCartesianMesh) -> int: def _run_flux(ctx: _Context) -> None: - scheme = CentralDifferenceDiffusion(ctx.boundary_conditions, _GAMMA) + scheme = CentralDifferenceDiffusion(ctx.boundary_conditions, {}, _GAMMA) try: ctx.flux = scheme.flux(ctx.scalar) except UnconfiguredBoundaryFaceError as exc: @@ -237,7 +237,7 @@ def _when_measure_convergence(ctx: _Context) -> None: "east": condition, "west": condition, } - scheme = CentralDifferenceDiffusion(boundary_conditions, _GAMMA) + scheme = CentralDifferenceDiffusion(boundary_conditions, {}, _GAMMA) discrete = accumulate_flux_to_cells(mesh, scheme.flux(field_values)) max_error = 0.0 @@ -261,7 +261,7 @@ def _when_advanced_many(ctx: _Context) -> None: def _advance(ctx: _Context, steps: int) -> None: assert ctx.dt is not None assert ctx.history is not None - scheme = CentralDifferenceDiffusion(ctx.boundary_conditions, _GAMMA) + scheme = CentralDifferenceDiffusion(ctx.boundary_conditions, {}, _GAMMA) scalar = ctx.scalar for _ in range(steps): flux = scheme.flux(scalar) diff --git a/tests/unit/test_configuration.py b/tests/unit/test_configuration.py index d408a88..129fa87 100644 --- a/tests/unit/test_configuration.py +++ b/tests/unit/test_configuration.py @@ -32,6 +32,9 @@ def test_defaults_are_valid() -> None: assert config.field_display.arrow_color == "#ffffff" assert config.field_display.arrow_scale == 0.3 assert config.field_display.show_legend is True + assert config.simulation.scalar_pattern is None + assert config.simulation.velocity_pattern is None + assert config.simulation.velocity == (1.0, 0.0) assert config.numerics.advection == "first_order_upwind" assert config.numerics.diffusion == "central_difference" assert config.numerics.diffusion_coefficient == 1.0 @@ -430,6 +433,50 @@ def test_load_config_rejects_a_non_boolean_show_legend(tmp_path: Path) -> None: load_config(config_file) +# -- SimulationConfig (TASK-030) ------------------------------------------ + + +def test_load_config_reads_simulation_section(tmp_path: Path) -> None: + config_file = tmp_path / "config.yaml" + config_file.write_text( + "simulation:\n" + " scalar_pattern: gaussian_blob\n" + " velocity_pattern: uniform\n" + " velocity: [2.0, -1.5]\n" + ) + + config = load_config(config_file) + + assert config.simulation.scalar_pattern == "gaussian_blob" + assert config.simulation.velocity_pattern == "uniform" + assert config.simulation.velocity == (2.0, -1.5) + assert isinstance(config.simulation.velocity, tuple) + + +def test_load_config_rejects_an_unknown_simulation_scalar_pattern(tmp_path: Path) -> None: + config_file = tmp_path / "config.yaml" + config_file.write_text("simulation:\n scalar_pattern: checkerboard\n") + + with pytest.raises(ValueError, match="scalar_pattern"): + load_config(config_file) + + +def test_load_config_rejects_an_unknown_simulation_velocity_pattern(tmp_path: Path) -> None: + config_file = tmp_path / "config.yaml" + config_file.write_text("simulation:\n velocity_pattern: spiral\n") + + with pytest.raises(ValueError, match="velocity_pattern"): + load_config(config_file) + + +def test_load_config_rejects_a_non_numeric_simulation_velocity(tmp_path: Path) -> None: + config_file = tmp_path / "config.yaml" + config_file.write_text("simulation:\n velocity: [not-a-number, 0.0]\n") + + with pytest.raises(ValueError, match="simulation.velocity"): + load_config(config_file) + + @pytest.mark.parametrize( ("yaml_text", "expected"), [ diff --git a/tests/unit/test_conjugate_gradient_solver.py b/tests/unit/test_conjugate_gradient_solver.py index 55ac912..86dc942 100644 --- a/tests/unit/test_conjugate_gradient_solver.py +++ b/tests/unit/test_conjugate_gradient_solver.py @@ -77,7 +77,7 @@ def _build_semidefinite_matrix(mesh: StructuredCartesianMesh) -> torch.Tensor: "east": condition, "west": condition, } - diffusion = CentralDifferenceDiffusion(boundary_conditions, diffusion_coefficient=1.0) + diffusion = CentralDifferenceDiffusion(boundary_conditions, {}, diffusion_coefficient=1.0) n = mesh.num_cells matrix = torch.zeros((n, n), dtype=torch.float64) diff --git a/tests/unit/test_dirichlet_boundary.py b/tests/unit/test_dirichlet_boundary.py index b7639af..f3b6691 100644 --- a/tests/unit/test_dirichlet_boundary.py +++ b/tests/unit/test_dirichlet_boundary.py @@ -151,13 +151,13 @@ def _given_real_dirichlet_for_diffusion(ctx: _Context) -> None: @when("the advective flux is computed using this condition") def _when_advective_flux(ctx: _Context) -> None: - scheme = FirstOrderUpwindAdvection(ctx.boundary_conditions) + scheme = FirstOrderUpwindAdvection(ctx.boundary_conditions, {}) ctx.flux = scheme.flux(ctx.scalar, ctx.velocity) @when("the diffusive flux is computed using this condition") def _when_diffusive_flux(ctx: _Context) -> None: - scheme = CentralDifferenceDiffusion(ctx.boundary_conditions, _GAMMA) + scheme = CentralDifferenceDiffusion(ctx.boundary_conditions, {}, _GAMMA) ctx.flux = scheme.flux(ctx.scalar) diff --git a/tests/unit/test_first_order_upwind_advection.py b/tests/unit/test_first_order_upwind_advection.py index 41c0861..cba30f8 100644 --- a/tests/unit/test_first_order_upwind_advection.py +++ b/tests/unit/test_first_order_upwind_advection.py @@ -98,7 +98,7 @@ def _zero_gradient_everywhere() -> dict[str, BoundaryCondition]: def _run_flux(ctx: _Context) -> None: - scheme = FirstOrderUpwindAdvection(ctx.boundary_conditions) + scheme = FirstOrderUpwindAdvection(ctx.boundary_conditions, {}) try: ctx.flux = scheme.flux(ctx.scalar, ctx.velocity) except UnconfiguredBoundaryFaceError as exc: @@ -269,7 +269,7 @@ def _when_advanced_many(ctx: _Context) -> None: def _advance(ctx: _Context, steps: int) -> None: assert ctx.dt is not None assert ctx.history is not None - scheme = FirstOrderUpwindAdvection(ctx.boundary_conditions) + scheme = FirstOrderUpwindAdvection(ctx.boundary_conditions, {}) scalar = ctx.scalar for _ in range(steps): flux = scheme.flux(scalar, ctx.velocity) diff --git a/tests/unit/test_generator.py b/tests/unit/test_generator.py index 1b1b992..6c16810 100644 --- a/tests/unit/test_generator.py +++ b/tests/unit/test_generator.py @@ -140,8 +140,8 @@ def test_top_level_key_order_matches_pyflowconfig_field_order() -> None: schema's own declared field order -- it only proves the dict iteration order was preserved. Check the parsed keys directly against `PyFlowConfig`'s declared order (`logging`, `rendering`, - `mesh`, `field_display`, `numerics`), not assumed from the dumper - flag. + `mesh`, `field_display`, `simulation`, `numerics`), not assumed from + the dumper flag. """ text = generate_config_yaml(PyFlowConfig()) @@ -152,5 +152,6 @@ def test_top_level_key_order_matches_pyflowconfig_field_order() -> None: "rendering", "mesh", "field_display", + "simulation", "numerics", ] diff --git a/tests/unit/test_main.py b/tests/unit/test_main.py index 7434e9d..8b7b6cd 100644 --- a/tests/unit/test_main.py +++ b/tests/unit/test_main.py @@ -97,6 +97,11 @@ def test_generate_config_with_no_output_prints_to_stdout( "arrow_scale": 0.3, "show_legend": True, }, + "simulation": { + "scalar_pattern": None, + "velocity_pattern": None, + "velocity": [1.0, 0.0], + }, "numerics": { "advection": "first_order_upwind", "diffusion": "central_difference", @@ -151,5 +156,12 @@ def test_generate_config_with_output_writes_file_and_prints_nothing( captured = capsys.readouterr() assert captured.out == "" written = yaml.safe_load(output_path.read_text()) - assert list(written.keys()) == ["logging", "rendering", "mesh", "field_display", "numerics"] + assert list(written.keys()) == [ + "logging", + "rendering", + "mesh", + "field_display", + "simulation", + "numerics", + ] assert written["mesh"]["extent"] == list(PyFlowConfig().mesh.extent) diff --git a/tests/unit/test_neumann_boundary.py b/tests/unit/test_neumann_boundary.py index 85f5c0f..343ad2d 100644 --- a/tests/unit/test_neumann_boundary.py +++ b/tests/unit/test_neumann_boundary.py @@ -146,13 +146,13 @@ def _given_inflow_boundary(ctx: _Context) -> None: @when("the advective flux is computed using this condition") def _when_advective_flux(ctx: _Context) -> None: - scheme = FirstOrderUpwindAdvection(ctx.boundary_conditions) + scheme = FirstOrderUpwindAdvection(ctx.boundary_conditions, {}) ctx.flux = scheme.flux(ctx.scalar, ctx.velocity) @when("the diffusive flux is computed using this condition") def _when_diffusive_flux(ctx: _Context) -> None: - scheme = CentralDifferenceDiffusion(ctx.boundary_conditions, _GAMMA) + scheme = CentralDifferenceDiffusion(ctx.boundary_conditions, {}, _GAMMA) ctx.flux = scheme.flux(ctx.scalar) diff --git a/tests/unit/test_periodic_boundary.py b/tests/unit/test_periodic_boundary.py new file mode 100644 index 0000000..1d1d3a2 --- /dev/null +++ b/tests/unit/test_periodic_boundary.py @@ -0,0 +1,329 @@ +"""Binds `tests/features/periodic_boundary.feature` (TASK-030) -- Stage +4's ninth and last task. Periodic bypasses `BoundaryCondition` entirely +(`docs/planning/roadmap.md` TASK-030's own Design decision): a real +Advection/Diffusion scheme wired with a `periodic_pairs` mapping must +consult the wrapped-neighbour cell directly, through the same formula it +already uses for a genuine interior neighbour, never a `BoundaryCondition` +-- so unlike `test_dirichlet_boundary.py`/`test_neumann_boundary.py`, +there is no condition class under test here at all, only the mesh's own +`wrapped_neighbour_cell` (`test_structured_cartesian_mesh.py`) and the two +real schemes' own wiring to it. + +Not a golden demo -- no config file under `examples/golden-demos/`, no +CLI run. Lives here, not under `tests/golden/`, per this directory's own +scope: isolated logic, no process boundary (`tests/unit/CLAUDE.md`). +""" + +from __future__ import annotations + +import math +from dataclasses import dataclass +from typing import Literal + +import torch +from pytest_bdd import given, scenarios, then, when + +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 +from pyflow.engine.numerics.boundary_condition import BoundaryCondition +from pyflow.engine.numerics.diffusion import CentralDifferenceDiffusion +from pyflow.engine.numerics.linear_solver import LinearSolver, LinearSolverResult +from pyflow.engine.numerics.pressure_coupling import PressureCoupling +from pyflow.engine.numerics.time_integrator import RK4Integrator +from pyflow.engine.scalar_field import ScalarField +from pyflow.engine.simulation import step as simulation_step +from pyflow.engine.vector_field import VectorField + +scenarios("periodic_boundary.feature") + +_GAMMA = 2.0 +"""Not `1.0` -- same reasoning as every other diffusion scenario's own +identically-named constant in this repository. +""" + + +class _FixedGradientCondition(BoundaryCondition): + """Zero-gradient, for the non-periodic boundary faces the round-trip + scenario's own mesh still needs *something* configured for (diffusion + has no inflow/outflow carve-out) -- never the mechanism under test. + """ + + def __init__(self, gradient: float = 0.0) -> None: + self._gradient = gradient + + @property + def kind(self) -> Literal["value", "gradient"]: + return "gradient" + + def evaluate(self, field: Field, face: int) -> float: + self._check_boundary_face(field, face) + return self._gradient + + +class _InertLinearSolver(LinearSolver): + """`AssembledNumerics` requires one to construct at all; `step` never + calls it -- same precedent as `test_simulation.py`'s own double. + """ + + def solve(self, matrix: torch.Tensor, rhs: torch.Tensor) -> LinearSolverResult: + return LinearSolverResult(solution=torch.zeros_like(rhs), converged=False, iterations=0) + + +class _InertPressureCoupling(PressureCoupling): + """Same reasoning as `_InertLinearSolver` above.""" + + def correct( + self, provisional_velocity: VectorField, dt: float + ) -> tuple[VectorField, ScalarField]: + del dt + return provisional_velocity.copy(), ScalarField(provisional_velocity.mesh, "pressure") + + +@dataclass +class _Context: + mesh: StructuredCartesianMesh + scalar: ScalarField + velocity: VectorField + boundary_conditions: dict[str, BoundaryCondition] + periodic_pairs: dict[str, str] + flux: torch.Tensor | None = None + target_face: int | None = None + wrapped_neighbour: int | None = None + velocity_x: float | None = None + round_trip_errors: tuple[float, float] | None = None + + +def _mesh() -> StructuredCartesianMesh: + # Non-"nice" origin/spacing and a non-square extent, matching every + # other contract suite's fixture in this repository. + return StructuredCartesianMesh(origin=(0.5, -1.0), spacing=(0.2, 0.3), extent=(4, 2)) + + +def _distinct_scalar(mesh: StructuredCartesianMesh) -> ScalarField: + # Every cell gets a genuinely different value -- a wrong wrapped-cell + # id (or a mirrored/clamped-to-owner fallback) would read a + # coincidentally-equal value for at most one cell, never for every + # face this file exercises. + return ScalarField(mesh, "tracer", initial_value=lambda x, y: 10 * x + 100 * y) + + +def _west_face(mesh: StructuredCartesianMesh) -> int: + return next(f for f in range(mesh.num_faces) if mesh.boundary_face_name(f) == "west") + + +def _face_normal_velocity(ctx: _Context, face: int, neighbour: int | None) -> float: + """Independently derived, not calling into either scheme under test -- + same reasoning `test_neumann_boundary.py`'s own identically-named + helper states. + """ + owner, _ = ctx.mesh.face_neighbours(face) + normal_x, normal_y = ctx.mesh.face_normal(face) + owner_x, owner_y = ctx.velocity.value_at(owner) + if neighbour is None: + v_x, v_y = owner_x, owner_y + else: + neighbour_x, neighbour_y = ctx.velocity.value_at(neighbour) + v_x, v_y = (owner_x + neighbour_x) / 2, (owner_y + neighbour_y) / 2 + return v_x * normal_x + v_y * normal_y + + +# -- Given ------------------------------------------------------------- + + +@given( + "a small, non-square, non-trivially-origined mesh whose cells each hold a distinct value", + target_fixture="ctx", +) +def _given_default_mesh() -> _Context: + mesh = _mesh() + velocity = VectorField(mesh, "velocity", num_components=2, initial_value=(0.0, 0.0)) + scalar = _distinct_scalar(mesh) + return _Context( + mesh=mesh, + scalar=scalar, + velocity=velocity, + boundary_conditions={ + "north": _FixedGradientCondition(), + "south": _FixedGradientCondition(), + "east": _FixedGradientCondition(), + }, + periodic_pairs={}, + ) + + +@given("a west boundary face configured periodic with its east partner") +def _given_periodic_west(ctx: _Context) -> None: + ctx.target_face = _west_face(ctx.mesh) + ctx.periodic_pairs = {"west": "east", "east": "west"} + ctx.wrapped_neighbour = ctx.mesh.wrapped_neighbour_cell(ctx.target_face) + + +@given("no boundary condition configured for that face at all") +def _given_no_condition_for_that_face(ctx: _Context) -> None: + ctx.boundary_conditions = { + name: condition for name, condition in ctx.boundary_conditions.items() if name != "west" + } + + +@given("a mesh whose east and west edges are both periodic") +def _given_fully_periodic_mesh(ctx: _Context) -> None: + ctx.periodic_pairs = {"west": "east", "east": "west"} + + +@given( + "a uniform velocity field that carries the domain's own scalar field once fully around it " + "in a whole number of real timesteps" +) +def _given_uniform_round_trip_velocity(ctx: _Context) -> None: + ctx.velocity_x = 1.0 + + +# -- When -------------------------------------------------------------- + + +@when("the advective flux is computed with inflow at that face") +def _when_advective_flux(ctx: _Context) -> None: + scheme = FirstOrderUpwindAdvection(ctx.boundary_conditions, ctx.periodic_pairs) + # west's canonical normal is (-1, 0); velocity (+1, 0) gives + # velocity_normal = 1*(-1) = -1 -- inflow, same reasoning as + # `test_neumann_boundary.py`'s own identically-shaped step. + ctx.velocity = VectorField(ctx.mesh, "velocity", num_components=2, initial_value=(1.0, 0.0)) + ctx.flux = scheme.flux(ctx.scalar, ctx.velocity) + + +@when("the diffusive flux is computed at that face") +def _when_diffusive_flux(ctx: _Context) -> None: + scheme = CentralDifferenceDiffusion(ctx.boundary_conditions, ctx.periodic_pairs, _GAMMA) + ctx.flux = scheme.flux(ctx.scalar) + + +def _round_trip_error( + nx: int, num_steps: int, velocity_x: float, periodic_pairs: dict[str, str] +) -> float: + """One lap around a periodic domain at resolution `nx`, `num_steps` + real RK4 timesteps (chosen equal to `nx`, i.e. Courant number ~1 -- + verified numerically that refining time alone barely moves this + error, since it is dominated by first-order upwind's own spatial + truncation, not time-integration error), returning the max absolute + difference between the field's final and starting values. + + A fresh sine-in-x pattern (period equal to the domain width, so it is + genuinely periodic-compatible -- unlike a plain linear ramp, which + creates an artificial discontinuity at the wrap seam that numerical + diffusion then smooths, confounding this specific claim, found by + running exactly that fixture first) plus a `100 * y` term, so every + row is checked, not only one. + """ + origin = (0.5, -1.0) + spacing = (0.8 / nx, 0.3) + mesh = StructuredCartesianMesh(origin=origin, spacing=spacing, extent=(nx, 2)) + dx, _ = spacing + domain_width = nx * dx + x0 = origin[0] + + def pattern(x: float, y: float) -> float: + return math.sin(2 * math.pi * (x - x0) / domain_width) + 100 * y + + scalar = ScalarField(mesh, "tracer", initial_value=pattern) + velocity = VectorField(mesh, "velocity", num_components=2, initial_value=(velocity_x, 0.0)) + solver = _InertLinearSolver() + numerics = AssembledNumerics( + advection=FirstOrderUpwindAdvection({}, periodic_pairs), + diffusion=_ZeroDiffusion(), + time_integration=RK4Integrator(), + linear_solver=solver, + pressure_coupling=_InertPressureCoupling(solver), + boundary_conditions={}, + names={}, + ) + dt = domain_width / (velocity_x * num_steps) + fields: dict[str, Field] = {"tracer": scalar} + for _ in range(num_steps): + fields = simulation_step(fields, velocity, numerics, dt) + tracer = fields["tracer"] + assert isinstance(tracer, ScalarField) + return float((tracer.values - scalar.values).abs().max()) + + +@when( + "the field is advected for exactly that many real timesteps, at two mesh resolutions " + "four times apart" +) +def _when_advected_at_two_resolutions(ctx: _Context) -> None: + assert ctx.velocity_x is not None + coarse = _round_trip_error(16, 16, ctx.velocity_x, ctx.periodic_pairs) + fine = _round_trip_error(64, 64, ctx.velocity_x, ctx.periodic_pairs) + ctx.round_trip_errors = (coarse, fine) + + +# -- Then ---------------------------------------------------------------- + + +@then( + "the inflow boundary face's implied value is the wrapped neighbour's own value, not the owner's" +) +def _then_advection_reads_wrapped_neighbour(ctx: _Context) -> None: + assert ctx.flux is not None + assert ctx.target_face is not None + assert ctx.wrapped_neighbour is not None + owner, _ = ctx.mesh.face_neighbours(ctx.target_face) + velocity_normal = _face_normal_velocity(ctx, ctx.target_face, ctx.wrapped_neighbour) + implied = float(ctx.flux[ctx.target_face]) / velocity_normal + assert implied == ctx.scalar.value_at(ctx.wrapped_neighbour) + assert implied != ctx.scalar.value_at(owner) + + +@then( + "that boundary face's flux equals the diffusion coefficient times the difference " + "between the wrapped neighbour's and the owner's own value, divided by one full " + "grid spacing" +) +def _then_diffusion_uses_wrapped_neighbour_at_full_spacing(ctx: _Context) -> None: + assert ctx.flux is not None + assert ctx.target_face is not None + assert ctx.wrapped_neighbour is not None + owner, _ = ctx.mesh.face_neighbours(ctx.target_face) + distance = 2 * ctx.mesh.face_centroid_distance(ctx.target_face) + expected = ( + _GAMMA + * (ctx.scalar.value_at(ctx.wrapped_neighbour) - ctx.scalar.value_at(owner)) + / distance + ) + assert math.isclose(float(ctx.flux[ctx.target_face]), expected, abs_tol=1e-9) + + +@then( + "the round-trip error at the finer resolution is well under two thirds of the " + "error at the coarser one" +) +def _then_error_shrinks_with_refinement(ctx: _Context) -> None: + assert ctx.round_trip_errors is not None + coarse, fine = ctx.round_trip_errors + assert fine < (2 / 3) * coarse, ( + f"fine-resolution round-trip error {fine} did not drop meaningfully below " + f"two thirds of the coarse-resolution error {coarse} -- a wrapped-neighbour " + "implementation should converge toward the starting distribution as the mesh " + "is refined; a mirrored/clamped one measurably does not (see this scenario's " + "own comment in periodic_boundary.feature)" + ) + + +# -- Local doubles (defined after use above is fine at import time; kept +# near the bottom since only the round-trip scenario needs it) ---------- + + +class _ZeroDiffusion(CentralDifferenceDiffusion): + """No diffusion at all -- the round-trip scenario is specifically an + *advection* claim (`periodic_boundary.feature`'s own comment): + diffusion's own periodic wiring is already checked directly by the + scenario above, and diffusion is irreversible, so including a real + one here would falsify this claim rather than test it. + """ + + def __init__(self) -> None: + super().__init__({}, {}, 1.0) + + def flux(self, field: Field) -> torch.Tensor: + return torch.zeros(field.mesh.num_faces, dtype=torch.float64) diff --git a/tests/unit/test_structured_cartesian_mesh.py b/tests/unit/test_structured_cartesian_mesh.py index d178d1d..20facab 100644 --- a/tests/unit/test_structured_cartesian_mesh.py +++ b/tests/unit/test_structured_cartesian_mesh.py @@ -16,7 +16,11 @@ from pyflow.configuration.schema import MeshConfig from pyflow.engine.coordinate_system import UniformVertexCoordinateSystem -from pyflow.engine.mesh import InvalidMeshEntityError, StructuredCartesianMesh +from pyflow.engine.mesh import ( + InvalidMeshEntityError, + NotABoundaryFaceError, + StructuredCartesianMesh, +) _ORIGIN = (1.5, -2.25) _SPACING = (0.1, 0.3) @@ -222,6 +226,56 @@ def test_interior_face_centroid_distance_equals_the_grid_spacing() -> None: assert math.isclose(mesh.face_centroid_distance(face), expected) +def test_wrapped_neighbour_cell_pairs_each_boundary_face_with_the_opposite_edge_cell() -> None: + # TASK-030: the west/east boundary faces of the same row wrap to each + # other's owner cell, and likewise north/south of the same column -- + # a periodic domain's own "the same relative position on the opposite + # edge" reading, checked directly rather than assumed. + mesh = _mesh() + nx, ny = _EXTENT + + for i, j in itertools.product(range(nx), range(ny)): + cell = mesh.cell_id(i, j) + west, east, south, north = mesh.cell_faces(cell) + + if i == 0: + assert mesh.wrapped_neighbour_cell(west) == mesh.cell_id(nx - 1, j) + if i == nx - 1: + assert mesh.wrapped_neighbour_cell(east) == mesh.cell_id(0, j) + if j == 0: + assert mesh.wrapped_neighbour_cell(south) == mesh.cell_id(i, ny - 1) + if j == ny - 1: + assert mesh.wrapped_neighbour_cell(north) == mesh.cell_id(i, 0) + + +def test_wrapped_neighbour_cell_rejects_an_interior_vertical_face() -> None: + mesh = _mesh() + nx, ny = _EXTENT + interior_cell = mesh.cell_id(nx // 2, ny // 2) + west, east, _south, _north = mesh.cell_faces(interior_cell) + interior_vertical_face = next(f for f in (west, east) if not mesh.is_boundary_face(f)) + + with pytest.raises(NotABoundaryFaceError): + mesh.wrapped_neighbour_cell(interior_vertical_face) + + +def test_wrapped_neighbour_cell_rejects_an_interior_horizontal_face() -> None: + mesh = _mesh() + nx, ny = _EXTENT + interior_cell = mesh.cell_id(nx // 2, ny // 2) + _west, _east, south, north = mesh.cell_faces(interior_cell) + interior_horizontal_face = next(f for f in (south, north) if not mesh.is_boundary_face(f)) + + with pytest.raises(NotABoundaryFaceError): + mesh.wrapped_neighbour_cell(interior_horizontal_face) + + +def test_wrapped_neighbour_cell_rejects_an_out_of_range_face() -> None: + mesh = _mesh() + with pytest.raises(InvalidMeshEntityError): + mesh.wrapped_neighbour_cell(mesh.num_faces) + + def test_boundary_face_centroid_distance_is_half_the_grid_spacing() -> None: # A boundary face sits exactly midway between where its owner's # centroid is and where a (non-existent) neighbour's would be, so