diff --git a/README.md b/README.md index d6b1be4..2c3140b 100644 --- a/README.md +++ b/README.md @@ -143,7 +143,15 @@ timestep against the stability limit, so refining a shipped demo's mesh blows up silently at step 17; and `BoundaryFaceConfig.velocity` is a validated, documented configuration field that no engine code reads. See `docs/planning/roadmap.md`'s own Stage 9 section for the criteria -and the measurements. It is placed before Better Numerics by dependency, +and the measurements. **TASK-052 has landed** -- a sealed box now +conserves a purely advected tracer exactly, where it lost 14.27% before, +and the lid-driven cavity's own error against Ghia, Ghia & Shin (1982) +*fell at every resolution* (9x9 0.1433 -> 0.1292, 13x13 0.0874 -> +0.0766, 17x17 0.0578 -> 0.0524), which is independent evidence the +change was physics rather than a re-fitted tolerance. Its own demo is +`uv run python -m pyflow run --demos sealed_box`. TASK-053 (a failed +frame failing the run) and TASK-054 (the timestep stability warning) +are drafted and not yet built. It is placed before Better Numerics by dependency, not preference -- Stage 10's own Rayleigh-Bénard criterion measures convection between heated walls, which is not meaningful while those walls leak. diff --git a/adr/ADR-009-pressure-coupling-dt.md b/adr/ADR-009-pressure-coupling-dt.md index 018d7f1..1cfdae7 100644 --- a/adr/ADR-009-pressure-coupling-dt.md +++ b/adr/ADR-009-pressure-coupling-dt.md @@ -10,9 +10,7 @@ TASK-021 (Stage 3, done 2026-08-23) gave `PressureCoupling` this abstract method: ```python -def correct( - self, provisional_velocity: VectorField -) -> tuple[VectorField, ScalarField]: ... +def correct(self, provisional_velocity: VectorField) -> tuple[VectorField, ScalarField]: ... ``` No timestep. That was sufficient through Stage 3, where no concrete diff --git a/docs/architecture/icds.md b/docs/architecture/icds.md index a8c41d8..370bc47 100644 --- a/docs/architecture/icds.md +++ b/docs/architecture/icds.md @@ -140,6 +140,17 @@ limit (`docs/handbook/numerical-methods/advection.md`, `time-integration.md`). Appropriate for MVP correctness validation, not for accuracy-sensitive production use. +**Boundary treatment, corrected 2026-09-12 (TASK-052, Stage 9):** a +genuine boundary face's *transporting* velocity is resolved from that +face's own boundary condition, through the same +`boundary_normal_velocity` function `GreenGaussDivergence` resolves it +through -- not from the cell inside the face, which is what this scheme +did until then, and which let material cross solid walls at up to 42% of +a flow's own peak speed. The transported field's own boundary value is a +separate question, answered the same way it always was (a Dirichlet +condition's prescribed value at inflow, zero-order extrapolation at +outflow or under Neumann). + **Limitations:** first-order accuracy only; smooths sharp gradients more than a user comparing against a higher-order reference might expect. The artificial diffusivity this amounts to is roughly $\rho |u| \Delta x / 2$ @@ -458,6 +469,22 @@ surface, which no Stage 5 task discharged. **Expected behaviour:** each condition type supplies the face value (Dirichlet), face gradient (Neumann), or wrapped-neighbour reference (periodic) the interior advection/diffusion schemes need at that face. +**Since TASK-052 (Stage 9, 2026-09-12) it also supplies the face's own +normal *velocity*** -- the number every operator that transports across +that face resolves through `boundary_normal_velocity` +(`src/pyflow/engine/numerics/boundary_condition.py`). A Dirichlet +condition with no per-field override for the velocity field prescribes +`0.0` there, a no-penetration wall; it used to fall through to +`scalar_value`, a transported scalar's boundary value, which is not a +velocity. + +**A fifth requirement, still open: `BoundaryFaceConfig.velocity` reaches +no scheme.** It is validated (mutual exclusivity with `pressure`, zero +net flux) and then read by nothing -- the per-component channel +`field_values["velocity.0"]`/`["velocity.1"]` is what a real +configuration uses and what the engine reads. Recorded 2026-09-12 as +Stage 9 Completion Criterion 3's remaining half rather than fixed in +the same change as the defect above. **Limitations:** limited to simple, axis-aligned domain edges -- internal boundaries and arbitrary-geometry surfaces are explicitly future work, diff --git a/docs/implementation/golden-demos.md b/docs/implementation/golden-demos.md index fc1360f..8680dad 100644 --- a/docs/implementation/golden-demos.md +++ b/docs/implementation/golden-demos.md @@ -591,6 +591,51 @@ scene (`tests/integration/test_playback_cli.py:: test_play_renders_both_arrows_and_a_panel_for_smoke_transport`), not only that the CLI exits cleanly. +## Sealed Box + +TASK-052's own golden demo, and Stage 9's -- the visible form of the +defect that opened that stage. A tracer stirred by a real, solved +lid-driven flow inside a closed no-slip cavity has nowhere to go, and +must still be entirely inside the cavity however long it is stirred. + +**The second demo, after Multi-Field Plume, that exists because an audit +found something wrong rather than because a task built a capability.** +Before TASK-052, `FirstOrderUpwindAdvection` took a boundary face's +transporting velocity from the cell inside it rather than from what the +configuration prescribed there, so material crossed solid walls at up to +42% of the flow's own peak speed -- 14.27% of a purely advected tracer +lost over 400 steps on `smoke_transport.yaml`'s own geometry, whose mesh, +timestep, viscosity and moving lid this demo shares. + +"Working" means, concretely: + +- the demo *is* `examples/golden-demos/sealed_box.yaml`, run via + `uv run python -m pyflow run --config examples/golden-demos/sealed_box.yaml` + (or `--demos sealed_box`); +- **every wall face carries exactly zero advective flux, while the cells + against those walls are still moving** -- both halves, checked through + `bootstrap()` on the committed config + (`tests/golden/test_sealed_box.py`). The second half is what stops the + first being satisfied by a flow that simply stopped; +- the tracer is genuinely carried by the recirculating flow, the same + "bootstrap at two frame counts, compare" shape + `test_smoke_transport.py` already established; +- it runs headlessly via `--backend offscreen`, same as every other demo. + +**This demo checks the wall, not the domain integral, and that is a real +constraint rather than a preference.** Diffusion to a zero-valued wall +removes tracer legitimately, so the integral is not exactly constant +here; a scenario asserting it was would be wrong. Switching diffusion off +entirely would make exact conservation the right claim, and +`FieldConfig` rejects a non-positive `diffusion_coefficient` -- a +rejection that is itself a bound Stage 6 acceptance criterion +(`tests/features/field_declaration.feature`). Relaxing it for one demo's +convenience would reopen another stage's criterion, so exact +conservation is proven at the engine level instead, against a purely +advective fixture: `tests/features/boundary_velocity.feature`'s own +sealed-domain scenario, which measures `+0.000000000000%` against a +pre-fix `-4.54%`. + ## Future Demos Add an entry here when a new capability is implemented, per diff --git a/docs/planning/roadmap.md b/docs/planning/roadmap.md index f5f48cc..f9580c1 100644 --- a/docs/planning/roadmap.md +++ b/docs/planning/roadmap.md @@ -381,7 +381,11 @@ 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): **1209 tests as of 2026-09-11**, up from 1200 on 2026-09-09 +(C1a/C1b): **1226 tests as of 2026-09-12**, up from 1209 on 2026-09-11 +(TASK-052, Stage 9: 9 in `tests/unit/test_boundary_velocity.py` and 3 in +`tests/golden/test_sealed_box.py` for the wall-permeability fix and its +own golden demo, plus 5 the fixtures those changed gained along the way). +Before that, 1209 on 2026-09-11 (the Stage 8 exit audit, closing the gap between three Completion Criteria and what actually checked them: 4 in `tests/integration/test_import_order.py` (Criterion 1's own "imports @@ -971,7 +975,7 @@ properties `PISO` (TASK-027, Stage 4) already computed but Stage 4's own criteria never had cause to check -- constant pressure for a divergence-free provisional field, the null-space remedy actually holding, `step` rejecting a `PressureField` -- against the real `PISO` -class throughout, no new pressure-solving mechanism. **144 of those 1003 +class throughout, no new pressure-solving mechanism. **156 of those 1003 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` @@ -9641,7 +9645,7 @@ Criterion 10 only ever asked for the first. | 3. Every field carries its own physical coefficients, through the mechanisms that already exist | **Met, with one qualification the criterion's own last bullet makes worth stating.** Two fields at two diffusivities, each decaying at its own analytic rate rather than merely differing: `tests/features/humidity_field.feature`'s first scenario -- the first time `CentralDifferenceDiffusion.coefficient_overrides` has carried **two** entries whose rates both had to be right (TASK-042 already gave it one non-momentum entry per declared field, so the single-field case was not the test; two fields at two rates, each measured against its own closed form, is what a shared-rate bug cannot pass). Two fields taking different values at the same wall through `BoundaryFaceConfig.field_values`: its second scenario, the first to reach that surface through real `load_config` → `assemble_numerics` rather than a hand-built condition. **The qualification: TASK-035 did add a new field-keyed dictionary** -- `buoyancy_couplings`, built in `bootstrap.py` and threaded to `BoussinesqBuoyancy` on the same "assemble_numerics stays field-name-agnostic" split `coefficient_overrides` established. Judged not a failure of this criterion, and the reasoning is recorded rather than assumed: the bullet forbids *replacing* the per-field coefficient mechanisms this stage exists to test, and a buoyancy reference/coefficient pair is a quantity neither `coefficient_overrides` nor `overrides` could carry. Diffusivity and wall values -- the two coefficients that already had mechanisms -- both went through those mechanisms unchanged. | | 4. Buoyancy is one coupling, not one per field | **Met after this audit; the one check that made it a *claim* rather than an implementation detail was overstated as it landed.** One `BoussinesqBuoyancy` serves both couplings -- TASK-036 added zero lines under `src/pyflow/` and reused the class with `c = +1/rho_0` where temperature uses `c = -beta`. Warm fluid rises and reversing a *configured* gravity reverses it (`temperature_field.feature`); a denser patch sinks (`density_field.feature`), which a sign error common to both couplings could not survive. The null case is exact: a uniform temperature field leaves velocity bit-identical to carrying no temperature field at all, `torch.equal`, not a tolerance. **What failed: `density_field.feature`'s "one instance, both couplings" scenario built the coupling map itself** -- a line-for-line copy of `bootstrap.py`'s loop -- and handed it to `assemble_numerics`, so it proved that `assemble_numerics` forwards a map the *test* wrote. Criterion 4's own words are "reached through the same configured seam in both", and a seam a test reimplements is not the seam. Rewritten in this audit to go through real `bootstrap()` and to probe the assembled instance with each declared field alone; verified to fail (`density's declared coupling drove nothing downward`) when `bootstrap.py` is made to drop the second declaration. **The sign was also re-derived independently against `docs/handbook/physics/buoyancy.md`** rather than taken from TASK-035's own derivation: that entry's per-volume form divided by `rho_0` gives `-beta (T - T_0) g` and `(rho - rho_0)/rho_0 g`, which is exactly `c * (phi - phi_0) * g` under the two coefficients the configuration supplies. | | 5. A passive tracer is exactly passive | **Met, both bullets, in the one scenario the criterion insisted they share.** The same configuration run with and without a declared tracer produces bit-identical velocity components (`torch.equal`), and in that same scenario the tracer itself is measurably different after several timesteps than after one -- so the exactness above cannot be passed by a tracer the engine ignores. Four tracers transported together are pairwise distinct and each identical to its own solo run. | -| 6. Physical correctness against a known answer, per case | **Met on all four cases; the criterion's own closing sentence had not been carried out, and is now.** Temperature: a sinusoidal mode's measured decay rate matches `Gamma * k^2` within 10%, through a real configured field, and `tests/features/heat_transport.feature` runs the identical check on the committed demo -- the agreement Stage 5's anonymous scalar and this stage's named field were required to show. Humidity and density: each field's own domain integral unchanged to 1e-9 under pure advection on a periodic domain with diffusion and source both zeroed, so only transport could leak. Buoyancy: a layer heated from below develops more than twice the vertical-velocity RMS of the same layer heated from above (Rayleigh number ~2.0e4 against a critical ~1708, so the unstable case is genuinely supercritical), which is the qualitative bar design question five settled. **What was missing: the criterion says `docs/planning/backlog.md`'s Rayleigh-Bénard item "is amended in the same change rather than left reading as though this stage owed it", and the commit that drafted it recorded "both documents that name it say so now". Only `docs/planning/implementation-plan.md` was amended.** Two documents name ~1708; the backlog was the other one, and its bullet still read as though Level 3 owed the number. Amended in this audit's own change. | +| 6. Physical correctness against a known answer, per case | **Met on all four cases; the criterion's own closing sentence had not been carried out, and is now.** Temperature: a sinusoidal mode's measured decay rate matches `Gamma * k^2` within 10%, through a real configured field, and `tests/features/heat_transport.feature` runs the identical check on the committed demo -- the agreement Stage 5's anonymous scalar and this stage's named field were required to show. Humidity and density: each field's own domain integral unchanged to 1e-9 under pure advection on a periodic domain with diffusion and source both zeroed, so only transport could leak. Buoyancy: a layer heated from below develops more than twice the vertical-velocity RMS of the same layer heated from above (Rayleigh number ~2.0e4 against a critical ~1708, so the unstable case is genuinely supercritical), which is the qualitative bar design question five settled. **That buoyancy sentence stopped being true on 2026-09-12, and the correction belongs here rather than only in the stage that made it.** TASK-052 (Stage 9) stopped `FirstOrderUpwindAdvection` transporting through solid walls, and this fixture's own top and bottom walls are solid -- before that fix heat and momentum crossed them freely, with the Dirichlet condition re-injecting boundary-temperature fluid, a spurious heat pump. The factor-of-two verdict recorded above was produced by that defect: on correct physics the two final RMS values are 0.1621 and 0.1267, a ratio of 1.28, and extending the run to 1200 steps only reaches 1.59 because the unstable case saturates. **The criterion is still met, on a re-derived check** (maintainer's call, 2026-09-12): the scenario now compares *growth* rather than final magnitude -- the layer heated from below is still growing a convective roll (0.0904 -> 0.1621, +79%) while the one heated from above has settled (0.1251 -> 0.1267, +1.3%), a separation of roughly sixty to one against bars of 25% and 5%, mutation-verified against an inverted buoyancy sign. That is what this fixture's own docstring had said it always meant to check, and what the 2026-08-31 audit found had never been implemented. **Recorded as a verdict a later stage changed, not rewritten as though it had always read this way** (root `CLAUDE.md`'s Integrity section). **What was missing: the criterion says `docs/planning/backlog.md`'s Rayleigh-Bénard item "is amended in the same change rather than left reading as though this stage owed it", and the commit that drafted it recorded "both documents that name it say so now". Only `docs/planning/implementation-plan.md` was amended.** Two documents name ~1708; the backlog was the other one, and its bullet still read as though Level 3 owed the number. Amended in this audit's own change. | | 7. Rejection paths are exercised against real bad input | **Met as to the six surfaces named in advance -- and a seventh existed, was reachable from an ordinary configuration file, and did nothing.** The six: the `simulation.scalar_pattern` migration break, a duplicate field name, a name colliding with a velocity component, a name colliding with the pressure field, a non-positive diffusivity and an unknown initial-condition pattern (all `tests/features/field_declaration.feature`), plus a buoyancy coupling in a run whose velocity is not solved (`temperature_field.feature`). Each fails at `load_config` with a message naming the field, and each rejection step definition calls `load_config` directly and never reaches `bootstrap()` -- which is the accessor half made structural rather than asserted. **The seventh, found by this audit: a field declaring a buoyancy coupling while `numerics.source_term` is left at its default `"none"`.** It loads cleanly, transports the field, and no force ever reaches momentum -- measured end to end through `bootstrap()` at maximum vertical velocity **0.0**, against **0.451** for the identical configuration with `boussinesq_buoyancy` selected. `tests/unit/test_density_field.py`'s own committed two-coupling configuration was one such file, and its comment (`"none"-equivalent`) shows the assumption that hid it. **Fixed rather than recorded**, at the maintainer's standing direction: `_validate_buoyancy_couplings` now rejects it with a message naming the field, `numerics.source_term`, and what to set it to; the scenario was written first and watched fail. The rule is phrased as "no source term selected" rather than "not `boussinesq_buoyancy`", so a future term does not have to be added to it. | | 8. Executable Gherkin criteria, `make check-scenarios` gating -- and the step-definition count reported as evidence for or against the stage's own goal | **Met, and the count is reported as a finding, as the criterion requires whatever else passes. It is large.** Stage 6's five tasks added **93 step definitions**, taking the repository from 241 to 334 -- 28% of its entire step vocabulary, in one stage -- across **8 new feature files** and **+36 scenarios** (95 to 131). Both figures are measured (`grep -cE '^@(given|when|then)'` per module, at this stage's base commit and at its last task) rather than estimated. **Reuse, measured the same way, is small and concentrated:** the three golden-demo modules (`test_heat_transport.py`, `test_thermal_buoyancy.py`, `test_smoke_transport.py`) define 2 definitions each and reuse `tests/golden/conftest.py`'s demo vocabulary for the rest -- 9 reused step usages across 3 files; the five mechanism-level modules reuse nothing and define 87 between them. **The honest reading, since the criterion asks for one:** the number does not support this stage's claim. It has a structural explanation -- `tests/features/CLAUDE.md`'s own convention is that a step only one feature could use lives in that feature's module, and five of the eight files are engine-level with no shared vocabulary to draw on -- but an explanation is not a defence, and 93 is the figure. What it is *not* evidence of is engine special-casing: the same stage added zero `src/pyflow/` lines in three of its five tasks. `make check-scenarios` gates and reports all 132 scenarios across 29 feature files bound and running; `check_references.py`'s `PLANNED` table is empty, every promised artifact having landed. Three feature files exist beyond the five the task entries promised (`heat_transport`, `thermal_buoyancy`, `smoke_transport`), each recorded in its own task's Status note when it was added. | | 9. Demonstrations: Heat Transport, Smoke Transport and Thermal Buoyancy run from committed configuration | **Met, and the three-list reconciliation the criterion names had genuinely been done when the criteria were drafted -- verified rather than taken on trust.** All three run: `examples/golden-demos/heat_transport.yaml`, `smoke_transport.yaml`, `thermal_buoyancy.yaml`, each with a `tests/golden/` module invoking the real CLI as a subprocess with the demo's own config file (the Definition of Done's own strongest clause) plus a second scenario checking something physical rather than exit status -- the named field's analytic decay rate, the smoke field genuinely carried by the recirculating flow, the warm patch's vertical velocity positive at its own warmest cell. Each has a section in `docs/implementation/golden-demos.md`, a row in `docs/planning/implementation-plan.md`'s Golden Demos table, an entity in `planning/data/demos.yaml` with a `validates` edge to Capability Level 3, and a name in Level 3's own Golden Demo list -- all four checked directly, in all three cases. | @@ -12766,12 +12770,6 @@ needs a Neumann velocity face plus pressure anchoring, neither of which this stage has explored -- and a demo whose own treatment is unexplored is a design session, not a demonstration. -Tasks include - -- TASK-052 — Prescribed Boundary Velocity Reaches The Schemes -- TASK-053 — A Failed Frame Fails The Run -- TASK-054 — Timestep Stability Warning - ### Completion Criteria **Written 2026-09-12, before this stage's first task**, per @@ -12860,6 +12858,284 @@ claim. | 6. Documentation matches the tree | Whichever task lands last | | 7. Re-baselined numbers recorded with their predecessors | TASK-052 | + +## TASK-052 — Prescribed Boundary Velocity Reaches The Schemes + +**Status: Done, 2026-09-12.** Discharges Completion Criteria 1, 2, 3 and 7. + +### Purpose + +Stop `FirstOrderUpwindAdvection` transporting material through solid +walls. It took a boundary face's transporting velocity from the cell +inside it rather than from what the configuration prescribed there, so a +sealed box lost **14.27%** of a purely advected tracer in 400 steps on +`examples/golden-demos/smoke_transport.yaml`'s own geometry, with the +wall-normal velocity in boundary cells reaching 42% of the flow's peak +speed. `GreenGaussDivergence` meanwhile treated the same walls as solid, +so the two operators disagreed about the domain they were solving in. + +This falsifies a use case **Stage 4 wrote down for itself** -- "prescribe +a wall's value or its gradient, or make a boundary periodic, and have the +interior schemes honour it". Found 2026-09-12 by an end-to-end audit, +behind a green `make ci`. + +### Dependencies + +`advection.py`/`divergence.py`/`boundary_condition.py` (Stage 3/4) only. +No configuration change, no registry change, no interface change. + +### Design decisions, recorded here + +1. **The resolution is per-face, through `BoundaryCondition.evaluate`, + not a per-edge number threaded from configuration.** A first + implementation built a third mapping in `assemble_numerics` + (`{face_name: BoundaryFaceConfig.velocity}`, alongside + `boundary_conditions`/`periodic_pairs`) and widened the advection and + pressure-coupling registries to carry it. **It was built, measured + working, and then abandoned** -- `tests/unit/numerics/ + test_divergence_contract.py`'s own linear-field exactness test cannot + be expressed through a per-edge scalar, because a linear velocity + field's normal component varies *along* an edge, and a future + non-uniform inlet (a parabolic channel profile) would have the same + problem. `evaluate` is already per-face and already open to a + user-supplied implementation, and it is what `GreenGaussDivergence` + was already calling -- so the shipped fix adds no new mechanism at + all. Recorded rather than presented as the first idea: the abandoned + version touched six files and two registries, the shipped one touches + three files and neither. +2. **One function, two callers, checked structurally.** + `boundary_condition.boundary_normal_velocity` is extracted from what + `GreenGaussDivergence._boundary_face_normal_velocity` already did, and + `FirstOrderUpwindAdvection.flux` now calls it too. The check that both + still do is a source-level assertion (`boundary_velocity.feature`'s + first scenario), not a behavioural one -- the two operators agreed for + every fixture in the repository at the moment they were written and + disagreed for every fixture with motion near a wall, so a behavioural + check has to guess which fixture exposes the difference. +3. **A vector field with no override of its own resolves to `0.0`, not to + `BoundaryFaceConfig.scalar_value`.** `DirichletBoundaryCondition. + evaluate` used to fall through to its plain prescribed value for any + field name it had no override for -- and for the velocity field that + value is a *transported scalar's* boundary value, which is not a + velocity. Every configuration this repository ships leaves + `scalar_value` at `0.0`, so walls came out impermeable by coincidence + of the defaults rather than by anything the configuration said: the + same shape of accident the pre-fix divergence operator already had. + **Found by a test, not by reading** -- `boundary_velocity.feature`'s + "The same boundary without a prescribed normal velocity transports + nothing through itself" failed against a fixture whose scalar boundary + value was `3.0`, because the wall then resolved to a normal velocity + of 3.0. +4. **`BoundaryFaceConfig.velocity` is still read by nothing, and that is + recorded rather than fixed here.** It is validated for mutual + exclusivity with `pressure` and for zero net flux, and reaches no + scheme; the per-component channel (`field_values["velocity.0"]`/ + `["velocity.1"]`, which is what a real configuration already uses for + the lid) is the one this fix reads through. Closing that properly is + Criterion 3's remaining half -- see Discharges below, which records + that criterion as partially closed rather than met. + +### Artifacts Produced + +- `src/pyflow/engine/numerics/boundary_condition.py` -- + `boundary_normal_velocity`, the shared resolver; and + `DirichletBoundaryCondition.evaluate`'s narrowed fallback (design + decision 3). +- `src/pyflow/engine/numerics/advection.py` -- `flux` resolves a genuine + boundary face's transporting velocity through it, instead of leaving + the owner cell's own in place. +- `src/pyflow/engine/numerics/divergence.py` -- the same call it always + made, now through the shared function; `_boundary_face_normal_velocity` + no longer takes the field only to hand it to `evaluate`. +- `examples/golden-demos/sealed_box.yaml` -- this stage's Golden Demo. +- `src/pyflow/configuration/golden_demos.py` -- `sealed_box` appended to + the curated registry, so `pyflow run --demos sealed_box` works. +- `tests/features/boundary_velocity.feature` (9 scenarios), + `tests/features/sealed_box.feature` (3). +- `tests/unit/_numerics.py` -- `prescribed_face_normal_velocity`, a + companion to `face_normal_velocity_toward` rather than a widening of + it, per `tests/unit/CLAUDE.md`. +- Tests: 9 in `tests/unit/test_boundary_velocity.py`, 3 in + `tests/golden/test_sealed_box.py`; fixtures updated in + `test_first_order_upwind_advection.py`, `test_dirichlet_boundary.py`, + `test_velocity_field_support.py` (each now states a boundary's normal + velocity separately from the scalar's own boundary value) and + `test_temperature_field.py` (see Discharges). + +### Acceptance Criteria + +`tests/features/boundary_velocity.feature` and +`tests/features/sealed_box.feature` are the criteria +(`adr/ADR-007-executable-acceptance-criteria.md`). Both are bound and +gated by `make check-scenarios`. Additionally, verified by hand against +the real CLI (root `CLAUDE.md`'s Feature Verification rule): `uv run +python -m pyflow run --config examples/golden-demos/sealed_box.yaml +--backend offscreen --max-frames 20` exits cleanly, and `uv run python -m +pyflow run --demos` lists `sealed_box` as demo 12. + +### Discharges + +- **Criterion 1 in full.** `boundary_velocity.feature`'s structural + scenario plus its sealed-domain conservation scenario, which measures + `+0.000000000000%` against a pre-fix `-4.54%` on the same fixture. + Mutation-verified: reverting the resolver to the owner cell's own + velocity fails it. +- **Criterion 2 in full.** The prescribed-inflow/outflow/no-velocity + trio, plus the divergence operator reading the same number. +- **Criterion 3 partially: the rejection half is not built.** The sweep + over `BoundaryFaceConfig`'s own fields is not here, and + `BoundaryFaceConfig.velocity` still reaches no scheme -- design + decision 4. What this task closes is the half that was actually wrong: + a scalar's boundary value is no longer read as a velocity. **Named here + rather than left for the exit audit to find**, per + `docs/practices.md`'s "A discharge line that names no artifact is a + rubber stamp". +- **Criterion 7 in full.** Both Ghia baselines are recorded, in + `tests/features/navier_stokes_timestep.feature` and in + `tests/unit/test_navier_stokes_timestep.py`: the error against Ghia's + profiles **fell at every resolution** (9x9 0.1433 -> 0.1292, 13x13 + 0.0874 -> 0.0766, 17x17 0.0578 -> 0.0524) and stayed monotonic, which + is independent evidence the change was physics. Couette, Taylor-Green + (both halves), kinetic-energy monotonicity and Stage 6's periodic + domain integrals and decay rates were all confirmed unmoved rather than + assumed so. + +**One check did move, and it belonged to another stage.** Stage 6's +Rayleigh-Bénard onset scenario compared the two orientations' final +vertical-velocity RMS and required a factor of two; it passed at 6.35 and +now measures 1.28. That verdict was produced by this defect -- heat and +momentum were crossing that fixture's own solid top and bottom walls, +with the Dirichlet condition re-injecting boundary-temperature fluid, a +spurious heat pump. Re-derived as a growth comparison (maintainer's +call), which is what that fixture's own docstring said it always meant to +check, and which separates the two by roughly sixty to one rather than +1.28. Mutation-verified against an inverted buoyancy sign. See Stage 6's +own status table. + + +## TASK-053 — A Failed Frame Fails The Run + +**Status: Not started, drafted 2026-09-12.** Will discharge Completion +Criterion 4. + +### Purpose + +Make `pyflow run` and `pyflow play` exit non-zero when the engine raises +inside a frame. Today they do not: `DivergenceDidNotConvergeError` raised +inside `RenderWindow._draw`'s own `on_frame` callback is swallowed by +`rendercanvas`'s own `with log_exception("Draw error")` block (in that +third-party package's `base.py`, not anything tracked here), so PyFlow never sees +it. Measured on the shipped cavity refined to 64x64 with its own +`numerics.timestep` untouched: 22 of 40 frames failed, and the CLI +printed `pyflow exited cleanly` and returned **0**. + +This falsifies the second of the two Stage 4 use cases this stage exists +for -- "solve a linear system, and be told when it did not converge +instead of receiving a plausible wrong answer". The headless `pyflow +record` path is already correct (exit 1, full divergence history), +because PyFlow owns that loop; the rendering path inverts control to +`rendercanvas` and loses the stack with it. + +**A second symptom, same cause**: on the `glfw` backend the raise also +skips `on_draw`'s own reschedule (`window.py`), so a `--max-frames` run +never reaches its frame budget and hangs. Observed at 300 s before being +killed. + +### Dependencies + +`rendering/window.py`, `bootstrap.py`. None on TASK-052. + +### Design decisions, recorded here + +1. **Catch in `RenderWindow._draw`, re-raise from `run`.** `rendercanvas` + exposes no error-handler API (checked directly: no `set_*error*`, + `error_handler`, `excepthook` or `on_error` anywhere in the package), + and `log_exception` de-duplicates by message hash, so a repeating + failure degrades to one-liners. The seam that works is PyFlow's own: + stash the exception, close the canvas, and re-raise once `run` + returns. That covers both backends and `play()` at once, since all + three go through the same `RenderWindow`. +2. **`frame_count` moves to after a successful `on_frame`.** It is + currently incremented before, so a frame that died in the simulation + still counts as drawn. +3. **`__main__.py` needs no change.** An exception out of `main()` + already gives exit 1 with the real traceback -- exactly how + `record`/`resume` propagate today. + +### Acceptance Criteria + +Prose bullets rather than a `.feature` file, the same scope judgement +TASK-044/045/049/050 already record: this is a mechanism claim, not a +physical one, and `adr/ADR-007-executable-acceptance-criteria.md`'s own +scope is "real simulation work". + +- A run that raises inside a frame exits non-zero and prints the engine's + own diagnostic, asserted as an exit code **and** a stderr substring, + per `tests/integration/`'s own convention -- an exit code alone does + not distinguish a real failure from argparse. +- Covered on both backends and from both `run` and `play`. +- A `--max-frames` `glfw` run that raises terminates rather than hanging. +- Mutation-verified: removing the re-raise fails the new test. + +--- + +## TASK-054 — Timestep Stability Warning + +**Status: Not started, drafted 2026-09-12.** Will discharge Completion +Criterion 5. + +### Purpose + +Tell a user their configured timestep is above the stability limit +*before* the run, rather than leaving them to infer it from the +explosion. `stable_timestep` (`engine/simulation.py`) already computes +that limit and is well-derived, but no live path calls it -- a gap Stage +5's own Criterion 12 verdict noted and filed rather than fixed. Measured +by refining the shipped cavity and leaving `numerics.timestep` alone: + +| mesh | configured dt | stable dt | ratio | outcome | +|------|---------------|-----------|-------|---------| +| 16x16 | 0.008 | 0.01562 | 0.51 | ran 150 steps | +| 32x32 | 0.008 | 0.00781 | 1.02 | ran 150 steps | +| 48x48 | 0.008 | 0.00521 | 1.54 | ran 150 steps | +| 64x64 | 0.008 | 0.00391 | 2.05 | diverged at step 17 | + +### Dependencies + +`simulation_run.py`, `engine/simulation.py`. None on TASK-052 or +TASK-053, though it pairs with TASK-053: together they turn a silent +explosion into a warning up front and a loud failure after. + +### Design decisions, recorded here + +1. **Warn, do not reject** (maintainer's call, 2026-09-12). + `_STABILITY_SAFETY_FACTOR` is `0.25` against a measured stable edge of + `0.3`, so a configured timestep above the derived limit is not + automatically unstable -- the 32x32 and 48x48 rows above demonstrably + run. Rejecting them would make a deliberately conservative heuristic + load-bearing. +2. **The check lives in `simulation_run.py`, not `bootstrap.py`**, so + `run`, `record` and `resume` all get it. `record` has the same defect + and is the path a long unattended run uses. +3. **`numerics.timestep: auto` is deliberately not built here.** Stage 10 + (Better Numerics)'s own Design Question Two already owns whether + adaptive timestepping means CFL-driven selection -- which is exactly + that -- or an embedded error-estimating integrator. See + `docs/planning/backlog.md` §15. + +### Acceptance Criteria + +Prose bullets, same scope judgement as TASK-053. + +- The warning names the configured timestep, the stable one, and their + ratio. A warning that says only "unstable" tells a user nothing they + can act on. +- Emitted on `run`, `record` and `resume` alike. +- **Absent** below the limit, checked as its own case -- a warning that + always fires is a warning nobody reads. +- Non-fatal: the run proceeds. +- Verified by hand against the real CLI at 64x64 and at 16x16. + --- # Stage 10 — Better Numerics diff --git a/docs/planning/status.md b/docs/planning/status.md index 9b3fcf1..1606f91 100644 --- a/docs/planning/status.md +++ b/docs/planning/status.md @@ -17,14 +17,14 @@ demand, not part of this file. ## Progress -**52/52 tasks complete (100%)** across 17 planned stages. For the full plan, including +**53/55 tasks complete (96%)** across 17 planned stages. For the full plan, including stages below not yet broken into tasks: [roadmap.md](roadmap.md). ```mermaid pie showData title "Tasks across the roadmap" - "Done" : 52 - "Not started" : 0 + "Done" : 53 + "Not started" : 2 ``` ### Milestones @@ -41,13 +41,13 @@ pie showData ### Up next -**Stage 9 -- Solver & Run Integrity** is next, and has not been broken into tasks yet. +**Stage 9 -- Solver & Run Integrity** is next, starting with TASK-053 (A Failed Frame Fails The Run), 1 more not yet started in this stage. ## Live repository facts - **49** `CLAUDE.md` files -- **1209** tests collected -- **144** Gherkin scenarios (`tests/features/*.feature`) +- **1226** tests collected +- **156** Gherkin scenarios (`tests/features/*.feature`) ## Stages @@ -168,7 +168,13 @@ pie showData ### Stage 9 -- Solver & Run Integrity -**no status recorded** -- not yet broken into tasks; 7 criteria defined, no status line yet +**no status recorded** -- `███░░░░░░░` 1/3 tasks; 7 criteria defined, no status line yet + +| Task | Status | Date | Artifact | +|------|--------|------|----------| +| TASK-052 -- Prescribed Boundary Velocity Reaches The Schemes | Done | 2026-09-12 | `examples/golden-demos/smoke_transport.yaml` | +| TASK-053 -- A Failed Frame Fails The Run | Not started | | | +| TASK-054 -- Timestep Stability Warning | Not started | | | ### Stage 10 -- Better Numerics diff --git a/docs/repository-inventory.md b/docs/repository-inventory.md index 611a253..a90b3c5 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. -**363 tracked files** across 49 directories; +**368 tracked files** across 49 directories; 2 are empty. ## (root) @@ -201,6 +201,7 @@ listing files. - `multi_field_plume.yaml` - `numerics_assembly.yaml` - `passive_scalar_transport.yaml` +- `sealed_box.yaml` - `smoke_transport.yaml` - `thermal_buoyancy.yaml` @@ -335,6 +336,7 @@ listing files. ## tests/features - `CLAUDE.md` +- `boundary_velocity.feature` - `central_difference_diffusion.feature` - `conjugate_gradient_solver.feature` - `density_field.feature` @@ -360,6 +362,7 @@ listing files. - `pressure_correction_loop.feature` - `pressure_field.feature` - `rk4_time_integration.feature` +- `sealed_box.feature` - `simulation_orchestrator.feature` - `smoke_transport.feature` - `temperature_field.feature` @@ -387,6 +390,7 @@ listing files. - `test_multi_field_plume.py` - `test_numerics_assembly.py` - `test_passive_scalar_transport.py` +- `test_sealed_box.py` - `test_smoke_transport.py` - `test_thermal_buoyancy.py` @@ -417,6 +421,7 @@ listing files. - `_numerics.py` - `test_benchmark_demos.py` - `test_bootstrap.py` +- `test_boundary_velocity.py` - `test_buoyancy.py` - `test_central_difference_diffusion.py` - `test_check_claims.py` diff --git a/docs/repository-manifest.md b/docs/repository-manifest.md index 4adae08..b0a9157 100644 --- a/docs/repository-manifest.md +++ b/docs/repository-manifest.md @@ -939,6 +939,19 @@ than after one, under the same recirculating lid-driven flow "passive" is proven separately, at the engine level, by `tests/unit/test_passive_tracers.py` (`tests/features/ passive_tracers.feature`). +`tests/golden/test_sealed_box.py` (TASK-052, Stage 9's own golden demo +module) binds `tests/features/sealed_box.feature` -- the required +CLI-subprocess scenario, plus two more: every wall face carrying exactly +zero advective flux *and* the cells against those walls still carrying +real motion (both halves, since the first alone is satisfied by a flow +that stopped), and the tracer genuinely being carried by the +recirculating flow. It is the one demo module here whose second scenario +checks something the engine does rather than something the frame shows, +reaching `window.assembled_numerics`/`window.simulation_fields` rather +than pixels -- a wall's flux is not a thing a frame can show. Exact +conservation is proven separately, at the engine level, by +`tests/unit/test_boundary_velocity.py` (`tests/features/ +boundary_velocity.feature`). `unit/` otherwise holds config/logging/rendering (D1/D2/D3), the tooling tests @@ -1100,8 +1113,13 @@ riding along a solved velocity with no coupling) and `golden-demos/multi_field_plume.yaml` (2026-09-04 -- Stage 6's own claim, four differently-named fields on one solved velocity, added because that stage's three demos each declare exactly one field against -a criterion asking for four in one run) are -the eleven demos so far: plain configuration files, no Python -- golden demos run +a criterion asking for four in one run) and +`golden-demos/sealed_box.yaml` (TASK-052, 2026-09-12 -- Stage 9's own, +and the second demo after Multi-Field Plume to exist because an audit +found a defect rather than because a task built a capability: a tracer +in a closed no-slip cavity, and every wall face carrying exactly zero +advective flux while the cells against those walls are still moving) are +the twelve 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 8c670b4..a8aa856 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. -Eleven demos live here as of 2026-09-04, one per stage that has produced +Twelve demos live here as of 2026-09-12, one per stage that has produced a visible capability, plus one that deliberately has nothing new to render and one added by an audit rather than by a task: @@ -108,6 +108,23 @@ render and one added by an audit rather than by a task: real CLI -- the second demo whose output is its point, after `numerics_assembly`. +- `sealed_box.yaml` (TASK-052, 2026-09-12), Stage 9's own -- and the + second demo after `multi_field_plume.yaml` to exist because an audit + found something wrong rather than because a task built a capability. + Identical mesh, timestep, viscosity and moving lid to + `smoke_transport.yaml`, plus one declared `tracer` field with a + deliberately small diffusivity. Its whole content is that a tracer in + a closed cavity stays in the cavity: before TASK-052, + `FirstOrderUpwindAdvection` took a boundary face's transporting + velocity from the cell inside it rather than from what the + configuration prescribed, and material crossed solid walls at up to + 42% of the flow's own peak speed. **Its own scenarios check the wall, + not the domain integral** -- diffusion to a zero-valued wall removes + tracer legitimately and `FieldConfig` rejects a non-positive + diffusivity, so exact conservation is proven at the engine level + instead (`tests/features/boundary_velocity.feature`). See the config + file's own header. + Every demo here should follow the same shape: - a plain YAML file setting only what makes that demo *that demo* -- diff --git a/examples/golden-demos/sealed_box.yaml b/examples/golden-demos/sealed_box.yaml new file mode 100644 index 0000000..d5182a9 --- /dev/null +++ b/examples/golden-demos/sealed_box.yaml @@ -0,0 +1,81 @@ +# Sealed Box golden demo (Stage 9, TASK-052) -- the visible form of the +# defect that opened Stage 9: a tracer in a closed, no-slip lid-driven +# cavity has nowhere to go, and must still be entirely inside the box +# however long it is stirred. +# +# Run it exactly the way any user would: +# +# uv run python -m pyflow run --config examples/golden-demos/sealed_box.yaml +# +# **Before TASK-052 it leaked.** `FirstOrderUpwindAdvection` took a +# boundary face's transporting velocity from the cell inside it rather +# than from what the configuration prescribed there, so material crossed +# solid walls at up to 42% of the flow's own peak speed -- measured on +# `smoke_transport.yaml`, whose geometry this demo shares, at 14.27% of a +# purely advected tracer lost over 400 steps. `GreenGaussDivergence` +# meanwhile treated the same walls as solid, so the two operators +# disagreed about the domain they were solving in. +# +# The same lid-driven-cavity shape `lid_driven_cavity.yaml` and +# `smoke_transport.yaml` already use (mesh, timestep, viscosity, moving +# north wall), with a deliberately small `diffusion_coefficient` so the +# tracer stays sharp enough to watch being carried around. +# +# **What this demo checks is the wall itself, not the domain integral.** +# Diffusion to a zero-valued wall removes tracer legitimately, so the +# integral is not exactly constant here and a scenario asserting that it +# was would be wrong. Exact conservation is the *engine*-level claim, and +# it is proven against a purely advective fixture in +# `tests/features/boundary_velocity.feature`. What this demo proves is +# the thing a viewer can see: every wall face carries exactly zero +# advective flux while the cells against those walls are still moving -- +# which is false on the pre-TASK-052 engine and true on this one. +# +# `fields.tracer.diffusion_coefficient` cannot be `0.0`: `FieldConfig` +# rejects a non-positive coefficient, and that rejection is a bound +# Stage 6 acceptance criterion (`tests/features/field_declaration.feature`). +# Relaxing it to allow a perfectly non-diffusive tracer is a real and +# defensible change, and deliberately not made here -- it would reopen +# another stage's criterion for one demo's convenience. + +mesh: + extent: [16, 16] + spacing: [0.0625, 0.0625] + +numerics: + timestep: 0.008 + boundary_conditions: + north: + type: dirichlet + field_values: + velocity.0: 1.0 + velocity.1: 0.0 + south: + type: dirichlet + east: + type: dirichlet + west: + type: dirichlet + +fields: + - name: tracer + initial_condition: gaussian_blob + # Small, not zero -- see the header for why zero is unavailable and + # why this demo's own claim does not need it. + diffusion_coefficient: 0.001 + +simulation: + velocity_solved: true + +fluid: + viscosity: 0.01 + +field_display: + low_color: "#0b1020" + high_color: "#ffd27f" + field_label: Tracer (model units) + vector_label: Velocity (model units) + +rendering: + title: Sealed Box + background_color: "#141422" diff --git a/planning/data/demos.yaml b/planning/data/demos.yaml index a91b610..7405261 100644 --- a/planning/data/demos.yaml +++ b/planning/data/demos.yaml @@ -305,6 +305,27 @@ entities: does not: every other level with a named demo had one here, and this one had a hole. + - id: demo-sealed-box + name: Sealed Box + description: >- + A tracer stirred by a solved lid-driven flow inside a closed + no-slip cavity, and still entirely inside it -- every wall face + carrying exactly zero advective flux while the cells against those + walls are still moving. + documented_in: docs/implementation/golden-demos.md + edges: + - type: validates + to: capability-level-2 + note: >- + The same Level Lid-Driven Cavity validates, and for a reason + that Level's own demo could not check: an incompressible + solver that lets material cross its own no-slip walls is not + solving the domain it was configured with. Recorded as a + `validates` edge rather than only a `demonstrates` one because + `demo-validates-something` is right to insist -- a demo + proving nothing named is mis-specified, and what this one + proves is that Level 2's boundary treatment holds. + # Level 10's Golden Demo is "To be defined as capabilities are # implemented" -- a real, deliberate absence rather than an oversight, so # there is no entity for it and none should be invented. diff --git a/planning/data/features.yaml b/planning/data/features.yaml index f87e333..85985c4 100644 --- a/planning/data/features.yaml +++ b/planning/data/features.yaml @@ -683,3 +683,33 @@ entities: to: stage-8 - type: depends_on to: task-047 + + - id: task-052 + name: "TASK-052 — Prescribed Boundary Velocity Reaches The Schemes" + documented_in: docs/planning/roadmap.md + must_appear_in: docs/planning/roadmap.md + edges: + - type: belongs_to + to: stage-9 + - type: depends_on + to: task-023 + + - id: task-053 + name: "TASK-053 — A Failed Frame Fails The Run" + documented_in: docs/planning/roadmap.md + must_appear_in: docs/planning/roadmap.md + edges: + - type: belongs_to + to: stage-9 + - type: depends_on + to: task-007 + + - id: task-054 + name: "TASK-054 — Timestep Stability Warning" + documented_in: docs/planning/roadmap.md + must_appear_in: docs/planning/roadmap.md + edges: + - type: belongs_to + to: stage-9 + - type: depends_on + to: task-034 diff --git a/src/pyflow/configuration/golden_demos.py b/src/pyflow/configuration/golden_demos.py index 80dcf3a..63a1df7 100644 --- a/src/pyflow/configuration/golden_demos.py +++ b/src/pyflow/configuration/golden_demos.py @@ -49,6 +49,7 @@ ("thermal_buoyancy", "thermal_buoyancy.yaml"), ("smoke_transport", "smoke_transport.yaml"), ("multi_field_plume", "multi_field_plume.yaml"), + ("sealed_box", "sealed_box.yaml"), ) _DEFAULT_GOLDEN_DEMOS_DIR = Path("examples/golden-demos") diff --git a/src/pyflow/engine/CLAUDE.md b/src/pyflow/engine/CLAUDE.md index c764284..430ebbd 100644 --- a/src/pyflow/engine/CLAUDE.md +++ b/src/pyflow/engine/CLAUDE.md @@ -457,13 +457,23 @@ concretely (`docs/handbook/numerical-methods/advection.md`). The face-normal velocity that decides which side is upstream is the average of owner and neighbour velocities for an interior face (exact on PyFlow's uniform MVP mesh, where both are equidistant from the shared -face) or the owner's own velocity alone for a boundary face, dotted with +face), dotted with `Mesh.face_normal`'s own canonical direction; `velocity_normal >= 0` means the owner is upstream (flow moving along the canonical direction, owner toward neighbour or outward at a boundary), matching how `accumulate_flux_to_cells`'s own sign convention reads that same direction (`simulation.py`'s `CLAUDE.md` entry, above). +**At a genuine boundary face it is the one that face's own condition +prescribes, resolved through `boundary_normal_velocity` (TASK-052, Stage +9, 2026-09-12) -- not the owner cell's own.** This sentence read "or the +owner's own velocity alone for a boundary face" until then, and that was +an accurate description of a defect: material crossed solid walls at up +to 42% of a flow's own peak speed, and a sealed box lost 14.27% of a +purely advected tracer over 400 steps. The owner's own velocity is still +what a *gradient* (Neumann) face extrapolates, which is how an outlet is +expressed and the one case where reading the interior is correct. + **At a boundary face with inflow, the exterior value comes from this scheme's own `boundary_conditions`, keyed by `StructuredCartesianMesh.boundary_face_name` (TASK-023's own Design @@ -662,6 +672,35 @@ either -- `bootstrap.py` decides, since that is the one place that legitimately knows a run's velocity field is conventionally named `"velocity"`. +**`boundary_condition.py` also holds `boundary_normal_velocity`, a free +function, since TASK-052 (Stage 9, 2026-09-12) -- the one place a +genuine boundary face's *transporting* velocity is resolved.** A +`"value"` face's own condition supplies it; a `"gradient"` face +extrapolates the owner cell's, which is how an outlet is expressed and +the one case where reading the interior is correct. It exists because +`FirstOrderUpwindAdvection` and `GreenGaussDivergence` disagreed about +it for sixteen days: divergence consulted the condition and got a solid +wall, advection used the owner cell's own velocity and transported +straight through the same wall -- 14.27% of a purely advected tracer +lost in 400 steps on a shipped demo. Neither operator decides it for +itself any more, and `tests/features/boundary_velocity.feature`'s first +scenario checks that structurally rather than behaviourally, since the +two agreed for every fixture in the repository at the moment they were +written and disagreed for every fixture with motion near a wall. + +**`DirichletBoundaryCondition.evaluate` no longer falls through to its +own prescribed value for a *vector* field with no override (TASK-052).** +That value is `BoundaryFaceConfig.scalar_value` -- a transported +scalar's boundary value, not a velocity -- and returning it as one made +a wall permeable at whatever a scalar happened to be pinned to there. +Every configuration this repository ships leaves it at `0.0`, so walls +came out impermeable by coincidence of the defaults rather than by +anything the configuration said. A vector field with no override now +resolves to `0.0`: a Dirichlet velocity boundary naming no normal +component is a no-penetration wall. **Found by a test, not by reading** +-- a new scenario failed against a fixture whose scalar boundary value +was `3.0`. + **`boundary_condition.py`** (TASK-019, done 2026-08-23) is `BoundaryCondition` -- two abstract members, not one: `evaluate(field, face) -> float` and a `kind: Literal["value", "gradient"]` property diff --git a/src/pyflow/engine/numerics/advection.py b/src/pyflow/engine/numerics/advection.py index 981568c..81dfe09 100644 --- a/src/pyflow/engine/numerics/advection.py +++ b/src/pyflow/engine/numerics/advection.py @@ -25,7 +25,10 @@ from pyflow.engine.collocated_field import CollocatedField from pyflow.engine.field import Field from pyflow.engine.mesh import StructuredCartesianMesh -from pyflow.engine.numerics.boundary_condition import BoundaryCondition +from pyflow.engine.numerics.boundary_condition import ( + BoundaryCondition, + boundary_normal_velocity, +) from pyflow.engine.vector_field import ( IncompatibleVelocityFieldError as IncompatibleVelocityFieldError, ) @@ -104,6 +107,18 @@ class FirstOrderUpwindAdvection(AdvectionScheme): keyed by named edge, and consults them itself rather than the orchestrator substituting a value into its output afterward. + **A genuine boundary face's own *transporting* velocity comes from + that face's condition too, not from the cell inside it (TASK-052, + Stage 9).** Resolved through `boundary_normal_velocity` + (`boundary_condition.py`), the same function `GreenGaussDivergence` + resolves it through -- which is the point: the two used to decide it + separately, and disagreed for any flow with motion near a wall. + Until then this scheme used the owner cell's own velocity there, and + transported straight through solid walls: 14.27% of a purely advected + tracer lost in 400 steps on a shipped demo. The *transported field's* + own boundary value is a separate question, answered exactly as + before (below). + **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 @@ -134,12 +149,39 @@ def flux(self, field: Field, velocity: VectorField) -> torch.Tensor: v_neighbour = velocity.values[geometry.neighbour_ids] # At a genuine boundary face, `neighbour_ids` is a placeholder # equal to `owner_ids` (built below), so this average reduces to - # `v_owner` exactly -- "a boundary face has only the owner's own - # velocity to draw on" falls out of the same formula, no separate - # branch needed. + # `v_owner` exactly. That used to be the answer at such a face; + # since TASK-052 (Stage 9) it is only the *input* to one -- the + # loop below replaces it with whatever that face's own condition + # prescribes, and the owner's own value survives only where a + # gradient face extrapolates it. v_avg = (v_owner + v_neighbour) / 2 velocity_normal = v_avg[:, 0] * geometry.normal_x + v_avg[:, 1] * geometry.normal_y + # A genuine boundary face's *transporting* velocity is whatever + # the configuration prescribes there, not whichever cell happens + # to be inside it (TASK-052, Stage 9). Resolved through + # `boundary_normal_velocity`, the same function + # `GreenGaussDivergence` resolves it through -- see that + # function's own docstring for why the two sharing it is the + # point rather than an economy. The vectorised value computed + # above is the owner's own normal velocity at these faces (the + # placeholder neighbour equals the owner), which is exactly what + # a gradient face wants extrapolated, so it is passed in rather + # than recomputed. + for face in geometry.boundary_faces: + condition = self._boundary_conditions.get(geometry.boundary_names[face]) + if condition is None: + # Unconfigured: left as the owner's own, and rejected + # below only if it turns out to be inflow -- the existing + # carve-out, deliberately unchanged (see this class's own + # docstring). Unreachable from a configuration file: + # `assemble_numerics` resolves a condition for every + # non-periodic face. + continue + velocity_normal[face] = boundary_normal_velocity( + condition, velocity, face, float(velocity_normal[face]) + ) + values = field.values owner_values = values[geometry.owner_ids] neighbour_values = values[geometry.neighbour_ids] diff --git a/src/pyflow/engine/numerics/boundary_condition.py b/src/pyflow/engine/numerics/boundary_condition.py index dbd331e..1431e38 100644 --- a/src/pyflow/engine/numerics/boundary_condition.py +++ b/src/pyflow/engine/numerics/boundary_condition.py @@ -26,6 +26,7 @@ from collections.abc import Mapping from typing import Literal +from pyflow.engine.collocated_field import CollocatedField from pyflow.engine.field import Field @@ -101,7 +102,31 @@ def kind(self) -> Literal["value", "gradient"]: def evaluate(self, field: Field, face: int) -> float: self._check_boundary_face(field, face) - return self._overrides.get(field.name, self._value) + override = self._overrides.get(field.name) + if override is not None: + return override + if isinstance(field, CollocatedField) and field.component_shape != (): + # A *vector* field with no override of its own prescribes + # nothing here, and `value` is `BoundaryFaceConfig. + # scalar_value` -- a transported scalar's boundary value, + # which is not a velocity and must not be returned as one + # (TASK-052, Stage 9). `0.0` is the honest answer: a + # Dirichlet velocity boundary that names no normal component + # is a no-penetration wall, which is what every + # configuration in this repository relies on. + # + # **Found by a test, not by reading.** The scenario + # `boundary_velocity.feature`'s "The same boundary without a + # prescribed normal velocity transports nothing through + # itself" failed against a fixture whose scalar boundary + # value was `3.0`, because the wall then resolved to a + # normal velocity of 3.0. Every shipped demo leaves + # `scalar_value` at `0.0`, so the walls were impermeable by + # coincidence of the defaults rather than by anything the + # configuration said -- the same shape of accident + # `GreenGaussDivergence`'s own pre-TASK-052 behaviour had. + return 0.0 + return self._value class NeumannBoundaryCondition(BoundaryCondition): @@ -126,3 +151,65 @@ def kind(self) -> Literal["value", "gradient"]: def evaluate(self, field: Field, face: int) -> float: self._check_boundary_face(field, face) return self._overrides.get(field.name, self._gradient) + + +def boundary_normal_velocity( + condition: BoundaryCondition, + velocity: Field, + face: int, + owner_normal_velocity: float, +) -> float: + """The face-normal velocity **transporting** material across a genuine + boundary face (TASK-052, Stage 9) -- the one source every operator + that needs that number resolves it through. + + `velocity` is the whole velocity field, and `face` the boundary face; + `owner_normal_velocity` is the owning cell's own velocity projected + onto the face normal, which is what a gradient face extrapolates. + + The rule, and it is the one `docs/handbook/numerical-methods/ + boundary-conditions.md` already states for every boundary quantity + ("the flux must instead be determined by the boundary condition + itself"): + + - A **value** (Dirichlet) face prescribes its own normal velocity, + and the condition supplies it. For every configuration this + repository ships that resolves to `0.0` -- a no-penetration wall. + - A **gradient** (Neumann) face prescribes none, by definition: the + normal velocity is whatever the interior brings to it, + extrapolated zero-order. This is how an outlet is expressed, and + it is the one case where reading the owner cell is correct. + + **Exists because two operators disagreed about it for sixteen days.** + This is exactly what `GreenGaussDivergence` already did, extracted + unchanged so that `FirstOrderUpwindAdvection` does it too. Advection + used the owner cell's own velocity instead and transported straight + through solid walls -- 14.27% of a purely advected tracer lost in + 400 steps on a shipped demo (`docs/planning/roadmap.md` Stage 9, + Completion Criterion 1). Neither operator decides this for itself any + more, which is what makes "they agree" structural rather than a + coincidence of whichever fixture is in front of them. + + **Deliberately per-face, not per-named-edge.** A first version of + this took a `Mapping[str, float | None]` built from + `BoundaryFaceConfig.velocity` -- one number per named edge -- and was + abandoned when `tests/unit/numerics/test_divergence_contract.py`'s + own linear-field exactness test could not be expressed through it: a + linear velocity field's normal component varies *along* an edge, so + a per-edge scalar cannot carry it, and a future non-uniform inlet (a + parabolic channel profile) would have the same problem. + `BoundaryCondition.evaluate` is already per-face and already open to + a user-supplied implementation, so it is the right shape; that it is + also what divergence was already calling means this change adds no + new mechanism at all. See `docs/planning/roadmap.md` TASK-052's own + Design decisions. + + Deliberately raises nothing: an unconfigured face is each calling + module's own exception vocabulary (`advection.py` and `divergence.py` + each own an `UnconfiguredBoundaryFaceError` of their own, for + reasons their docstrings record), so a caller checks for a missing + condition before reaching here. + """ + if condition.kind == "gradient": + return owner_normal_velocity + return condition.evaluate(velocity, face) diff --git a/src/pyflow/engine/numerics/divergence.py b/src/pyflow/engine/numerics/divergence.py index 6804e0b..45aeb38 100644 --- a/src/pyflow/engine/numerics/divergence.py +++ b/src/pyflow/engine/numerics/divergence.py @@ -21,7 +21,10 @@ from pyflow.engine.collocated_field import CollocatedField from pyflow.engine.field import Field from pyflow.engine.mesh import StructuredCartesianMesh -from pyflow.engine.numerics.boundary_condition import BoundaryCondition +from pyflow.engine.numerics.boundary_condition import ( + BoundaryCondition, + boundary_normal_velocity, +) from pyflow.engine.simulation import accumulate_flux_to_cells @@ -210,14 +213,30 @@ def _boundary_face_normal_velocity( normal_x: float, normal_y: float, ) -> float: + """The face-normal velocity at a genuine boundary face. + + **No longer takes the `field` (TASK-052, Stage 9).** It used to, + so it could call `condition.evaluate(field, face)` and read + whatever number that returned as a wall-normal velocity -- which + for a `VectorField` named `"velocity"` resolved to + `BoundaryFaceConfig.scalar_value`, a *transported scalar's* + boundary value, not a velocity at all. Every shipped + configuration left both at `0.0`, so the wall came out + impermeable by coincidence of the defaults rather than because + anything read the velocity the configuration prescribed. It now + reads that velocity, through the shared + `boundary_normal_velocity` resolver, which is what makes this + class's own docstring claim about + `BoundaryFaceConfig.velocity`'s convention true. + """ condition = self._boundary_conditions.get(boundary_name) if condition is None: raise UnconfiguredBoundaryFaceError( f"face {face} (boundary {boundary_name!r}) has no BoundaryCondition configured" ) - if condition.kind == "value": - return condition.evaluate(field, face) - return owner_x * normal_x + owner_y * normal_y + return boundary_normal_velocity( + condition, field, face, owner_x * normal_x + owner_y * normal_y + ) @dataclass(frozen=True) diff --git a/tests/features/boundary_velocity.feature b/tests/features/boundary_velocity.feature new file mode 100644 index 0000000..90b59a0 --- /dev/null +++ b/tests/features/boundary_velocity.feature @@ -0,0 +1,109 @@ +# The acceptance criteria for Prescribed Boundary Velocity Reaches The +# Schemes (TASK-052, Stage 9's first task). Not a golden demo -- no +# config file under `examples/golden-demos/`, no CLI run; +# `tests/unit/test_boundary_velocity.py` binds these scenarios directly, +# per `tests/unit/CLAUDE.md`'s own scope. The stage's visible +# demonstration is Sealed Box (`tests/features/sealed_box.feature`), +# which shows the same claim through the public CLI. +# +# **Every number below was measured against the pre-fix engine before +# these scenarios were written**, so each is a claim about a real +# difference rather than a bound chosen to pass. The headline: the +# sealed-box scenario loses 14.27% of the field on the shipped +# `smoke_transport.yaml` geometry before the fix, against the exact +# conservation pure advection in a closed domain must give. + +Feature: Prescribed Boundary Velocity + Stage 9 Completion Criteria 1 and 2. A boundary face's *transporting* + velocity comes from what the configuration prescribes there, not from + whichever cell happens to be inside it -- and every operator that needs + that number gets it from the same place, so two of them cannot disagree + about whether a wall is solid. + + # -- Criterion 1: one source, checked structurally. + # + # This is the scenario that would have caught the original defect, and + # it is deliberately a source-level check rather than a behavioural + # one. `FirstOrderUpwindAdvection` and `GreenGaussDivergence` agreed + # about walls for every fixture in the repository at the moment they + # were written, and disagreed for every fixture with motion near one -- + # so a behavioural check has to guess which fixture exposes the + # difference. This one cannot be passed by two implementations that + # happen to coincide today. + + Scenario: Both operators resolve a boundary's normal velocity through the same function + Given the advection scheme's source and the divergence scheme's source + Then both call boundary_normal_velocity, and neither decides a boundary's normal velocity for itself + + # -- Criterion 1: the physical claim, on a sealed domain. + # + # Pure advection in a closed domain must conserve exactly: with zero + # normal velocity on every wall, every boundary face's flux is zero, + # and interior faces cancel pairwise inside `accumulate_flux_to_cells`. + # "Exactly" means to floating-point tolerance, not to a loose bound -- + # a bound is what let the original defect hide for sixteen days. + + Scenario: A sealed domain conserves a purely advected field exactly + Given a closed no-slip cavity carrying a tracer, with the tracer's diffusion switched off + When the simulation is advanced for many timesteps + Then the tracer's domain integral is unchanged to floating-point tolerance + + # The negative control for the scenario above. Without it, "the + # integral did not change" is also satisfied by a tracer that never + # moved at all -- which is the shape of vacuous pass the Stage 4 exit + # audit found in this claim's own predecessor + # (`first_order_upwind_advection.feature`'s closed-domain scenario; see + # `docs/practices.md`'s "When a test is found weak, ask why its fixture + # had to be that shape"). + + Scenario: The same sealed domain genuinely moves the tracer it conserves + Given a closed no-slip cavity carrying a tracer, with the tracer's diffusion switched off + When the simulation is advanced for many timesteps + Then the tracer's own spatial distribution has measurably changed + + # -- Criterion 1: the wall itself, measured rather than inferred. + # The second Then is what stops the first being satisfied by a flow + # that simply stopped: the cells behind those walls still move. + + Scenario: No advective flux crosses a wall whose prescribed normal velocity is zero + Given a closed no-slip cavity carrying a tracer, with the tracer's diffusion switched off + When the simulation is advanced for many timesteps + Then every wall face's advective flux is exactly zero + And the cells behind those walls still carry real motion + + # -- Criterion 2: a prescribed velocity reaches the solver, and the + # scalar's own boundary value is not mistaken for one. + + Scenario: A prescribed inflow carries the boundary's own value into the domain + Given a boundary prescribing an inward normal velocity, and a different boundary value for the transported field + When the advective flux is computed + Then that boundary's advective flux is the prescribed velocity times the prescribed field value + + Scenario: The same boundary without a prescribed normal velocity transports nothing through itself + Given a boundary prescribing no normal velocity, and a boundary value for the transported field + When the advective flux is computed + Then that boundary's advective flux is exactly zero + + Scenario: A prescribed outflow carries the interior cell's own value out of the domain + Given a boundary prescribing an outward normal velocity, and a different boundary value for the transported field + When the advective flux is computed + Then that boundary's advective flux uses the interior cell's own value, not the boundary's + + # -- Criterion 2: the divergence operator reads the same number, so a + # prescribed inflow is a real divergence rather than a wall. + + Scenario: The divergence operator sees a prescribed inflow as a real divergence + Given a boundary prescribing an inward normal velocity, and a different boundary value for the transported field + When the velocity field's divergence is computed + Then the cell behind that boundary reports the prescribed inflow, not its own interior velocity + + # -- Criterion 2: a gradient face still extrapolates, and says so. + # A zero-gradient velocity boundary is how an outlet is expressed: the + # normal velocity is whatever the interior brings to it. This is the + # one case where reading the owner cell is correct, and it has to + # survive the fix rather than be swept up by it. + + Scenario: A boundary prescribing a velocity gradient extrapolates from the interior + Given a boundary prescribing a velocity gradient rather than a value + When the advective flux is computed + Then that boundary's normal velocity is the interior cell's own diff --git a/tests/features/navier_stokes_timestep.feature b/tests/features/navier_stokes_timestep.feature index 1c8a9a6..df9d435 100644 --- a/tests/features/navier_stokes_timestep.feature +++ b/tests/features/navier_stokes_timestep.feature @@ -83,14 +83,25 @@ Feature: Navier-Stokes Timestep # scenario made no absolute accuracy claim at all -- errors of 10, 5 # and 2 would have satisfied monotonic decrease exactly as well as the # real ones do. Measured on real runs at this exact origin, spacing and - # steadiness criterion: 0.1433 at 9x9, 0.0874 at 13x13, 0.0578 at - # 17x17, so the bound keeps roughly 38% margin at the finest while + # steadiness criterion: 0.1292 at 9x9, 0.0766 at 13x13, 0.0524 at + # 17x17, so the bound keeps roughly 34% margin at the finest while # sitting well below what the coarsest scores. A velocity field of # zeros -- the cheapest "solved nothing" failure -- scores 0.3366 - # against these same 34 tabulated points, nearly six times the bound. + # against these same 34 tabulated points, over six times the bound. # The convergence claim above it is still the gating one; this is what # stops the trend being a trend towards nothing in particular. # + # **Those three numbers were 0.1433, 0.0874 and 0.0578 until + # 2026-09-12, and they moved because the physics did** (TASK-052, + # Stage 9, which stopped advection transporting momentum through this + # cavity's own solid walls). The error fell at every resolution and + # stayed monotonic. Both sets are recorded rather than the old ones + # overwritten, per Stage 9 Completion Criterion 7 -- an error moving + # toward a published reference at every resolution is evidence the + # change was physics rather than a re-fitted tolerance, and that + # evidence needs its predecessor to compare against. The 0.08 bound is + # deliberately unchanged. + # # **This fixture takes two deliberate exceptions to Criterion 7's # degenerate-fixture rule, and they are forced by the reference, not # chosen.** That rule asks every fixture for a non-square mesh, a diff --git a/tests/features/sealed_box.feature b/tests/features/sealed_box.feature new file mode 100644 index 0000000..fa76332 --- /dev/null +++ b/tests/features/sealed_box.feature @@ -0,0 +1,42 @@ +# The acceptance criteria for the Sealed Box golden demo (TASK-052, +# Stage 9). `tests/golden/test_sealed_box.py` binds these scenarios. +# +# Stage 9's own Golden Demo, and the visible form of the defect that +# opened the stage: a tracer in a closed, no-slip cavity has nowhere to +# go. Before TASK-052 it left through the walls, because +# `FirstOrderUpwindAdvection` took a boundary face's transporting +# velocity from the cell inside it rather than from what the +# configuration prescribed there. +# +# **This demo checks the wall, not the domain integral** -- see the +# config file's own header for why. Diffusion to a zero-valued wall +# removes tracer legitimately, and `FieldConfig` will not accept a +# perfectly non-diffusive one, so exact conservation is proven at the +# engine level instead (`tests/features/boundary_velocity.feature`) +# against a purely advective fixture. What is checkable here, on a +# committed config through the public API, is that no tracer crosses a +# wall by advection while the flow against those walls is still moving. + +Feature: Sealed Box + A tracer stirred by a real, solved lid-driven flow inside a closed + no-slip cavity -- and still entirely inside it. + + Background: + Given the golden demo "sealed_box" + + Scenario: A user can run it with the documented command + When it is run through the public CLI, headless + Then the command exits cleanly + + # The second scenario `tests/golden/CLAUDE.md` requires of any demo + # whose output is its point: a clean exit alone would pass against a + # box that leaked, and against one whose flow never started. + + Scenario: No tracer crosses a wall, and the flow against those walls is still moving + When it is bootstrapped for several real timesteps + Then every wall face carries exactly zero advective flux + And the cells against those walls still carry real motion + + Scenario: The tracer is genuinely carried by the recirculating flow + When it is bootstrapped for several real timesteps + Then the tracer field has measurably changed from where it started diff --git a/tests/features/temperature_field.feature b/tests/features/temperature_field.feature index 7b823af..def3f7e 100644 --- a/tests/features/temperature_field.feature +++ b/tests/features/temperature_field.feature @@ -69,13 +69,30 @@ Feature: Temperature # -- Criterion: convection onset, the qualitative bar design question # five settled -- rolls form heated from below, not heated from above. # The critical Rayleigh number (~1708, rigid-rigid) is explicitly not - # this stage's bar; the quantitative comparison is deferred to Stage 8. - - Scenario: A fluid layer heated from below convects; the same layer heated from above does not + # this stage's bar; the quantitative comparison is deferred to Stage 10 + # (Better Numerics). **That pointer said "Stage 8" until 2026-09-12** -- + # correct when written, stale from the 2026-09-07 renumbering onward, + # and missed by that sweep because a `.feature` file is prose nothing + # greps for a stage number. + # + # **Re-derived 2026-09-12 (TASK-052, Stage 9): this compares growth, + # not final magnitude.** It compared magnitude until then, and passed + # at a ratio of 6.35 against a bar of 2 -- a verdict produced by the + # defect TASK-052 fixed, since heat and momentum were crossing this + # fixture's own solid top and bottom walls and the Dirichlet condition + # was re-injecting boundary-temperature fluid. On correct physics the + # final magnitudes are 0.1621 and 0.1267 and no magnitude bar + # separates them. Growth does, by roughly sixty to one, and is the + # physically meaningful statement: an unstable layer grows a + # convective roll, a stable one settles to a steady forced response to + # its own initial perturbation. See `tests/unit/test_temperature_field.py`'s + # own `_rayleigh_benard_rms` for the measurements. + + Scenario: A fluid layer heated from below keeps convecting; the same layer heated from above settles Given a closed, no-slip fluid layer heated from below And the same layer heated from above instead When both are advanced for many Navier-Stokes timesteps - Then the layer heated from below develops a substantially larger vertical velocity than the one heated from above + Then the layer heated from below keeps growing a vertical velocity while the one heated from above has settled # -- Criterion: rejection paths exercised against real bad input -- # the sixth named surface, belonging to this task since only it knows diff --git a/tests/golden/CLAUDE.md b/tests/golden/CLAUDE.md index 764176e..085a5fd 100644 --- a/tests/golden/CLAUDE.md +++ b/tests/golden/CLAUDE.md @@ -36,6 +36,18 @@ the real CLI", which reads the report back out of the subprocess's stderr. Ask of any new demo whether exit-code-zero really covers what it claims to show. +**`test_sealed_box.py` (TASK-052, Stage 9, added 2026-09-12) is the one +demo module here whose second scenario checks something the *engine* +does rather than something the frame shows** -- every wall face carries +exactly zero advective flux, while the cells against those walls are +still moving. Both halves are needed: the first alone is satisfied by a +flow that stopped. It reaches `window.assembled_numerics` and +`window.simulation_fields`, the same two report-back attributes +`test_numerics_assembly.py` and the Stage 6 demo modules already use, +rather than inspecting pixels -- a wall's flux is not a thing a frame +can show, and asserting on the rendered image instead would be a weaker +claim dressed as a stronger one. + **When adding a demo:** write the feature file first, bind it from a module here, and add a step only when the shared vocabulary genuinely cannot express the criterion -- not to make a scenario easier to write. diff --git a/tests/golden/test_sealed_box.py b/tests/golden/test_sealed_box.py new file mode 100644 index 0000000..5e77e78 --- /dev/null +++ b/tests/golden/test_sealed_box.py @@ -0,0 +1,114 @@ +"""Sealed Box golden demo (TASK-052, Stage 9). + +The acceptance criteria are `tests/features/sealed_box.feature` +(`adr/ADR-007-executable-acceptance-criteria.md`). Stage 9's own Golden +Demo, and the visible form of the defect that opened the stage: a tracer +in a closed, no-slip lid-driven cavity has nowhere to go, and before +TASK-052 it left through the walls anyway. + +**This module checks the wall, not the domain integral**, and the +config file's own header says why: diffusion to a zero-valued wall +removes tracer legitimately, and `FieldConfig` rejects a perfectly +non-diffusive one, so exact conservation is proven at the engine level +instead (`tests/features/boundary_velocity.feature`) against a purely +advective fixture. What this demo proves, on a committed config through +the public API, is that no tracer crosses a wall by advection while the +flow against those walls is still moving -- which is false on the +pre-TASK-052 engine and true on this one. +""" + +from __future__ import annotations + +import math + +from pytest_bdd import scenarios, then, when + +from pyflow.bootstrap import bootstrap +from pyflow.engine.mesh import StructuredCartesianMesh +from pyflow.engine.scalar_field import ScalarField +from pyflow.engine.vector_field import VectorField +from pyflow.rendering.window import RenderWindow + +from ._demo import DemoRun + +scenarios("sealed_box.feature") + +_FRAMES = 12 +"""Measured directly against this demo's own config before choosing it, +the same discipline `test_smoke_transport.py`'s own `_FRAMES` uses: by +twelve real timesteps the lid-driven circulation has reached the walls +(the fastest wall-adjacent cell is moving at ~0.29, well clear of the +bound below) and the tracer field has visibly changed from its initial +blob. Fewer frames and the flow has not yet reached the walls, which +would make a zero wall flux prove nothing. +""" + +_MOVING_WALL_CELL_SPEED = 0.05 +"""A cell against a wall must be moving at least this fast for "no flux +crosses the wall" to mean the wall is impermeable rather than the flow +dead. Measured ~0.29 at `_FRAMES`, so this keeps a factor of five in +hand. +""" + + +def _velocity(window: RenderWindow) -> VectorField: + assert window.simulation_fields is not None + components = [ + window.simulation_fields[VectorField.component_name("velocity", index)] + for index in range(2) + ] + return VectorField.assemble(components, "velocity") # type: ignore[arg-type] + + +def _boundary_faces(mesh: StructuredCartesianMesh) -> list[int]: + return [face for face in range(mesh.num_faces) if mesh.face_neighbours(face)[1] is None] + + +@when("it is bootstrapped for several real timesteps", target_fixture="window") +def _when_bootstrapped(demo: DemoRun) -> RenderWindow: + return bootstrap(demo.config_path, backend="offscreen", max_frames=_FRAMES) + + +@then("every wall face carries exactly zero advective flux") +def _then_no_wall_flux(window: RenderWindow) -> None: + assert window.simulation_fields is not None + assert window.assembled_numerics is not None + tracer = window.simulation_fields["tracer"] + assert isinstance(tracer, ScalarField) + mesh = tracer.mesh + assert isinstance(mesh, StructuredCartesianMesh) + flux = window.assembled_numerics.advection.flux(tracer, _velocity(window)) + for face in _boundary_faces(mesh): + assert float(flux[face]) == 0.0, ( + f"wall face {face} ({mesh.boundary_face_name(face)}) carries advective flux " + f"{float(flux[face])}; a sealed box must carry none through any wall" + ) + + +@then("the cells against those walls still carry real motion") +def _then_wall_cells_move(window: RenderWindow) -> None: + assert window.simulation_fields is not None + velocity = _velocity(window) + mesh = velocity.mesh + assert isinstance(mesh, StructuredCartesianMesh) + owners = {mesh.face_neighbours(face)[0] for face in _boundary_faces(mesh)} + fastest = max(math.hypot(*velocity.value_at(cell)) for cell in owners) + assert fastest > _MOVING_WALL_CELL_SPEED, ( + f"the cells against the walls must still be moving for a zero wall flux to mean the wall " + f"is impermeable rather than the flow dead; fastest was {fastest}" + ) + + +@then("the tracer field has measurably changed from where it started") +def _then_tracer_transported(demo: DemoRun, window: RenderWindow) -> None: + initial = bootstrap(demo.config_path, backend="offscreen", max_frames=1) + assert initial.simulation_fields is not None + assert window.simulation_fields is not None + early = initial.simulation_fields["tracer"] + late = window.simulation_fields["tracer"] + assert isinstance(early, ScalarField) + assert isinstance(late, ScalarField) + assert not early.values.equal(late.values), ( + "expected the tracer to keep being carried by the recirculating flow, not freeze after " + "one step" + ) diff --git a/tests/unit/CLAUDE.md b/tests/unit/CLAUDE.md index dca8ef3..f1954b0 100644 --- a/tests/unit/CLAUDE.md +++ b/tests/unit/CLAUDE.md @@ -397,6 +397,38 @@ shared one**: `test_periodic_boundary.py` needs face has no mesh-reported neighbour, and that difference is visible at the call site precisely because it was not hidden behind a default. +**`test_boundary_velocity.py` (TASK-052, added 2026-09-12) is the +fourteenth, and Stage 9's first** -- binding +`tests/features/boundary_velocity.feature`'s nine scenarios: a boundary +face's transporting velocity comes from what the configuration +prescribes there, and both operators that need that number resolve it +through one shared function. Same shape as every module before it: its +own `_Context` dataclass, no golden-demo config file or CLI run. + +**Two things about it are new to this directory, and both are +deliberate.** Its sealed-cavity fixture goes through +`assemble_numerics` rather than constructing the six schemes by hand -- +the defect it guards was a *disagreement between two operators a real +run assembles together*, so a fixture wiring them consistently itself +would prove nothing about what a configuration produces. And its first +scenario is a source-level assertion (`inspect.getsource` over +`advection.py` and `divergence.py`) rather than a behavioural one: the +two operators agreed for every fixture in the repository at the moment +they were written and disagreed for every fixture with motion near a +wall, so a behavioural check has to guess which fixture exposes the +difference. The same reasoning `test_velocity_field_support.py`'s own +no-special-casing check already uses. + +**`_numerics.py` gained `prescribed_face_normal_velocity` in the same +task -- a companion to `face_normal_velocity_toward`, not a widening of +it**, per this file's own "when a fixture detail genuinely differs, copy +it rather than adding a parameter to the shared one". That function +encodes what every scheme did at a boundary face before TASK-052 +(extrapolate the owner's own velocity), and several scenarios still +legitimately want exactly that -- a gradient boundary, an unconfigured +one. The new one encodes the rule for a face whose condition prescribes +a value. + **`test_golden_demo_annotations.py` (added 2026-09-03, Stage 7 (Rendering Annotations) exit audit) is in none of the lineages above, and the distinction is worth keeping.** It binds no feature file, diff --git a/tests/unit/_numerics.py b/tests/unit/_numerics.py index f3b7bce..cde3c55 100644 --- a/tests/unit/_numerics.py +++ b/tests/unit/_numerics.py @@ -40,6 +40,7 @@ from __future__ import annotations +from collections.abc import Mapping from typing import Literal from pyflow.engine.field import Field @@ -184,3 +185,35 @@ def face_normal_velocity(mesh: StructuredCartesianMesh, velocity: VectorField, f """ _owner, neighbour = mesh.face_neighbours(face) return face_normal_velocity_toward(mesh, velocity, face, neighbour) + + +def prescribed_face_normal_velocity( + mesh: StructuredCartesianMesh, + boundary_conditions: Mapping[str, BoundaryCondition], + velocity: VectorField, + face: int, +) -> float: + """The face-normal velocity a scheme should use at the genuine + boundary face `face` (TASK-052, Stage 9) -- derived here + independently of `boundary_condition.boundary_normal_velocity`, for + the same reason `face_normal_velocity_toward` above derives its own: + a test's notion of the right answer must not be the implementation's. + + **A companion to `face_normal_velocity_toward`, not a replacement + for it**, per `tests/unit/CLAUDE.md`'s "when a fixture detail + genuinely differs, copy it rather than adding a parameter to the + shared one". That function encodes what *every* scheme did at a + boundary face before TASK-052 -- extrapolate the owner's own + velocity -- and several scenarios still legitimately want exactly + that (a gradient boundary, an unconfigured one). This one encodes + the rule for a face whose condition prescribes a value. + + A condition with no entry for the velocity field's own name falls + through to its plain prescribed value, which is + `DirichletBoundaryCondition`'s documented behaviour and, for every + configuration this repository ships, `0.0` -- a no-penetration wall. + """ + condition = boundary_conditions.get(mesh.boundary_face_name(face) or "") + if condition is None or condition.kind == "gradient": + return face_normal_velocity(mesh, velocity, face) + return float(condition.evaluate(velocity, face)) diff --git a/tests/unit/test_boundary_velocity.py b/tests/unit/test_boundary_velocity.py new file mode 100644 index 0000000..e024a54 --- /dev/null +++ b/tests/unit/test_boundary_velocity.py @@ -0,0 +1,464 @@ +"""Binds `tests/features/boundary_velocity.feature` (TASK-052, Stage 9's +first task) -- Stage 9 Completion Criteria 1 and 2: a boundary face's +transporting velocity comes from what the configuration prescribes +there, and every operator that needs that number resolves it through one +shared function. + +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`). The +stage's visible demonstration is Sealed Box +(`tests/golden/test_sealed_box.py`), which makes the same claim through +the public CLI. + +**The sealed-cavity fixture goes through `assemble_numerics` rather than +constructing the six schemes directly**, unlike most modules here. That +is deliberate: the defect this task fixes was a disagreement between two +operators that a real run assembles together, so a fixture that builds +them separately could wire them consistently by hand and prove nothing +about what a configuration actually produces. +""" + +from __future__ import annotations + +import inspect +import math +from dataclasses import dataclass + +import torch +from pytest_bdd import given, scenarios, then, when + +from pyflow.configuration.schema import ( + BoundaryConditionsConfig, + BoundaryFaceConfig, + NumericsConfig, +) +from pyflow.engine.field import Field +from pyflow.engine.mesh import StructuredCartesianMesh +from pyflow.engine.numerics import advection as advection_module +from pyflow.engine.numerics import divergence as divergence_module +from pyflow.engine.numerics.advection import FirstOrderUpwindAdvection +from pyflow.engine.numerics.assembly import AssembledNumerics, assemble_numerics +from pyflow.engine.numerics.boundary_condition import ( + BoundaryCondition, + DirichletBoundaryCondition, + NeumannBoundaryCondition, +) +from pyflow.engine.numerics.divergence import GreenGaussDivergence +from pyflow.engine.scalar_field import ScalarField +from pyflow.engine.simulation import navier_stokes_step, stable_timestep +from pyflow.engine.vector_field import VectorField + +from ._numerics import default_mesh, west_face + +scenarios("boundary_velocity.feature") + + +# -- The sealed cavity ------------------------------------------------------- + +_CAVITY_EXTENT = (10, 10) +_CAVITY_SPACING = 0.1 +_CAVITY_ORIGIN = (0.35, -0.2) +"""Non-trivial, per Stage 5 Completion Criterion 7's degenerate-fixture +rule -- nothing here compares against a reference that fixes the frame, +so there is no reason to sit at the origin. +""" + +_CAVITY_VISCOSITY = 0.01 +_CAVITY_STEPS = 150 +"""Measured, not guessed: at this size and step count the tracer's own +peak-to-peak change is 0.75 (so it genuinely moves) and the pre-fix +engine loses **4.54%** of it through the walls, against the exact +conservation below. A shorter run leaks less and proves less. +""" + +_CONSERVATION_TOLERANCE = 1e-12 +"""Relative. "Exactly" is the claim -- with zero normal velocity on every +wall each boundary face's flux is exactly `0.0 * phi`, and interior faces +cancel pairwise inside `accumulate_flux_to_cells`, so the only drift +available is floating-point summation order. Measured drift is `0.0` to +every digit printed; this bound is a floor, not a fitted tolerance. +""" + + +def _cavity_numerics() -> AssembledNumerics: + """A lid-driven cavity's own numerics, assembled the way a real run + assembles them -- the tracer's diffusion coefficient is `0.0` so that + advection is the only transport left, which is what makes exact + conservation the right claim rather than an approximate one. + """ + lid = BoundaryFaceConfig(type="dirichlet", field_values={"velocity.0": 1.0, "velocity.1": 0.0}) + wall = BoundaryFaceConfig(type="dirichlet") + return assemble_numerics( + NumericsConfig( + boundary_conditions=BoundaryConditionsConfig( + north=lid, south=wall, east=wall, west=wall + ), + pressure_correction_tolerance=1e-8, + pressure_correction_max_iterations=200, + ), + 0.0, + {"velocity.0": _CAVITY_VISCOSITY, "velocity.1": _CAVITY_VISCOSITY}, + ) + + +def _cavity_mesh() -> StructuredCartesianMesh: + return StructuredCartesianMesh( + origin=_CAVITY_ORIGIN, + spacing=(_CAVITY_SPACING, _CAVITY_SPACING), + extent=_CAVITY_EXTENT, + ) + + +def _tracer(mesh: StructuredCartesianMesh) -> ScalarField: + """A blob offset below the lid, so the cavity's own circulation + carries it against three of the four walls over the run rather than + sitting still in the middle. + """ + width = _CAVITY_EXTENT[0] * _CAVITY_SPACING + center_x = _CAVITY_ORIGIN[0] + 0.5 * width + center_y = _CAVITY_ORIGIN[1] + 0.35 * width + sigma = 0.18 * width + return ScalarField( + mesh, + "tracer", + initial_value=lambda x, y: math.exp( + -((x - center_x) ** 2 + (y - center_y) ** 2) / (2 * sigma**2) + ), + ) + + +@dataclass +class _Context: + mesh: StructuredCartesianMesh + numerics: AssembledNumerics | None = None + fields: dict[str, Field] | None = None + initial_tracer: torch.Tensor | None = None + flux: torch.Tensor | None = None + divergence: torch.Tensor | None = None + boundary_conditions: dict[str, BoundaryCondition] | None = None + velocity: VectorField | None = None + scalar: ScalarField | None = None + target_face: int | None = None + prescribed_velocity: float | None = None + prescribed_value: float | None = None + + +# -- Given ------------------------------------------------------------------- + + +@given("the advection scheme's source and the divergence scheme's source", target_fixture="ctx") +def _given_both_sources() -> _Context: + return _Context(mesh=default_mesh()) + + +@given( + "a closed no-slip cavity carrying a tracer, with the tracer's diffusion switched off", + target_fixture="ctx", +) +def _given_sealed_cavity() -> _Context: + mesh = _cavity_mesh() + velocity = VectorField(mesh, "velocity", num_components=2, initial_value=(0.0, 0.0)) + fields: dict[str, Field] = {c.name: c for c in velocity.decompose()} + fields["tracer"] = _tracer(mesh) + return _Context(mesh=mesh, numerics=_cavity_numerics(), fields=fields) + + +def _boundary_fixture(prescribed_velocity: float | None, gradient: bool = False) -> _Context: + """One west-boundary fixture for every Criterion 2 scenario. + + `prescribed_value` (`3.0`) is the *scalar's* boundary value and + `prescribed_velocity` the boundary's own normal component -- two + numbers, deliberately different from each other and from every + interior value, so no scenario below can pass by confusing one for + another (`docs/practices.md`, "Verify a conversion where its factors + are distinct"). West's canonical normal is `(-1, 0)` and the + convention is positive outward, so a negative prescribed velocity is + inflow. + """ + mesh = default_mesh(extent=(3, 2)) + velocity = VectorField(mesh, "velocity", num_components=2, initial_value=(0.5, 0.0)) + scalar = ScalarField(mesh, "tracer", initial_value=1.0) + condition: BoundaryCondition + if gradient: + condition = NeumannBoundaryCondition(0.0) + else: + overrides = {} if prescribed_velocity is None else {"velocity": prescribed_velocity} + condition = DirichletBoundaryCondition(3.0, overrides) + # Only west carries the scenario's own condition; the other three + # extrapolate, so the west face is the only thing that can move any + # assertion below. With a uniform horizontal velocity the north and + # south faces then contribute exactly zero, which is what makes the + # divergence scenario's hand-derivation a two-term sum. + conditions: dict[str, BoundaryCondition] = { + name: NeumannBoundaryCondition(0.0) for name in ("north", "south", "east") + } + conditions["west"] = condition + return _Context( + mesh=mesh, + velocity=velocity, + scalar=scalar, + boundary_conditions=conditions, + target_face=west_face(mesh), + prescribed_velocity=prescribed_velocity, + prescribed_value=3.0, + ) + + +@given( + "a boundary prescribing an inward normal velocity, and a different boundary value for the " + "transported field", + target_fixture="ctx", +) +def _given_inflow_boundary() -> _Context: + return _boundary_fixture(prescribed_velocity=-2.0) + + +@given( + "a boundary prescribing an outward normal velocity, and a different boundary value for the " + "transported field", + target_fixture="ctx", +) +def _given_outflow_boundary() -> _Context: + return _boundary_fixture(prescribed_velocity=2.0) + + +@given( + "a boundary prescribing no normal velocity, and a boundary value for the transported field", + target_fixture="ctx", +) +def _given_no_prescribed_velocity() -> _Context: + # No `"velocity"` override at all, so the condition falls through to + # its own prescribed value's default of `0.0` -- the no-penetration + # wall every configuration in this repository ships. + return _boundary_fixture(prescribed_velocity=None) + + +@given("a boundary prescribing a velocity gradient rather than a value", target_fixture="ctx") +def _given_gradient_boundary() -> _Context: + return _boundary_fixture(prescribed_velocity=None, gradient=True) + + +# -- When -------------------------------------------------------------------- + + +@when("the simulation is advanced for many timesteps") +def _when_advanced(ctx: _Context) -> None: + assert ctx.numerics is not None + assert ctx.fields is not None + tracer = ctx.fields["tracer"] + assert isinstance(tracer, ScalarField) + ctx.initial_tracer = tracer.values.clone() + dt = stable_timestep(ctx.mesh, _CAVITY_VISCOSITY, velocity_scale=1.0) + fields = ctx.fields + for _ in range(_CAVITY_STEPS): + fields = navier_stokes_step(fields, "velocity", ctx.numerics, dt).fields + ctx.fields = fields + + +@when("the advective flux is computed") +def _when_flux_computed(ctx: _Context) -> None: + assert ctx.boundary_conditions is not None + assert ctx.velocity is not None + assert ctx.scalar is not None + scheme = FirstOrderUpwindAdvection(ctx.boundary_conditions, {}) + ctx.flux = scheme.flux(ctx.scalar, ctx.velocity) + + +@when("the velocity field's divergence is computed") +def _when_divergence_computed(ctx: _Context) -> None: + assert ctx.boundary_conditions is not None + assert ctx.velocity is not None + scheme = GreenGaussDivergence(ctx.boundary_conditions, {}) + ctx.divergence = scheme.divergence(ctx.velocity) + + +# -- Then -------------------------------------------------------------------- + + +@then( + "both call boundary_normal_velocity, and neither decides a boundary's normal velocity for " + "itself" +) +def _then_one_shared_resolver(ctx: _Context) -> None: + del ctx + advection_source = inspect.getsource(advection_module) + divergence_source = inspect.getsource(divergence_module) + for name, source in ( + ("advection", advection_source), + ("divergence", divergence_source), + ): + assert "boundary_normal_velocity(" in source, ( + f"{name}.py no longer calls the shared resolver; two operators deciding a wall's " + "normal velocity separately is the defect Stage 9 Criterion 1 exists to prevent" + ) + assert "def boundary_normal_velocity" not in source, ( + f"{name}.py defines its own boundary-normal-velocity resolver; there must be exactly " + "one, in boundary_condition.py" + ) + + +@then("the tracer's domain integral is unchanged to floating-point tolerance") +def _then_tracer_conserved(ctx: _Context) -> None: + assert ctx.fields is not None + assert ctx.initial_tracer is not None + tracer = ctx.fields["tracer"] + assert isinstance(tracer, ScalarField) + before = float(ctx.initial_tracer.sum()) + after = float(tracer.values.sum()) + assert abs(after - before) <= _CONSERVATION_TOLERANCE * abs(before), ( + f"a sealed domain must conserve a purely advected tracer exactly; its integral went " + f"{before} -> {after}, a relative change of {(after - before) / before}" + ) + + +@then("the tracer's own spatial distribution has measurably changed") +def _then_tracer_moved(ctx: _Context) -> None: + assert ctx.fields is not None + assert ctx.initial_tracer is not None + tracer = ctx.fields["tracer"] + assert isinstance(tracer, ScalarField) + moved = float((tracer.values - ctx.initial_tracer).abs().max()) + # Measured 0.75 on this fixture. A tracer that never moved would + # conserve its integral trivially, which is the vacuous pass the + # scenario above needs this one to rule out. + assert moved > 0.1, ( + f"the tracer must actually be transported for its conservation to mean anything; its " + f"largest per-cell change over the run was {moved}" + ) + + +@then("every wall face's advective flux is exactly zero") +def _then_no_wall_flux(ctx: _Context) -> None: + assert ctx.fields is not None + assert ctx.numerics is not None + tracer = ctx.fields["tracer"] + assert isinstance(tracer, ScalarField) + components = [ctx.fields[VectorField.component_name("velocity", i)] for i in range(2)] + velocity = VectorField.assemble(components, "velocity") # type: ignore[arg-type] + flux = ctx.numerics.advection.flux(tracer, velocity) + for face in range(ctx.mesh.num_faces): + if ctx.mesh.face_neighbours(face)[1] is None: + assert float(flux[face]) == 0.0, ( + f"wall face {face} ({ctx.mesh.boundary_face_name(face)}) carries advective flux " + f"{float(flux[face])}; a wall prescribing zero normal velocity must carry none" + ) + + +@then("the cells behind those walls still carry real motion") +def _then_wall_cells_move(ctx: _Context) -> None: + assert ctx.fields is not None + components = [ctx.fields[VectorField.component_name("velocity", i)] for i in range(2)] + velocity = VectorField.assemble(components, "velocity") # type: ignore[arg-type] + boundary_owners = { + ctx.mesh.face_neighbours(face)[0] + for face in range(ctx.mesh.num_faces) + if ctx.mesh.face_neighbours(face)[1] is None + } + fastest = max(math.hypot(*velocity.value_at(cell)) for cell in boundary_owners) + # Without this, "no flux crosses the wall" is also satisfied by a + # flow that stopped everywhere. Measured ~0.3 on this fixture. + assert fastest > 0.05, ( + f"the cells against the walls must still be moving for a zero wall flux to mean the wall " + f"is impermeable rather than the flow dead; fastest was {fastest}" + ) + + +def _boundary_flux(ctx: _Context) -> float: + assert ctx.flux is not None + assert ctx.target_face is not None + return float(ctx.flux[ctx.target_face]) + + +@then("that boundary's advective flux is the prescribed velocity times the prescribed field value") +def _then_inflow_carries_prescribed_value(ctx: _Context) -> None: + assert ctx.prescribed_velocity is not None + assert ctx.prescribed_value is not None + expected = ctx.prescribed_velocity * ctx.prescribed_value + assert _boundary_flux(ctx) == expected, ( + f"expected the prescribed inflow to carry the prescribed boundary value in: " + f"{ctx.prescribed_velocity} * {ctx.prescribed_value} = {expected}, " + f"got {_boundary_flux(ctx)}" + ) + + +@then("that boundary's advective flux is exactly zero") +def _then_no_transport(ctx: _Context) -> None: + assert _boundary_flux(ctx) == 0.0, ( + f"a boundary prescribing no normal velocity must transport nothing through itself; got " + f"{_boundary_flux(ctx)}" + ) + + +@then("that boundary's advective flux uses the interior cell's own value, not the boundary's") +def _then_outflow_carries_owner_value(ctx: _Context) -> None: + assert ctx.prescribed_velocity is not None + assert ctx.prescribed_value is not None + assert ctx.scalar is not None + assert ctx.target_face is not None + owner, _ = ctx.mesh.face_neighbours(ctx.target_face) + expected = ctx.prescribed_velocity * ctx.scalar.value_at(owner) + assert _boundary_flux(ctx) == expected, ( + f"expected outflow to carry the interior cell's own value out: got " + f"{_boundary_flux(ctx)}, expected {expected}" + ) + assert _boundary_flux(ctx) != ctx.prescribed_velocity * ctx.prescribed_value + + +@then("the cell behind that boundary reports the prescribed inflow, not its own interior velocity") +def _then_divergence_sees_inflow(ctx: _Context) -> None: + assert ctx.divergence is not None + assert ctx.target_face is not None + assert ctx.prescribed_velocity is not None + assert ctx.velocity is not None + owner, _ = ctx.mesh.face_neighbours(ctx.target_face) + + # Hand-derived on this fixture's own mesh (origin (0.5, -1.0), + # spacing (0.2, 0.3), extent (3, 2), uniform velocity (0.5, 0)). + # The west owner cell has four faces, and only two contribute: + # west boundary, normal (-1, 0), area dy: prescribed * dy + # interior east face, normal (+1, 0), area dy: 0.5 * dy + # Its north and south faces see a normal velocity of exactly zero, + # because the flow is purely horizontal and both extrapolate. + # Divergence is that sum over the cell volume. + area = ctx.mesh.face_area(ctx.target_face) + volume = ctx.mesh.cell_volume(owner) + interior_face_velocity = 0.5 + expected = (ctx.prescribed_velocity * area + interior_face_velocity * area) / volume + + # What the same cell would report if the west face used the owner's + # own velocity instead -- the pre-TASK-052 answer, and exactly zero + # here, since a uniform flow entering and leaving one cell is + # divergence-free. The two differ by a factor of infinity, which is + # the strongest form this comparison can take. + owner_x, owner_y = ctx.velocity.value_at(owner) + normal_x, normal_y = ctx.mesh.face_normal(ctx.target_face) + interior_answer = ( + (owner_x * normal_x + owner_y * normal_y) * area + interior_face_velocity * area + ) / volume + assert interior_answer == 0.0, "fixture error: the interior-velocity answer should be zero" + + assert math.isclose(float(ctx.divergence[owner]), expected, rel_tol=1e-12, abs_tol=1e-15), ( + f"expected the cell behind a prescribed inflow to report it: {expected}, got " + f"{float(ctx.divergence[owner])}" + ) + assert float(ctx.divergence[owner]) != interior_answer + + +@then("that boundary's normal velocity is the interior cell's own") +def _then_gradient_extrapolates(ctx: _Context) -> None: + assert ctx.scalar is not None + assert ctx.velocity is not None + assert ctx.target_face is not None + owner, _ = ctx.mesh.face_neighbours(ctx.target_face) + owner_x, owner_y = ctx.velocity.value_at(owner) + normal_x, normal_y = ctx.mesh.face_normal(ctx.target_face) + expected_normal = owner_x * normal_x + owner_y * normal_y + # A gradient face extrapolates zero-order for the transported field + # too, so the face value is the owner's own and the flux is the + # product of the two. + expected = expected_normal * ctx.scalar.value_at(owner) + assert math.isclose(_boundary_flux(ctx), expected, rel_tol=1e-12, abs_tol=1e-15), ( + f"a gradient boundary must extrapolate the interior's own normal velocity " + f"({expected_normal}); got a flux of {_boundary_flux(ctx)} against {expected}" + ) diff --git a/tests/unit/test_dirichlet_boundary.py b/tests/unit/test_dirichlet_boundary.py index 8a6158e..683c64a 100644 --- a/tests/unit/test_dirichlet_boundary.py +++ b/tests/unit/test_dirichlet_boundary.py @@ -29,7 +29,7 @@ from ._numerics import ( default_mesh, - face_normal_velocity, + prescribed_face_normal_velocity, west_face, zero_gradient_everywhere, ) @@ -81,8 +81,17 @@ def _given_inflow_boundary(ctx: _Context) -> None: "interior cell's own value" ) def _given_real_dirichlet_for_advection(ctx: _Context) -> None: + # The `"velocity"` override is the boundary's own normal velocity, + # separate from `42.0`, which is the scalar value that inflow carries + # in (TASK-052, Stage 9 -- west's canonical normal is `(-1, 0)` and + # the convention is positive outward, so `-1.0` is inflow). This is + # still the real `DirichletBoundaryCondition`, which is what this + # scenario's own claim is about; only the fixture now says which of + # the two numbers it prescribes is which. ctx.prescribed_value = 42.0 - ctx.boundary_conditions = {"west": DirichletBoundaryCondition(ctx.prescribed_value)} + ctx.boundary_conditions = { + "west": DirichletBoundaryCondition(ctx.prescribed_value, {"velocity": -1.0}) + } @given("a boundary face") @@ -125,7 +134,9 @@ def _then_inflow_uses_prescribed_value(ctx: _Context) -> None: assert ctx.flux is not None assert ctx.target_face is not None assert ctx.prescribed_value is not None - velocity_normal = face_normal_velocity(ctx.mesh, ctx.velocity, ctx.target_face) + velocity_normal = prescribed_face_normal_velocity( + ctx.mesh, ctx.boundary_conditions, ctx.velocity, ctx.target_face + ) implied = float(ctx.flux[ctx.target_face]) / velocity_normal assert implied == ctx.prescribed_value diff --git a/tests/unit/test_first_order_upwind_advection.py b/tests/unit/test_first_order_upwind_advection.py index f1ea0e7..b31928d 100644 --- a/tests/unit/test_first_order_upwind_advection.py +++ b/tests/unit/test_first_order_upwind_advection.py @@ -27,16 +27,19 @@ FirstOrderUpwindAdvection, UnconfiguredBoundaryFaceError, ) -from pyflow.engine.numerics.boundary_condition import BoundaryCondition +from pyflow.engine.numerics.boundary_condition import ( + BoundaryCondition, + DirichletBoundaryCondition, +) from pyflow.engine.scalar_field import ScalarField from pyflow.engine.simulation import accumulate_flux_to_cells from pyflow.engine.vector_field import VectorField from ._numerics import ( FixedGradientCondition, - FixedValueCondition, default_mesh, face_normal_velocity, + prescribed_face_normal_velocity, zero_gradient_everywhere, ) @@ -104,8 +107,16 @@ def _given_outflow_boundary(ctx: _Context) -> None: @given("that boundary's own condition prescribes a value the interior cell does not have") def _given_outflow_condition(ctx: _Context) -> None: + # Two separate numbers on one condition since TASK-052 (Stage 9): + # `99.0` is the *scalar's* prescribed boundary value, and the + # `"velocity"` override is the boundary's own outward normal + # velocity, which is what now decides inflow from outflow. Before + # that task a scheme read the owner cell's own velocity here, so one + # number served both and the two could not be told apart -- which is + # exactly how a wall came to be permeable (`docs/planning/roadmap.md` + # Stage 9, Completion Criterion 1). ctx.prescribed_value = 99.0 - condition = FixedValueCondition(ctx.prescribed_value) + condition = DirichletBoundaryCondition(ctx.prescribed_value, {"velocity": 2.0}) ctx.boundary_conditions = { "north": condition, "south": condition, @@ -127,8 +138,15 @@ def _given_inflow_boundary(ctx: _Context) -> None: @given("that boundary's own condition prescribes a fixed value") def _given_inflow_value_condition(ctx: _Context) -> None: + # `-1.0` is the boundary's own normal velocity (west's canonical + # normal is `(-1, 0)`, and the convention is positive outward, so a + # negative value is inflow); `7.5` is the scalar value that inflow + # should carry in. See `_given_outflow_condition` above for why the + # two are separate numbers. ctx.prescribed_value = 7.5 - ctx.boundary_conditions = {"west": FixedValueCondition(ctx.prescribed_value)} + ctx.boundary_conditions = { + "west": DirichletBoundaryCondition(ctx.prescribed_value, {"velocity": -1.0}) + } @given("that boundary's own condition prescribes a gradient instead of a value") @@ -307,7 +325,9 @@ def _then_outflow_uses_owner(ctx: _Context) -> None: assert ctx.flux is not None assert ctx.target_face is not None owner, _ = ctx.mesh.face_neighbours(ctx.target_face) - velocity_normal = face_normal_velocity(ctx.mesh, ctx.velocity, ctx.target_face) + velocity_normal = prescribed_face_normal_velocity( + ctx.mesh, ctx.boundary_conditions, ctx.velocity, ctx.target_face + ) implied = float(ctx.flux[ctx.target_face]) / velocity_normal assert isinstance(ctx.scalar, ScalarField) assert implied == ctx.scalar.value_at(owner) @@ -320,7 +340,9 @@ def _then_inflow_uses_prescribed_value(ctx: _Context) -> None: assert ctx.flux is not None assert ctx.target_face is not None assert ctx.prescribed_value is not None - velocity_normal = face_normal_velocity(ctx.mesh, ctx.velocity, ctx.target_face) + velocity_normal = prescribed_face_normal_velocity( + ctx.mesh, ctx.boundary_conditions, ctx.velocity, ctx.target_face + ) implied = float(ctx.flux[ctx.target_face]) / velocity_normal assert implied == ctx.prescribed_value diff --git a/tests/unit/test_navier_stokes_timestep.py b/tests/unit/test_navier_stokes_timestep.py index c680faa..82a5b9e 100644 --- a/tests/unit/test_navier_stokes_timestep.py +++ b/tests/unit/test_navier_stokes_timestep.py @@ -738,13 +738,28 @@ def _then_taylor_green_mismatches(tg_result: _TaylorGreenContext) -> None: ones do. **Defended against three measured numbers, not chosen for comfort.** On -the finest mesh here (17x17) the real run scores 0.0578, so this bound -keeps roughly 38% margin. The coarsest (9x9) scores 0.1433, so this is +the finest mesh here (17x17) the real run scores 0.0524, so this bound +keeps roughly 34% margin. The coarsest (9x9) scores 0.1292, so this is genuinely a claim about the finest mesh rather than one any resolution would pass. And a velocity field of zeros everywhere -- the cheapest possible "solved nothing" failure -- scores 0.3366 against these same -34 tabulated points, so the bound sits nearly six times tighter than -doing nothing at all. It is not tight enough to call first-order upwind +34 tabulated points, so the bound sits over six times tighter than +doing nothing at all. + +**Those three numbers were 0.1433, 0.0874 and 0.0578 until 2026-09-12, +and they moved because the physics did** (TASK-052, Stage 9). That task +stopped `FirstOrderUpwindAdvection` transporting momentum through the +cavity's own solid walls; the error against Ghia's profiles fell at +every resolution -- 9x9 0.1433 -> 0.1292, 13x13 0.0874 -> 0.0766, 17x17 +0.0578 -> 0.0524 -- and stayed monotonic. **Both sets are recorded +rather than the old ones overwritten**, per Stage 9 Completion Criterion +7: a solver whose error moves toward a published reference at every +resolution is evidence the change was physics rather than a tolerance +edit, and that evidence only exists if the predecessor is still here to +compare against. The bound itself is deliberately unchanged at 0.08 -- +tightening it to the new margin was considered and not done, since this +criterion's gating claim is convergence and a bound re-fitted to each +improvement stops being an independent check. It is not tight enough to call first-order upwind accurate at this resolution, and is not meant to be: Criterion 5's gating claim is the convergence one, and `docs/implementation/ upgrade-paths.md` is where a less diffusive scheme lands. diff --git a/tests/unit/test_temperature_field.py b/tests/unit/test_temperature_field.py index ac5fc55..5fc7963 100644 --- a/tests/unit/test_temperature_field.py +++ b/tests/unit/test_temperature_field.py @@ -239,8 +239,8 @@ class _Context: without_field_state: dict[str, Field] | None = None marker_result: dict[str, Field] | None = None doubled_marker_result: dict[str, Field] | None = None - below_rms: float | None = None - above_rms: float | None = None + below_rms: tuple[float, float] | None = None + above_rms: tuple[float, float] | None = None # -- Given ----------------------------------------------------------------- @@ -411,6 +411,7 @@ def _given_boussinesq_selected_with_no_couplings(ctx: _Context) -> None: _RB_COEFFICIENT = -0.0102 _RB_DELTA_T = 20.0 _RB_STEPS = 300 +_RB_MIDPOINT = 150 def _rayleigh_benard_numerics(heated_from_below: bool) -> tuple[AssembledNumerics, float]: @@ -446,9 +447,9 @@ def _rayleigh_benard_numerics(heated_from_below: bool) -> tuple[AssembledNumeric return numerics, dt -def _rayleigh_benard_rms(heated_from_below: bool) -> float: - """The vertical velocity's own RMS after `_RB_STEPS` steps -- the one - number the onset scenario compares between the two configurations. +def _rayleigh_benard_rms(heated_from_below: bool) -> tuple[float, float]: + """The vertical velocity's own RMS at `_RB_MIDPOINT` and at + `_RB_STEPS` -- the pair the onset scenario compares. **This said "at four checkpoints -- a monotone rise for the unstable case, essentially flat for the stable one" until 2026-08-31, when @@ -456,10 +457,37 @@ def _rayleigh_benard_rms(heated_from_below: bool) -> float: checkpoint, not four, reached through a vestigial `[step for step in (_RB_STEPS,)]` comprehension left over from a multi-checkpoint draft, and only its last (and only) entry was ever - read. Nothing measured monotonicity or flatness. The description is - now what the function does; the qualitative bar design question five - settled is a comparison of the two final values, which is what the - scenario's own `Then` states. + read. Nothing measured monotonicity or flatness. + + **It returns two checkpoints since 2026-09-12 (TASK-052, Stage 9), + and the scenario compares growth rather than final magnitude -- + which is what that 2026-08-31 note said this was always meant to + do.** The reason it changed is not a preference: TASK-052 stopped + `FirstOrderUpwindAdvection` transporting through solid walls, and + this fixture's own top and bottom walls are solid. Before that fix + heat and momentum crossed them freely, with the Dirichlet condition + re-injecting boundary-temperature fluid -- a spurious heat pump that + drove the unstable case to an RMS of 0.4379 and suppressed the + stable one to 0.0690, a ratio of 6.35 against a bar of 2. On correct + physics the two final magnitudes are 0.1621 and 0.1267, a ratio of + 1.28, and extending the run to 1200 steps only reaches 1.59 because + the unstable case saturates. **The old verdict was produced by the + defect this fixture was meant to be independent of.** + + Growth separates them decisively where magnitude does not, and is + the physically meaningful statement in any case -- an unstable layer + *grows* a convective roll, a stable one settles to a steady forced + response to its own initial perturbation and stays there: + + | steps | heated from below | heated from above | + |-------|-------------------|-------------------| + | 150 | 0.0904 | 0.1251 | + | 300 | 0.1621 | 0.1267 | + + That is +79% against +1.3% over the same interval, a separation of + roughly sixty to one, on a bar of 25% and 5%. See + `docs/planning/roadmap.md`'s Stage 6 status table for the verdict + this changed, and Stage 9's Completion Criterion 7. """ numerics, dt = _rayleigh_benard_numerics(heated_from_below) mesh = StructuredCartesianMesh(origin=(0.0, 0.0), spacing=_RB_SPACING, extent=_RB_EXTENT) @@ -479,11 +507,17 @@ def temperature_ic(x: float, y: float) -> float: state = _initial_velocity_state(mesh) state["temperature"] = ScalarField(mesh, "temperature", initial_value=temperature_ic) - for _ in range(_RB_STEPS): + def vertical_rms() -> float: + v = state["velocity.1"] + assert isinstance(v, ScalarField) + return math.sqrt(sum(v.value_at(c) ** 2 for c in range(mesh.num_cells)) / mesh.num_cells) + + midpoint = 0.0 + for step in range(1, _RB_STEPS + 1): state = navier_stokes_step(state, "velocity", numerics, dt).fields - v = state["velocity.1"] - assert isinstance(v, ScalarField) - return math.sqrt(sum(v.value_at(c) ** 2 for c in range(mesh.num_cells)) / mesh.num_cells) + if step == _RB_MIDPOINT: + midpoint = vertical_rms() + return midpoint, vertical_rms() @given("a closed, no-slip fluid layer heated from below", target_fixture="ctx") @@ -687,18 +721,39 @@ def _then_runs_identical(ctx: _Context) -> None: ) +_RB_UNSTABLE_GROWTH = 1.25 +"""The unstable layer's vertical-velocity RMS must grow by at least this +factor over the run's second half. Measured 1.79 (0.0904 -> 0.1621), so +the bar keeps real margin below what a correct solver actually does. +""" + +_RB_STABLE_GROWTH = 1.05 +"""And the stable layer's must not. Measured 1.013 (0.1251 -> 0.1267) -- +a settled forced response, not a growing one. The two bars do not +overlap, which is the claim: 1.25 is above 1.013 and 1.05 is below 1.79, +so neither configuration can satisfy the other's half by drifting a +little. +""" + + @then( - "the layer heated from below develops a substantially larger vertical velocity than the " - "one heated from above" + "the layer heated from below keeps growing a vertical velocity while the one heated from " + "above has settled" ) def _then_below_convects_above_does_not(ctx: _Context) -> None: assert ctx.below_rms is not None assert ctx.above_rms is not None - below_final = ctx.below_rms - above_final = ctx.above_rms - assert below_final > 2 * above_final, ( - f"expected heated-from-below's own RMS vertical velocity ({below_final}) to be at " - f"least twice heated-from-above's ({above_final})" + below_mid, below_final = ctx.below_rms + above_mid, above_final = ctx.above_rms + assert below_final > _RB_UNSTABLE_GROWTH * below_mid, ( + f"expected the layer heated from below to still be growing a convective roll: its " + f"vertical-velocity RMS went {below_mid} -> {below_final}, a factor of " + f"{below_final / below_mid}, against a bar of {_RB_UNSTABLE_GROWTH}" + ) + assert above_final < _RB_STABLE_GROWTH * above_mid, ( + f"expected the layer heated from above to have settled: its vertical-velocity RMS went " + f"{above_mid} -> {above_final}, a factor of {above_final / above_mid}, against a bar of " + f"{_RB_STABLE_GROWTH}" ) diff --git a/tests/unit/test_velocity_field_support.py b/tests/unit/test_velocity_field_support.py index b268132..1530414 100644 --- a/tests/unit/test_velocity_field_support.py +++ b/tests/unit/test_velocity_field_support.py @@ -277,8 +277,20 @@ def _given_self_advected_velocity(ctx: _Context) -> _Context: # `step`-driving scenario in this file uses. u_name = VectorField.component_name("velocity", 0) v_name = VectorField.component_name("velocity", 1) + # The `"velocity"` entry is west's own prescribed *normal* velocity + # (negative is inward, since west's canonical normal is `(-1, 0)`), + # added by TASK-052 (Stage 9). It reads `-2.0` because that is the + # value this fixture always relied on -- before that task a boundary + # face's transporting velocity was the owner cell's own, and cell 0's + # `u` is `2.0`, so west was an inflow by accident of the interior + # rather than by anything the fixture said. Stating it keeps the + # hand-derivation below unchanged and makes the boundary's own + # intent explicit; dropping it makes west impermeable and `u` + # becomes `[1.4, 1.4]`, which is the correct answer to a different + # question. ctx.advection = FirstOrderUpwindAdvection( - {"west": DirichletBoundaryCondition(0.0, {u_name: 3.0, v_name: 0.0})}, {} + {"west": DirichletBoundaryCondition(0.0, {"velocity": -2.0, u_name: 3.0, v_name: 0.0})}, + {}, ) return ctx