rtx-cfd: OversetPisoSolver::momentum_residual — the background predictor's own staggered stencil (u_rhs/v_rhs, factored out of the predictor bit-identically) evaluated on every face of a NaN-masked field; solved faces read rounding, the active–fringe interface reads the composite's pressure level offset δ·h (cancels in the sum), prescribed fringe–fringe / fringe–hole faces read the stamping's momentum injection; hole ghosts (p, u, v) from a band widened three rows into the hole make every ring face evaluable; overset_cfd1 prints the buckets, the ring x-bands and δ at the settled state; pin: residual vanishes on the solved faces

Co-Authored-By: Claude Fable 5.1 <[email protected]>
Claude-Session: https://claude.ai/code/session_01X2GmJXeQ2njUecEKiJZ1G2

crates/specialized/rtx-cfd/src/solvers/incompressible/embedded/mod.rs

@@ -41,7 +41,7 @@ mod projection;
 use super::ale::{AleBoundaries, SideBoundary};
 use super::embedded_body::{EmbeddedBody, EmbeddedMask, FaceKind};
 use super::poisson::{
-    MgPrecision, MultigridParameters, PoissonProblem, PoissonSolverKind, solve_multigrid_pcg,
+    solve_multigrid_pcg, MgPrecision, MultigridParameters, PoissonProblem, PoissonSolverKind,
 };
 use super::simple::ConvectionScheme;
 use super::{FlowField, SolverResult};
@@ -480,147 +480,280 @@ impl EmbeddedPisoSolver {
         }
     }
 
+    /// The predictor's right-hand side on the u face `(j, i)`, `i = 1..nx`,
+    /// from `field.u_old`, `field.v_old` and `field.p`: `−conv + diff −
+    /// ∇p/ρ + f/ρ`, expression for expression the fixed-grid PISO's. The
+    /// predictor writes `u_old + dt · rhs` on the fluid faces; the overset
+    /// momentum-residual diagnostic (P4 option B) evaluates the same
+    /// operator on the prescribed faces, so "the solver's own stencil" is
+    /// this function by construction.
+    #[allow(clippy::too_many_lines)]
+    pub(crate) fn u_rhs(&self, field: &FlowField, j: usize, i: usize, t_old: f64) -> f64 {
+        let (nx, ny, dx, dy) = field.grid_info();
+        let rho = self.config.density;
+        let nu = self.config.viscosity / rho;
+        let b = self.parameters.boundaries;
+        let velocity = SideBoundary::Velocity;
+        let uo = &field.u_old;
+        let vo = &field.v_old;
+        let u_p = uo[(j, i)];
+
+        let ue_face = 0.5 * (uo[(j, i)] + uo[(j, i + 1)]);
+        let uw_face = 0.5 * (uo[(j, i - 1)] + uo[(j, i)]);
+
+        let south_is_wall = j == 0;
+        let north_is_wall = j + 1 == ny;
+
+        // Transverse face velocities from the stored v faces — on a
+        // domain side these are the prescribed boundary normals
+        // (zero on a wall, the outflow on an outlet). The fixed-grid
+        // PISO zeroes them on its walls, which is the same number on
+        // a wall and wrong on an outlet: the outgoing mass flux must
+        // carry momentum out, or the last row accumulates it.
+        let vn_face = 0.5 * (vo[(j + 1, i - 1)] + vo[(j + 1, i)]);
+        let vs_face = 0.5 * (vo[(j, i - 1)] + vo[(j, i)]);
+
+        // Upwind value across a domain side: the boundary function's
+        // tangential value on a Velocity side, the interior value
+        // otherwise (zero-gradient).
+        let beyond_north = if b.top == velocity {
+            self.boundary(i as f64 * dx, ny as f64 * dy, t_old).0
+        } else {
+            u_p
+        };
+        let beyond_south = if b.bottom == velocity {
+            self.boundary(i as f64 * dx, 0.0, t_old).0
+        } else {
+            u_p
+        };
+
+        let conv_x = (ue_face * Self::upwind(ue_face, uo[(j, i)], uo[(j, i + 1)])
+            - uw_face * Self::upwind(uw_face, uo[(j, i - 1)], uo[(j, i)]))
+            / dx;
+        let conv_y = (vn_face
+            * if north_is_wall {
+                Self::upwind(vn_face, u_p, beyond_north)
+            } else {
+                Self::upwind(vn_face, uo[(j, i)], uo[(j + 1, i)])
+            }
+            - vs_face
+                * if south_is_wall {
+                    Self::upwind(vs_face, beyond_south, u_p)
+                } else {
+                    Self::upwind(vs_face, uo[(j - 1, i)], uo[(j, i)])
+                })
+            / dy;
+
+        // Limited (TVD) corrections to the four convective face
+        // values; exactly zero-cost on the default upwind scheme.
+        let scheme = self.parameters.convection_scheme;
+        let mut conv_x = conv_x;
+        let mut conv_y = conv_y;
+        if scheme != ConvectionScheme::Upwind {
+            let delta_e = if ue_face >= 0.0 {
+                scheme.face_correction(Some(uo[(j, i - 1)]), uo[(j, i)], uo[(j, i + 1)])
+            } else {
+                let far = (i + 2 <= nx).then(|| uo[(j, i + 2)]);
+                scheme.face_correction(far, uo[(j, i + 1)], uo[(j, i)])
+            };
+            let delta_w = if uw_face >= 0.0 {
+                let far = (i >= 2).then(|| uo[(j, i - 2)]);
+                scheme.face_correction(far, uo[(j, i - 1)], uo[(j, i)])
+            } else {
+                scheme.face_correction(Some(uo[(j, i + 1)]), uo[(j, i)], uo[(j, i - 1)])
+            };
+            let delta_n = if north_is_wall {
+                0.0
+            } else if vn_face >= 0.0 {
+                let far = (j >= 1).then(|| uo[(j - 1, i)]);
+                scheme.face_correction(far, uo[(j, i)], uo[(j + 1, i)])
+            } else {
+                let far = (j + 2 < ny).then(|| uo[(j + 2, i)]);
+                scheme.face_correction(far, uo[(j + 1, i)], uo[(j, i)])
+            };
+            let delta_s = if south_is_wall {
+                0.0
+            } else if vs_face >= 0.0 {
+                let far = (j >= 2).then(|| uo[(j - 2, i)]);
+                scheme.face_correction(far, uo[(j - 1, i)], uo[(j, i)])
+            } else {
+                let far = (j + 1 < ny).then(|| uo[(j + 1, i)]);
+                scheme.face_correction(far, uo[(j, i)], uo[(j - 1, i)])
+            };
+            conv_x += (ue_face * delta_e - uw_face * delta_w) / dx;
+            conv_y += (vn_face * delta_n - vs_face * delta_s) / 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 on a
+        // Velocity side; a slip wall or outlet carries no shear.
+        let flux_north = if north_is_wall {
+            if b.top == velocity {
+                let u_wall = self.boundary(i as f64 * dx, ny as f64 * dy, t_old).0;
+                nu * (u_wall - u_p) / (0.5 * dy)
+            } else {
+                0.0
+            }
+        } else {
+            nu * (uo[(j + 1, i)] - u_p) / dy
+        };
+        let flux_south = if south_is_wall {
+            if b.bottom == velocity {
+                let u_wall = self.boundary(i as f64 * dx, 0.0, t_old).0;
+                nu * (u_p - u_wall) / (0.5 * dy)
+            } else {
+                0.0
+            }
+        } else {
+            nu * (u_p - uo[(j - 1, i)]) / dy
+        };
+        let diff_y = (flux_north - flux_south) / dy;
+
+        let pressure_gradient = -(field.p[(j, i)] - field.p[(j, i - 1)]) / (rho * dx);
+
+        let body_force = self.momentum_source.as_ref().map_or(0.0, |f| {
+            f(i as f64 * dx, (j as f64 + 0.5) * dy, t_old).0 / rho
+        });
+
+        -conv_x - conv_y + diff_x + diff_y + pressure_gradient + body_force
+    }
+
+    /// See [`Self::u_rhs`]: the v face `(j, i)`, `j = 1..ny`.
+    #[allow(clippy::too_many_lines)]
+    pub(crate) fn v_rhs(&self, field: &FlowField, j: usize, i: usize, t_old: f64) -> f64 {
+        let (nx, ny, dx, dy) = field.grid_info();
+        let rho = self.config.density;
+        let nu = self.config.viscosity / rho;
+        let b = self.parameters.boundaries;
+        let velocity = SideBoundary::Velocity;
+        let uo = &field.u_old;
+        let vo = &field.v_old;
+        let v_p = vo[(j, i)];
+
+        let vn_face = 0.5 * (vo[(j, i)] + vo[(j + 1, i)]);
+        let vs_face = 0.5 * (vo[(j - 1, i)] + vo[(j, i)]);
+
+        let west_is_wall = i == 0;
+        let east_is_wall = i + 1 == nx;
+
+        let ue_face = 0.5 * (uo[(j - 1, i + 1)] + uo[(j, i + 1)]);
+        let uw_face = 0.5 * (uo[(j - 1, i)] + uo[(j, i)]);
+        let beyond_east = if b.right == velocity {
+            self.boundary(nx as f64 * dx, j as f64 * dy, t_old).1
+        } else {
+            v_p
+        };
+        let beyond_west = if b.left == velocity {
+            self.boundary(0.0, j as f64 * dy, t_old).1
+        } else {
+            v_p
+        };
+
+        let conv_y = (vn_face * Self::upwind(vn_face, vo[(j, i)], vo[(j + 1, i)])
+            - vs_face * Self::upwind(vs_face, vo[(j - 1, i)], vo[(j, i)]))
+            / dy;
+        let conv_x = (ue_face
+            * if east_is_wall {
+                Self::upwind(ue_face, v_p, beyond_east)
+            } else {
+                Self::upwind(ue_face, vo[(j, i)], vo[(j, i + 1)])
+            }
+            - uw_face
+                * if west_is_wall {
+                    Self::upwind(uw_face, beyond_west, v_p)
+                } else {
+                    Self::upwind(uw_face, vo[(j, i - 1)], vo[(j, i)])
+                })
+            / dx;
+
+        let scheme = self.parameters.convection_scheme;
+        let mut conv_x = conv_x;
+        let mut conv_y = conv_y;
+        if scheme != ConvectionScheme::Upwind {
+            let delta_n = if vn_face >= 0.0 {
+                scheme.face_correction(Some(vo[(j - 1, i)]), vo[(j, i)], vo[(j + 1, i)])
+            } else {
+                let far = (j + 2 <= ny).then(|| vo[(j + 2, i)]);
+                scheme.face_correction(far, vo[(j + 1, i)], vo[(j, i)])
+            };
+            let delta_s = if vs_face >= 0.0 {
+                let far = (j >= 2).then(|| vo[(j - 2, i)]);
+                scheme.face_correction(far, vo[(j - 1, i)], vo[(j, i)])
+            } else {
+                scheme.face_correction(Some(vo[(j + 1, i)]), vo[(j, i)], vo[(j - 1, i)])
+            };
+            let delta_e = if east_is_wall {
+                0.0
+            } else if ue_face >= 0.0 {
+                let far = (i >= 1).then(|| vo[(j, i - 1)]);
+                scheme.face_correction(far, vo[(j, i)], vo[(j, i + 1)])
+            } else {
+                let far = (i + 2 < nx).then(|| vo[(j, i + 2)]);
+                scheme.face_correction(far, vo[(j, i + 1)], vo[(j, i)])
+            };
+            let delta_w = if west_is_wall {
+                0.0
+            } else if uw_face >= 0.0 {
+                let far = (i >= 2).then(|| vo[(j, i - 2)]);
+                scheme.face_correction(far, vo[(j, i - 1)], vo[(j, i)])
+            } else {
+                let far = (i + 1 < nx).then(|| vo[(j, i + 1)]);
+                scheme.face_correction(far, vo[(j, i)], vo[(j, i - 1)])
+            };
+            conv_y += (vn_face * delta_n - vs_face * delta_s) / dy;
+            conv_x += (ue_face * delta_e - uw_face * delta_w) / 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 {
+            if b.right == velocity {
+                let v_wall = self.boundary(nx as f64 * dx, j as f64 * dy, t_old).1;
+                nu * (v_wall - v_p) / (0.5 * dx)
+            } else {
+                0.0
+            }
+        } else {
+            nu * (vo[(j, i + 1)] - v_p) / dx
+        };
+        let flux_west = if west_is_wall {
+            if b.left == velocity {
+                let v_wall = self.boundary(0.0, j as f64 * dy, t_old).1;
+                nu * (v_p - v_wall) / (0.5 * dx)
+            } else {
+                0.0
+            }
+        } else {
+            nu * (v_p - vo[(j, i - 1)]) / dx
+        };
+        let diff_x = (flux_east - flux_west) / dx;
+
+        let pressure_gradient = -(field.p[(j, i)] - field.p[(j - 1, i)]) / (rho * dy);
+
+        let body_force = self.momentum_source.as_ref().map_or(0.0, |f| {
+            f((i as f64 + 0.5) * dx, j as f64 * dy, t_old).1 / rho
+        });
+
+        -conv_x - conv_y + diff_x + diff_y + pressure_gradient + body_force
+    }
+
     /// Explicit momentum predictor on the fluid faces, expression for
     /// expression the fixed-grid PISO's (so the no-body case is identical
     /// to the bit), plus the slip-wall / outlet arms of the ALE solver on
     /// the domain sides. Non-fluid faces keep their prescribed values.
     #[allow(clippy::too_many_lines)]
     fn momentum_predictor(&self, field: &mut FlowField, dt: f64, t_old: f64) -> CfdResult<()> {
-        let (nx, ny, dx, dy) = field.grid_info();
-        let rho = self.config.density;
-        let nu = self.config.viscosity / rho;
+        let (nx, ny, _, _) = field.grid_info();
         let b = self.parameters.boundaries;
-        let velocity = SideBoundary::Velocity;
 
         for j in 0..ny {
             for i in 1..nx {
                 if !self.u_is_fluid(j, i) {
                     continue;
                 }
-                let uo = &field.u_old;
-                let vo = &field.v_old;
-                let u_p = uo[(j, i)];
-
-                let ue_face = 0.5 * (uo[(j, i)] + uo[(j, i + 1)]);
-                let uw_face = 0.5 * (uo[(j, i - 1)] + uo[(j, i)]);
-
-                let south_is_wall = j == 0;
-                let north_is_wall = j + 1 == ny;
-
-                // Transverse face velocities from the stored v faces — on a
-                // domain side these are the prescribed boundary normals
-                // (zero on a wall, the outflow on an outlet). The fixed-grid
-                // PISO zeroes them on its walls, which is the same number on
-                // a wall and wrong on an outlet: the outgoing mass flux must
-                // carry momentum out, or the last row accumulates it.
-                let vn_face = 0.5 * (vo[(j + 1, i - 1)] + vo[(j + 1, i)]);
-                let vs_face = 0.5 * (vo[(j, i - 1)] + vo[(j, i)]);
-
-                // Upwind value across a domain side: the boundary function's
-                // tangential value on a Velocity side, the interior value
-                // otherwise (zero-gradient).
-                let beyond_north = if b.top == velocity {
-                    self.boundary(i as f64 * dx, ny as f64 * dy, t_old).0
-                } else {
-                    u_p
-                };
-                let beyond_south = if b.bottom == velocity {
-                    self.boundary(i as f64 * dx, 0.0, t_old).0
-                } else {
-                    u_p
-                };
-
-                let conv_x = (ue_face * Self::upwind(ue_face, uo[(j, i)], uo[(j, i + 1)])
-                    - uw_face * Self::upwind(uw_face, uo[(j, i - 1)], uo[(j, i)]))
-                    / dx;
-                let conv_y = (vn_face
-                    * if north_is_wall {
-                        Self::upwind(vn_face, u_p, beyond_north)
-                    } else {
-                        Self::upwind(vn_face, uo[(j, i)], uo[(j + 1, i)])
-                    }
-                    - vs_face
-                        * if south_is_wall {
-                            Self::upwind(vs_face, beyond_south, u_p)
-                        } else {
-                            Self::upwind(vs_face, uo[(j - 1, i)], uo[(j, i)])
-                        })
-                    / dy;
-
-                // Limited (TVD) corrections to the four convective face
-                // values; exactly zero-cost on the default upwind scheme.
-                let scheme = self.parameters.convection_scheme;
-                let mut conv_x = conv_x;
-                let mut conv_y = conv_y;
-                if scheme != ConvectionScheme::Upwind {
-                    let delta_e = if ue_face >= 0.0 {
-                        scheme.face_correction(Some(uo[(j, i - 1)]), uo[(j, i)], uo[(j, i + 1)])
-                    } else {
-                        let far = (i + 2 <= nx).then(|| uo[(j, i + 2)]);
-                        scheme.face_correction(far, uo[(j, i + 1)], uo[(j, i)])
-                    };
-                    let delta_w = if uw_face >= 0.0 {
-                        let far = (i >= 2).then(|| uo[(j, i - 2)]);
-                        scheme.face_correction(far, uo[(j, i - 1)], uo[(j, i)])
-                    } else {
-                        scheme.face_correction(Some(uo[(j, i + 1)]), uo[(j, i)], uo[(j, i - 1)])
-                    };
-                    let delta_n = if north_is_wall {
-                        0.0
-                    } else if vn_face >= 0.0 {
-                        let far = (j >= 1).then(|| uo[(j - 1, i)]);
-                        scheme.face_correction(far, uo[(j, i)], uo[(j + 1, i)])
-                    } else {
-                        let far = (j + 2 < ny).then(|| uo[(j + 2, i)]);
-                        scheme.face_correction(far, uo[(j + 1, i)], uo[(j, i)])
-                    };
-                    let delta_s = if south_is_wall {
-                        0.0
-                    } else if vs_face >= 0.0 {
-                        let far = (j >= 2).then(|| uo[(j - 2, i)]);
-                        scheme.face_correction(far, uo[(j - 1, i)], uo[(j, i)])
-                    } else {
-                        let far = (j + 1 < ny).then(|| uo[(j + 1, i)]);
-                        scheme.face_correction(far, uo[(j, i)], uo[(j - 1, i)])
-                    };
-                    conv_x += (ue_face * delta_e - uw_face * delta_w) / dx;
-                    conv_y += (vn_face * delta_n - vs_face * delta_s) / 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 on a
-                // Velocity side; a slip wall or outlet carries no shear.
-                let flux_north = if north_is_wall {
-                    if b.top == velocity {
-                        let u_wall = self.boundary(i as f64 * dx, ny as f64 * dy, t_old).0;
-                        nu * (u_wall - u_p) / (0.5 * dy)
-                    } else {
-                        0.0
-                    }
-                } else {
-                    nu * (uo[(j + 1, i)] - u_p) / dy
-                };
-                let flux_south = if south_is_wall {
-                    if b.bottom == velocity {
-                        let u_wall = self.boundary(i as f64 * dx, 0.0, t_old).0;
-                        nu * (u_p - u_wall) / (0.5 * dy)
-                    } else {
-                        0.0
-                    }
-                } else {
-                    nu * (u_p - uo[(j - 1, i)]) / dy
-                };
-                let diff_y = (flux_north - flux_south) / dy;
-
-                let pressure_gradient = -(field.p[(j, i)] - field.p[(j, i - 1)]) / (rho * dx);
-
-                let body_force = self.momentum_source.as_ref().map_or(0.0, |f| {
-                    f(i as f64 * dx, (j as f64 + 0.5) * dy, t_old).0 / rho
-                });
-
-                field.u[(j, i)] = u_p
-                    + dt * (-conv_x - conv_y + diff_x + diff_y + pressure_gradient + body_force);
+                let rhs = self.u_rhs(field, j, i, t_old);
+                field.u[(j, i)] = field.u_old[(j, i)] + dt * rhs;
             }
         }
 
@@ -629,116 +762,8 @@ impl EmbeddedPisoSolver {
                 if !self.v_is_fluid(j, i) {
                     continue;
                 }
-                let uo = &field.u_old;
-                let vo = &field.v_old;
-                let v_p = vo[(j, i)];
-
-                let vn_face = 0.5 * (vo[(j, i)] + vo[(j + 1, i)]);
-                let vs_face = 0.5 * (vo[(j - 1, i)] + vo[(j, i)]);
-
-                let west_is_wall = i == 0;
-                let east_is_wall = i + 1 == nx;
-
-                let ue_face = 0.5 * (uo[(j - 1, i + 1)] + uo[(j, i + 1)]);
-                let uw_face = 0.5 * (uo[(j - 1, i)] + uo[(j, i)]);
-                let beyond_east = if b.right == velocity {
-                    self.boundary(nx as f64 * dx, j as f64 * dy, t_old).1
-                } else {
-                    v_p
-                };
-                let beyond_west = if b.left == velocity {
-                    self.boundary(0.0, j as f64 * dy, t_old).1
-                } else {
-                    v_p
-                };
-
-                let conv_y = (vn_face * Self::upwind(vn_face, vo[(j, i)], vo[(j + 1, i)])
-                    - vs_face * Self::upwind(vs_face, vo[(j - 1, i)], vo[(j, i)]))
-                    / dy;
-                let conv_x = (ue_face
-                    * if east_is_wall {
-                        Self::upwind(ue_face, v_p, beyond_east)
-                    } else {
-                        Self::upwind(ue_face, vo[(j, i)], vo[(j, i + 1)])
-                    }
-                    - uw_face
-                        * if west_is_wall {
-                            Self::upwind(uw_face, beyond_west, v_p)
-                        } else {
-                            Self::upwind(uw_face, vo[(j, i - 1)], vo[(j, i)])
-                        })
-                    / dx;
-
-                let scheme = self.parameters.convection_scheme;
-                let mut conv_x = conv_x;
-                let mut conv_y = conv_y;
-                if scheme != ConvectionScheme::Upwind {
-                    let delta_n = if vn_face >= 0.0 {
-                        scheme.face_correction(Some(vo[(j - 1, i)]), vo[(j, i)], vo[(j + 1, i)])
-                    } else {
-                        let far = (j + 2 <= ny).then(|| vo[(j + 2, i)]);
-                        scheme.face_correction(far, vo[(j + 1, i)], vo[(j, i)])
-                    };
-                    let delta_s = if vs_face >= 0.0 {
-                        let far = (j >= 2).then(|| vo[(j - 2, i)]);
-                        scheme.face_correction(far, vo[(j - 1, i)], vo[(j, i)])
-                    } else {
-                        scheme.face_correction(Some(vo[(j + 1, i)]), vo[(j, i)], vo[(j - 1, i)])
-                    };
-                    let delta_e = if east_is_wall {
-                        0.0
-                    } else if ue_face >= 0.0 {
-                        let far = (i >= 1).then(|| vo[(j, i - 1)]);
-                        scheme.face_correction(far, vo[(j, i)], vo[(j, i + 1)])
-                    } else {
-                        let far = (i + 2 < nx).then(|| vo[(j, i + 2)]);
-                        scheme.face_correction(far, vo[(j, i + 1)], vo[(j, i)])
-                    };
-                    let delta_w = if west_is_wall {
-                        0.0
-                    } else if uw_face >= 0.0 {
-                        let far = (i >= 2).then(|| vo[(j, i - 2)]);
-                        scheme.face_correction(far, vo[(j, i - 1)], vo[(j, i)])
-                    } else {
-                        let far = (i + 1 < nx).then(|| vo[(j, i + 1)]);
-                        scheme.face_correction(far, vo[(j, i)], vo[(j, i - 1)])
-                    };
-                    conv_y += (vn_face * delta_n - vs_face * delta_s) / dy;
-                    conv_x += (ue_face * delta_e - uw_face * delta_w) / 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 {
-                    if b.right == velocity {
-                        let v_wall = self.boundary(nx as f64 * dx, j as f64 * dy, t_old).1;
-                        nu * (v_wall - v_p) / (0.5 * dx)
-                    } else {
-                        0.0
-                    }
-                } else {
-                    nu * (vo[(j, i + 1)] - v_p) / dx
-                };
-                let flux_west = if west_is_wall {
-                    if b.left == velocity {
-                        let v_wall = self.boundary(0.0, j as f64 * dy, t_old).1;
-                        nu * (v_p - v_wall) / (0.5 * dx)
-                    } else {
-                        0.0
-                    }
-                } else {
-                    nu * (v_p - vo[(j, i - 1)]) / dx
-                };
-                let diff_x = (flux_east - flux_west) / dx;
-
-                let pressure_gradient = -(field.p[(j, i)] - field.p[(j - 1, i)]) / (rho * dy);
-
-                let body_force = self.momentum_source.as_ref().map_or(0.0, |f| {
-                    f((i as f64 + 0.5) * dx, j as f64 * dy, t_old).1 / rho
-                });
-
-                field.v[(j, i)] = v_p
-                    + dt * (-conv_x - conv_y + diff_x + diff_y + pressure_gradient + body_force);
+                let rhs = self.v_rhs(field, j, i, t_old);
+                field.v[(j, i)] = field.v_old[(j, i)] + dt * rhs;
             }
         }