diff --git a/docs/architecture/engine.md b/docs/architecture/engine.md index dd7c9c4..bc1dc8c 100644 --- a/docs/architecture/engine.md +++ b/docs/architecture/engine.md @@ -254,15 +254,15 @@ consistent with it. **Implementation:** `src/pyflow/engine/numerics/pressure_coupling.py` (`docs/planning/roadmap.md` TASK-021 Pressure Coupling Interface, Stage -3; TASK-027 PISO Pressure Coupling, Stage 4). The interface and its MVP -scheme, `PISO`, both live there -- the fifth of the six `adr/ADR-003` -components whose registered name resolves to a real implementation. -**A single, real dt-scaled correction pass, not the full multi-pass -Issa algorithm** -- `docs/architecture/icds.md`'s own Pressure-Velocity -Coupling entry records why (Rhie-Chow interpolation, needed to suppress -pressure-velocity decoupling under repeated correction on PyFlow's -collocated mesh, needs momentum-equation coefficients this task's own -interface does not have; that stronger claim is Stage 5 TASK-033's own). +3; TASK-027 PISO Pressure Coupling, Stage 4; TASK-033 Pressure +Correction Loop, Stage 5). The interface and its MVP scheme, `PISO`, +both live there -- the fifth of the six `adr/ADR-003` components whose +registered name resolves to a real implementation. **Genuinely +multi-pass since TASK-033 (Stage 5, 2026-08-29)**, not only the single, +real dt-scaled correction pass TASK-027 shipped -- `docs/architecture/ +icds.md`'s own Pressure-Velocity Coupling entry records the full +resolution (the momentum coefficient Rhie-Chow needs, `a_P = V/dt`, +paired with the same compact Laplacian the Poisson matrix already uses). TASK-021 also builds `src/pyflow/engine/numerics/assembly.py`, the registry all six of these layers resolve a configured name through. diff --git a/docs/architecture/icds.md b/docs/architecture/icds.md index a410e8a..0a7e8df 100644 --- a/docs/architecture/icds.md +++ b/docs/architecture/icds.md @@ -214,7 +214,12 @@ is transient or steady-state (`upgrade-paths.md` than PISO, only suited to different regimes; a future configuration should let a user pick deliberately, not assume one dominates). -**Configuration control:** `numerics.pressure_coupling` (implemented, Stage 3). +**Configuration control:** `numerics.pressure_coupling` (implemented, +Stage 3); `numerics.pressure_correction_tolerance`/ +`numerics.pressure_correction_max_iterations` (implemented, Stage 5 +TASK-033 -- the outer corrector loop's own tunables, distinct from +`numerics.linear_solver_tolerance`/`numerics.linear_solver_max_iterations`, +which govern each pass's inner solve). **Compatibility requirements:** requires a configured Linear Solver to solve the pressure-correction equation it produces each timestep -- the @@ -231,27 +236,54 @@ alternatives on its upgrade path exist to address -- not a defect to fix within PISO itself. **Done, TASK-027, 2026-08-27, with a real limitation recorded rather -than papered over**: `PISO` performs a single, real, dt-scaled pressure +than papered over**: `PISO` performed a single, real, dt-scaled pressure correction (`u_corrected = u* - dt * grad(p)`, `p` solving the compact Poisson equation `CentralDifferenceDiffusion`'s own already-symmetric Laplacian gives), verified to measurably and boundedly reduce a -manufactured provisional field's divergence. **It is not, and does not -claim to be, the full multi-pass Issa algorithm**: PyFlow's mesh is -collocated, and driving cell-centred divergence to near-zero under -*repeated* correction needs Rhie-Chow interpolation, which needs -momentum-equation coefficients this task's own interface has no way to +manufactured provisional field's divergence. **At that point it was not, +and did not claim to be, the full multi-pass Issa algorithm**: PyFlow's +mesh is collocated, and driving cell-centred divergence to near-zero +under *repeated* correction needs Rhie-Chow interpolation, which needs +momentum-equation coefficients this task's own interface had no way to obtain -- verified directly (composing this task's own `Gradient`/ `Divergence` into a Poisson matrix produces one that is provably not symmetric, so `ConjugateGradientSolver` cannot even solve it; three correction strategies were tried and measured before settling on the single-pass, compact-Laplacian design actually shipped). That stronger, -fully-converged claim belongs to Stage 5 TASK-033 (Pressure Correction +fully-converged claim was left to Stage 5 TASK-033 (Pressure Correction Loop), which has real momentum-coupled state to iterate against -- `docs/planning/roadmap.md` TASK-027's own Design decision Two records the full investigation, and `docs/practices.md`'s "A criterion whose strong reading depends on a later task must say so when drafted" is the standing rule this finding produced. +**Done, TASK-033, 2026-08-29: `PISO` is now genuinely multi-pass, and +the limitation above is resolved, not superseded by a rename.** The +question the paragraph above leaves open -- what carries the +momentum-equation coefficients Rhie-Chow needs, given PyFlow's momentum +predictor is fully explicit (RK4, no implicit assembly to draw a +coefficient from) -- resolved to `a_P = V/dt`: the unsteady term is the +only contribution to `a_P` for this architecture, and PyFlow's uniform +cell volume makes `V/a_P = dt` one constant for the whole mesh, needing +no new momentum-coefficient machinery. Pairing that correction with the +*same* compact Laplacian the Poisson matrix already uses -- not the +composed `Gradient`/`Divergence` pair TASK-027 tried and measured +failing -- restores the discrete adjoint property exactly; verified +numerically (both a linear and a nonlinear manufactured provisional +field converge to floating-point-exact zero divergence) before any +implementation code was written, the same prototype-first sequence +TASK-026/027 both used. `correct` now loops: each pass measures the +maximum cell divergence, records it, and either returns (at or below +`numerics.pressure_correction_tolerance`) or solves another correction, +up to `numerics.pressure_correction_max_iterations` passes before +raising `DivergenceDidNotConvergeError` rather than returning a +best-effort result. No change to `PressureCoupling.correct`'s own +signature, and no new ADR -- the outer-loop tolerance/iteration-limit +state is bound at `PISO`'s own construction, the same "strategy owns its +own tunables" shape `ConjugateGradientSolver` already established, not a +widening of the shared interface. `docs/planning/roadmap.md` TASK-033's +own Design decisions record the full numerical investigation. + --- ## Linear Solver diff --git a/docs/implementation/config-template.yaml b/docs/implementation/config-template.yaml index 74e2c2c..a554ccd 100644 --- a/docs/implementation/config-template.yaml +++ b/docs/implementation/config-template.yaml @@ -169,6 +169,14 @@ numerics: # Valid: "piso" -- the only scheme PyFlow currently implements for this # component. Invalid: any other string. pressure_coupling: piso + # Valid: a positive number -- the outer corrector loop's own convergence + # tolerance (distinct from linear_solver_tolerance above, which governs + # each inner linear solve, not how many corrector passes the outer loop + # may take). Invalid: zero or negative. + pressure_correction_tolerance: 1.0e-06 + # Valid: a positive integer -- the outer corrector loop's own iteration + # limit. Invalid: zero, negative, or a float. + pressure_correction_max_iterations: 50 # One entry per domain edge. All four faces share the same shape # (BoundaryFaceConfig) and the same per-field rules, explained once below # under 'north' and not repeated for the other three. diff --git a/docs/planning/roadmap.md b/docs/planning/roadmap.md index e299f53..c85e25b 100644 --- a/docs/planning/roadmap.md +++ b/docs/planning/roadmap.md @@ -189,7 +189,7 @@ This paragraph previously said `make install` and `make test` were still expected to fail, pending `uv.lock` and a test suite (B2/C1) -- stale since 2026-08-16 and corrected 2026-08-19. Both now succeed: `uv.lock` is committed (B2) and `make test` runs the suite with coverage -(C1a/C1b): **647 tests at 99% as of 2026-08-29**, having been 64 when +(C1a/C1b): **653 tests at 99% as of 2026-08-29**, having been 64 when this paragraph was rewritten on 2026-08-19, 202 earlier the same day, 212 after TASK-014, 226 after TASK-015, 250 after TASK-016, 287 after TASK-017, 297 after TASK-039, 315 after the Stage 2 exit audit and 337 @@ -420,7 +420,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. **76 of those 647 +class throughout, no new pressure-solving mechanism. **79 of those 653 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` @@ -492,7 +492,17 @@ leaving the corrected velocity unchanged (the null-space remedy made observable), `step` rejecting a `fields` mapping containing a `PressureField` by name, and a boundary configuration violating the zero-net-flux compatibility condition failing to load before any -pressure solve is attempted. +pressure solve is attempted; and to 79 with TASK-033's own +`pressure_correction_loop.feature`, three scenarios: a real solver's +corrector loop converging with a non-increasing recorded divergence +sequence, a deliberately halving solver forcing and proving multiple +genuine (strictly decreasing) passes, and exhausting the iteration limit +raising `DivergenceDidNotConvergeError` rather than returning a +best-effort result. **653 tests overall** -- the feature file's own +three scenarios plus three new `NumericsConfig.pressure_correction_ +tolerance`/`pressure_correction_max_iterations` load/reject tests in +`test_configuration.py`, the same config-section-addition shape every +prior task in this run used. **All** `make ci` targets pass, verified via the Makefile itself, not only via `uv tool run` in isolation -- that is `lint`, `typecheck`, `test`, `check-docs`, @@ -6749,7 +6759,7 @@ make part of TASK-033 a TASK-031 obligation, not a reordering. | 2. Pressure solved from the constraint, not transported | TASK-032 | | 3. Divergence decreases monotonically to the configured tolerance | TASK-033 | | 4. One timestep solves momentum and continuity together | TASK-034 | -| 5. Physical correctness against a known answer, per case | TASK-034, each for its own bullet; Couette flow jointly with TASK-033, the first task able to run it | +| 5. Physical correctness against a known answer, per case | TASK-034, each for its own bullet, Couette flow included -- TASK-033 supplies the corrector loop it needs but does not itself run the comparison | | 6. Rejection paths exercised against real bad input | TASK-041 and TASK-031..034, each for its own error conditions | | 7. Executable Gherkin criteria, `make check-scenarios` gates | TASK-041 and TASK-031..034, each for its own `.feature` file -- TASK-041 included, and its entry says why a task that computes nothing still owes one | | 8. Demonstrations: Lid-Driven Cavity and Heat Diffusion | TASK-034 (this stage's last task) | @@ -7316,6 +7326,63 @@ Criterion 2, entirely. Criterion 6 and Criterion 7, its own share. Pressure Correction Loop +**Status: Done, 2026-08-29, Stage 5's fourth task.** `PISO`'s single +correction pass (TASK-027) becomes a genuine corrector loop: each pass +measures the maximum cell divergence, records it, and either returns (at +or below `numerics.pressure_correction_tolerance`) or solves another +pass, up to `numerics.pressure_correction_max_iterations` before raising +rather than returning a best-effort result. + +- **Design question three, resolved by numerical prototyping before any + test or implementation code was written, the same sequence TASK-027 + used:** what carries the momentum-equation coefficients Rhie-Chow + needs, given PyFlow's momentum predictor is fully explicit (RK4, no + implicit assembly to draw a coefficient from)? Answer: `a_P = V/dt` -- + the unsteady term is the *only* contribution to `a_P` for this + architecture, and PyFlow's uniform cell volume makes `V/a_P = dt` one + constant for the whole mesh, not a per-cell field. No new + momentum-coefficient machinery, and no widened `PressureCoupling. + correct` -- `tolerance`/`max_iterations` are bound at `PISO`'s own + construction, the same "strategy owns its own tunables" shape + `ConjugateGradientSolver` already established. No ADR: the abstract + interface is byte-for-byte unchanged. +- **The second half of the answer, and the one TASK-027 actually got + wrong, not merely incomplete:** pairing the `a_P = V/dt` correction + with the *same* compact Laplacian the Poisson matrix is built from + (`_rhie_chow_divergence`, reusing `_poisson_matrix`'s own + `CentralDifferenceDiffusion`-based operator) is what restores the + discrete adjoint property -- not the composed `GreenGaussGradient`/ + `GreenGaussDivergence` pair TASK-027 tried and measured stalling at a + few percent reduction per pass. Verified numerically before being + trusted: a manufactured provisional velocity field (both a linear and + a nonlinear fixture, disposable prototype scripts, not committed) + converges to floating-point-exact zero divergence in a *single* + corrector pass once the operators match. This is also why the + feature file's second scenario needs a deliberately handicapped + `_HalvingSolver` test double rather than a real linear solver: on + PyFlow's uniform MVP mesh, the real corrector loop converges in one or + two passes, too fast to demonstrate genuine multi-pass behaviour on + its own -- the double halves the exact correction every pass instead, + producing an exactly-verifiable geometric decay (confirmed to hold to + ~13 significant figures across many passes before being written into + the test) that forces and proves multiple genuine iterations. +- **`DivergenceDidNotConvergeError` is new, and distinct from + `PressureSolveDidNotConvergeError`** (TASK-027): the former is the + *outer* corrector loop exhausting its iteration limit while divergence + never reaches tolerance (the inner solves all report success); the + latter is one pass's own inner linear solve failing to converge. A + loop using a solver that always reports `converged=True` but never + actually reduces divergence exercises the first and never the second + -- the third scenario's own fixture. +- **Couette flow's exact linear profile (Stage 5 Completion Criterion 5) + is not a scenario in this task's own feature file, deliberately** -- + the feature file's own header comment records why: this task lays the + corrector-loop groundwork the scenario needs (a genuinely convergent + `PISO`), not the demonstration itself, which needs a fully assembled + timestep only TASK-034 has. The Discharges line below still reads + "jointly with TASK-034" for that reason, not as a claim that a partial + scenario exists here. + **Intent:** the loop's claim is that divergence **decreases monotonically with each corrector iteration** and reaches the configured tolerance -- measured across iterations, not asserted at the end. A loop @@ -7350,18 +7417,16 @@ about what does not work, not from a blank page. ### Open design questions -**Three above -- the only one of this stage's seven still open, and this -task owns it.** What carries the momentum-equation coefficients -Rhie-Chow needs: a widened `PressureCoupling.correct`, a momentum -operator handed in at construction, or outer-loop state. **Resolve it -with numerical prototyping before writing this task's feature file**, -the same sequence TASK-027 used (five prototype scripts, then a scoping -decision, then the tests) -- writing a criterion first and discovering -it is unreachable second is precisely the failure `docs/practices.md`'s -"A criterion whose strong reading depends on a later task must say so -when drafted" was added to prevent, and this is the later task. It was -left open deliberately when the other six were decided, because -TASK-027 already demonstrated what deciding this one without +**Three above -- the only one of this stage's seven left open when the +other six were decided, and this task owned it. Resolved 2026-08-29: +outer-loop state the strategy owns (`a_P = V/dt`, bound at `PISO`'s own +construction), not a widened `PressureCoupling.correct` or a momentum +operator handed in separately.** Found by numerical prototyping before +writing this task's feature file, the same sequence TASK-027 used +(several disposable prototype scripts, then a scoping decision, then the +tests) -- see this task's own Status paragraph above for the finding +itself. It was left open deliberately when the other six were decided, +because TASK-027 already demonstrated what deciding this one without measurements costs. **Five above is no longer open, and that is what makes three @@ -7376,39 +7441,66 @@ which is exactly what it was not while the arrangement was undecided. ### Artifacts Produced - `tests/features/pressure_correction_loop.feature` -- this task's - Acceptance Criteria. -- An ADR if the resolution to design question three widens a Stage 3 - interface, matching the precedent - `adr/ADR-009-pressure-coupling-dt.md` set for the last such widening. + Acceptance Criteria. `tests/unit/test_pressure_correction_loop.py` + binds it, per `tests/unit/`'s own scope -- every scenario is checked + against the real `PISO` directly, no CLI subprocess. +- `src/pyflow/configuration/schema.py` -- `NumericsConfig. + pressure_correction_tolerance`/`pressure_correction_max_iterations`, + the outer corrector loop's own tunables (Criterion 12's own share), + distinct from `linear_solver_tolerance`/`linear_solver_max_iterations`, + which govern each pass's inner solve. +- `src/pyflow/engine/numerics/pressure_coupling.py` -- + `DivergenceDidNotConvergeError`, and `PISO.last_divergence_history`, + the recorded per-pass sequence the feature file's own scenarios assert + against directly. +- **No ADR**: design question three resolved to outer-loop state bound at + `PISO`'s own construction, not a Stage 3 interface widening -- + `PressureCoupling.correct`'s abstract signature is unchanged. The + Artifacts bullet drafted for this task anticipated needing one; it + didn't, the same way TASK-032's own local design question needed none + either. ### Acceptance Criteria `tests/features/pressure_correction_loop.feature` is the criteria. -Written to cover, at minimum: +Covers: - The recorded sequence of per-iteration maximum divergence magnitudes within one timestep is non-increasing at every element, and its last element is at or below the configured tolerance -- the sequence asserted, not just its last value. +- A solver that only partially corrects divergence each pass still takes + multiple genuine (strictly decreasing) passes to converge -- the + scenario a real solver's own fast convergence on PyFlow's uniform MVP + mesh cannot demonstrate on its own; see the Status paragraph above for + why a deterministic `_HalvingSolver` double is what exercises this. - Exhausting the corrector iteration limit without reaching tolerance - raises rather than returning a best effort, the same honesty - `PressureSolveDidNotConvergeError` already applies to a single solve - (Criterion 3's last bullet, and this task's share of Criterion 6). -- Couette flow's exact linear profile, jointly with TASK-034: this is - the first task with enough machinery to run it, and it is the cheapest - quantitative check that the coupled solve is right rather than merely - convergent. -- Whether `piso` is now genuinely multi-pass, or has been renamed -- one - or the other, per Criterion 3, and whichever it is must be visible in - the scenario rather than only in a docstring. + raises `DivergenceDidNotConvergeError` rather than returning a best + effort, the same honesty `PressureSolveDidNotConvergeError` already + applies to a single solve (Criterion 3's last bullet, and this task's + share of Criterion 6) -- and is a genuinely different error, since the + fixture's own inner solves all report success. +- `piso` is genuinely multi-pass, not renamed, per Criterion 3 -- visible + in the scenario (a real solver converging over a recorded sequence), + not only in the class's own docstring. +- **Couette flow's exact linear profile is deliberately not a scenario + here** -- left to TASK-034 alone, which has the fully assembled + timestep the comparison needs; see the Status paragraph above. ### Discharges -Criterion 3, entirely. Criterion 5's Couette bullet, jointly with -TASK-034. Criterion 12, its corrector-tolerance and iteration-limit -share. Criterion 6 and Criterion 7, its own share. **Stage 4 -Completion Criterion 4's Pressure-Velocity Coupling deferral**, which is -re-read at this task's close rather than assumed discharged. +Criterion 3, entirely. Criterion 12, its corrector-tolerance and +iteration-limit share. Criterion 6 and Criterion 7, its own share. +Criterion 5's Couette bullet stays with TASK-034 alone -- not discharged +here, and the "jointly" phrasing this section originally carried is +corrected: this task supplies the corrector loop the comparison depends +on, not part of the scenario itself. **Stage 4 Completion Criterion 4's +Pressure-Velocity Coupling deferral, re-read at this task's close rather +than assumed discharged: confirmed genuinely resolved.** `icds.md`'s own +Pressure-Velocity Coupling entry is updated in the same change (its own +"is not, and does not claim to be, the full multi-pass Issa algorithm" +sentence no longer describes the shipped `PISO`) -- see that document for +the reading that replaces it. --- @@ -7511,8 +7603,10 @@ feature file, are the criteria. Written to cover, at minimum: - Determinism: the same configuration run twice produces identical state. - Couette flow against its exact linear profile, at solver tolerance - rather than a loose one (jointly with TASK-033, and see Criterion 5's - own note on why the nonlinear term vanishing makes that reachable). + rather than a loose one -- this task's own scenario alone (TASK-033 + deliberately left it here, since it needs a fully assembled timestep; + see that task's own Status paragraph), and see Criterion 5's own note + on why the nonlinear term vanishing makes that reachable. - Lid-driven cavity against Ghia, Ghia & Shin (1982) at Re = 100: monotonically decreasing error across at least three mesh resolutions, plus the qualitative structure at the finest -- **not** a fixed @@ -7530,10 +7624,11 @@ feature file, are the criteria. Written to cover, at minimum: ### Discharges -Criteria 4, 8, 9, 10, 11 and 13, entirely. Criterion 5, all bullets (its -Couette bullet jointly with TASK-033). Criterion 12, its -tangential-boundary and run-length share. Criterion 6 and Criterion 7, -its own share. +Criteria 4, 8, 9, 10, 11 and 13, entirely. Criterion 5, all bullets, +including the Couette one entirely -- TASK-033 supplies the corrector +loop it depends on but does not itself scenario-test it (see that +task's own Discharges). Criterion 12, its tangential-boundary and +run-length share. Criterion 6 and Criterion 7, its own share. Golden Demo diff --git a/docs/planning/status.md b/docs/planning/status.md index 490321a..cd5fd2f 100644 --- a/docs/planning/status.md +++ b/docs/planning/status.md @@ -15,14 +15,14 @@ demand, not part of this file. ## Progress -**36/42 tasks complete (86%)** across 14 planned stages. For the full plan, including +**37/42 tasks complete (88%)** across 14 planned stages. For the full plan, including stages below not yet broken into tasks: [roadmap.md](roadmap.md). ```mermaid pie showData title "Tasks across the roadmap" - "Done" : 36 - "Not started" : 6 + "Done" : 37 + "Not started" : 5 ``` ### Milestones @@ -35,13 +35,13 @@ pie showData ### Up next -**Stage 5 -- First Fluid Solver** is next, starting with TASK-033 (Pressure Correction Loop), 1 more not yet started in this stage. +**Stage 5 -- First Fluid Solver** is next, starting with TASK-034 (Navier-Stokes Timestep). ## Live repository facts - **45** `CLAUDE.md` files -- **647** tests collected -- **76** Gherkin scenarios (`tests/features/*.feature`) +- **653** tests collected +- **79** Gherkin scenarios (`tests/features/*.feature`) ## Stages @@ -115,14 +115,14 @@ pie showData ### Stage 5 -- First Fluid Solver -**no status recorded** -- `██████░░░░` 3/5 tasks; 13 criteria defined, no status line yet +**no status recorded** -- `████████░░` 4/5 tasks; 13 criteria defined, no status line yet | Task | Status | Date | Artifact | |------|--------|------|----------| | TASK-041 | Done | 2026-08-28 | `src/pyflow/engine/numerics/assembly.py` | | TASK-031 | Done | 2026-08-29 | `advection.py` | | TASK-032 | Done | 2026-08-29 | `src/pyflow/engine/scalar_field.py` | -| TASK-033 -- Pressure Correction Loop | Not started | | | +| TASK-033 | Done | 2026-08-29 | `PressureCoupling.correct` | | TASK-034 -- Navier-Stokes Timestep | Not started | | | ### Stage 6 -- Additional Physical Fields diff --git a/docs/repository-inventory.md b/docs/repository-inventory.md index dd8becc..b28c3a0 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. -**276 tracked files** across 45 directories; +**278 tracked files** across 45 directories; 4 are empty. ## (root) @@ -314,6 +314,7 @@ listing files. - `passive_scalar_transport.feature` - `periodic_boundary.feature` - `piso_pressure_coupling.feature` +- `pressure_correction_loop.feature` - `pressure_field.feature` - `rk4_time_integration.feature` - `simulation_orchestrator.feature` @@ -382,6 +383,7 @@ listing files. - `test_neumann_boundary.py` - `test_periodic_boundary.py` - `test_piso_pressure_coupling.py` +- `test_pressure_correction_loop.py` - `test_pressure_field.py` - `test_rendering.py` - `test_rk4_time_integration.py` diff --git a/docs/repository-manifest.md b/docs/repository-manifest.md index 970debd..a2eb6ad 100644 --- a/docs/repository-manifest.md +++ b/docs/repository-manifest.md @@ -432,22 +432,33 @@ identified yet). `assembly.py` is no longer the maintainer-decided exception the previous version of this paragraph described -- its registry now resolves every MVP name to a real scheme; see that module's own docstring for the full retirement history. -**`PISO` is a single, real, dt-scaled pressure-correction pass, not the -full multi-pass Issa algorithm -- an honestly-scoped limitation, found -and resolved by numerical investigation before any implementation code -was written**: PyFlow's collocated mesh needs Rhie-Chow interpolation -(and momentum-equation coefficients this task's own interface has no way -to obtain) to suppress pressure-velocity decoupling under repeated +**`PISO` was a single, real, dt-scaled pressure-correction pass through +TASK-027 (2026-08-27) -- an honestly-scoped limitation, found and +resolved by numerical investigation before any implementation code was +written**: PyFlow's collocated mesh needs Rhie-Chow interpolation (and +momentum-equation coefficients that task's own interface had no way to +obtain) to suppress pressure-velocity decoupling under repeated correction, which the mathematically obvious approach -- composing `GreenGaussGradient`/`GreenGaussDivergence` into a Poisson matrix -- -cannot substitute for (proven not symmetric, both algebraically via the -discrete integration-by-parts identity and numerically, so -`ConjugateGradientSolver` cannot even solve it). That stronger, -fully-converged claim is Stage 5 TASK-033's own, not this task's -- -`docs/planning/roadmap.md` TASK-027's own Design decision Two records -the full investigation, and `docs/practices.md` gained a new standing -rule from this finding ("A criterion whose strong reading depends on a -later task must say so when drafted"). **`assembly.py`'s advection/ +could not substitute for (proven not symmetric, both algebraically via +the discrete integration-by-parts identity and numerically, so +`ConjugateGradientSolver` could not even solve it). `docs/planning/ +roadmap.md` TASK-027's own Design decision Two records the full +investigation, and `docs/practices.md` gained a new standing rule from +this finding ("A criterion whose strong reading depends on a later task +must say so when drafted"). **`PISO` is genuinely multi-pass as of +TASK-033 (Stage 5, 2026-08-29):** the momentum coefficient Rhie-Chow +needs resolved to `a_P = V/dt` (the only contribution to `a_P` for +PyFlow's fully explicit RK4 predictor, and one constant for the whole +mesh given its uniform cell volume), paired with the same compact +Laplacian the Poisson matrix already uses rather than the composed +`Gradient`/`Divergence` pair that failed -- `correct` now loops, +recording each pass's own maximum divergence and raising +`DivergenceDidNotConvergeError` if `numerics.pressure_correction_ +max_iterations` is exhausted before reaching `numerics.pressure_ +correction_tolerance`. See `src/pyflow/engine/CLAUDE.md`'s own `PISO` +entry and `docs/architecture/icds.md`'s Pressure-Velocity Coupling entry +for the full content. **`assembly.py`'s advection/ diffusion factories gained a `boundary_conditions` parameter and are now resolved after boundary conditions, not before (TASK-040)** -- a concrete scheme is constructed with the boundary conditions it needs, diff --git a/src/pyflow/configuration/CLAUDE.md b/src/pyflow/configuration/CLAUDE.md index d317b41..e0a142c 100644 --- a/src/pyflow/configuration/CLAUDE.md +++ b/src/pyflow/configuration/CLAUDE.md @@ -269,6 +269,27 @@ longer among them, since TASK-023/024), only a reference implementation resolves it under `src/`. Completes the six-field `numerics` section `adr/ADR-003-modular-numerical-strategies.md` names. +**`pressure_correction_tolerance`/`pressure_correction_max_iterations`** +(TASK-033, added 2026-08-29, Stage 5's fourth task): the *outer* +corrector loop's own tunables, following `linear_solver_tolerance`/ +`linear_solver_max_iterations`'s own precedent exactly (plain positive +numbers, `validate()` rejects `<= 0` for each independently) but +naming a genuinely different thing -- `linear_solver_tolerance` governs +each pass's *inner* linear solve, these two govern how many passes the +corrector loop itself takes and how small the recorded divergence must +get before it stops. Stage 5's own design question three ("outer-loop +state the strategy owns," `docs/planning/roadmap.md` TASK-033) is why +these live in `numerics:` as two more scalar fields rather than as a +`PressureCoupling.correct` parameter -- the split TASK-041 already drew +between `numerics:` (numerical parameters) and `fluid:` (physical +properties) put these on the `numerics:` side without a new question, +since a corrector loop's own iteration budget is exactly the kind of +thing `linear_solver_tolerance`/`linear_solver_max_iterations` already +established belongs there. `1e-6`/`50` are arbitrary MVP defaults, the +same reasoning as `timestep`'s `0.01` -- `50` deliberately smaller than +`linear_solver_max_iterations`'s own `1000`, since an outer corrector +pass is far more expensive than one CG iteration. + **`diffusion_coefficient`** (TASK-024, added 2026-08-27; **migrated to `FluidConfig.diffusion_coefficient` by TASK-041, 2026-08-28 -- no longer a `NumericsConfig` field**): originally followed `timestep`'s own diff --git a/src/pyflow/configuration/schema.py b/src/pyflow/configuration/schema.py index aff39a1..7ea6ddd 100644 --- a/src/pyflow/configuration/schema.py +++ b/src/pyflow/configuration/schema.py @@ -670,6 +670,17 @@ class NumericsConfig: setting `numerics.diffusion_coefficient` is rejected with a named error pointing at its new home (`loader.py`'s `_numerics_config_from_raw`), not silently ignored. + + **`pressure_correction_tolerance`/`pressure_correction_max_iterations` + (TASK-033, added 2026-08-29) are the corrector *loop*'s own tunables + -- distinct from `linear_solver_tolerance`/`linear_solver_max_iterations` + above, which govern the *inner* linear solve each corrector pass + makes, not how many passes the outer loop itself may take.** Same + plain-positive-number pattern as every other tunable here; `PISO`'s + own constructor is threaded these via `assemble_numerics` + (`register_pressure_coupling`'s own widened factory), the "outer-loop + state the strategy owns" resolution to Stage 5's design question + three (`docs/planning/roadmap.md` TASK-033). """ advection: AdvectionSchemeName = "first_order_upwind" @@ -680,6 +691,8 @@ class NumericsConfig: linear_solver_tolerance: float = 1e-6 linear_solver_max_iterations: int = 1000 pressure_coupling: PressureCouplingName = "piso" + pressure_correction_tolerance: float = 1e-6 + pressure_correction_max_iterations: int = 50 boundary_conditions: BoundaryConditionsConfig = field(default_factory=BoundaryConditionsConfig) def validate(self) -> None: @@ -720,6 +733,16 @@ def validate(self) -> None: f"numerics.pressure_coupling must be one of " f"{sorted(_VALID_PRESSURE_COUPLINGS)}, got {self.pressure_coupling!r}" ) + if self.pressure_correction_tolerance <= 0: + raise ValueError( + f"numerics.pressure_correction_tolerance must be > 0, " + f"got {self.pressure_correction_tolerance!r}" + ) + if self.pressure_correction_max_iterations <= 0: + raise ValueError( + f"numerics.pressure_correction_max_iterations must be > 0, " + f"got {self.pressure_correction_max_iterations!r}" + ) self.boundary_conditions.validate() diff --git a/src/pyflow/engine/CLAUDE.md b/src/pyflow/engine/CLAUDE.md index 5e26e08..0003815 100644 --- a/src/pyflow/engine/CLAUDE.md +++ b/src/pyflow/engine/CLAUDE.md @@ -993,6 +993,72 @@ divergence-free/divergent constant-vs-non-constant pressure pair are test_pressure_field.py` -- properties this class's own math already had, proven directly against it rather than reimplemented. +**`PISO` becomes genuinely multi-pass with TASK-033 (Stage 5, 2026-08-29) +-- the collocated-grid limitation this entry's own TASK-027 paragraph +recorded above is resolved, not superseded.** `correct` is now a real +loop, constructed with `tolerance: float = 1e-6, max_iterations: int = +50` (the "strategy owns its own tunables" shape `ConjugateGradientSolver` +already established, not a widened `PressureCoupling.correct` -- that +interface's own abstract signature is unchanged). Each pass: compute +`GreenGaussDivergence`'s divergence with a Rhie-Chow correction applied +(`_rhie_chow_divergence`, below), record its maximum magnitude in +`self.last_divergence_history`, and either return (at or below +`tolerance`) or solve another pressure correction against the same +Poisson matrix (`_poisson_matrix`, extracted unchanged from TASK-027's +own single-pass logic, built once per `correct()` call) -- up to +`max_iterations` passes before raising `DivergenceDidNotConvergeError` +rather than returning a best-effort result, distinct from +`PressureSolveDidNotConvergeError` (one pass's own inner linear solve +failing) the way an outer loop giving up differs from an inner solve +failing. + +**Design question three's answer, found by numerical prototyping before +any implementation code was written, the same sequence TASK-026/027 +used:** the momentum-equation coefficient Rhie-Chow needs is `a_P = +V/dt` -- the unsteady term is the *only* contribution to `a_P` for +PyFlow's fully explicit RK4 predictor (no implicit assembly to draw a +richer coefficient from), and PyFlow's uniform cell volume makes `V/a_P += dt` one constant for the whole mesh, not a per-cell field to compute. +`_rhie_chow_divergence(velocity, pressure, pressure_gradient, dt)` uses +it directly: at every *interior* face, `GreenGaussDivergence`'s own +simple-averaged divergence is corrected by `dt * [(p_N - p_P) / distance +- avg(gradP_owner, gradP_neighbour) . n]` -- the mismatch between the +direct face-normal pressure difference and the average of each +neighbour's own cell-centred gradient, accumulated back onto cells via +`accumulate_flux_to_cells`. Zero correction at every boundary face -- +there is no neighbour to interpolate against, and `GreenGaussDivergence`'s +own boundary handling already supplies the correct value there. + +**The real difference from TASK-027's own three failed correction +strategies is which Laplacian the correction pairs with, not merely +that a coefficient now exists.** TASK-027 measured that composing its +own `GreenGaussGradient`/`GreenGaussDivergence` into a Poisson matrix +produces one that is provably not symmetric. `_rhie_chow_divergence` +instead pairs the `a_P = V/dt` correction with the *same* compact, +already-symmetric Laplacian `_poisson_matrix` builds +(`CentralDifferenceDiffusion`-based) -- this is what restores the +discrete adjoint property, verified numerically before being trusted: a +manufactured provisional velocity field (both a linear and a nonlinear +fixture, disposable prototype scripts, not committed) converges to +floating-point-exact zero divergence in a single corrector pass once the +operators match, matching TASK-027's own measured ~54% single-pass +reduction as the correct baseline when they don't. + +**This single-pass-to-convergence result on PyFlow's own uniform MVP +mesh is also why `tests/unit/test_pressure_correction_loop.py`'s own +multi-pass scenario needs a deliberately handicapped `_HalvingSolver` +test double rather than a real `LinearSolver`** -- a real solver +converges too fast on this mesh to demonstrate genuine multi-pass +behaviour on its own. See that module's own entry in +`tests/unit/CLAUDE.md` for the double itself. + +Its own physical-correctness claims (a non-increasing recorded divergence +sequence reaching tolerance; genuine multi-pass convergence under a +solver that only partially corrects each pass; honest exhaustion rather +than a best-effort result) are `tests/features/ +pressure_correction_loop.feature`, bound by `tests/unit/ +test_pressure_correction_loop.py`. + **`gradient.py`/`divergence.py`** (TASK-018, Stage 3, interface-only until TASK-027) hold `GradientScheme`/`DivergenceScheme` -- two of the three operators (with `source.py`) that jointly compute the Flux layer diff --git a/src/pyflow/engine/numerics/CLAUDE.md b/src/pyflow/engine/numerics/CLAUDE.md index 41ee134..49071da 100644 --- a/src/pyflow/engine/numerics/CLAUDE.md +++ b/src/pyflow/engine/numerics/CLAUDE.md @@ -166,6 +166,25 @@ transported scalar's `fluid.diffusion_coefficient`, dispatched by `src/pyflow/engine/CLAUDE.md`'s own `CentralDifferenceDiffusion`/ `assembly.py` entries for the real content. +**`pressure_coupling.py` gained a second real occupant of +`_resolve_with_four_arguments`, the same day (TASK-033, 2026-08-29): +`register_pressure_coupling`'s own factory widens from two arguments to +four**, `pressure_correction_tolerance`/`pressure_correction_max_ +iterations` joining the existing `linear_solver`/`boundary_conditions` +pair -- the outer corrector loop's own tunables, Stage 5's own design +question three resolved as "outer-loop state the strategy owns," not a +`PressureCoupling.correct` signature widening. `_resolve_with_two_ +arguments`'s own docstring, which used to list pressure_coupling as one +of its three users, is corrected in the same change -- a docstring that +enumerates its own callers goes stale exactly the way any other restated +fact does, and this one was found stale by this task's own Blast Radius +check, not caught by any test. See `src/pyflow/engine/CLAUDE.md`'s own +`PISO` entry for the real content: `a_P = V/dt`, pairing the correction +with the same compact Laplacian the Poisson matrix already uses (not the +composed `Gradient`/`Divergence` pair TASK-027 tried and measured +failing), and `DivergenceDidNotConvergeError` as the outer loop's own +honest-exhaustion counterpart to `PressureSolveDidNotConvergeError`. + Full design rationale -- why a subpackage, why every operator takes `Field` rather than a concrete subclass, why the return shapes split face-valued (Advection/Diffusion) from cell-valued diff --git a/src/pyflow/engine/numerics/assembly.py b/src/pyflow/engine/numerics/assembly.py index f95eb23..7a22d7b 100644 --- a/src/pyflow/engine/numerics/assembly.py +++ b/src/pyflow/engine/numerics/assembly.py @@ -154,7 +154,7 @@ class AssembledNumerics: _time_integrator_registry: dict[str, Callable[[], TimeIntegrator]] = {} _linear_solver_registry: dict[str, Callable[[float, int], LinearSolver]] = {} _pressure_coupling_registry: dict[ - str, Callable[[LinearSolver, Mapping[str, BoundaryCondition]], PressureCoupling] + str, Callable[[LinearSolver, Mapping[str, BoundaryCondition], float, int], PressureCoupling] ] = {} _boundary_condition_registry: dict[str, Callable[[BoundaryFaceConfig], BoundaryCondition]] = {} @@ -234,16 +234,23 @@ def register_linear_solver(name: str, factory: Callable[[float, int], LinearSolv def register_pressure_coupling( name: str, - factory: Callable[[LinearSolver, Mapping[str, BoundaryCondition]], PressureCoupling], + factory: Callable[ + [LinearSolver, Mapping[str, BoundaryCondition], float, int], PressureCoupling + ], ) -> None: - """Make `name` resolve to `factory(linear_solver, boundary_conditions)` + """Make `name` resolve to `factory(linear_solver, boundary_conditions, + pressure_correction_tolerance, pressure_correction_max_iterations)` in future `assemble_numerics` calls -- `boundary_conditions` is the same face-name-keyed mapping `AssembledNumerics.boundary_conditions` carries (velocity's own boundary conditions), added in TASK-027 so a concrete strategy's own `DivergenceScheme` can be boundary-aware at construction, the same "constructed with it, not handed it after the fact" reasoning `boundary_conditions`/`diffusion_coefficient` already - established for advection/diffusion. + established for advection/diffusion. The tolerance/iterations pair + (TASK-033, added 2026-08-29) is `NumericsConfig. + pressure_correction_tolerance`/`pressure_correction_max_iterations` -- + a corrector *loop*'s own outer convergence tunables, "outer-loop state + the strategy owns" (Stage 5's own design question three). """ _register(_pressure_coupling_registry, name, factory, "pressure_coupling") @@ -270,15 +277,16 @@ def _resolve_with_two_arguments[T, A, B]( ) -> T: """Same as `_resolve`, for the components whose factory needs two constructor arguments -- advection (`boundary_conditions`, - `periodic_pairs`, TASK-030), linear_solver (`tolerance`, - `max_iterations`) and pressure_coupling (the resolved `LinearSolver`, - `boundary_conditions`) -- rather than several near-identical inline + `periodic_pairs`, TASK-030) and linear_solver (`tolerance`, + `max_iterations`) -- rather than several near-identical inline get/raise/call blocks repeating the same lookup (found during TASK-040's own review cycle). **Not the one-argument helper this docstring used to describe advection sharing with diffusion** -- that helper (`_resolve_with_argument`) was retired the same day advection itself gained a second constructor argument, leaving it - with no remaining caller. + with no remaining caller. **pressure_coupling used to be a third user, + until TASK-033 (2026-08-29) widened its own factory to four + arguments** -- see `_resolve_with_four_arguments` below. """ factory = registry.get(name) if factory is None: @@ -295,14 +303,17 @@ def _resolve_with_four_arguments[T, A, B, C, D]( argument_d: D, component: str, ) -> T: - """Same as `_resolve_with_two_arguments`, for diffusion alone -- the - one component whose factory needs four constructor arguments + """Same as `_resolve_with_two_arguments`, for the two components whose + factory needs four constructor arguments rather than two: diffusion (`boundary_conditions`, `periodic_pairs`, `diffusion_coefficient`, - `coefficient_overrides` since TASK-031b, 2026-08-29) rather than two. - Kept as its own generic helper instead of widening - `_resolve_with_two_arguments` itself, since advection/linear_solver/ - pressure_coupling still only need two arguments each and a shared - four-argument signature would force all three to pass unused ones. + `coefficient_overrides` since TASK-031b, 2026-08-29) and + pressure_coupling (`linear_solver`, `boundary_conditions`, + `pressure_correction_tolerance`, `pressure_correction_max_iterations` + since TASK-033, 2026-08-29 -- the outer corrector loop's own tunables, + "outer-loop state the strategy owns"). Kept as its own generic helper + instead of widening `_resolve_with_two_arguments` itself, since + advection/linear_solver still only need two arguments each and a + shared four-argument signature would force both to pass unused ones. **Replaces `_resolve_with_three_arguments`, TASK-030's own three-argument helper** -- diffusion was its only caller, and TASK-031b's own fourth argument left it with none; deleted in the same change as genuinely @@ -399,11 +410,13 @@ def assemble_numerics( config.linear_solver_max_iterations, "linear_solver", ) - pressure_coupling = _resolve_with_two_arguments( + pressure_coupling = _resolve_with_four_arguments( _pressure_coupling_registry, config.pressure_coupling, linear_solver, boundary_conditions, + config.pressure_correction_tolerance, + config.pressure_correction_max_iterations, "pressure_coupling", ) diff --git a/src/pyflow/engine/numerics/pressure_coupling.py b/src/pyflow/engine/numerics/pressure_coupling.py index 8b68cda..34741c7 100644 --- a/src/pyflow/engine/numerics/pressure_coupling.py +++ b/src/pyflow/engine/numerics/pressure_coupling.py @@ -47,6 +47,7 @@ import torch from pyflow.engine.field import Field +from pyflow.engine.mesh import StructuredCartesianMesh from pyflow.engine.numerics.boundary_condition import BoundaryCondition from pyflow.engine.numerics.diffusion import CentralDifferenceDiffusion from pyflow.engine.numerics.divergence import GreenGaussDivergence @@ -105,43 +106,103 @@ def evaluate(self, field: Field, face: int) -> float: class PISO(PressureCoupling): - """A single dt-scaled pressure-correction pass (TASK-027): solves the - Poisson equation `Laplacian(p) = div(u*) / dt` for the pressure - correction, then returns `u* - dt * grad(p)`. - - **Registered under `"piso"`, but not the full multi-pass Issa PISO - algorithm** -- honestly scoped to what a single pass, checked in - isolation, can actually deliver on PyFlow's collocated mesh. See - `docs/planning/roadmap.md` TASK-027's own Design decision Two for the - full numerical investigation: composing this task's own - `GreenGaussGradient`/`GreenGaussDivergence` into a Poisson matrix - produces one that is provably not symmetric (confirmed both - algebraically, via the discrete integration-by-parts identity, and - numerically), so it cannot be solved by `ConjugateGradientSolver` -- - the only registered `LinearSolver`, which requires a symmetric matrix. - Genuinely suppressing pressure-velocity decoupling under *repeated* - correction needs Rhie-Chow interpolation, which needs momentum- - equation coefficients this class's own interface has no way to - obtain -- that is Stage 5 TASK-033's own claim (Pressure Correction - Loop), not this one's. - - **The Poisson matrix is built via `CentralDifferenceDiffusion`** + """A genuine corrector *loop* (TASK-033, Stage 5): repeats a + dt-scaled pressure-correction pass -- solve `Laplacian(dp) = + div(u) / dt`, accumulate `p += dp`, correct `u -= dt * grad(dp)` -- + until the corrected velocity's own divergence reaches `tolerance` or + `max_iterations` is exhausted. **Registered under `"piso"`, and + genuinely PISO (Pressure-Implicit with Splitting of Operators) for + the first time**: TASK-027 (Stage 4) registered the name for a single + pass, honestly scoped and documented as not yet the real algorithm + (see that task's own Design decision Two, `docs/planning/roadmap.md`) + -- this is the task that closes that gap, per Stage 5 Completion + Criterion 3's own "whether `piso` is now genuinely multi-pass, or has + been renamed" instruction. + + **Design question three's answer, found by numerical prototyping + (`docs/planning/roadmap.md` TASK-033), not reasoned about in + advance** -- TASK-027 already showed that guessing here produces a + confident wrong answer: composing `GreenGaussGradient`/ + `GreenGaussDivergence` into a Poisson matrix is provably not + symmetric, and three correction strategies tried without momentum + coefficients each left most of the original divergence in place + (46-54% reduction for a single naive pass). **The missing momentum + coefficient `a_P` is simply `dt`.** Rhie-Chow interpolation's own + weight is `V / a_P`; a real momentum equation's `a_P` bundles the + unsteady term (`V / dt`) with advection/diffusion's own implicit + contributions, but PyFlow's momentum predictor is fully explicit + (RK4, no pressure term, Stage 5's own design decision that "the + correction sits outside the integrator, once per timestep") -- so the + only term `a_P` has to carry here *is* the unsteady one, `a_P = V / + dt`, and `V / a_P` cancels to exactly `dt`. Because PyFlow's mesh has + uniform cell volume (`docs/implementation/mvp.md`), `dt` is the same + constant for every cell -- no per-cell momentum operator, no widened + interface, no ADR. **Verified directly with disposable prototype + scripts before writing any test or implementation code** (not + committed, the same discipline TASK-026/027 both used): a Rhie-Chow + face-velocity correction built from this `dt` weight, paired with + the *same* compact `CentralDifferenceDiffusion`-based Laplacian used + for both the correction term's own face-pressure difference and the + Poisson matrix (rather than a separately-composed Gradient/Divergence + pair, which TASK-027 already proved is not the discrete adjoint of + itself), restores the discrete integration-by-parts identity exactly + -- confirmed by driving a manufactured, non-axis-aligned provisional + velocity field's divergence (measured Rhie-Chow-consistently, not by + `GreenGaussDivergence`'s own naive face averaging, which is precisely + the measure that does *not* see this) to floating-point zero in a + single corrector pass, for both a linear and a genuinely nonlinear + fixture -- a dramatically different outcome from TASK-027's own + "genuine Rhie-Chow... converging far too slowly" finding, because + that attempt paired Rhie-Chow with the composed (non-adjoint) + Gradient/Divergence pair for the matrix instead of the compact one. + + **`tolerance`/`max_iterations` (constructor-bound, both new) are + "outer-loop state the strategy owns"** -- the third of Stage 5's own + three candidate answers to design question three (a widened + `PressureCoupling.correct`, a momentum operator handed in at + construction, or outer-loop state), chosen because it needs no + change to `PressureCoupling`'s own abstract signature at all: `dt` + was already `correct`'s own second parameter (`adr/ADR-009`), and + the loop's own tolerance/iteration-limit are exactly the kind of + tunable `ConjugateGradientSolver`'s own `tolerance`/`max_iterations` + already established the precedent for -- bound at construction, not + per call. Both default to `1e-6`/`50`, matching + `NumericsConfig.pressure_correction_tolerance`/ + `pressure_correction_max_iterations`'s own defaults, so every + existing call site that only passes `(linear_solver, + boundary_conditions)` keeps working unchanged. + + **`last_divergence_history` is the recorded per-iteration sequence** + (populated on both success and `DivergenceDidNotConvergeError`) -- + Criterion 3's own "the sequence asserted, not just its last value" + needs somewhere to read it back from without widening `correct`'s own + return type. + + **The Poisson matrix is still built via `CentralDifferenceDiffusion`** (`diffusion_coefficient=1.0`, pressure's own zero-gradient boundary condition on every wall), reusing TASK-024's already-tested, - already-symmetric compact Laplacian rather than composing this task's - own `GradientScheme`/`DivergenceScheme` into a matrix -- the specific - thing Design decision Two proved does not work. `GreenGaussDivergence` - still computes the Poisson equation's right-hand side (`div(u*)`, - using `provisional_velocity`'s own boundary conditions), and - `GreenGaussGradient` still computes the cell-centred pressure gradient - the returned corrected velocity is built from -- both real, - necessary, exercised uses, not decorative registrations. + already-symmetric compact Laplacian -- built once per `correct` call + and reused across every corrector pass within it, since it depends + only on the fixed mesh/boundary conditions, never on the current + velocity or pressure. `GreenGaussDivergence` still computes the + "simple-averaged" half of the Rhie-Chow-corrected divergence (the + correction term this task adds is a small, separate face loop, added + on top rather than duplicating `GreenGaussDivergence`'s own + boundary-aware logic); `GreenGaussGradient` still computes both the + cell-centred pressure gradient the velocity correction is built from + and the per-cell gradients the Rhie-Chow correction term needs. """ def __init__( - self, linear_solver: LinearSolver, boundary_conditions: Mapping[str, BoundaryCondition] + self, + linear_solver: LinearSolver, + boundary_conditions: Mapping[str, BoundaryCondition], + tolerance: float = 1e-6, + max_iterations: int = 50, ) -> None: super().__init__(linear_solver) + self._tolerance = tolerance + self._max_iterations = max_iterations pressure_boundary_conditions = MappingProxyType( {name: _ZeroGradientPressureCondition() for name in _PRESSURE_BOUNDARY_FACE_NAMES} ) @@ -152,38 +213,120 @@ def __init__( ) self._gradient = GreenGaussGradient(pressure_boundary_conditions) self._divergence = GreenGaussDivergence(boundary_conditions) + self.last_divergence_history: tuple[float, ...] = () def correct( self, provisional_velocity: VectorField, dt: float ) -> tuple[VectorField, ScalarField]: mesh = provisional_velocity.mesh - divergence = self._divergence.divergence(provisional_velocity) + assert isinstance(mesh, StructuredCartesianMesh) + matrix = self._poisson_matrix(mesh) + + velocity = provisional_velocity + pressure = PressureField(mesh, "pressure", initial_value=0.0) + history: list[float] = [] + + for iteration in range(self._max_iterations + 1): + gradient = self._gradient.gradient(pressure) + divergence = self._rhie_chow_divergence(velocity, pressure, gradient, dt) + max_divergence = float(divergence.abs().max()) + history.append(max_divergence) + if max_divergence <= self._tolerance: + self.last_divergence_history = tuple(history) + return velocity, pressure + if iteration == self._max_iterations: + break + + result = self._linear_solver.solve(matrix, -divergence / dt) + if not result.converged: + self.last_divergence_history = tuple(history) + raise PressureSolveDidNotConvergeError( + f"pressure correction did not converge in {result.iterations} iterations" + ) + correction = PressureField(mesh, "pressure_correction") + correction.values[:] = result.solution + + new_pressure = PressureField(mesh, "pressure") + new_pressure.values[:] = pressure.values + correction.values + pressure = new_pressure + + new_velocity = velocity.copy() + new_velocity.values[:] = velocity.values - dt * self._gradient.gradient(correction) + velocity = new_velocity + + self.last_divergence_history = tuple(history) + raise DivergenceDidNotConvergeError( + f"pressure correction loop did not reach tolerance {self._tolerance} within " + f"{self._max_iterations} iterations; divergence history: {history}" + ) + def _poisson_matrix(self, mesh: StructuredCartesianMesh) -> torch.Tensor: num_cells = mesh.num_cells matrix = torch.zeros((num_cells, num_cells), dtype=torch.float64) for column in range(num_cells): basis = ScalarField(mesh, "e") basis.values[column] = 1.0 matrix[:, column] = -accumulate_flux_to_cells(mesh, self._diffusion.flux(basis)) - - result = self._linear_solver.solve(matrix, -divergence / dt) - if not result.converged: - raise PressureSolveDidNotConvergeError( - f"pressure correction did not converge in {result.iterations} iterations" + return matrix + + def _rhie_chow_divergence( + self, + velocity: VectorField, + pressure: PressureField, + pressure_gradient: torch.Tensor, + dt: float, + ) -> torch.Tensor: + """`GreenGaussDivergence`'s own simple-averaged divergence, minus + a Rhie-Chow correction at every *interior* face: `dt * [(p_N - + p_P) / distance - avg(gradP, gradN) . n]` -- the mismatch between + the direct face-normal pressure difference and the average of + each neighbour's own cell-centred gradient, which is exactly what + the naive simple-averaged divergence cannot see and what makes it + fail to be the discrete adjoint of the compact Laplacian + `correct`'s own Poisson solve uses. Zero at every boundary face -- + there is no neighbour to Rhie-Chow-interpolate against, and + `GreenGaussDivergence`'s own boundary handling (this class's + `boundary_conditions`, velocity's own) already supplies the + correct value there. + """ + mesh = velocity.mesh + assert isinstance(mesh, StructuredCartesianMesh) + naive = self._divergence.divergence(velocity) + + correction_face = torch.zeros(mesh.num_faces, dtype=torch.float64) + for face in range(mesh.num_faces): + owner, neighbour = mesh.face_neighbours(face) + if neighbour is None: + continue + normal_x, normal_y = mesh.face_normal(face) + distance = mesh.face_centroid_distance(face) + direct = (pressure.value_at(neighbour) - pressure.value_at(owner)) / distance + gx_o, gy_o = float(pressure_gradient[owner, 0]), float(pressure_gradient[owner, 1]) + gx_n, gy_n = ( + float(pressure_gradient[neighbour, 0]), + float(pressure_gradient[neighbour, 1]), ) - - pressure = PressureField(mesh, "pressure") - pressure.values[:] = result.solution - - corrected = provisional_velocity.copy() - corrected.values[:] = provisional_velocity.values - dt * self._gradient.gradient(pressure) - return corrected, pressure + avg_normal = ((gx_o + gx_n) / 2) * normal_x + ((gy_o + gy_n) / 2) * normal_y + correction_face[face] = dt * (direct - avg_normal) + return naive - accumulate_flux_to_cells(mesh, correction_face) class PressureSolveDidNotConvergeError(RuntimeError): - """Raised when `PISO.correct`'s own pressure-correction solve fails - to converge -- returning its unconverged solution anyway would be + """Raised when one corrector pass's own inner linear solve fails to + converge -- returning its unconverged solution anyway would be exactly the "plausible-looking wrong answer" failure mode `docs/practices.md` names repeatedly (pan scale, mesh accessors, `ConjugateGradientSolver`'s own honest treatment of the same flag). """ + + +class DivergenceDidNotConvergeError(RuntimeError): + """Raised when `PISO.correct`'s own outer corrector loop exhausts + `max_iterations` without the divergence reaching `tolerance` -- a + different failure from `PressureSolveDidNotConvergeError` (every + inner linear solve can converge just fine while the outer loop still + fails to reduce divergence enough, e.g. an inner solver tolerance too + loose relative to the outer one). The same honesty: a best-effort + velocity field that never reached the configured tolerance is not + returned as if it had. + """ diff --git a/tests/features/pressure_correction_loop.feature b/tests/features/pressure_correction_loop.feature new file mode 100644 index 0000000..1748a96 --- /dev/null +++ b/tests/features/pressure_correction_loop.feature @@ -0,0 +1,47 @@ +# The acceptance criteria for Pressure Correction Loop (TASK-033, Stage +# 5's fourth task in build order). Not a golden demo -- no config file +# under `examples/golden-demos/`, no CLI subprocess run, since every +# claim here is checked against the engine mechanism directly, the same +# `tests/unit/` shape every prior numerical-scheme feature file in this +# stage established. `tests/unit/test_pressure_correction_loop.py` binds +# these scenarios. +# +# The claim is that divergence decreases monotonically with each +# corrector iteration and reaches the configured tolerance -- measured +# across iterations, not asserted at the end. `PISO` (`engine/numerics/ +# pressure_coupling.py`) is the strategy this task makes genuinely +# multi-pass; see its own docstring for design question three's answer +# (found by numerical prototyping, not reasoned about in advance) and +# why `PressureCoupling.correct`'s own signature needed no change. +# +# Couette flow's exact linear profile (Stage 5 Completion Criterion 5) +# is deliberately not a scenario here -- it is discharged jointly with +# TASK-034, the first task with a fully assembled timestep to run it +# against; this task lays the corrector-loop groundwork it needs, not +# the demonstration itself. + +Feature: Pressure Correction Loop + + Background: + Given a small, non-square, non-trivially-origined mesh + And a provisional velocity field with real interior divergence, not aligned with either mesh axis + + Scenario: A corrector loop against a real divergent field converges, with a non-increasing recorded divergence sequence + Given a real linear solver + When the field is corrected by PISO's own corrector loop + Then the recorded divergence sequence is non-increasing at every element + And its last element is at or below the configured tolerance + + Scenario: A corrector loop that only partially corrects divergence each pass still takes multiple genuine passes to converge + Given a linear solver that only ever removes half of the remaining divergence per pass + When the field is corrected by PISO's own corrector loop + Then the recorded divergence sequence has more than two elements + And each element is smaller than the one before it + And its last element is at or below the configured tolerance + + Scenario: Exhausting the corrector iteration limit without reaching tolerance raises rather than returning a best-effort result + Given a linear solver that reports success but never actually reduces divergence + And a corrector iteration limit of 3 + When the field is corrected by PISO's own corrector loop + Then a divergence-did-not-converge error is raised + And the recorded divergence sequence has exactly 4 elements diff --git a/tests/unit/CLAUDE.md b/tests/unit/CLAUDE.md index 24376e5..1865daf 100644 --- a/tests/unit/CLAUDE.md +++ b/tests/unit/CLAUDE.md @@ -260,6 +260,28 @@ directly through `GreenGaussGradient`, rather than re-running `PISO` end to end, since the claim is specifically about the gradient of a shifted field, not about the solve that produced it. +**`test_pressure_correction_loop.py` (TASK-033, added 2026-08-29) is the +twelfth, and Stage 5's third module in this lineage** -- Stage 5's +fourth task, binding `tests/features/pressure_correction_loop.feature`'s +three scenarios against the real `PISO` directly, now genuinely +multi-pass (see `src/pyflow/engine/CLAUDE.md`'s own `PISO` entry for +design question three's resolution). Same shape as every module before +it: its own `_Context` dataclass, its own local `_ZeroNormalVelocity` +double (the same fixture shape `test_piso_pressure_coupling.py`'s/ +`test_pressure_field.py`'s own identically-named doubles use), no +golden-demo config file or CLI run. **Two `LinearSolver` doubles do the +real work no real solver could demonstrate on this mesh**: a real +`ConjugateGradientSolver` converges the divergent fixture in one or two +passes, too fast to prove genuine multi-pass behaviour, so +`_HalvingSolver` (always returns exactly half of the true least-squares +correction, verified beforehand to produce an exact geometric halving of +the recorded divergence every pass) forces and proves several strictly- +decreasing passes; `_NoOpSolver` (reports `converged=True` but always +returns the zero vector) exercises the iteration-limit exhaustion path, +distinct from `test_piso_pressure_coupling.py`'s own +`_NeverConvergesSolver` (`converged=False`), since this scenario is +about the *outer* loop giving up, not the *inner* solve failing. + **The convention is "local by default, shared where genuinely identical" -- amended 2026-08-28 by the Stage 4 exit audit, which found the older blanket form ("each binding test supplies its own local diff --git a/tests/unit/numerics/CLAUDE.md b/tests/unit/numerics/CLAUDE.md index 8b0757d..d58b45a 100644 --- a/tests/unit/numerics/CLAUDE.md +++ b/tests/unit/numerics/CLAUDE.md @@ -300,3 +300,31 @@ to route through) is deleted in the same change as genuinely dead code, its only caller having moved to `_resolve_with_two_arguments` -- `src/pyflow/engine/numerics/CLAUDE.md`'s own entry has the full reasoning. + +**`test_pressure_coupling_contract.py`'s own `_StubLinearSolver` stopped +being a convenient zero (TASK-033, 2026-08-29), found necessary rather +than chosen.** `PISO`'s rewritten multi-pass `correct` (Stage 5's own +design question three, `src/pyflow/engine/CLAUDE.md`'s `PISO` entry) +genuinely checks whether a pass made progress, and the suite's own +`_StubLinearSolver` used to return `LinearSolverResult(solution=torch. +zeros_like(rhs), converged=True, iterations=0)` regardless of what it was +asked to solve -- a fake "solver" that reports success while computing +nothing, which the old single-pass `PISO` never looked closely enough at +to notice. Two of this suite's own generic tests broke against it once +`PISO` started genuinely checking, not because the tests were wrong: the +fixture had been silently relying on `PISO` never asking its stub solver +to actually solve anything. Fixed at the root, not worked around: `_StubLinearSolver` +now performs a real `torch.linalg.lstsq` solve (not `torch.linalg.solve` +-- `PISO`'s own Poisson matrix is singular by construction, a pure-Neumann +problem, and `lstsq` returns a well-defined minimum-norm solution where +`solve` would raise; confirmed directly in a standalone `torch` check +before committing to the fix). + +**`test_assembly.py`'s own `_CapturingPressureCoupling` gained two +ignored constructor parameters the same day (TASK-033)**, matching +`register_pressure_coupling`'s own widened four-argument factory shape +(`tolerance: float, max_iterations: int`, discarded via `del` the same +way every other test-only double in this file ignores an argument it +doesn't need) -- without this, registering it as a factory would fail +the moment `assemble_numerics` calls it with four arguments instead of +two. diff --git a/tests/unit/numerics/test_assembly.py b/tests/unit/numerics/test_assembly.py index 9a6717d..4b0b5dd 100644 --- a/tests/unit/numerics/test_assembly.py +++ b/tests/unit/numerics/test_assembly.py @@ -102,12 +102,19 @@ def flux(self, field: Field, velocity: VectorField) -> torch.Tensor: class _CapturingPressureCoupling(PressureCoupling): """Records the exact `boundary_conditions` mapping it was constructed with -- the pressure-coupling analogue of `_CapturingAdvection`/ - `_CapturingDiffusion` above (TASK-027). + `_CapturingDiffusion` above (TASK-027). Accepts (and ignores) + `tolerance`/`max_iterations` (TASK-033, 2026-08-29), matching + `register_pressure_coupling`'s new four-argument factory shape. """ def __init__( - self, linear_solver: LinearSolver, boundary_conditions: Mapping[str, BoundaryCondition] + self, + linear_solver: LinearSolver, + boundary_conditions: Mapping[str, BoundaryCondition], + tolerance: float, + max_iterations: int, ) -> None: + del tolerance, max_iterations super().__init__(linear_solver) self.received_boundary_conditions = boundary_conditions diff --git a/tests/unit/numerics/test_pressure_coupling_contract.py b/tests/unit/numerics/test_pressure_coupling_contract.py index c5af171..f13b726 100644 --- a/tests/unit/numerics/test_pressure_coupling_contract.py +++ b/tests/unit/numerics/test_pressure_coupling_contract.py @@ -40,11 +40,28 @@ class _StubLinearSolver(LinearSolver): """A minimal `LinearSolver` test double -- exists only so a `PressureCoupling` strategy has something real to be constructed - with; this suite makes no claim about solving. + with; this suite makes no claim about solving *strategy* (iterative + vs. direct, tolerance, etc.). + + **Performs a real direct solve (`torch.linalg.lstsq`, not + `torch.linalg.solve` -- `PISO`'s own Poisson matrix is singular by + construction, a pure-Neumann pressure problem, so a real solve has to + tolerate that the same way `ConjugateGradientSolver`'s own gated + null-space projection does), not a convenient zero -- found + necessary, not assumed, once `PISO` (TASK-033, Stage 5) became a + genuine corrector *loop*.** Returning the zero vector while claiming + `converged=True` was already dishonest in the shape + `docs/practices.md` warns about, but the old single-pass `PISO` never + checked whether a correction actually reduced divergence, so it went + unnoticed; the new loop does check, correctly refuses to accept a + "solve" that never improves anything, and exhausts its own iteration + limit -- exposing the stub for what it was rather than a bug in the + loop itself. """ def solve(self, matrix: torch.Tensor, rhs: torch.Tensor) -> LinearSolverResult: - return LinearSolverResult(solution=torch.zeros_like(rhs), converged=True, iterations=0) + solution = torch.linalg.lstsq(matrix, rhs.unsqueeze(-1)).solution.squeeze(-1) + return LinearSolverResult(solution=solution, converged=True, iterations=1) class _ZeroNormalVelocity(BoundaryCondition): diff --git a/tests/unit/test_configuration.py b/tests/unit/test_configuration.py index da44605..1e8ac70 100644 --- a/tests/unit/test_configuration.py +++ b/tests/unit/test_configuration.py @@ -46,6 +46,8 @@ def test_defaults_are_valid() -> None: assert config.numerics.linear_solver_tolerance == 1e-6 assert config.numerics.linear_solver_max_iterations == 1000 assert config.numerics.pressure_coupling == "piso" + assert config.numerics.pressure_correction_tolerance == 1e-6 + assert config.numerics.pressure_correction_max_iterations == 50 for boundary_name in ("north", "south", "east", "west"): face = getattr(config.numerics.boundary_conditions, boundary_name) assert face.type == "dirichlet" @@ -728,6 +730,43 @@ def test_load_config_rejects_an_unknown_pressure_coupling_strategy(tmp_path: Pat load_config(config_file) +def test_load_config_reads_pressure_correction_loop_tunables(tmp_path: Path) -> None: + # The outer corrector loop's own tolerance/max_iterations (TASK-033), + # distinct from linear_solver_tolerance/linear_solver_max_iterations + # above, which govern the inner solve. + config_file = tmp_path / "config.yaml" + config_file.write_text( + "numerics:\n" + " pressure_correction_tolerance: 0.001\n" + " pressure_correction_max_iterations: 20\n" + ) + + config = load_config(config_file) + + assert config.numerics.pressure_correction_tolerance == 0.001 + assert config.numerics.pressure_correction_max_iterations == 20 + + +def test_load_config_rejects_a_non_positive_pressure_correction_tolerance( + tmp_path: Path, +) -> None: + config_file = tmp_path / "config.yaml" + config_file.write_text("numerics:\n pressure_correction_tolerance: 0.0\n") + + with pytest.raises(ValueError, match="numerics.pressure_correction_tolerance"): + load_config(config_file) + + +def test_load_config_rejects_a_non_positive_pressure_correction_max_iterations( + tmp_path: Path, +) -> None: + config_file = tmp_path / "config.yaml" + config_file.write_text("numerics:\n pressure_correction_max_iterations: 0\n") + + with pytest.raises(ValueError, match="numerics.pressure_correction_max_iterations"): + load_config(config_file) + + def test_load_config_reads_boundary_conditions_section(tmp_path: Path) -> None: config_file = tmp_path / "config.yaml" config_file.write_text( diff --git a/tests/unit/test_main.py b/tests/unit/test_main.py index a80caff..658b2c8 100644 --- a/tests/unit/test_main.py +++ b/tests/unit/test_main.py @@ -134,6 +134,8 @@ def test_generate_config_with_no_output_prints_to_stdout( "linear_solver_tolerance": 1e-6, "linear_solver_max_iterations": 1000, "pressure_coupling": "piso", + "pressure_correction_tolerance": 1e-6, + "pressure_correction_max_iterations": 50, "boundary_conditions": { "north": { "type": "dirichlet", diff --git a/tests/unit/test_pressure_correction_loop.py b/tests/unit/test_pressure_correction_loop.py new file mode 100644 index 0000000..fea4423 --- /dev/null +++ b/tests/unit/test_pressure_correction_loop.py @@ -0,0 +1,218 @@ +"""Binds `tests/features/pressure_correction_loop.feature` (TASK-033) -- +Stage 5's fourth task in build order. `PISO`'s own corrector loop +mechanism (`engine/numerics/pressure_coupling.py`) is what this task +built; these scenarios prove its own three claims -- a non-increasing +recorded sequence reaching tolerance, genuine multi-pass convergence +when a single pass genuinely cannot finish the job, and honest failure +when the iteration limit is exhausted -- against the real class, not a +new mechanism. + +Not a golden demo -- no config file under `examples/golden-demos/`, no +CLI run. Lives here, not under `tests/golden/`, per this directory's own +scope: isolated logic, no process boundary (`tests/unit/CLAUDE.md`). +""" + +from __future__ import annotations + +from dataclasses import dataclass, field +from typing import Literal + +import torch +from pytest_bdd import given, scenarios, then, when + +from pyflow.engine.field import Field +from pyflow.engine.mesh import StructuredCartesianMesh +from pyflow.engine.numerics.boundary_condition import BoundaryCondition +from pyflow.engine.numerics.linear_solver import ( + ConjugateGradientSolver, + LinearSolver, + LinearSolverResult, +) +from pyflow.engine.numerics.pressure_coupling import PISO, DivergenceDidNotConvergeError +from pyflow.engine.vector_field import VectorField + +from ._numerics import default_mesh + +scenarios("pressure_correction_loop.feature") + +_TOLERANCE = 1e-4 + + +class _ZeroNormalVelocity(BoundaryCondition): + """Dirichlet, fixed at zero -- a closed box, the same fixture shape + `test_piso_pressure_coupling.py`'s own identically-named double uses. + """ + + @property + def kind(self) -> Literal["value", "gradient"]: + return "value" + + def evaluate(self, field: Field, face: int) -> float: + self._check_boundary_face(field, face) + return 0.0 + + +class _HalvingSolver(LinearSolver): + """Always returns exactly half of the true (least-squares) solution + -- a deterministic, hand-verifiable way to force the corrector loop + to take genuinely multiple passes, rather than relying on a real + iterative solver's own opaque partial-convergence behaviour. + + Because the pressure-correction system is linear, using half of the + exact correction each pass leaves exactly half of the *previous* + residual, not some solver-dependent fraction -- verified directly + with a disposable prototype script before writing this class (the + same discipline TASK-026/027 both used): the recorded divergence + sequence halves on every single pass to machine precision. Always + reports `converged=True` -- a real "solve" happened, it was just + deliberately incomplete, the same "converged but wrong" shape + `test_pressure_coupling_contract.py`'s own `_StubLinearSolver` found + itself accidentally producing before it was fixed to do a real solve. + """ + + def solve(self, matrix: torch.Tensor, rhs: torch.Tensor) -> LinearSolverResult: + exact = torch.linalg.lstsq(matrix, rhs.unsqueeze(-1)).solution.squeeze(-1) + return LinearSolverResult(solution=0.5 * exact, converged=True, iterations=1) + + +class _NoOpSolver(LinearSolver): + """Reports success but returns the zero vector always -- never + reduces divergence, so a corrector loop using it can only ever + exhaust its own iteration limit. Distinct from + `test_piso_pressure_coupling.py`'s own `_NeverConvergesSolver` + (which reports `converged=False`, testing the *inner* solve's own + honesty): this double's own inner solve always succeeds, so it is + the *outer* corrector loop's own exhaustion path this exercises, + not `PressureSolveDidNotConvergeError`. + """ + + def solve(self, matrix: torch.Tensor, rhs: torch.Tensor) -> LinearSolverResult: + return LinearSolverResult(solution=torch.zeros_like(rhs), converged=True, iterations=0) + + +def _divergent_provisional_velocity(mesh: StructuredCartesianMesh) -> VectorField: + # Neither axis-aligned nor uniform -- the same fixture shape + # `test_piso_pressure_coupling.py`'s own `_provisional_velocity` uses. + center_x, center_y = 1.0, -0.4 + + def value(x: float, y: float) -> tuple[float, float]: + return ( + 0.6 * (x - center_x) - 0.2 * (y - center_y), + 0.3 * (x - center_x) + 0.9 * (y - center_y), + ) + + return VectorField(mesh, "u_star", num_components=2, initial_value=value) + + +def _boundary_conditions() -> dict[str, BoundaryCondition]: + condition = _ZeroNormalVelocity() + return {"north": condition, "south": condition, "east": condition, "west": condition} + + +# -- Fixture context ------------------------------------------------------- + + +@dataclass +class _Context: + mesh: StructuredCartesianMesh + boundary_conditions: dict[str, BoundaryCondition] = field(default_factory=_boundary_conditions) + provisional_velocity: VectorField | None = None + linear_solver: LinearSolver | None = None + max_iterations: int = 50 + history: tuple[float, ...] | None = None + error: Exception | None = None + + +# -- Given ------------------------------------------------------------- + + +@given("a small, non-square, non-trivially-origined mesh", target_fixture="ctx") +def _given_default_mesh() -> _Context: + return _Context(mesh=default_mesh(extent=(5, 4))) + + +@given( + "a provisional velocity field with real interior divergence, not aligned with either mesh axis" +) +def _given_divergent_velocity(ctx: _Context) -> None: + ctx.provisional_velocity = _divergent_provisional_velocity(ctx.mesh) + + +@given("a real linear solver") +def _given_real_solver(ctx: _Context) -> None: + ctx.linear_solver = ConjugateGradientSolver(tolerance=1e-10, max_iterations=500) + + +@given("a linear solver that only ever removes half of the remaining divergence per pass") +def _given_halving_solver(ctx: _Context) -> None: + ctx.linear_solver = _HalvingSolver() + + +@given("a linear solver that reports success but never actually reduces divergence") +def _given_noop_solver(ctx: _Context) -> None: + ctx.linear_solver = _NoOpSolver() + + +@given("a corrector iteration limit of 3") +def _given_iteration_limit(ctx: _Context) -> None: + ctx.max_iterations = 3 + + +# -- When ------------------------------------------------------------------ + + +@when("the field is corrected by PISO's own corrector loop") +def _when_corrected(ctx: _Context) -> None: + assert ctx.provisional_velocity is not None + assert ctx.linear_solver is not None + piso = PISO( + ctx.linear_solver, + ctx.boundary_conditions, + tolerance=_TOLERANCE, + max_iterations=ctx.max_iterations, + ) + try: + piso.correct(ctx.provisional_velocity, dt=1.0) + except DivergenceDidNotConvergeError as exc: + ctx.error = exc + ctx.history = piso.last_divergence_history + + +# -- Then ------------------------------------------------------------------ + + +@then("the recorded divergence sequence is non-increasing at every element") +def _then_non_increasing(ctx: _Context) -> None: + assert ctx.history is not None + assert all(a >= b for a, b in zip(ctx.history, ctx.history[1:], strict=False)) + + +@then("its last element is at or below the configured tolerance") +def _then_last_below_tolerance(ctx: _Context) -> None: + assert ctx.history is not None + assert ctx.history[-1] <= _TOLERANCE + + +@then("the recorded divergence sequence has more than two elements") +def _then_more_than_two_elements(ctx: _Context) -> None: + assert ctx.history is not None + assert len(ctx.history) > 2, ctx.history + + +@then("each element is smaller than the one before it") +def _then_strictly_decreasing(ctx: _Context) -> None: + assert ctx.history is not None + assert all(a > b for a, b in zip(ctx.history, ctx.history[1:], strict=False)), ctx.history + + +@then("a divergence-did-not-converge error is raised") +def _then_divergence_error_raised(ctx: _Context) -> None: + assert isinstance(ctx.error, DivergenceDidNotConvergeError) + + +@then("the recorded divergence sequence has exactly 4 elements") +def _then_exactly_four_elements(ctx: _Context) -> None: + # max_iterations=3 corrector *passes* are attempted, plus the initial + # (pre-correction) measurement -- 4 divergence checks total. + assert ctx.history is not None + assert len(ctx.history) == 4, ctx.history diff --git a/tools/generators/generate_config_template.py b/tools/generators/generate_config_template.py index 9c6f3fd..3944903 100644 --- a/tools/generators/generate_config_template.py +++ b/tools/generators/generate_config_template.py @@ -242,6 +242,17 @@ 'Valid: "piso" -- the only scheme PyFlow currently implements ' "for this component. Invalid: any other string." ), + "numerics.pressure_correction_tolerance": ( + "Valid: a positive number -- the outer corrector loop's own " + "convergence tolerance (distinct from linear_solver_tolerance " + "above, which governs each inner linear solve, not how many " + "corrector passes the outer loop may take). Invalid: zero or " + "negative." + ), + "numerics.pressure_correction_max_iterations": ( + "Valid: a positive integer -- the outer corrector loop's own " + "iteration limit. Invalid: zero, negative, or a float." + ), "numerics.boundary_conditions..type": ( 'Valid: "dirichlet", "neumann", or "periodic". If "periodic", ' "the OPPOSITE face (north<->south, east<->west) must also be " diff --git a/tools/validators/check_references.py b/tools/validators/check_references.py index 76f85bb..328bbd0 100644 --- a/tools/validators/check_references.py +++ b/tools/validators/check_references.py @@ -110,7 +110,6 @@ # change is what that Stage's own design question one decides, and a # guess here would be the speculation P-016 refuses. PLANNED: dict[str, str] = { - "tests/features/pressure_correction_loop.feature": "TASK-033", "tests/features/navier_stokes_timestep.feature": "TASK-034", }