diff --git a/docs/planning/roadmap.md b/docs/planning/roadmap.md
index edef501..e299f53 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): **642 tests at 99% as of 2026-08-29**, having been 64 when
+(C1a/C1b): **647 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
@@ -413,8 +413,15 @@ fields across two sections rather than one -- and one new
 decomposes/steps/reassembles velocity, not only `simulation.step()`
 called directly (found necessary by its own coverage report: the new
 `velocity_solved` branches in `_add_passive_scalar_transport` were
-otherwise unexercised by anything in this run). **71 of those 642 are
-Gherkin scenarios rather than pytest functions**
+otherwise unexercised by anything in this run). 647 after TASK-032
+(Pressure Field, Stage 5's third task): five new Gherkin scenarios in
+`tests/unit/test_pressure_field.py` (`pressure_field.feature`), proving
+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
+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`
 joining, TASK-021; to 24 with TASK-040's own
@@ -477,7 +484,15 @@ unchanged by whether the velocity carrying it was solved or prescribed,
 self-advection matching a hand-derived result, the existing
 `IncompatibleVelocityFieldError` check surviving the new path, and the
 orchestrator's own source still carrying no field-name-specific
-branching for velocity).
+branching for velocity); and to 76 with TASK-032's own
+`pressure_field.feature`, five scenarios: a divergence-free provisional
+field yielding pressure constant to solver tolerance, a divergent one
+yielding non-constant pressure, adding a constant to the solved pressure
+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.
 **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`,
@@ -7196,6 +7211,32 @@ is TASK-041's**. Criterion 6 and Criterion 7, its own share.
 
 Pressure Field
 
+**Status: Done, 2026-08-29, Stage 5's third task.** `PISO` (TASK-027,
+Stage 4) already performs the solve this task's own criteria describe;
+this task adds no new pressure-solving mechanism, only the properties
+Stage 4's own criteria never had cause to check and the one genuinely
+missing API-level guard:
+
+- **The local design question resolved: pressure gets a type of its
+  own, `PressureField(ScalarField)`** (`src/pyflow/engine/scalar_field.py`)
+  -- a marker subclass with no behaviour of its own, `PISO.correct` now
+  constructs one instead of a plain `ScalarField`. `simulation.step`
+  gains a real `isinstance` check, raising `PressureFieldTransportError`
+  (a new class, `simulation.py`) if `fields` contains one -- the "stated
+  at the API level" shape Criterion 2 asks for, checked directly against
+  the object handed in rather than by name.
+- **`PressureField` lives in `scalar_field.py`, not
+  `engine/numerics/pressure_coupling.py`** (its one real producer) --
+  `simulation.py` needs to import it and cannot import from
+  `pressure_coupling.py` without a circular import (that module already
+  imports `accumulate_flux_to_cells` from `simulation.py`). The same
+  circular-import reasoning TASK-031a's own `IncompatibleVelocityFieldError`
+  move used.
+- **`PressureCoupling.correct`'s own abstract signature is unchanged**
+  (`-> tuple[VectorField, ScalarField]`) -- no Stage 3 interface change;
+  a `PressureField` instance already satisfies it (valid covariant
+  return narrowing under `mypy --strict`).
+
 **Intent:** pressure is *not* transported -- it is solved for, from the
 incompressibility constraint. A criterion that treats it as another
 advected scalar has misunderstood the task. See
@@ -7228,13 +7269,19 @@ empty finds an answer next to it.
 ### Artifacts Produced
 
 - `tests/features/pressure_field.feature` -- this task's Acceptance
-  Criteria.
-- Source artifacts named when this task is drafted. **TASK-031's own
-  design questions are no longer what this waits on** -- all three were
-  answered on 2026-08-28 -- but one local question genuinely remains
-  this task's to decide: whether pressure needs a type of its own, or is
-  a `ScalarField` the coupling owns. Nothing before it has cause to
-  choose, so it is not escalated as a stage-level design question.
+  Criteria. `tests/unit/test_pressure_field.py` binds them, per
+  `tests/unit/`'s own scope (isolated logic, no process boundary) --
+  every scenario is checked against the real `PISO`/`GreenGaussGradient`
+  machinery directly, none needs a CLI subprocess.
+- `src/pyflow/engine/scalar_field.py` -- `PressureField(ScalarField)`,
+  resolving this task's one local design question (a type of its own,
+  not a plain `ScalarField` the coupling owns).
+- `src/pyflow/engine/simulation.py` -- `step` gains an `isinstance`
+  guard against `PressureField`, raising a new
+  `PressureFieldTransportError`.
+- `src/pyflow/engine/numerics/pressure_coupling.py` -- `PISO.correct`
+  constructs a `PressureField` instead of a plain `ScalarField`; no
+  signature change.
 
 ### Acceptance Criteria
 
diff --git a/docs/planning/status.md b/docs/planning/status.md
index 1459959..490321a 100644
--- a/docs/planning/status.md
+++ b/docs/planning/status.md
@@ -15,14 +15,14 @@ demand, not part of this file.
 
 ## Progress
 
-**35/42 tasks complete (83%)** across 14 planned stages. For the full plan, including
+**36/42 tasks complete (86%)** 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" : 35
-    "Not started" : 7
+    "Done" : 36
+    "Not started" : 6
 ```
 
 ### Milestones
@@ -35,13 +35,13 @@ pie showData
 
 ### Up next
 
-**Stage 5 -- First Fluid Solver** is next, starting with TASK-032 (Pressure Field), 2 more not yet started in this stage.
+**Stage 5 -- First Fluid Solver** is next, starting with TASK-033 (Pressure Correction Loop), 1 more not yet started in this stage.
 
 ## Live repository facts
 
 - **45** `CLAUDE.md` files
-- **642** tests collected
-- **71** Gherkin scenarios (`tests/features/*.feature`)
+- **647** tests collected
+- **76** Gherkin scenarios (`tests/features/*.feature`)
 
 ## Stages
 
@@ -115,13 +115,13 @@ pie showData
 
 ### Stage 5 -- First Fluid Solver
 
-**no status recorded** -- `████░░░░░░` 2/5 tasks; 13 criteria defined, no status line yet
+**no status recorded** -- `██████░░░░` 3/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 -- Pressure Field | Not started |  |  |
+| TASK-032 | Done | 2026-08-29 | `src/pyflow/engine/scalar_field.py` |
 | TASK-033 -- Pressure Correction Loop | Not started |  |  |
 | TASK-034 -- Navier-Stokes Timestep | Not started |  |  |
 
diff --git a/docs/repository-inventory.md b/docs/repository-inventory.md
index 2a2b31b..dd8becc 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.
 
-**274 tracked files** across 45 directories;
+**276 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_field.feature`
 - `rk4_time_integration.feature`
 - `simulation_orchestrator.feature`
 - `velocity_field_support.feature`
@@ -381,6 +382,7 @@ listing files.
 - `test_neumann_boundary.py`
 - `test_periodic_boundary.py`
 - `test_piso_pressure_coupling.py`
+- `test_pressure_field.py`
 - `test_rendering.py`
 - `test_rk4_time_integration.py`
 - `test_scalar_field.py`
diff --git a/src/pyflow/engine/CLAUDE.md b/src/pyflow/engine/CLAUDE.md
index b1ebd03..5e26e08 100644
--- a/src/pyflow/engine/CLAUDE.md
+++ b/src/pyflow/engine/CLAUDE.md
@@ -263,6 +263,22 @@ implementation without per-implementation casts.
 claims (plain `float` return, exact formula, copy independence for its
 own storage).
 
+**`PressureField(ScalarField)` (TASK-032, added 2026-08-29) is a thin
+marker subclass, no new behaviour at all.** Stage 5 Completion
+Criterion 2's own "pressure is not among the fields `step` advances"
+needed pressure to be identifiable at the API level -- a real
+`isinstance` check, since there is no configuration surface that names
+which fields are transported. `PISO.correct`
+(`engine/numerics/pressure_coupling.py`) constructs one instead of a
+plain `ScalarField`; `simulation.step` raises
+`PressureFieldTransportError` if `fields` contains one. Lives here, not
+in `pressure_coupling.py` (its one real producer), because
+`simulation.py` needs to import it and cannot import from
+`pressure_coupling.py` without a circular import -- the identical
+reasoning `IncompatibleVelocityFieldError`'s own TASK-031a move to
+`vector_field.py` used, applied to the same shape of problem a second
+time.
+
 **There are two contract suites here, one per interface, and the split
 is load-bearing** (2026-08-22, Stage 2 exit audit -- `docs/planning/
 roadmap.md`). `field.py` deliberately carries no storage so that a
@@ -962,6 +978,21 @@ boundedly reduces divergence; non-convergence is reported, not returned
 as a plausible answer) is `tests/features/piso_pressure_coupling.feature`,
 bound by `tests/unit/test_piso_pressure_coupling.py`.
 
+**`correct` constructs a `PressureField`, not a plain `ScalarField`,
+since TASK-032 (Stage 5, 2026-08-29).** No signature change --
+`PressureCoupling.correct`'s own abstract return type stays
+`tuple[VectorField, ScalarField]`, and a `PressureField` instance
+already satisfies it (covariant return narrowing). This is what makes
+`simulation.step`'s own new `PressureFieldTransportError` guard
+meaningful: `PISO`'s real output is now identifiable as pressure at the
+API level, not just by the string `"pressure"` its constructor happens
+to name it. Its own null-space property (adding a constant to the
+solved pressure leaves the corrected velocity unchanged) and the
+divergence-free/divergent constant-vs-non-constant pressure pair are
+`tests/features/pressure_field.feature`, bound by `tests/unit/
+test_pressure_field.py` -- properties this class's own math already had,
+proven directly against it rather than reimplemented.
+
 **`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
@@ -1292,6 +1323,16 @@ the buildable reading chosen instead, stated explicitly per root
 `CLAUDE.md`'s Integrity section -- see TASK-040's own entry in
 `docs/planning/roadmap.md` for the full reasoning.
 
+**`PressureFieldTransportError` (TASK-032, added 2026-08-29) is
+`step`'s second rejection, checked in the same per-field loop as
+`MismatchedMeshError`, before either.** Raised if `fields` contains a
+`PressureField` (`scalar_field.py`, above) -- Stage 5 Completion
+Criterion 2's own claim that pressure is solved, not transported, made
+real: a config-level check was not possible (`SimulationConfig` names no
+surface for "which fields are transported" at all), so this is a real
+`isinstance` check against the object itself, the same "stated at the
+API level" shape the criterion asks for.
+
 **`simulation_orchestrator.feature` is not a golden demo** -- no config
 file under `examples/golden-demos/`, no CLI subprocess run, since this
 is the mechanism a future demo (TASK-030) is built on top of, not a demo
diff --git a/src/pyflow/engine/numerics/pressure_coupling.py b/src/pyflow/engine/numerics/pressure_coupling.py
index fbc4781..8b68cda 100644
--- a/src/pyflow/engine/numerics/pressure_coupling.py
+++ b/src/pyflow/engine/numerics/pressure_coupling.py
@@ -52,7 +52,7 @@
 from pyflow.engine.numerics.divergence import GreenGaussDivergence
 from pyflow.engine.numerics.gradient import GreenGaussGradient
 from pyflow.engine.numerics.linear_solver import LinearSolver
-from pyflow.engine.scalar_field import ScalarField
+from pyflow.engine.scalar_field import PressureField, ScalarField
 from pyflow.engine.simulation import accumulate_flux_to_cells
 from pyflow.engine.vector_field import VectorField
 
@@ -172,7 +172,7 @@ def correct(
                 f"pressure correction did not converge in {result.iterations} iterations"
             )
 
-        pressure = ScalarField(mesh, "pressure")
+        pressure = PressureField(mesh, "pressure")
         pressure.values[:] = result.solution
 
         corrected = provisional_velocity.copy()
diff --git a/src/pyflow/engine/scalar_field.py b/src/pyflow/engine/scalar_field.py
index b33f4f3..692798b 100644
--- a/src/pyflow/engine/scalar_field.py
+++ b/src/pyflow/engine/scalar_field.py
@@ -1,6 +1,10 @@
 """ScalarField (TASK-015): a single value per cell -- the collocated
 scalar leaf of the `Field` hierarchy. See `collocated_field.py` for the
 shared storage/initialisation/access logic this builds on.
+
+**`PressureField` (TASK-032, added 2026-08-29)** is a thin marker
+subclass -- see its own docstring below for why it exists and why it
+lives here rather than in `engine/numerics/pressure_coupling.py`.
 """
 
 from __future__ import annotations
@@ -33,3 +37,34 @@ def copy(self) -> ScalarField:
         clone = ScalarField(self.mesh, self.name)
         clone._values = self._values.clone()
         return clone
+
+
+class PressureField(ScalarField):
+    """A `ScalarField` that is specifically pressure -- solved from the
+    incompressibility constraint, never transported (Stage 5 Completion
+    Criterion 2, `docs/planning/roadmap.md` TASK-032: "pressure is not
+    among the fields `step` advances, and handing `step` a `fields`
+    mapping that contains the pressure field raises a named error").
+
+    Carries no behaviour of its own beyond `ScalarField`'s -- purely a
+    type marker `simulation.step` can `isinstance`-check, the "stated at
+    the API level deliberately, not the configuration level" shape
+    Criterion 2 asks for: there is no configuration surface that names
+    which fields are transported, so the guard has to live where a real
+    `Field` object is actually in hand.
+
+    **Lives here, not in `engine/numerics/pressure_coupling.py` (this
+    class's one real producer, `PISO.correct`), because `simulation.py`
+    needs to import it and cannot import from `pressure_coupling.py`
+    without a circular import** -- that module already imports
+    `accumulate_flux_to_cells` from `simulation.py`. The same reasoning
+    `IncompatibleVelocityFieldError`'s own TASK-031a move to
+    `vector_field.py` used (`src/pyflow/engine/CLAUDE.md`'s
+    `vector_field.py` entry).
+
+    `PressureCoupling.correct`'s own abstract signature is unchanged
+    (`-> tuple[VectorField, ScalarField]`) -- no Stage 3 interface
+    change, since a `PressureField` instance already satisfies it
+    (covariant return narrowing, valid under `mypy --strict`). Only
+    `PISO`'s own concrete implementation constructs one.
+    """
diff --git a/src/pyflow/engine/simulation.py b/src/pyflow/engine/simulation.py
index 92d9aed..3dc2a3a 100644
--- a/src/pyflow/engine/simulation.py
+++ b/src/pyflow/engine/simulation.py
@@ -21,6 +21,7 @@
 
 from pyflow.engine.field import Field
 from pyflow.engine.mesh import Mesh
+from pyflow.engine.scalar_field import PressureField
 from pyflow.engine.vector_field import VectorField
 
 if TYPE_CHECKING:
@@ -45,6 +46,18 @@ class MismatchedMeshError(ValueError):
     """
 
 
+class PressureFieldTransportError(ValueError):
+    """Raised when `step`'s own `fields` mapping contains a
+    `PressureField` (TASK-032, Stage 5 Completion Criterion 2,
+    `docs/planning/roadmap.md`): pressure is solved from the
+    incompressibility constraint, not transported, so quietly advecting
+    it would silently discard exactly the property that makes it
+    meaningful. Stated at the API level -- a real `isinstance` check
+    against the field object itself, since there is no configuration
+    surface today that names which fields are transported.
+    """
+
+
 def accumulate_flux_to_cells(mesh: Mesh, face_values: torch.Tensor) -> torch.Tensor:
     """Reduce a face-valued array to a cell-valued one via the discrete
     Gauss theorem: `sum(value * area * outward_normal_sign) / volume`,
@@ -122,10 +135,16 @@ def step(
     decisions for the recorded version of this paragraph.
 
     Raises `MismatchedMeshError` if any field in `fields` is not defined
-    over the same mesh as `velocity`.
+    over the same mesh as `velocity`, and `PressureFieldTransportError`
+    (TASK-032) if `fields` contains a `PressureField` -- pressure is
+    solved from the incompressibility constraint, never transported.
     """
     mesh = velocity.mesh
     for name, field in fields.items():
+        if isinstance(field, PressureField):
+            raise PressureFieldTransportError(
+                f"field {name!r} is a PressureField -- pressure is solved, not transported"
+            )
         if field.mesh is not mesh:
             raise MismatchedMeshError(
                 f"field {name!r} is defined over a different mesh than the velocity field"
diff --git a/tests/features/pressure_field.feature b/tests/features/pressure_field.feature
new file mode 100644
index 0000000..b821f72
--- /dev/null
+++ b/tests/features/pressure_field.feature
@@ -0,0 +1,48 @@
+# The acceptance criteria for Pressure Field (TASK-032, Stage 5's third
+# 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_field.py` binds these
+# scenarios.
+#
+# Pressure is *not* transported -- it is solved for, from the
+# incompressibility constraint. A criterion that treats it as another
+# advected scalar has misunderstood the task
+# (`docs/handbook/numerical-methods/pressure-velocity-coupling.md`).
+# `PISO` (TASK-027, Stage 4) already performs the solve this task's own
+# criteria describe -- these scenarios prove the properties Stage 4's
+# own criteria never had cause to check (constant pressure for a
+# divergence-free input, the null-space remedy actually holding), not a
+# new pressure-solving mechanism.
+
+Feature: Pressure Field
+
+  Background:
+    Given a small, non-square, non-trivially-origined mesh
+
+  Scenario: A divergence-free provisional velocity field yields a pressure field constant to solver tolerance
+    Given a uniform provisional velocity field with zero-gradient boundaries on every wall
+    When the field is corrected by one PISO pass
+    Then the solved pressure field is constant to solver tolerance
+
+  Scenario: A provisional velocity field with known nonzero divergence yields a pressure field that is not constant
+    Given a provisional velocity field with real interior divergence, not aligned with either mesh axis
+    When the field is corrected by one PISO pass
+    Then the solved pressure field is not constant
+
+  Scenario: Adding a constant to the pressure field leaves the corrected velocity unchanged
+    Given a provisional velocity field with real interior divergence, not aligned with either mesh axis
+    When the field is corrected by one PISO pass
+    And a nonzero constant is added to the solved pressure field everywhere
+    Then the velocity correction computed from the shifted pressure field matches the original exactly
+
+  Scenario: A fields mapping containing a pressure field is rejected by step, not silently transported
+    Given a pressure field produced by a real PISO correction
+    When the simulation is stepped with that pressure field included among the transported fields
+    Then a named error says pressure is not transported
+
+  Scenario: A boundary configuration whose prescribed velocities violate the zero-net-flux compatibility condition fails to load
+    Given a configuration prescribing a nonzero net velocity flux across all four boundaries
+    When the configuration is loaded
+    Then loading is rejected before any pressure solve could be attempted
diff --git a/tests/unit/CLAUDE.md b/tests/unit/CLAUDE.md
index 7c4d243..24376e5 100644
--- a/tests/unit/CLAUDE.md
+++ b/tests/unit/CLAUDE.md
@@ -237,6 +237,29 @@ being trusted, the same discipline every other hand-derived scenario in
 this repository uses -- see this module's own commit message for the
 full per-face derivation.
 
+**`test_pressure_field.py` (TASK-032, added 2026-08-29) is the
+eleventh, and Stage 5's second module in this lineage** -- Stage 5's
+third task, binding `tests/features/pressure_field.feature`'s five
+scenarios against the real `PISO` (TASK-027, Stage 4) directly, not a
+new pressure-solving mechanism: this task's own job was proving
+properties Stage 4's own criteria never had cause to check. 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 own identically-named double uses),
+no golden-demo config file or CLI run. **Its divergence-free fixture is
+a uniform velocity field with zero-gradient boundaries, not a
+hand-crafted Dirichlet one** -- a uniform field's own divergence is
+exactly zero by the discrete Gauss theorem's closure identity (the sum
+of a closed cell's own face-area-weighted outward normals is the zero
+vector) regardless of the field's value, and zero-gradient boundaries
+extrapolate the interior's own uniform value back at every wall with no
+artificial jump, unlike a Dirichlet condition that would have to
+happen to match it exactly. The null-space scenario shifts the real
+solved pressure by a nonzero constant and recomputes the correction
+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.
+
 **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/test_pressure_field.py b/tests/unit/test_pressure_field.py
new file mode 100644
index 0000000..d5742ba
--- /dev/null
+++ b/tests/unit/test_pressure_field.py
@@ -0,0 +1,270 @@
+"""Binds `tests/features/pressure_field.feature` (TASK-032) -- Stage 5's
+third task in build order. `PISO` (TASK-027, Stage 4) already performs
+the solve this task's own criteria describe; these scenarios prove the
+properties Stage 4's own criteria never had cause to check (constant
+pressure for a divergence-free input, the null-space remedy actually
+holding, pressure rejected by `step`), not a new pressure-solving
+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
+from pathlib import Path
+from typing import Literal
+
+import pytest
+import torch
+from pytest_bdd import given, scenarios, then, when
+
+from pyflow.configuration import load_config
+from pyflow.engine import simulation
+from pyflow.engine.field import Field
+from pyflow.engine.mesh import StructuredCartesianMesh
+from pyflow.engine.numerics.boundary_condition import BoundaryCondition
+from pyflow.engine.numerics.gradient import GreenGaussGradient
+from pyflow.engine.numerics.linear_solver import ConjugateGradientSolver
+from pyflow.engine.numerics.pressure_coupling import PISO
+from pyflow.engine.scalar_field import PressureField
+from pyflow.engine.vector_field import VectorField
+
+from ._numerics import FixedGradientCondition, default_mesh
+
+scenarios("pressure_field.feature")
+
+
+class _ZeroNormalVelocity(BoundaryCondition):
+    """Dirichlet, fixed at zero -- a closed box: every boundary
+    prescribes zero normal velocity, the same fixture shape
+    `test_piso_pressure_coupling.py`'s own identically-named double
+    uses, for the divergent scenario's own provisional field.
+    """
+
+    @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
+
+
+def _divergent_provisional_velocity(mesh: StructuredCartesianMesh) -> VectorField:
+    # Neither axis-aligned nor uniform -- a real divergence field, not a
+    # degenerate one a wrong implementation could satisfy by luck
+    # (docs/practices.md, "distinct factors").
+    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 _closed_box_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]
+    provisional_velocity: VectorField
+    dt: float = 1.0
+    pressure: PressureField | None = None
+    corrected_velocity: VectorField | None = None
+    shifted_correction: torch.Tensor | None = None
+    error: Exception | None = None
+    config_error: ValueError | None = None
+
+
+def _piso(boundary_conditions: dict[str, BoundaryCondition]) -> PISO:
+    return PISO(ConjugateGradientSolver(tolerance=1e-10, max_iterations=500), boundary_conditions)
+
+
+# -- Given -------------------------------------------------------------
+
+
+@given("a small, non-square, non-trivially-origined mesh", target_fixture="ctx")
+def _given_default_mesh() -> _Context:
+    mesh = default_mesh(extent=(5, 4))
+    return _Context(
+        mesh=mesh,
+        boundary_conditions=_closed_box_boundary_conditions(),
+        provisional_velocity=VectorField(
+            mesh, "u_star", num_components=2, initial_value=(0.0, 0.0)
+        ),
+    )
+
+
+@given("a uniform provisional velocity field with zero-gradient boundaries on every wall")
+def _given_uniform_divergence_free_velocity(ctx: _Context) -> None:
+    # A uniform vector field's own divergence is exactly zero by the
+    # discrete Gauss theorem's own closure identity (sum of face
+    # area-weighted outward normals over any closed cell is the zero
+    # vector), regardless of the field's own value -- verified directly
+    # in this scenario's own `Then` step, not just assumed.
+    ctx.provisional_velocity = VectorField(
+        ctx.mesh, "u_star", num_components=2, initial_value=(1.3, -0.7)
+    )
+    condition = FixedGradientCondition(0.0)
+    ctx.boundary_conditions = {
+        "north": condition,
+        "south": condition,
+        "east": condition,
+        "west": condition,
+    }
+
+
+@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)
+    ctx.boundary_conditions = _closed_box_boundary_conditions()
+
+
+@given("a pressure field produced by a real PISO correction", target_fixture="ctx")
+def _given_real_pressure_field() -> _Context:
+    mesh = default_mesh(extent=(5, 4))
+    boundary_conditions = _closed_box_boundary_conditions()
+    provisional = _divergent_provisional_velocity(mesh)
+    _corrected, pressure = _piso(boundary_conditions).correct(provisional, dt=1.0)
+    # `PressureCoupling.correct`'s own abstract signature stays
+    # `ScalarField` (no Stage 3 interface change) -- this assertion is
+    # what actually proves `PISO`'s concrete implementation constructs
+    # the narrower `PressureField`, which is what makes the "rejected by
+    # step" scenario below meaningful at all.
+    assert isinstance(pressure, PressureField)
+    ctx = _Context(
+        mesh=mesh, boundary_conditions=boundary_conditions, provisional_velocity=provisional
+    )
+    ctx.pressure = pressure
+    return ctx
+
+
+@given("a configuration prescribing a nonzero net velocity flux across all four boundaries")
+def _given_incompatible_boundary_config(ctx: _Context) -> None:
+    # nx=4, ny=2, dx=dy=1: north/south length 4, east/west length 2.
+    # 1*4 + 0*4 + (-2)*2 + (-1)*2 = 4 - 4 - 2 = -2, nonzero -- the same
+    # fixture shape `test_configuration.py`'s own net-flux rejection test
+    # uses, since this is exactly the existing check TASK-019 already
+    # built (`_validate_boundary_conditions_jointly`), not reimplemented
+    # here.
+    del ctx
+
+
+# -- When ------------------------------------------------------------------
+
+
+@when("the field is corrected by one PISO pass")
+def _when_corrected(ctx: _Context) -> None:
+    piso = PISO(
+        ConjugateGradientSolver(tolerance=1e-10, max_iterations=500), ctx.boundary_conditions
+    )
+    corrected, pressure = piso.correct(ctx.provisional_velocity, dt=ctx.dt)
+    assert isinstance(pressure, PressureField)
+    ctx.corrected_velocity = corrected
+    ctx.pressure = pressure
+
+
+@when("a nonzero constant is added to the solved pressure field everywhere")
+def _when_constant_added(ctx: _Context) -> None:
+    assert ctx.pressure is not None
+    shifted = PressureField(ctx.mesh, "pressure_shifted")
+    shifted.values[:] = ctx.pressure.values + 7.5
+    pressure_boundary_conditions = {
+        name: FixedGradientCondition(0.0) for name in ("north", "south", "east", "west")
+    }
+    gradient_scheme = GreenGaussGradient(pressure_boundary_conditions)
+    ctx.shifted_correction = ctx.provisional_velocity.values - ctx.dt * gradient_scheme.gradient(
+        shifted
+    )
+
+
+@when("the simulation is stepped with that pressure field included among the transported fields")
+def _when_stepped_with_pressure(ctx: _Context) -> None:
+    assert ctx.pressure is not None
+    velocity = VectorField(ctx.mesh, "velocity", num_components=2, initial_value=(0.0, 0.0))
+    fields: dict[str, Field] = {"pressure": ctx.pressure}
+    try:
+        simulation.step(fields, velocity, _StubNumerics(), 0.1)  # type: ignore[arg-type]
+    except simulation.PressureFieldTransportError as exc:
+        ctx.error = exc
+
+
+class _StubNumerics:
+    """Never reached -- `step`'s own pressure guard raises before
+    `numerics` is used for anything, so this only needs to exist to
+    satisfy the parameter.
+    """
+
+
+@when("the configuration is loaded")
+def _when_config_loaded(ctx: _Context, tmp_path: Path) -> None:
+    config_file = tmp_path / "config.yaml"
+    config_file.write_text(
+        "mesh:\n"
+        "  extent: [4, 2]\n"
+        "  spacing: [1.0, 1.0]\n"
+        "numerics:\n"
+        "  boundary_conditions:\n"
+        "    north:\n      velocity: 1.0\n"
+        "    south:\n      velocity: 0.0\n"
+        "    east:\n      velocity: -2.0\n"
+        "    west:\n      velocity: -1.0\n"
+    )
+    try:
+        load_config(config_file)
+    except ValueError as exc:
+        ctx.config_error = exc
+
+
+# -- Then ------------------------------------------------------------------
+
+
+@then("the solved pressure field is constant to solver tolerance")
+def _then_pressure_constant(ctx: _Context) -> None:
+    assert ctx.pressure is not None
+    spread = float(ctx.pressure.values.max() - ctx.pressure.values.min())
+    assert spread == pytest.approx(0.0, abs=1e-6), f"pressure spread {spread} is not ~0"
+
+
+@then("the solved pressure field is not constant")
+def _then_pressure_not_constant(ctx: _Context) -> None:
+    assert ctx.pressure is not None
+    spread = float(ctx.pressure.values.max() - ctx.pressure.values.min())
+    assert spread > 1e-3, f"pressure spread {spread} is too close to constant"
+
+
+@then(
+    "the velocity correction computed from the shifted pressure field matches the original exactly"
+)
+def _then_shifted_correction_matches(ctx: _Context) -> None:
+    assert ctx.corrected_velocity is not None
+    assert ctx.shifted_correction is not None
+    torch.testing.assert_close(
+        ctx.shifted_correction, ctx.corrected_velocity.values, rtol=1e-9, atol=1e-9
+    )
+
+
+@then("a named error says pressure is not transported")
+def _then_pressure_transport_error(ctx: _Context) -> None:
+    assert isinstance(ctx.error, simulation.PressureFieldTransportError)
+
+
+@then("loading is rejected before any pressure solve could be attempted")
+def _then_config_rejected(ctx: _Context) -> None:
+    assert ctx.config_error is not None
+    assert "net flux" in str(ctx.config_error)
diff --git a/tools/validators/check_references.py b/tools/validators/check_references.py
index 635b1c7..76f85bb 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_field.feature": "TASK-032",
     "tests/features/pressure_correction_loop.feature": "TASK-033",
     "tests/features/navier_stokes_timestep.feature": "TASK-034",
 }