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    pub mesh: Arc<FvMesh>,
45    pub ldu: LduMatrix,
46    pub source: Field<Vector3>,
47}
48
49impl FvVectorMatrix {
50    pub fn new(mesh: Arc<FvMesh>) -> Self {
51        let n_cells = mesh.n_cells;
52        let owner = mesh.owner[..mesh.n_internal_faces].to_vec();
53        let neighbour = mesh.neighbour.clone();
54        Self {
55            ldu: LduMatrix::new(n_cells, owner, neighbour),
56            source: Field::from_fn(n_cells, |_| Vector3::ZERO),
57            mesh,
58        }
59    }
60
61    pub fn add_to_diag(&mut self, coeff: &Field<f64>) {
62        for c in 0..self.mesh.n_cells {
63            self.ldu.diag[c] += coeff[c];
64        }
65    }
66
67    pub fn add_to_source(&mut self, term: &Field<Vector3>) {
68        for c in 0..self.mesh.n_cells {
69            self.source[c] = self.source[c] + term[c];
70        }
71    }
72
73    /// Pin one cell's velocity to a fixed value (reference cell for closed domains).
74    pub fn set_reference(&mut self, cell: usize, value: Vector3) {
75        self.ldu.diag[cell] += 1e30;
76        self.source[cell] = self.source[cell] + value * 1e30;
77    }
78
79    /// Diagonal coefficient per cell: `A[c] = diag[c]`.
80    ///
81    /// Used in PISO: `rAU = 1 / UEqn.a_field()`.
82    pub fn a_field(&self) -> VolScalarField {
83        let mesh = self.mesh.clone();
84        let boundary = mesh
85            .patches
86            .iter()
87            .map(|p| PatchField::zero_gradient(p.size))
88            .collect();
89        VolScalarField::new(
90            "A",
91            mesh.clone(),
92            Field::new(self.ldu.diag.clone()),
93            boundary,
94        )
95    }
96
97    /// Off-diagonal + source residual: `H[c] = source[c] − Σ off-diag · U`.
98    ///
99    /// For a zero field x this returns `source[c]` directly.
100    /// Used in PISO: `HbyA = rAU * UEqn.h_field(U)`.
101    pub fn h_field(&self, u: &VolVectorField) -> VolVectorField {
102        let mesh = self.mesh.clone();
103        let n = mesh.n_cells;
104        let mut h = vec![Vector3::ZERO; n];
105        for c in 0..n {
106            h[c] = self.source[c];
107        }
108        for f in 0..mesh.n_internal_faces {
109            let o = self.ldu.owner[f];
110            let nb = self.ldu.neighbour[f];
111            h[o] = h[o] - u.internal[nb] * self.ldu.upper[f];
112            h[nb] = h[nb] - u.internal[o] * self.ldu.lower[f];
113        }
114        let boundary = mesh
115            .patches
116            .iter()
117            .map(|p| PatchField::zero_gradient_vec(p.size))
118            .collect();
119        VolVectorField::new("H", mesh, Field::new(h), boundary)
120    }
121
122    /// Solve each component (x, y, z) as an independent scalar Gauss-Seidel problem.
123    pub fn solve(
124        &self,
125        name: &str,
126        settings: SolverSettings,
127    ) -> (VolVectorField, SolverPerformance) {
128        let mesh = self.mesh.clone();
129        let n = mesh.n_cells;
130
131        let bx: Vec<f64> = (0..n).map(|c| self.source[c].x).collect();
132        let by: Vec<f64> = (0..n).map(|c| self.source[c].y).collect();
133        let bz: Vec<f64> = (0..n).map(|c| self.source[c].z).collect();
134
135        let mut xx = vec![0.0_f64; n];
136        let mut xy = vec![0.0_f64; n];
137        let mut xz = vec![0.0_f64; n];
138
139        let (ix, rx) = gs(
140            &self.ldu,
141            &bx,
142            &mut xx,
143            settings.tolerance,
144            settings.max_iter,
145        );
146        let (iy, ry) = gs(
147            &self.ldu,
148            &by,
149            &mut xy,
150            settings.tolerance,
151            settings.max_iter,
152        );
153        let (iz, rz) = gs(
154            &self.ldu,
155            &bz,
156            &mut xz,
157            settings.tolerance,
158            settings.max_iter,
159        );
160
161        let internal = Field::from_fn(n, |c| Vector3::new(xx[c], xy[c], xz[c]));
162        let boundary = mesh
163            .patches
164            .iter()
165            .map(|p| PatchField::zero_gradient_vec(p.size))
166            .collect();
167        let u = VolVectorField::new(name, mesh, internal, boundary);
168
169        let n_iters = ix.max(iy).max(iz);
170        let final_res = rx.max(ry).max(rz);
171        let perf = SolverPerformance {
172            n_iterations: n_iters,
173            final_residual: final_res,
174            converged: final_res < settings.tolerance,
175        };
176        (u, perf)
177    }
178}
179
180// ── Arithmetic ────────────────────────────────────────────────────────────────
181
182impl Add for FvVectorMatrix {
183    type Output = Self;
184    fn add(mut self, rhs: Self) -> Self {
185        for (a, b) in self.ldu.diag.iter_mut().zip(&rhs.ldu.diag) {
186            *a += b;
187        }
188        for (a, b) in self.ldu.lower.iter_mut().zip(&rhs.ldu.lower) {
189            *a += b;
190        }
191        for (a, b) in self.ldu.upper.iter_mut().zip(&rhs.ldu.upper) {
192            *a += b;
193        }
194        for c in 0..self.source.len() {
195            self.source[c] = self.source[c] + rhs.source[c];
196        }
197        self
198    }
199}
200
201impl Sub for FvVectorMatrix {
202    type Output = Self;
203    fn sub(mut self, rhs: Self) -> Self {
204        for (a, b) in self.ldu.diag.iter_mut().zip(&rhs.ldu.diag) {
205            *a -= b;
206        }
207        for (a, b) in self.ldu.lower.iter_mut().zip(&rhs.ldu.lower) {
208            *a -= b;
209        }
210        for (a, b) in self.ldu.upper.iter_mut().zip(&rhs.ldu.upper) {
211            *a -= b;
212        }
213        for c in 0..self.source.len() {
214            self.source[c] = self.source[c] - rhs.source[c];
215        }
216        self
217    }
218}
219
220impl Neg for FvVectorMatrix {
221    type Output = Self;
222    fn neg(mut self) -> Self {
223        for x in self.ldu.diag.iter_mut() {
224            *x = -*x;
225        }
226        for x in self.ldu.lower.iter_mut() {
227            *x = -*x;
228        }
229        for x in self.ldu.upper.iter_mut() {
230            *x = -*x;
231        }
232        for c in 0..self.source.len() {
233            self.source[c] = -self.source[c];
234        }
235        self
236    }
237}
238
239impl AddAssign for FvVectorMatrix {
240    fn add_assign(&mut self, rhs: Self) {
241        for (a, b) in self.ldu.diag.iter_mut().zip(&rhs.ldu.diag) {
242            *a += b;
243        }
244        for (a, b) in self.ldu.lower.iter_mut().zip(&rhs.ldu.lower) {
245            *a += b;
246        }
247        for (a, b) in self.ldu.upper.iter_mut().zip(&rhs.ldu.upper) {
248            *a += b;
249        }
250        for c in 0..self.source.len() {
251            self.source[c] = self.source[c] + rhs.source[c];
252        }
253    }
254}
255
256impl SubAssign for FvVectorMatrix {
257    fn sub_assign(&mut self, rhs: Self) {
258        for (a, b) in self.ldu.diag.iter_mut().zip(&rhs.ldu.diag) {
259            *a -= b;
260        }
261        for (a, b) in self.ldu.lower.iter_mut().zip(&rhs.ldu.lower) {
262            *a -= b;
263        }
264        for (a, b) in self.ldu.upper.iter_mut().zip(&rhs.ldu.upper) {
265            *a -= b;
266        }
267        for c in 0..self.source.len() {
268            self.source[c] = self.source[c] - rhs.source[c];
269        }
270    }
271}
272
273#[cfg(test)]
274mod tests {
275    use super::*;
276    use crate::mesh::fv_mesh::{BoundaryPatch, FvMeshBuilder, PatchKind};
277
278    fn unit_mesh() -> Arc<FvMesh> {
279        Arc::new(
280            FvMeshBuilder::new()
281                .n_cells(2)
282                .n_internal_faces(1)
283                .owner(vec![0, 1, 0])
284                .neighbour(vec![1])
285                .patches(vec![
286                    BoundaryPatch::new("right", 1, 1, PatchKind::Wall),
287                    BoundaryPatch::new("left", 2, 1, PatchKind::Wall),
288                ])
289                .cell_volumes(vec![1.0, 1.0])
290                .cell_centres(vec![
291                    Vector3::new(0.25, 0.0, 0.0),
292                    Vector3::new(0.75, 0.0, 0.0),
293                ])
294                .face_area_vectors(vec![
295                    Vector3::new(1.0, 0.0, 0.0),
296                    Vector3::new(1.0, 0.0, 0.0),
297                    Vector3::new(-1.0, 0.0, 0.0),
298                ])
299                .face_centres(vec![
300                    Vector3::new(0.5, 0.0, 0.0),
301                    Vector3::new(1.0, 0.0, 0.0),
302                    Vector3::new(0.0, 0.0, 0.0),
303                ])
304                .build()
305                .unwrap(),
306        )
307    }
308
309    #[test]
310    fn diagonal_system_solves_vector() {
311        let m = unit_mesh();
312        let mut mat = FvVectorMatrix::new(m.clone());
313        mat.ldu.diag[0] = 2.0;
314        mat.ldu.diag[1] = 3.0;
315        mat.source[0] = Vector3::new(4.0, 6.0, 8.0);
316        mat.source[1] = Vector3::new(6.0, 9.0, 12.0);
317        let (u, perf) = mat.solve("U", SolverSettings::default());
318        assert!(perf.converged, "residual = {}", perf.final_residual);
319        assert!((u.internal[0].x - 2.0).abs() < 1e-8);
320        assert!((u.internal[0].y - 3.0).abs() < 1e-8);
321        assert!((u.internal[1].x - 2.0).abs() < 1e-8);
322    }
323
324    #[test]
325    fn a_field_returns_diagonal() {
326        let m = unit_mesh();
327        let mut mat = FvVectorMatrix::new(m.clone());
328        mat.ldu.diag[0] = 5.0;
329        mat.ldu.diag[1] = 7.0;
330        let a = mat.a_field();
331        assert!((a.internal[0] - 5.0).abs() < 1e-12);
332        assert!((a.internal[1] - 7.0).abs() < 1e-12);
333    }
334}
outram_foam_basic_lib/ldu_matrix/fv_vector_matrix.rs

outram_foam_basic_lib/ldu_matrix/
fv_vector_matrix.rs
