rtx-cfd: PISO validated by manufactured solution — after fixing the inverted projection

PisoSolver was the only major solver in the workspace with no verification
of any kind. Writing the MMS harness for it (tests/mms_piso.rs) and
inspecting the implementation found the census's defect species again:

- The pressure correction had its SIGN inverted: it solved
  -lap(p') = +rho div(u*)/dt and then corrected with u = u* - (dt/rho)
  grad(p'), so each projection DOUBLED the divergence instead of removing
  it.
- The momentum sweeps froze the near-wall lines (1..ny-1) and the pressure
  correction skipped the outer ring of cells (1..nx-1) — both exactly the
  defects repaired in SIMPLE.
- The "explicit" predictor read neighbours the same sweep had already
  overwritten, so the step depended on sweep order.
- The pressure gradient was dropped entirely on the last interior face.

Rewritten as a genuinely explicit predictor plus anchored-Neumann
projection on the staggered grid, with the conventions SIMPLE now embodies:
near-wall lines are unknowns with half-cell wall diffusion, continuity on
every cell, boundary faces are prescribed data. Momentum-source and
wall-velocity hooks added so the manufactured solution can reach it.

Measured (16 -> 32 -> 64): L2 velocity 3.516214e-2, 1.953750e-2,
1.037512e-2 — orders 0.85 and 0.91, first-order upwind's rate — with
max |div u| ~ 1e-9 in every cell. The errors agree with SIMPLE's on the
same meshes to six or seven significant figures: an implicit under-relaxed
outer iteration and an explicit time-marching projection land on the same
discrete steady solution, which is what sharing a spatial discretisation
must produce and is very hard for two independently wrong solvers to fake.

285 tests, 0 failing.

Co-Authored-By: Claude Fable 5 <[email protected]>

crates/specialized/rtx-cfd/src/solvers/incompressible/piso.rs

@@ -1,8 +1,41 @@
 //! PISO (Pressure-Implicit with Splitting of Operators) algorithm
 //!
-//! The PISO algorithm is a non-iterative pressure-velocity coupling algorithm
+//! A transient pressure-velocity coupling method: one explicit momentum
-//! particularly well-suited for transient flow problems. It consists of one
+//! predictor per time step, followed by pressure-correction (projection)
-//! predictor step followed by two or more corrector steps.
+//! steps that make the velocity field divergence-free. Marching it in time
+//! with a steady forcing converges to the steady discrete solution, which is
+//! how `tests/mms_piso.rs` verifies it against a manufactured solution.
+//!
+//! # Grid convention
+//!
+//! The staggered layout is the one `FlowField` and the SIMPLE solver define:
+//! `u[(j, i)]` lives at `(i dx, (j + 0.5) dy)` for `i = 0..=nx`, `v[(j, i)]`
+//! at `((i + 0.5) dx, j dy)` for `j = 0..=ny`, `p[(j, i)]` at cell centres.
+//! The only velocity components on a domain boundary are the normal ones —
+//! u faces `i = 0`, `i = nx` and v faces `j = 0`, `j = ny`. Everything else,
+//! including the near-wall lines, is an unknown and is updated every step.
+//!
+//! # History
+//!
+//! The previous implementation had never had a test of any kind, and
+//! inspection plus the manufactured-solution harness found the same defect
+//! species the SIMPLE census recorded:
+//!
+//! - **The pressure correction had its sign inverted.** It solved
+//!   `-lap(p') = +rho div(u*) / dt` and then corrected with
+//!   `u = u* - (dt/rho) grad(p')`, so each projection *doubled* the
+//!   divergence instead of removing it.
+//! - The momentum sweeps froze the near-wall lines (`1..ny-1`), imposing the
+//!   wall half a cell inside the domain, and the pressure correction skipped
+//!   the outer ring of cells (`1..nx-1`), so ring cells had no continuity
+//!   equation — both exactly as in SIMPLE before its repair.
+//! - The predictor read neighbours that the same sweep had already
+//!   overwritten, so the "explicit" step mixed old and new values in sweep
+//!   order.
+//! - The pressure gradient was dropped entirely on the last interior face
+//!   (`if i < nx - 1 { ... } else { 0.0 }`).
+//! - Convective face fluxes fell back to the centre value at the sweep edges
+//!   instead of using the prescribed boundary faces that exist there.
 
 use super::{BoundaryConditions, FlowField, IncompressibleSolver, SolverResult};
 use crate::{CfdConfig, CfdResult};
@@ -11,11 +44,15 @@ use async_trait::async_trait;
 /// Parameters for PISO algorithm
 #[derive(Debug, Clone)]
 pub struct PisoParameters {
-    /// Number of corrector steps (typically 2-3)
+    /// Number of corrector steps (typically 2-3). With an explicit predictor
+    /// the first projection already removes the divergence; the extra
+    /// correctors are cheap no-ops kept for the algorithm's shape.
     pub corrector_steps: usize,
-    /// Time step size
+    /// Time step size. The predictor is explicit, so stability requires
+    /// `dt < dx^2 / (4 nu)` and `dt < dx / |u|_max`.
     pub time_step: f64,
-    /// Convergence tolerance
+    /// Convergence tolerance on the normalised mass imbalance after
+    /// correction.
     pub tolerance: f64,
 }
 
@@ -42,6 +79,18 @@ pub struct PisoResult {
 pub struct PisoSolver {
     config: CfdConfig,
     parameters: PisoParameters,
+    /// Optional volumetric momentum source `f(x, y) -> (f_x, f_y)`, per unit
+    /// volume — the hook a manufactured solution enters through, exactly as
+    /// on [`super::SimpleSolver`].
+    #[allow(clippy::type_complexity)]
+    momentum_source: Option<Box<dyn Fn(f64, f64) -> (f64, f64) + Send + Sync>>,
+    /// Optional wall velocity `f(x, y) -> (u_wall, v_wall)`, sampled at the
+    /// wall face position. The near-wall control volumes need the tangential
+    /// wall velocity for their half-cell diffusion term, and on this
+    /// staggered layout there is nowhere to store it. Falls back to the value
+    /// on the near-wall line itself when unset.
+    #[allow(clippy::type_complexity)]
+    wall_velocity: Option<Box<dyn Fn(f64, f64) -> (f64, f64) + Send + Sync>>,
 }
 
 impl PisoSolver {
@@ -49,442 +98,336 @@ impl PisoSolver {
     pub fn new(config: CfdConfig, parameters: PisoParameters) -> CfdResult<Self> {
         config.validate()?;
 
-        Ok(Self { config, parameters })
+        Ok(Self {
+            config,
+            parameters,
+            momentum_source: None,
+            wall_velocity: None,
+        })
     }
 
-    /// Solve momentum predictor step
+    /// Set a volumetric momentum source. See [`Self::momentum_source`].
-    /// Discretize: ∂u/∂t + ∇·(u⊗u) = -∇p^n/ρ + ν∇²u
+    pub fn set_momentum_source<F>(&mut self, source: F)
-    fn solve_momentum_predictor(
+    where
-        &self,
+        F: Fn(f64, f64) -> (f64, f64) + Send + Sync + 'static,
-        flow_field: &mut FlowField,
+    {
-        dt: f64,
+        self.momentum_source = Some(Box::new(source));
-        rho: f64,
-        nu: f64,
-    ) -> CfdResult<()> {
-        let (_nx, _ny, dx, dy) = flow_field.grid_info();
-
-        // Solve u-momentum equation
-        self.solve_u_momentum(flow_field, dt, rho, nu, dx, dy)?;
-
-        // Solve v-momentum equation
-        self.solve_v_momentum(flow_field, dt, rho, nu, dx, dy)?;
-
-        // Store predicted velocities
-        flow_field.copy_to_starred();
-
-        Ok(())
     }
 
-    /// Solve u-momentum equation using finite volume method
+    /// Set the wall velocity as a function of position. See
-    fn solve_u_momentum(
+    /// [`Self::wall_velocity`].
-        &self,
+    pub fn set_wall_velocity<F>(&mut self, f: F)
-        flow_field: &mut FlowField,
+    where
-        dt: f64,
+        F: Fn(f64, f64) -> (f64, f64) + Send + Sync + 'static,
-        rho: f64,
+    {
-        nu: f64,
+        self.wall_velocity = Some(Box::new(f));
-        dx: f64,
+    }
-        dy: f64,
-    ) -> CfdResult<()> {
-        let (nx, ny, _, _) = flow_field.grid_info();
 
-        // For each u-velocity control volume (i+1/2, j)
+    fn u_wall(&self, flow_field: &FlowField, i: usize, j: usize, y_wall: f64, dx: f64) -> f64 {
-        for j in 1..(ny - 1) {
+        self.wall_velocity
+            .as_ref()
+            .map_or(flow_field.u_old[(j, i)], |f| f(i as f64 * dx, y_wall).0)
+    }
+
+    fn v_wall(&self, flow_field: &FlowField, i: usize, j: usize, x_wall: f64, dy: f64) -> f64 {
+        self.wall_velocity
+            .as_ref()
+            .map_or(flow_field.v_old[(j, i)], |f| f(x_wall, j as f64 * dy).1)
+    }
+
+    /// Upwind face value: the value carried across the face is the one from
+    /// the side the flow comes from.
+    fn upwind(face_velocity: f64, upstream: f64, downstream: f64) -> f64 {
+        if face_velocity >= 0.0 {
+            upstream
+        } else {
+            downstream
+        }
+    }
+
+    /// Explicit momentum predictor:
+    /// `u* = u_old + dt (-conv + nu lap(u) - grad(p)/rho + f/rho)`,
+    /// every term evaluated from `u_old`/`v_old`, so the step is genuinely
+    /// explicit and independent of sweep order.
+    fn momentum_predictor(&self, flow_field: &mut FlowField, dt: f64) -> CfdResult<()> {
+        let (nx, ny, dx, dy) = flow_field.grid_info();
+        let rho = self.config.density;
+        let nu = self.config.viscosity / rho;
+
+        // u faces: every row is an unknown; only i = 0 and i = nx are
+        // boundary data.
+        for j in 0..ny {
             for i in 1..nx {
-                // Time derivative term: ∂u/∂t ≈ (u_new - u_old)/dt
+                let uo = &flow_field.u_old;
-                let time_coeff = 1.0 / dt;
+                let vo = &flow_field.v_old;
-                let time_source = flow_field.u_old[(j, i)] / dt;
+                let u_p = uo[(j, i)];
 
-                // Convective terms: ∇·(u⊗u)
+                // Cell-centre velocities on the east/west faces of the u
-                // Face velocities for convection (interpolated)
+                // control volume. The neighbours i-1 and i+1 always exist:
-                let u_east = if i < nx - 1 {
+                // they are boundary faces at the sweep edges, which hold
-                    0.5 * (flow_field.u[(j, i)] + flow_field.u[(j, i + 1)])
+                // prescribed data rather than needing a fallback.
-                } else {
+                let ue_face = 0.5 * (uo[(j, i)] + uo[(j, i + 1)]);
-                    flow_field.u[(j, i)]
+                let uw_face = 0.5 * (uo[(j, i - 1)] + uo[(j, i)]);
-                };
-                let u_west = if i > 1 {
-                    0.5 * (flow_field.u[(j, i - 1)] + flow_field.u[(j, i)])
-                } else {
-                    flow_field.u[(j, i)]
-                };
-                let _u_north = if j < ny - 1 {
-                    0.5 * (flow_field.u[(j, i)] + flow_field.u[(j + 1, i)])
-                } else {
-                    flow_field.u[(j, i)]
-                };
-                let _u_south = if j > 1 {
-                    0.5 * (flow_field.u[(j - 1, i)] + flow_field.u[(j, i)])
-                } else {
-                    flow_field.u[(j, i)]
-                };
 
-                // Transverse velocities
+                let south_is_wall = j == 0;
-                let v_north = if i > 0 && i < nx && j < ny {
+                let north_is_wall = j + 1 == ny;
-                    0.5 * (flow_field.v[(j + 1, i - 1)] + flow_field.v[(j + 1, i)])
+
-                } else {
+                // Transverse face velocities; a solid wall passes no mass.
+                let vn_face = if north_is_wall {
                     0.0
-                };
-                let v_south = if i > 0 && i < nx && j > 0 {
-                    0.5 * (flow_field.v[(j, i - 1)] + flow_field.v[(j, i)])
                 } else {
+                    0.5 * (vo[(j + 1, i - 1)] + vo[(j + 1, i)])
+                };
+                let vs_face = if south_is_wall {
                     0.0
+                } else {
+                    0.5 * (vo[(j, i - 1)] + vo[(j, i)])
                 };
 
-                // Convective fluxes (upwind scheme)
+                let conv_x = (ue_face * Self::upwind(ue_face, uo[(j, i)], uo[(j, i + 1)])
-                let conv_east = u_east
+                    - uw_face * Self::upwind(uw_face, uo[(j, i - 1)], uo[(j, i)]))
-                    * if u_east > 0.0 {
+                    / dx;
-                        flow_field.u[(j, i)]
+                let conv_y = (vn_face
-                    } else if i < nx - 1 {
+                    * if north_is_wall {
-                        flow_field.u[(j, i + 1)]
-                    } else {
-                        flow_field.u[(j, i)]
-                    };
-                let conv_west = u_west
-                    * if u_west > 0.0 {
-                        if i > 1 {
-                            flow_field.u[(j, i - 1)]
-                        } else {
-                            flow_field.u[(j, i)]
-                        }
-                    } else {
-                        flow_field.u[(j, i)]
-                    };
-                let conv_north = v_north
-                    * if v_north > 0.0 {
-                        flow_field.u[(j, i)]
-                    } else if j < ny - 1 {
-                        flow_field.u[(j + 1, i)]
-                    } else {
-                        flow_field.u[(j, i)]
-                    };
-                let conv_south = v_south
-                    * if v_south > 0.0 {
-                        if j > 1 {
-                            flow_field.u[(j - 1, i)]
-                        } else {
-                            flow_field.u[(j, i)]
-                        }
-                    } else {
-                        flow_field.u[(j, i)]
-                    };
-
-                let convection = (conv_east - conv_west) / dx + (conv_north - conv_south) / dy;
-
-                // Diffusive terms: ν∇²u
-                let u_center = flow_field.u[(j, i)];
-                let u_east_diff = if i < nx - 1 {
-                    flow_field.u[(j, i + 1)]
-                } else {
-                    u_center
-                };
-                let u_west_diff = if i > 1 {
-                    flow_field.u[(j, i - 1)]
-                } else {
-                    u_center
-                };
-                let u_north_diff = if j < ny - 1 {
-                    flow_field.u[(j + 1, i)]
-                } else {
-                    u_center
-                };
-                let u_south_diff = if j > 1 {
-                    flow_field.u[(j - 1, i)]
-                } else {
-                    u_center
-                };
-
-                let diffusion_x = (u_east_diff - 2.0 * u_center + u_west_diff) / (dx * dx);
-                let diffusion_y = (u_north_diff - 2.0 * u_center + u_south_diff) / (dy * dy);
-                let diffusion = nu * (diffusion_x + diffusion_y);
-
-                // Pressure gradient: -∂p/∂x / ρ (using pressure from previous time step)
-                let pressure_grad = if i < nx - 1 {
-                    -(flow_field.p[(j, i)] - flow_field.p[(j, i - 1)]) / (rho * dx)
-                } else {
                         0.0
+                    } else {
+                        Self::upwind(vn_face, uo[(j, i)], uo[(j + 1, i)])
+                    }
+                    - vs_face
+                        * if south_is_wall {
+                            0.0
+                        } else {
+                            Self::upwind(vs_face, uo[(j - 1, i)], uo[(j, i)])
+                        })
+                    / dy;
+
+                let diff_x = nu * (uo[(j, i + 1)] - 2.0 * u_p + uo[(j, i - 1)]) / (dx * dx);
+
+                // Wall-adjacent diffusive fluxes act over half a cell: the
+                // node beyond the wall face is the wall itself, dy/2 away.
+                let flux_north = if north_is_wall {
+                    nu * (self.u_wall(flow_field, i, j, ny as f64 * dy, dx) - u_p) / (0.5 * dy)
+                } else {
+                    nu * (uo[(j + 1, i)] - u_p) / dy
                 };
+                let flux_south = if south_is_wall {
+                    nu * (u_p - self.u_wall(flow_field, i, j, 0.0, dx)) / (0.5 * dy)
+                } else {
+                    nu * (u_p - uo[(j - 1, i)]) / dy
+                };
+                let diff_y = (flux_north - flux_south) / dy;
 
-                // Source term
+                // The pressure gradient acts on every unknown face — dropping
-                let source = time_source + diffusion + pressure_grad;
+                // it anywhere solves a different equation there.
+                let pressure_gradient =
+                    -(flow_field.p[(j, i)] - flow_field.p[(j, i - 1)]) / (rho * dx);
 
-                // Solve: (1/dt + convection_coeff) * u_new = source
+                let body_force = self
-                let total_coeff = time_coeff;
+                    .momentum_source
-                flow_field.u[(j, i)] = (source - convection) / total_coeff;
+                    .as_ref()
+                    .map_or(0.0, |f| f(i as f64 * dx, (j as f64 + 0.5) * dy).0 / rho);
+
+                flow_field.u[(j, i)] = u_p
+                    + dt * (-conv_x - conv_y + diff_x + diff_y + pressure_gradient + body_force);
             }
         }
 
-        Ok(())
+        // v faces, mirrored.
-    }
-
-    /// Solve v-momentum equation using finite volume method
-    fn solve_v_momentum(
-        &self,
-        flow_field: &mut FlowField,
-        dt: f64,
-        rho: f64,
-        nu: f64,
-        dx: f64,
-        dy: f64,
-    ) -> CfdResult<()> {
-        let (nx, ny, _, _) = flow_field.grid_info();
-
-        // For each v-velocity control volume (i, j+1/2)
         for j in 1..ny {
-            for i in 1..(nx - 1) {
+            for i in 0..nx {
-                // Time derivative term
+                let uo = &flow_field.u_old;
-                let time_coeff = 1.0 / dt;
+                let vo = &flow_field.v_old;
-                let time_source = flow_field.v_old[(j, i)] / dt;
+                let v_p = vo[(j, i)];
 
-                // Convective terms
+                let vn_face = 0.5 * (vo[(j, i)] + vo[(j + 1, i)]);
-                let _v_east = if i < nx - 1 {
+                let vs_face = 0.5 * (vo[(j - 1, i)] + vo[(j, i)]);
-                    0.5 * (flow_field.v[(j, i)] + flow_field.v[(j, i + 1)])
-                } else {
-                    flow_field.v[(j, i)]
-                };
-                let _v_west = if i > 1 {
-                    0.5 * (flow_field.v[(j, i - 1)] + flow_field.v[(j, i)])
-                } else {
-                    flow_field.v[(j, i)]
-                };
-                let v_north = if j < ny - 1 {
-                    0.5 * (flow_field.v[(j, i)] + flow_field.v[(j + 1, i)])
-                } else {
-                    flow_field.v[(j, i)]
-                };
-                let v_south = if j > 1 {
-                    0.5 * (flow_field.v[(j - 1, i)] + flow_field.v[(j, i)])
-                } else {
-                    flow_field.v[(j, i)]
-                };
 
-                // Transverse velocities
+                let west_is_wall = i == 0;
-                let u_east = if j > 0 && j < ny && i < nx - 1 {
+                let east_is_wall = i + 1 == nx;
-                    0.5 * (flow_field.u[(j - 1, i + 1)] + flow_field.u[(j, i + 1)])
+
-                } else {
+                let ue_face = if east_is_wall {
                     0.0
-                };
-                let u_west = if j > 0 && j < ny && i > 0 {
-                    0.5 * (flow_field.u[(j - 1, i)] + flow_field.u[(j, i)])
                 } else {
+                    0.5 * (uo[(j - 1, i + 1)] + uo[(j, i + 1)])
+                };
+                let uw_face = if west_is_wall {
                     0.0
+                } else {
+                    0.5 * (uo[(j - 1, i)] + uo[(j, i)])
                 };
 
-                // Convective fluxes (upwind)
+                let conv_y = (vn_face * Self::upwind(vn_face, vo[(j, i)], vo[(j + 1, i)])
-                let conv_east = u_east
+                    - vs_face * Self::upwind(vs_face, vo[(j - 1, i)], vo[(j, i)]))
-                    * if u_east > 0.0 {
+                    / dy;
-                        flow_field.v[(j, i)]
+                let conv_x = (ue_face
-                    } else if i < nx - 1 {
+                    * if east_is_wall {
-                        flow_field.v[(j, i + 1)]
-                    } else {
-                        flow_field.v[(j, i)]
-                    };
-                let conv_west = u_west
-                    * if u_west > 0.0 {
-                        if i > 1 {
-                            flow_field.v[(j, i - 1)]
-                        } else {
-                            flow_field.v[(j, i)]
-                        }
-                    } else {
-                        flow_field.v[(j, i)]
-                    };
-                let conv_north = v_north
-                    * if v_north > 0.0 {
-                        flow_field.v[(j, i)]
-                    } else if j < ny - 1 {
-                        flow_field.v[(j + 1, i)]
-                    } else {
-                        flow_field.v[(j, i)]
-                    };
-                let conv_south = v_south
-                    * if v_south > 0.0 {
-                        if j > 1 {
-                            flow_field.v[(j - 1, i)]
-                        } else {
-                            flow_field.v[(j, i)]
-                        }
-                    } else {
-                        flow_field.v[(j, i)]
-                    };
-
-                let convection = (conv_east - conv_west) / dx + (conv_north - conv_south) / dy;
-
-                // Diffusive terms
-                let v_center = flow_field.v[(j, i)];
-                let v_east_diff = if i < nx - 1 {
-                    flow_field.v[(j, i + 1)]
-                } else {
-                    v_center
-                };
-                let v_west_diff = if i > 1 {
-                    flow_field.v[(j, i - 1)]
-                } else {
-                    v_center
-                };
-                let v_north_diff = if j < ny - 1 {
-                    flow_field.v[(j + 1, i)]
-                } else {
-                    v_center
-                };
-                let v_south_diff = if j > 1 {
-                    flow_field.v[(j - 1, i)]
-                } else {
-                    v_center
-                };
-
-                let diffusion_x = (v_east_diff - 2.0 * v_center + v_west_diff) / (dx * dx);
-                let diffusion_y = (v_north_diff - 2.0 * v_center + v_south_diff) / (dy * dy);
-                let diffusion = nu * (diffusion_x + diffusion_y);
-
-                // Pressure gradient: -∂p/∂y / ρ
-                let pressure_grad = if j < ny - 1 {
-                    -(flow_field.p[(j, i)] - flow_field.p[(j - 1, i)]) / (rho * dy)
-                } else {
                         0.0
+                    } else {
+                        Self::upwind(ue_face, vo[(j, i)], vo[(j, i + 1)])
+                    }
+                    - uw_face
+                        * if west_is_wall {
+                            0.0
+                        } else {
+                            Self::upwind(uw_face, vo[(j, i - 1)], vo[(j, i)])
+                        })
+                    / dx;
+
+                let diff_y = nu * (vo[(j + 1, i)] - 2.0 * v_p + vo[(j - 1, i)]) / (dy * dy);
+
+                let flux_east = if east_is_wall {
+                    nu * (self.v_wall(flow_field, i, j, nx as f64 * dx, dy) - v_p) / (0.5 * dx)
+                } else {
+                    nu * (vo[(j, i + 1)] - v_p) / dx
                 };
+                let flux_west = if west_is_wall {
+                    nu * (v_p - self.v_wall(flow_field, i, j, 0.0, dy)) / (0.5 * dx)
+                } else {
+                    nu * (v_p - vo[(j, i - 1)]) / dx
+                };
+                let diff_x = (flux_east - flux_west) / dx;
 
-                let source = time_source + diffusion + pressure_grad;
+                let pressure_gradient =
+                    -(flow_field.p[(j, i)] - flow_field.p[(j - 1, i)]) / (rho * dy);
 
-                let total_coeff = time_coeff;
+                let body_force = self
-                flow_field.v[(j, i)] = (source - convection) / total_coeff;
+                    .momentum_source
+                    .as_ref()
+                    .map_or(0.0, |f| f((i as f64 + 0.5) * dx, j as f64 * dy).1 / rho);
+
+                flow_field.v[(j, i)] = v_p
+                    + dt * (-conv_x - conv_y + diff_x + diff_y + pressure_gradient + body_force);
             }
         }
 
+        flow_field.copy_to_starred();
         Ok(())
     }
 
-    /// Solve pressure correction equation
+    /// One projection: solve the pressure-correction Poisson equation and
-    /// ∇²p' = ρ∇·u*/dt
+    /// subtract `(dt/rho) grad(p')` from the predicted velocities, so the
-    fn solve_pressure_correction(
+    /// corrected field is discretely divergence-free.
-        &self,
+    ///
-        flow_field: &mut FlowField,
+    /// Continuity is enforced on every cell. A coefficient is zero exactly
-        dt: f64,
+    /// when its face is a domain boundary, where the normal velocity is
-        rho: f64,
+    /// prescribed and not correctable. With velocity prescribed on the whole
-        dx: f64,
+    /// boundary the system is pure Neumann; one cell is anchored to fix the
-        dy: f64,
+    /// level, which is legitimate because the source telescopes to the net
-    ) -> CfdResult<f64> {
+    /// boundary flux — zero for a closed box — so exactly one equation is
-        let (nx, ny, _, _) = flow_field.grid_info();
+    /// redundant.
+    ///
+    /// Returns the normalised mass imbalance of the *corrected* field — what
+    /// the projection failed to remove, which is the inner solver's
+    /// truncation and is the step's honest convergence measure.
+    fn project(&self, flow_field: &mut FlowField, dt: f64) -> CfdResult<f64> {
+        let (nx, ny, dx, dy) = flow_field.grid_info();
+        let rho = self.config.density;
 
-        // Reset pressure correction
         flow_field.p_prime.fill(0.0);
 
-        // Iterative solution using Gauss-Seidel
+        // Mass imbalance of the predicted field, per cell, as a flux.
-        let mut max_residual = 0.0;
+        for j in 0..ny {
+            for i in 0..nx {
+                let divergence_flux = rho
+                    * ((flow_field.u_star[(j, i + 1)] - flow_field.u_star[(j, i)]) * dy
+                        + (flow_field.v_star[(j + 1, i)] - flow_field.v_star[(j, i)]) * dx);
+                flow_field.sp[(j, i)] = -divergence_flux;
+            }
+        }
 
-        for _iter in 0..100 {
+        // With the correction `u = u* - (dt/rho) (p'_P - p'_W)/dx`, continuity
-            // Inner iterations for pressure correction
+        // of the corrected field gives neighbour coefficients
-            let mut residual: f64 = 0.0;
+        // `rho (dt/rho) A / delta = dt A / delta`.
+        let ae_interior = dt * dy / dx;
+        let an_interior = dt * dx / dy;
 
-            for j in 1..(ny - 1) {
+        for _sweep in 0..400 {
-                for i in 1..(nx - 1) {
+            let mut residual = 0.0;
-                    // Compute mass imbalance (divergence of velocity)
+            for j in 0..ny {
-                    let mass_imbalance =
+                for i in 0..nx {
-                        ((flow_field.u_star[(j, i + 1)] - flow_field.u_star[(j, i)]) / dx
+                    if i == 1 && j == 1 {
-                            + (flow_field.v_star[(j + 1, i)] - flow_field.v_star[(j, i)]) / dy)
+                        flow_field.p_prime[(j, i)] = 0.0;
-                            * rho
+                        continue;
-                            / dt;
+                    }
 
-                    // Coefficients for pressure correction equation
+                    let ae = if i + 1 == nx { 0.0 } else { ae_interior };
-                    let ae = 1.0 / (dx * dx);
+                    let aw = if i == 0 { 0.0 } else { ae_interior };
-                    let aw = 1.0 / (dx * dx);
+                    let an = if j + 1 == ny { 0.0 } else { an_interior };
-                    let an = 1.0 / (dy * dy);
+                    let as_ = if j == 0 { 0.0 } else { an_interior };
-                    let as_ = 1.0 / (dy * dy);
                     let ap = ae + aw + an + as_;
 
-                    // Neighboring pressure corrections
+                    let east = if i + 1 < nx {
-                    let p_east = if i < nx - 2 {
+                        ae * flow_field.p_prime[(j, i + 1)]
-                        flow_field.p_prime[(j, i + 1)]
                     } else {
                         0.0
                     };
-                    let p_west = if i > 1 {
+                    let west = if i > 0 {
-                        flow_field.p_prime[(j, i - 1)]
+                        aw * flow_field.p_prime[(j, i - 1)]
                     } else {
                         0.0
                     };
-                    let p_north = if j < ny - 2 {
+                    let north = if j + 1 < ny {
-                        flow_field.p_prime[(j + 1, i)]
+                        an * flow_field.p_prime[(j + 1, i)]
                     } else {
                         0.0
                     };
-                    let p_south = if j > 1 {
+                    let south = if j > 0 {
-                        flow_field.p_prime[(j - 1, i)]
+                        as_ * flow_field.p_prime[(j - 1, i)]
                     } else {
                         0.0
                     };
 
-                    // Gauss-Seidel update
+                    let p_new = (flow_field.sp[(j, i)] + east + west + north + south) / ap;
-                    let p_prime_new =
+                    let correction = p_new - flow_field.p_prime[(j, i)];
-                        (ae * p_east + aw * p_west + an * p_north + as_ * p_south + mass_imbalance)
+                    residual += correction * correction;
-                            / ap;
+                    flow_field.p_prime[(j, i)] = p_new;
-
-                    let correction_residual = (p_prime_new - flow_field.p_prime[(j, i)]).abs();
-                    residual = residual.max(correction_residual);
-
-                    flow_field.p_prime[(j, i)] = p_prime_new;
                 }
             }
-
+            if residual.sqrt() < 1e-12 {
-            max_residual = residual;
-
-            // Check inner convergence
-            if residual < 1e-8 {
                 break;
             }
         }
 
-        // Update pressure: p^(n+1) = p^n + p'
+        // Correct exactly the faces the equation above treated as
+        // correctable: every interior face.
+        for j in 0..ny {
+            for i in 1..nx {
+                let dp_dx = (flow_field.p_prime[(j, i)] - flow_field.p_prime[(j, i - 1)]) / dx;
+                flow_field.u[(j, i)] = flow_field.u_star[(j, i)] - (dt / rho) * dp_dx;
+            }
+        }
+        for j in 1..ny {
+            for i in 0..nx {
+                let dp_dy = (flow_field.p_prime[(j, i)] - flow_field.p_prime[(j - 1, i)]) / dy;
+                flow_field.v[(j, i)] = flow_field.v_star[(j, i)] - (dt / rho) * dp_dy;
+            }
+        }
+
+        // Fold the correction into the pressure. No under-relaxation: PISO
+        // corrects rather than iterates within the step.
         for j in 0..ny {
             for i in 0..nx {
                 flow_field.p[(j, i)] += flow_field.p_prime[(j, i)];
             }
         }
 
-        Ok(max_residual)
+        // What is left after the correction.
+        let mut mass_imbalance = 0.0;
+        for j in 0..ny {
+            for i in 0..nx {
+                let divergence_flux = rho
+                    * ((flow_field.u[(j, i + 1)] - flow_field.u[(j, i)]) * dy
+                        + (flow_field.v[(j + 1, i)] - flow_field.v[(j, i)]) * dx);
+                mass_imbalance += divergence_flux.abs();
+            }
         }
 
-    /// Correct velocities based on pressure correction
+        let reference = rho * self.config.reference_velocity * self.config.reference_length;
-    /// u^(n+1) = u* - (dt/ρ)∇p'
+        Ok(if reference > 0.0 {
-    fn correct_velocities(
+            mass_imbalance / reference
-        &self,
-        flow_field: &mut FlowField,
-        dt: f64,
-        rho: f64,
-        dx: f64,
-        dy: f64,
-    ) -> CfdResult<()> {
-        let (nx, ny, _, _) = flow_field.grid_info();
-
-        // Correct u-velocities
-        for j in 1..(ny - 1) {
-            for i in 1..nx {
-                let dp_dx = if i > 0 && i < nx {
-                    (flow_field.p_prime[(j, i.min(nx - 1))]
-                        - flow_field.p_prime[(j, (i - 1).max(0))])
-                        / dx
         } else {
-                    0.0
+            mass_imbalance
-                };
+        })
-
-                flow_field.u[(j, i)] = flow_field.u_star[(j, i)] - (dt / rho) * dp_dx;
-            }
-        }
-
-        // Correct v-velocities
-        for j in 1..ny {
-            for i in 1..(nx - 1) {
-                let dp_dy = if j > 0 && j < ny {
-                    (flow_field.p_prime[(j.min(ny - 1), i)]
-                        - flow_field.p_prime[((j - 1).max(0), i)])
-                        / dy
-                } else {
-                    0.0
-                };
-
-                flow_field.v[(j, i)] = flow_field.v_star[(j, i)] - (dt / rho) * dp_dy;
-            }
-        }
-
-        Ok(())
     }
 }
 
@@ -505,54 +448,35 @@ impl IncompressibleSolver for PisoSolver {
     ) -> CfdResult<Self::Result> {
         let start_time = std::time::Instant::now();
         let mut residual_history = Vec::new();
-        let (_nx, _ny, dx, dy) = flow_field.grid_info();
 
-        // Physical properties from config
+        // The state at the start of the step is what the explicit predictor
-        let rho = self.config.density;
+        // differentiates.
-        let nu = self.config.viscosity / rho; // kinematic viscosity
+        flow_field.apply_boundary_conditions(boundary_conditions)?;
-
-        // Store old values for time derivative
         flow_field.update_old_values();
 
-        // Apply boundary conditions
+        self.momentum_predictor(flow_field, dt)?;
-        flow_field.apply_boundary_conditions(boundary_conditions)?;
-
-        // STEP 1: MOMENTUM PREDICTOR
-        // Solve momentum equations with pressure from previous time step
-        // ∂u/∂t + ∇·(u⊗u) = -∇p/ρ + ν∇²u
-        self.solve_momentum_predictor(flow_field, dt, rho, nu)?;
 
         let mut total_corrector_steps = 0;
-
+        let mut final_residual = f64::INFINITY;
-        // PRESSURE-VELOCITY CORRECTION LOOP
+        for _corrector in 0..self.parameters.corrector_steps.max(1) {
-        for _corrector in 0..self.parameters.corrector_steps {
+            let mass_residual = self.project(flow_field, dt)?;
-            // STEP 2: PRESSURE CORRECTION
+            residual_history.push(mass_residual);
-            // Solve pressure Poisson equation: ∇²p' = ρ∇·u*/dt
+            final_residual = mass_residual;
-            let pressure_residual = self.solve_pressure_correction(flow_field, dt, rho, dx, dy)?;
-            residual_history.push(pressure_residual);
-
-            // STEP 3: VELOCITY CORRECTION
-            // Update velocities: u = u* - (dt/ρ)∇p'
-            self.correct_velocities(flow_field, dt, rho, dx, dy)?;
-
-            // Apply boundary conditions after correction
-            flow_field.apply_boundary_conditions(boundary_conditions)?;
-
             total_corrector_steps += 1;
 
-            // Check convergence
+            if mass_residual < self.parameters.tolerance {
-            if pressure_residual < self.parameters.tolerance {
                 break;
             }
+
+            // Re-project from the corrected field: with an explicit predictor
+            // the second pass mops up the inner solver's truncation.
+            flow_field.copy_to_starred();
         }
 
         let solve_time = start_time.elapsed();
-        let final_residual = residual_history.last().copied().unwrap_or(0.0);
-        let converged = final_residual < self.parameters.tolerance;
-
         Ok(PisoResult {
             solver_result: SolverResult {
-                converged,
+                converged: final_residual < self.parameters.tolerance,
                 iterations: total_corrector_steps,
                 final_residual,
                 residual_history,

crates/specialized/rtx-cfd/tests/mms_piso.rs

@@ -0,0 +1,216 @@
+//! Code verification of the PISO solver by manufactured solution.
+//!
+//! The manufactured field, its momentum source and the grid convention are
+//! exactly those of `tests/mms_navier_stokes.rs` — see that file for the
+//! derivation. PISO is a transient stepper, so instead of iterating an outer
+//! loop it is marched in time under the steady forcing until the field stops
+//! changing; the steady state it lands on satisfies the same spatial
+//! discretisation (first-order upwind convection, second-order diffusion,
+//! half-cell wall treatment), so the observed order should match SIMPLE's:
+//! approaching 1, limited by upwind's `O(h)` numerical viscosity.
+//!
+//! Until this file existed PISO had no verification of any kind — not a unit
+//! test, not a benchmark. The first run of this measurement, against the old
+//! implementation, is what confirmed the inverted pressure-correction sign
+//! and the frozen near-wall lines recorded in `piso.rs`'s module docs.
+
+use rtx_cfd::solvers::incompressible::{
+    BoundaryConditions, FlowField, IncompressibleSolver, PisoParameters, PisoSolver,
+};
+use rtx_cfd::{CfdConfig, CfdResult};
+use std::f64::consts::PI;
+
+const RHO: f64 = 1.0;
+const MU: f64 = 0.05;
+
+fn u_exact(x: f64, y: f64) -> f64 {
+    (PI * x).sin() * (PI * y).cos()
+}
+
+fn v_exact(x: f64, y: f64) -> f64 {
+    -(PI * x).cos() * (PI * y).sin()
+}
+
+fn source(x: f64, y: f64) -> (f64, f64) {
+    let fx = RHO * 0.5 * PI * (2.0 * PI * x).sin()
+        + 2.0 * PI * PI * MU * u_exact(x, y)
+        + PI * (PI * x).cos() * (PI * y).sin();
+    let fy = RHO * 0.5 * PI * (2.0 * PI * y).sin()
+        + 2.0 * PI * PI * MU * v_exact(x, y)
+        + PI * (PI * x).sin() * (PI * y).cos();
+    (fx, fy)
+}
+
+struct Measurement {
+    l2_velocity: f64,
+    max_div: f64,
+}
+
+/// March the manufactured problem on an `n` by `n` grid to steady state.
+async fn measure(n: usize) -> CfdResult<Measurement> {
+    let dx = 1.0 / n as f64;
+    let dy = dx;
+
+    // Explicit predictor: dt must respect the diffusion limit `dx^2 / (4 nu)`
+    // (the binding one here, with nu = 0.05 and |u| <= 1).
+    let nu = MU / RHO;
+    let dt = 0.4 * (dx * dx / (4.0 * nu)).min(dx);
+
+    let config = CfdConfig::new()
+        .with_density(RHO)
+        .with_viscosity(MU)
+        .with_reference_velocity(1.0)
+        .with_reference_length(1.0);
+    let params = PisoParameters {
+        corrector_steps: 2,
+        time_step: dt,
+        tolerance: 1e-8,
+    };
+    let mut solver = PisoSolver::new(config, params)?;
+    solver.set_momentum_source(source);
+    solver.set_wall_velocity(|x, y| (u_exact(x, y), v_exact(x, y)));
+
+    let mut field = FlowField::new(n, n, dx, dy)?;
+    for j in 0..n {
+        let y = (j as f64 + 0.5) * dy;
+        field.u[(j, 0)] = u_exact(0.0, y);
+        field.u[(j, n)] = u_exact(1.0, y);
+    }
+    for i in 0..n {
+        let x = (i as f64 + 0.5) * dx;
+        field.v[(0, i)] = v_exact(x, 0.0);
+        field.v[(n, i)] = v_exact(x, 1.0);
+    }
+
+    // March to steady state: stop when the field stops moving, measured as
+    // `max |u^{n+1} - u^n| / dt`, the discrete time derivative.
+    let empty = BoundaryConditions::new();
+    let mut steady_residual = f64::INFINITY;
+    for _step in 0..200_000 {
+        let u_before = field.u.clone();
+        let v_before = field.v.clone();
+
+        solver.solve_time_step(&mut field, &empty, dt).await?;
+
+        let mut max_change: f64 = 0.0;
+        for (a, b) in field.u.iter().zip(u_before.iter()) {
+            max_change = max_change.max((a - b).abs());
+        }
+        for (a, b) in field.v.iter().zip(v_before.iter()) {
+            max_change = max_change.max((a - b).abs());
+        }
+        steady_residual = max_change / dt;
+        if steady_residual < 1e-7 {
+            break;
+        }
+    }
+    assert!(
+        steady_residual < 1e-7,
+        "PISO did not reach a steady state: |du/dt| = {steady_residual:.3e}"
+    );
+
+    let mut squared = 0.0;
+    let mut volume = 0.0;
+    for j in 0..n {
+        for i in 1..n {
+            let e = field.u[(j, i)] - u_exact(i as f64 * dx, (j as f64 + 0.5) * dy);
+            squared += e * e * dx * dy;
+            volume += dx * dy;
+        }
+    }
+    for j in 1..n {
+        for i in 0..n {
+            let e = field.v[(j, i)] - v_exact((i as f64 + 0.5) * dx, j as f64 * dy);
+            squared += e * e * dx * dy;
+            volume += dx * dy;
+        }
+    }
+
+    let mut max_div: f64 = 0.0;
+    for j in 0..n {
+        for i in 0..n {
+            let div = (field.u[(j, i + 1)] - field.u[(j, i)]) / dx
+                + (field.v[(j + 1, i)] - field.v[(j, i)]) / dy;
+            max_div = max_div.max(div.abs());
+        }
+    }
+
+    Ok(Measurement {
+        l2_velocity: squared.sqrt() / volume.sqrt(),
+        max_div,
+    })
+}
+
+/// The steady state PISO marches to must converge to the exact solution at
+/// the rate the spatial discretisation dictates — order approaching 1 for
+/// first-order upwind — and must be divergence-free in every cell.
+///
+/// Measured (16 -> 32 -> 64): L2 velocity 3.516214e-2, 1.953750e-2,
+/// 1.037512e-2, orders 0.85 and 0.91, max |div u| ~ 1e-9 everywhere. The
+/// errors agree with SIMPLE's on the same meshes (3.516212e-2, 1.953751e-2,
+/// 1.037523e-2) to six or seven significant figures: two different
+/// algorithms — implicit under-relaxed outer iteration against explicit time
+/// marching with projection — land on the same discrete steady solution,
+/// which is exactly what sharing a spatial discretisation must produce and
+/// is very hard for two independently wrong solvers to fake.
+#[tokio::test]
+async fn piso_observed_order_matches_the_convection_scheme() -> CfdResult<()> {
+    let resolutions = [16usize, 32, 64];
+    let mut measurements = Vec::new();
+    for &n in &resolutions {
+        measurements.push(measure(n).await?);
+    }
+    let errors: Vec<f64> = measurements.iter().map(|m| m.l2_velocity).collect();
+    let rates: Vec<f64> = errors
+        .windows(2)
+        .map(|pair| (pair[0] / pair[1]).log2())
+        .collect();
+
+    for (i, &n) in resolutions.iter().enumerate() {
+        let rate = if i == 0 {
+            String::from("    -")
+        } else {
+            format!("{:5.2}", rates[i - 1])
+        };
+        println!(
+            "  n = {n:3}   L2 velocity error = {:.6e}   observed order = {rate}   \
+             max |div u| = {:.6e}",
+            errors[i], measurements[i].max_div
+        );
+    }
+
+    assert!(
+        errors.windows(2).all(|pair| pair[1] < pair[0]),
+        "the error must fall under refinement; got {errors:?}"
+    );
+
+    for (i, &rate) in rates.iter().enumerate() {
+        assert!(
+            rate > 0.75,
+            "refinement {} -> {}: observed order {rate:.3}, below the order 1 \
+             first-order upwind must deliver. Errors: {errors:?}",
+            resolutions[i],
+            resolutions[i + 1]
+        );
+        assert!(
+            rate < 2.3,
+            "refinement {} -> {}: observed order {rate:.3}, above what this \
+             scheme can deliver — suspect the error measure. Errors: {errors:?}",
+            resolutions[i],
+            resolutions[i + 1]
+        );
+    }
+
+    // Every cell, outer ring included, must satisfy continuity: the
+    // projection exists for no other reason.
+    for (m, &n) in measurements.iter().zip(&resolutions) {
+        assert!(
+            m.max_div < 1e-5,
+            "max |div u| = {:.3e} at n = {n}: the projection is not removing \
+             the divergence",
+            m.max_div
+        );
+    }
+
+    Ok(())
+}