outram_foam_basic_lib/ldu_matrix/
fv_vector_matrix.rs

1// SPDX-License-Identifier: GPL-3.0-only
2// Copyright (C) 2026 OUTRAM PARK contributors
3// Derived from OpenFOAM (www.openfoam.com)
4// Copyright (C) 2004-2023 OpenFOAM Foundation
5// Copyright (C) 2016-2023 OpenCFD Ltd.
6//
7// This file is part of OUTRAM PARK.
8//
9// OUTRAM PARK is free software: you can redistribute it and/or modify it
10// under the terms of the GNU General Public License as published by the
11// Free Software Foundation, either version 3 of the License, or (at your
12// option) any later version.
13//
14// OUTRAM PARK is distributed in the hope that it will be useful, but
15// WITHOUT ANY WARRANTY; without even the implied warranty of
16// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
17// General Public License for more details.
18//
19// You should have received a copy of the GNU General Public License along
20// with OUTRAM PARK.  If not, see <https://www.gnu.org/licenses/>.
21
22use std::ops::{Add, AddAssign, Neg, Sub, SubAssign};
23use std::sync::Arc;
24
25use super::fv_matrix::{SolverPerformance, SolverSettings};
26use super::ldu_matrix::LduMatrix;
27use super::solvers::gauss_seidel::gauss_seidel as gs;
28use crate::fields::boundary::bc::PatchField;
29use crate::fields::field::Field;
30use crate::fields::vol_field::{VolScalarField, VolVectorField};
31use crate::mesh::fv_mesh::FvMesh;
32use crate::primitives::Vector3;
33
34/// Implicit vector equation `A·U = b` for a `VolVectorField`.
35///
36/// Mirrors `Foam::fvVectorMatrix` (`fvMatrix<vector>`).
37///
38/// The LDU coefficients are **scalar** — they multiply the entire velocity
39/// vector equally in all three directions.  The source vector is a
40/// `Field<Vector3>`.  Solving decomposes into three independent scalar
41/// Gauss-Seidel solves (one per component).
42#[derive(Debug, Clone)]
43pub struct FvVectorMatrix {
44    /// Mesh the equation is defined on (shares the face addressing).
45    pub mesh: Arc<FvMesh>,
46    /// Scalar LDU coefficients of the operator `A` (shared by all 3 components).
47    pub ldu: LduMatrix,
48    /// Right-hand-side vector source per cell, length `n_cells`.
49    pub source: Field<Vector3>,
50}
51
52impl FvVectorMatrix {
53    /// Allocate a zero-initialised vector matrix for `mesh` (zero coefficients,
54    /// zero source).
55    pub fn new(mesh: Arc<FvMesh>) -> Self {
56        let n_cells = mesh.n_cells;
57        let owner = mesh.owner[..mesh.n_internal_faces].to_vec();
58        let neighbour = mesh.neighbour.clone();
59        Self {
60            ldu: LduMatrix::new(n_cells, owner, neighbour),
61            source: Field::from_fn(n_cells, |_| Vector3::ZERO),
62            mesh,
63        }
64    }
65
66    /// Add `coeff[c]` to the diagonal of cell `c` (e.g. a time-derivative term).
67    pub fn add_to_diag(&mut self, coeff: &Field<f64>) {
68        for c in 0..self.mesh.n_cells {
69            self.ldu.diag[c] += coeff[c];
70        }
71    }
72
73    /// Add `term[c]` to the vector source of cell `c`.
74    pub fn add_to_source(&mut self, term: &Field<Vector3>) {
75        for c in 0..self.mesh.n_cells {
76            self.source[c] = self.source[c] + term[c];
77        }
78    }
79
80    /// Pin one cell's velocity to a fixed value (reference cell for closed domains).
81    pub fn set_reference(&mut self, cell: usize, value: Vector3) {
82        self.ldu.diag[cell] += 1e30;
83        self.source[cell] = self.source[cell] + value * 1e30;
84    }
85
86    /// Diagonal coefficient per cell: `A[c] = diag[c]`.
87    ///
88    /// Used in PISO: `rAU = 1 / UEqn.a_field()`.
89    pub fn a_field(&self) -> VolScalarField {
90        let mesh = self.mesh.clone();
91        let boundary = mesh
92            .patches
93            .iter()
94            .map(|p| PatchField::zero_gradient(p.size))
95            .collect();
96        VolScalarField::new(
97            "A",
98            mesh.clone(),
99            Field::new(self.ldu.diag.clone()),
100            boundary,
101        )
102    }
103
104    /// Off-diagonal + source residual: `H[c] = source[c] − Σ off-diag · U`.
105    ///
106    /// For a zero field x this returns `source[c]` directly.
107    /// Used in PISO: `HbyA = rAU * UEqn.h_field(U)`.
108    pub fn h_field(&self, u: &VolVectorField) -> VolVectorField {
109        let mesh = self.mesh.clone();
110        let n = mesh.n_cells;
111        let mut h = vec![Vector3::ZERO; n];
112        for c in 0..n {
113            h[c] = self.source[c];
114        }
115        for f in 0..mesh.n_internal_faces {
116            let o = self.ldu.owner[f];
117            let nb = self.ldu.neighbour[f];
118            h[o] = h[o] - u.internal[nb] * self.ldu.upper[f];
119            h[nb] = h[nb] - u.internal[o] * self.ldu.lower[f];
120        }
121        let boundary = mesh
122            .patches
123            .iter()
124            .map(|p| PatchField::zero_gradient_vec(p.size))
125            .collect();
126        VolVectorField::new("H", mesh, Field::new(h), boundary)
127    }
128
129    /// Solve each component (x, y, z) as an independent scalar Gauss-Seidel problem.
130    pub fn solve(
131        &self,
132        name: &str,
133        settings: SolverSettings,
134    ) -> (VolVectorField, SolverPerformance) {
135        let mesh = self.mesh.clone();
136        let n = mesh.n_cells;
137
138        let bx: Vec<f64> = (0..n).map(|c| self.source[c].x).collect();
139        let by: Vec<f64> = (0..n).map(|c| self.source[c].y).collect();
140        let bz: Vec<f64> = (0..n).map(|c| self.source[c].z).collect();
141
142        let mut xx = vec![0.0_f64; n];
143        let mut xy = vec![0.0_f64; n];
144        let mut xz = vec![0.0_f64; n];
145
146        let (ix, rx) = gs(
147            &self.ldu,
148            &bx,
149            &mut xx,
150            settings.tolerance,
151            settings.max_iter,
152        );
153        let (iy, ry) = gs(
154            &self.ldu,
155            &by,
156            &mut xy,
157            settings.tolerance,
158            settings.max_iter,
159        );
160        let (iz, rz) = gs(
161            &self.ldu,
162            &bz,
163            &mut xz,
164            settings.tolerance,
165            settings.max_iter,
166        );
167
168        let internal = Field::from_fn(n, |c| Vector3::new(xx[c], xy[c], xz[c]));
169        let boundary = mesh
170            .patches
171            .iter()
172            .map(|p| PatchField::zero_gradient_vec(p.size))
173            .collect();
174        let u = VolVectorField::new(name, mesh, internal, boundary);
175
176        let n_iters = ix.max(iy).max(iz);
177        let final_res = rx.max(ry).max(rz);
178        let perf = SolverPerformance {
179            n_iterations: n_iters,
180            final_residual: final_res,
181            converged: final_res < settings.tolerance,
182        };
183        (u, perf)
184    }
185}
186
187// ── Arithmetic ────────────────────────────────────────────────────────────────
188
189impl Add for FvVectorMatrix {
190    type Output = Self;
191    fn add(mut self, rhs: Self) -> Self {
192        for (a, b) in self.ldu.diag.iter_mut().zip(&rhs.ldu.diag) {
193            *a += b;
194        }
195        for (a, b) in self.ldu.lower.iter_mut().zip(&rhs.ldu.lower) {
196            *a += b;
197        }
198        for (a, b) in self.ldu.upper.iter_mut().zip(&rhs.ldu.upper) {
199            *a += b;
200        }
201        for c in 0..self.source.len() {
202            self.source[c] = self.source[c] + rhs.source[c];
203        }
204        self
205    }
206}
207
208impl Sub for FvVectorMatrix {
209    type Output = Self;
210    fn sub(mut self, rhs: Self) -> Self {
211        for (a, b) in self.ldu.diag.iter_mut().zip(&rhs.ldu.diag) {
212            *a -= b;
213        }
214        for (a, b) in self.ldu.lower.iter_mut().zip(&rhs.ldu.lower) {
215            *a -= b;
216        }
217        for (a, b) in self.ldu.upper.iter_mut().zip(&rhs.ldu.upper) {
218            *a -= b;
219        }
220        for c in 0..self.source.len() {
221            self.source[c] = self.source[c] - rhs.source[c];
222        }
223        self
224    }
225}
226
227impl Neg for FvVectorMatrix {
228    type Output = Self;
229    fn neg(mut self) -> Self {
230        for x in self.ldu.diag.iter_mut() {
231            *x = -*x;
232        }
233        for x in self.ldu.lower.iter_mut() {
234            *x = -*x;
235        }
236        for x in self.ldu.upper.iter_mut() {
237            *x = -*x;
238        }
239        for c in 0..self.source.len() {
240            self.source[c] = -self.source[c];
241        }
242        self
243    }
244}
245
246impl AddAssign for FvVectorMatrix {
247    fn add_assign(&mut self, rhs: Self) {
248        for (a, b) in self.ldu.diag.iter_mut().zip(&rhs.ldu.diag) {
249            *a += b;
250        }
251        for (a, b) in self.ldu.lower.iter_mut().zip(&rhs.ldu.lower) {
252            *a += b;
253        }
254        for (a, b) in self.ldu.upper.iter_mut().zip(&rhs.ldu.upper) {
255            *a += b;
256        }
257        for c in 0..self.source.len() {
258            self.source[c] = self.source[c] + rhs.source[c];
259        }
260    }
261}
262
263impl SubAssign for FvVectorMatrix {
264    fn sub_assign(&mut self, rhs: Self) {
265        for (a, b) in self.ldu.diag.iter_mut().zip(&rhs.ldu.diag) {
266            *a -= b;
267        }
268        for (a, b) in self.ldu.lower.iter_mut().zip(&rhs.ldu.lower) {
269            *a -= b;
270        }
271        for (a, b) in self.ldu.upper.iter_mut().zip(&rhs.ldu.upper) {
272            *a -= b;
273        }
274        for c in 0..self.source.len() {
275            self.source[c] = self.source[c] - rhs.source[c];
276        }
277    }
278}
279
280#[cfg(test)]
281mod tests {
282    use super::*;
283    use crate::mesh::fv_mesh::{BoundaryPatch, FvMeshBuilder, PatchKind};
284
285    fn unit_mesh() -> Arc<FvMesh> {
286        Arc::new(
287            FvMeshBuilder::new()
288                .n_cells(2)
289                .n_internal_faces(1)
290                .owner(vec![0, 1, 0])
291                .neighbour(vec![1])
292                .patches(vec![
293                    BoundaryPatch::new("right", 1, 1, PatchKind::Wall),
294                    BoundaryPatch::new("left", 2, 1, PatchKind::Wall),
295                ])
296                .cell_volumes(vec![1.0, 1.0])
297                .cell_centres(vec![
298                    Vector3::new(0.25, 0.0, 0.0),
299                    Vector3::new(0.75, 0.0, 0.0),
300                ])
301                .face_area_vectors(vec![
302                    Vector3::new(1.0, 0.0, 0.0),
303                    Vector3::new(1.0, 0.0, 0.0),
304                    Vector3::new(-1.0, 0.0, 0.0),
305                ])
306                .face_centres(vec![
307                    Vector3::new(0.5, 0.0, 0.0),
308                    Vector3::new(1.0, 0.0, 0.0),
309                    Vector3::new(0.0, 0.0, 0.0),
310                ])
311                .build()
312                .unwrap(),
313        )
314    }
315
316    #[test]
317    fn diagonal_system_solves_vector() {
318        let m = unit_mesh();
319        let mut mat = FvVectorMatrix::new(m.clone());
320        mat.ldu.diag[0] = 2.0;
321        mat.ldu.diag[1] = 3.0;
322        mat.source[0] = Vector3::new(4.0, 6.0, 8.0);
323        mat.source[1] = Vector3::new(6.0, 9.0, 12.0);
324        let (u, perf) = mat.solve("U", SolverSettings::default());
325        assert!(perf.converged, "residual = {}", perf.final_residual);
326        assert!((u.internal[0].x - 2.0).abs() < 1e-8);
327        assert!((u.internal[0].y - 3.0).abs() < 1e-8);
328        assert!((u.internal[1].x - 2.0).abs() < 1e-8);
329    }
330
331    #[test]
332    fn a_field_returns_diagonal() {
333        let m = unit_mesh();
334        let mut mat = FvVectorMatrix::new(m.clone());
335        mat.ldu.diag[0] = 5.0;
336        mat.ldu.diag[1] = 7.0;
337        let a = mat.a_field();
338        assert!((a.internal[0] - 5.0).abs() < 1e-12);
339        assert!((a.internal[1] - 7.0).abs() < 1e-12);
340    }
341}
outram_foam_basic_lib/ldu_matrix/fv_vector_matrix.rs

outram_foam_basic_lib/ldu_matrix/
fv_vector_matrix.rs
