outram-foam-basic-lib

⚠️ Unverified until validated. All code in this workspace is unverified and untrusted unless a specific verification & validation (V&V) case demonstrates otherwise. V&V cases are human-reviewed and are intended for journal / arXiv publication — that is the trust workflow. See the workspace VERIFICATION_AND_VALIDATION.md and RESPONSIBLE_USE.md. Not for nuclear facility operation, reactor control, safety-critical, or licensing decisions.

This is OUTRAM PARK's independent Rust translation of selected OpenFOAM® algorithms. It is not the official OpenFOAM® software and is not affiliated with, endorsed by, or sanctioned by OpenCFD Ltd. or the ESI Group. OpenFOAM® is a registered trademark of OpenCFD Limited — see TRADEMARKS.md for the full attribution and non-affiliation notice. Translated from OpenFOAM/OpenFOAM-dev, master branch — no commit is pinned (translation was done by reading the C++ source directly, not from an ongoing codegen-from-clone pipeline); see upstream_source/README.md for the full provenance record.

Pure-Rust translation of the OpenFOAM primitive and finite-volume library layer — tensor algebra, polynomial solvers, ODE solvers, interpolation, thermophysics, fields, mesh, FV operators, and fluid/solid thermodynamics needed to build compressible and conjugate-heat-transfer CFD solvers.

Quick start

[dependencies]
outram-foam-basic-lib = "0.1.0"

use outram_foam_basic_lib::prelude::*;

// Tensor algebra
let u = Vector3::new(1.0, 0.0, 0.0);
let v = Vector3::new(0.0, 1.0, 0.0);
let cross = u.cross(v);              // z-axis unit vector

// Polynomial root finding
let cubic = CubicEqn::new(1.0, -6.0, 11.0, -6.0); // (x-1)(x-2)(x-3)
let roots = cubic.roots();

// Dense LU solve (no external BLAS required)
let mut m = SquareMatrix::new(3);
m.set(0, 0, 2.0); m.set(0, 1, 1.0); m.set(0, 2, 0.0);
m.set(1, 0, 1.0); m.set(1, 1, 3.0); m.set(1, 2, 1.0);
m.set(2, 0, 0.0); m.set(2, 1, 1.0); m.set(2, 2, 4.0);
let x = m.solve(&[7.0, 10.0, 15.0]);  // returns Vec<f64>

// FV operators
use outram_foam_basic_lib::fv_operators::{fvm, fvc};
// fvm::ddt, fvm::laplacian, fvm::div, fvm::laplacian_vec, fvm::div_vec,
// fvm::ddt_coeff, fvm::ddt_coeff_vec
// fvc::grad, fvc::div, fvc::flux, fvc::reconstruct, fvc::buoyancy_flux, ...

What's implemented

Layers 1a–1h — Primitives and thermophysics

Layer 2 — Fields and mesh

Layer 3 — Finite-volume operators

Layer 4 — Field-level fluid thermodynamics

Limitations

Read this before depending on the crate. Every item below is grounded in the current source; nothing here is aspirational. See also the unverified-until-validated banner at the top.

Scope boundaries

• Layers 1–4 only — no solver loops. This crate provides the mathematical building blocks (tensor algebra, primitives, FV operators, thermophysics kernels, field thermodynamics). The Layer 5 logic — PISO/PIMPLE loops, multi-region coupling drivers, turbulence-model registries — lives in separate downstream crates (openfoam-icof, openfoam-cht, openfoam-rho), which are not yet published. On their own the operators here do not run a time-marching CFD simulation.
• No turbulence models. LaminarModel / kOmegaSST and divDevRhoReff are planned for the CHT solver crate, not present here. Only laminar terms can be assembled.
• No OpenFOAM case / polyMesh file I/O. Meshes are constructed in code via FvMeshBuilder (you supply owner/neighbour connectivity and the geometric Vecs) or via interface::one_dimensional_meshing::create_one_d_mesh for 1-D system-code meshes. Reading an OpenFOAM constant/polyMesh directory is a downstream concern, not in this crate. (The DIC/warm-start discussion below refers to a read_poly_mesh that lives in the appbuilder crate.)
• Serial only. No MPI / domain decomposition; there are no processor boundary patches. The linear solvers and operators run single-threaded on one mesh.

Finite-volume numerics

• Implicit convection is first-order upwind only. fvm::div / fvm::div_vec assemble a first-order upwind matrix (src/fv_operators/fvm/div.rs). There is no implicit higher-order/TVD convection (linearUpwind, limitedLinear, vanLeer, etc.). Second-order TVD reconstruction exists only explicitly, via fvc::muscl (Upwind / Linear / limiter variants), intended for density-based central-upwind fluxes — the Kurganov–Tadmor flux assembly itself is Layer-5 work and is not wired up here.
• Laplacian is Gauss-orthogonal only — no non-orthogonal / skewness correction. fvm::laplacian and fvc::sn_grad use the orthogonal face-area-over-centre-distance coefficient (src/fv_operators/fvm/laplacian.rs), ignoring the angle between the face-area vector and the owner→neighbour vector. Diffusion is therefore accurate on orthogonal meshes (hex / 1-D) but under-resolved on non-orthogonal or skewed meshes; there is no over-relaxed/minimum-correction non-orthogonal term.
• Time integration is first-order. Only implicit Euler fvm::ddt (plus fvm::ddt_coeff, the vector forms, and a d2dt2) is provided. No higher-order/backward/Crank–Nicolson time schemes.
• Boundary conditions are a small closed set. BoundaryCondition (src/fields/boundary/bc.rs) supports FixedValue, FixedField, ZeroGradient, Symmetry, Empty, and Calculated only. There is no inletOutlet, non-zero fixedGradient, mixed/Robin, wall-function, cyclic/periodic, or processor boundary condition.

Linear solvers

• Krylov / multigrid solvers assume symmetric positive-definite systems. conjugate_gradient (DIC-preconditioned) and gamg target the symmetric pressure system. Asymmetric systems fall back to gauss_seidel; there is no BiCGStab/GMRES.
• GAMG rebuilds its hierarchy on every solve and is slower than DIC-PCG below ~10^4–105 cells (see the performance section below). It is mesh-independent in V-cycle count but not yet cached across time steps or usable as a CG preconditioner. PCG remains the default.
• Dense matrices provide LU only. SquareMatrix implements LU with scaled partial pivoting (solve returns Result<Vec<f64>, MatrixError>). The QR / Cholesky / SVD factorisations named in the design doc are not implemented.

Thermophysics accuracy (known-failing / unverified)

• Peng–Robinson EOS is inaccurate above the critical pressure. Density errors of 7–26 % vs NIST are observed for Pr > 1 / near-critical states; the corresponding NIST cross-checks are #[ignore]-d pending a suspected Z-root-selection or α-function fix (src/thermophysics/eos/peng_robinson.rs, docs/porting-roadmap.md). Prefer PerfectGas, RhoConst, or IcoPolynomial where density accuracy matters.
• JANAF T(H) inversion can fail to converge from a far initial guess. The Newton solve stalls crossing the Tcommon = 1000 K coefficient discontinuity when started far away (e.g. t0 = 100 K targeting 3000 K); one such test is #[ignore]-d (src/thermophysics/thermo/janaf.rs). Seed the iteration with a reasonable temperature and handle Err(NonConvergent).

ODE and polynomial restrictions

• The stiff Rosenbrock23 solver needs a user Jacobian. OdeSystem::jacobian has a default that panics (unimplemented!, src/ode/mod.rs); only the explicit Euler and Rkf45 solvers work without one.
• Polynomial<N>::integral() (degree-raising form) is not implemented. The type-level definite integral that returns a Polynomial<{N+1}> needs nightly generic_const_exprs; use the scalar integral(x1, x2) -> f64 / integral_minus1 forms instead (src/polynomial/polynomial.rs).

Validation status

• Unverified until validated. This crate's own test suite is unit / verification tests against analytic and reference values, not full-solver validation. The only end-to-end physics validation cited (lid-driven cavity, Ghia Re = 100) is run in the downstream outram-foam-appbuilder-lib crate, not here. Treat all outputs as untrusted draft until a specific V&V case covers your usage. Not for nuclear facility operation, reactor control, safety-critical, or licensing decisions.

Platform / build

• Pure Rust, no BLAS required by the library. Runtime dependencies are only thiserror, uom, and ndarray, so the library builds headless (including Android and other no-system-BLAS targets). The LAPACK comparison benchmark (tests/matrix_bench.rs) does need system OpenBLAS (Linux) / Intel MKL (Windows/macOS); it is a dev-dependency, target-gated off Android, and is not part of the library build.

SquareMatrix vs LAPACK benchmark

SquareMatrix::solve (pure-Rust LU, no external BLAS) compared to ndarray-linalg / OpenBLAS Array2::solve (LAPACK DGESV) — release mode, Linux x86-64, 2026-06-24:

SquareMatrix is faster for n ≤ 10 because OpenBLAS DGESV has ~300–400 ns of per-call FFI overhead that dominates at small sizes. The crossover is around n ≈ 20–50. For typical finite-volume networks (10–50 unknowns per implicit system), SquareMatrix eliminates the system BLAS dependency with no performance penalty. Reproduce with:

cargo test -p outram-foam-basic-lib --test matrix_bench --release -- --nocapture

Pressure linear-solver performance (DIC preconditioner + warm start)

The pressure Poisson solve dominates the cost of an incompressible transient run, so the PCG path here mirrors OpenFOAM's two key efficiency choices. The reasoning, measured on the pimple_foam_cavity Ghia Re=100 case (outram-foam-appbuilder-lib, 41×41 mesh, 6000 steps to t=12 s):

The diagnosis. Instrumenting every pressure solve showed 12 000 solves (6000 steps × 2 PISO correctors), each running a flat 272 PCG iterations — the same count on the first step and the last — for ≈ 3.26 million sparse matrix-vector products. That was essentially the whole wall-clock time. Two things were wrong:

1. Cold start every solve. conjugate_gradient began from x = 0 each call, so even near steady state — where the pressure barely changes between steps — it paid full convergence from scratch.
2. Jacobi preconditioner. The preconditioner was diagonal only (M⁻¹ = 1/diag), whose PCG iteration count scales with the mesh (∝ √κ ≈ O(Nₓ)).

The fix (both taken straight from the OpenFOAM source).

• Warm start. conjugate_gradient now accepts an initial guess x0, and FvMatrix::solve_cg_with_guess seeds the solve with the previous time step's field. This is exactly what Foam::PCG::scalarSolve does — it computes the initial residual as source − A·psi from the incoming psi, i.e. the field passed in is the initial guess. A transient solver near steady state then converges in a few iterations (the solver returns immediately, 0 iterations, if the guess already meets the tolerance).
• DIC preconditioner. The Jacobi preconditioner was replaced with DIC (Diagonal-based Incomplete Cholesky), a faithful port of Foam::DICPreconditioner: a one-time reciprocal-diagonal factorisation (calcReciprocalD) plus a forward/backward face sweep (precondition). It requires faces in upper-triangular order (owner[f] < neighbour[f]), which is how read_poly_mesh loads OpenFOAM polyMesh files.

Result. Per-solve iterations dropped from a flat 272 to ~73 (cold, DIC alone) and ~40 near steady state (DIC + warm start) — 6.6× fewer total CG iterations (3.26M → 0.50M) and the cavity test roughly halved in wall time, with the solution field bit-for-bit unchanged (Ghia RMS 0.0113).

GAMG (algebraic multigrid) — implemented; mesh-independent but not a win at small sizes. gamg is a serial port of OpenFOAM's GAMGSolver with algebraicPairGAMGAgglomeration: pairwise agglomeration on the matrix coefficients (|upper|), Galerkin coarse operators (Pᵀ A P), a recursive correction-scheme V-cycle with Gauss-Seidel pre/post smoothing and an optimal line-search correction scale, and a dense-LU coarsest solve. Its V-cycle count stays flat as the mesh is refined (verified in gamg::tests::gamg_cycle_count_is_mesh_independent) — the property PCG lacks.

On the 1681-cell cavity, however, GAMG reproduces the PCG solution bit-for-bit (field diff 2.6e-9) but runs slower (42.8 s vs 12.4 s): at this scale DIC-PCG already converges in ~40 cheap iterations per corrector, below the crossover (~10⁴–10⁵ cells) where multigrid's mesh-independence pays off, and the free gamg function rebuilds the hierarchy on every solve. So PCG stays the default; GAMG is the right tool once larger meshes are in play. Two follow-ups would make GAMG competitive at scale: cache the agglomeration across time steps (re-restrict coefficients only), and use GAMG as a CG preconditioner for Krylov robustness.

Prelude

use outram_foam_basic_lib::prelude::*;

Includes all tensor types, polynomial solvers, math functions, SquareMatrix, ODE solvers, interpolation, all thermophysics types, all field and mesh types, all LDU matrix types, FV operator modules (fvc, fvm), adjust_phi, and all fluid/solid thermo types.

Running tests

# Library unit tests (no external BLAS required)
cargo test -p outram-foam-basic-lib --lib --tests

# Matrix benchmark (release mode for meaningful numbers)
cargo test -p outram-foam-basic-lib --test matrix_bench --release -- --nocapture

Layer roadmap

✅ Layer 1a — Tensor algebra      (Vector3, Tensor, SymmTensor, SphericalTensor)
✅ Layer 1b — Dense matrices      (SquareMatrix: LU with scaled partial pivoting)
✅ Layer 1c — Polynomial eqns     (LinearEqn, QuadraticEqn, CubicEqn, Roots<N>)
✅ Layer 1d — Polynomial eval     (Polynomial<N>: Horner, derivative, integral)
✅ Layer 1e — ODE solvers         (Euler, Rkf45, Rosenbrock23; OdeSystem trait)
✅ Layer 1f — Interpolation       (interpolate_xy, interpolate_spline_xy)
✅ Layer 1g — Math functions      (erf_inv, inc_gamma_*, inv_inc_gamma)
✅ Layer 1h — Thermophysics       (EOS, Thermo, Transport traits + 4 impls each)
✅ Layer 2  — Fields + Mesh       (VolField, SurfaceField, FvMesh, LduMatrix,
                                   FvMatrix, FvVectorMatrix, RegionInterface,
                                   Gauss-Seidel + CG solvers)
✅ Layer 3  — FV operators        (fvm: ddt, ddt_coeff, ddt_vec, ddt_coeff_vec,
                                        laplacian, laplacian_vec, div, div_vec
                                   fvc: grad, div, interpolate, sn_grad, flux,
                                        reconstruct, ddt_corr, buoyancy_flux
                                   adjust_phi)
✅ Layer 4  — Field thermodynamics (FluidThermo, PsiThermo, RhoThermo,
                                    SolidThermo, ConstSolidThermo)
⬜ Layer 5  — Solver logic         (icoFoam PISO loop → openfoam-icof;
                                    chtMultiRegionFoam → openfoam-cht;
                                    rhoPimpleFoam → openfoam-rho)

Changelog

Unreleased

• Added GAMG (algebraic multigrid) solver (ldu_matrix::gamg, FvMatrix::solve_gamg{,_with_guess}). Serial port of Foam::GAMGSolver + algebraicPairGAMGAgglomeration: pairwise agglomeration, Galerkin coarse operators, recursive V-cycle (GS pre/post-smoothing, line-search correction scale, dense-LU coarsest). V-cycle count is mesh-independent. On the 1681-cell cavity it reproduces the PCG field exactly but is slower at that scale (see "Pressure linear-solver performance"); PCG remains the default.

• Pressure PCG: DIC preconditioner + warm start (≈6.6× fewer iterations). conjugate_gradient now uses a DIC preconditioner (Foam::DICPreconditioner port) instead of Jacobi and takes an optional initial guess x0; FvMatrix::solve_cg_with_guess warm-starts from the previous field. On the pimple_foam_cavity Ghia Re=100 case (41×41) this cut per-solve PCG iterations from a flat 272 to ~40 and total CG iterations 3.26M → 0.50M, with the solution unchanged. See "Pressure linear-solver performance" above. Note: conjugate_gradient gained a parameter (x0) — update callers from conjugate_gradient(&ldu, &b, &settings) to conjugate_gradient(&ldu, &b, None, &settings).

• Fixed an exponential-memory bug in field arithmetic (the "24 GB cavity"). VolField/SurfaceField's Add/Sub/Neg rebuilt the field's name string on every call (self.name = format!("({} + {})", self.name, rhs.name)). In a solver loop where a persistent field is reassigned from an expression containing itself — e.g. rho = rho + div(phi), with phi's name embedding interpolate(rho) — each operation glued two copies of the previous name together, so the name String doubled every timestep (2^step). The rhoPimpleFoam compressible_lid_cavity test grew to ~1.3 GB by step 24 and was killed by OOM at step 25, with each step running ~2× slower than the last.

  The trap: the field data was completely healthy — internal/boundary Vecs stayed the correct size and the physics (|U|, |p|, |phi|) stayed bounded — so it presented as a "leak"/hang rather than a numerical blow-up. It was localised by printing field.name.len() per step: rho/phi names doubled while freshly-solve()d fields (U, p, he) stayed constant.

  Fix: arithmetic operators now leave self.name as the left operand's name. This also matches OpenFOAM, where a GeometricField keeps its fixed registered IOobject name and solvers print residual labels rather than accumulating an expression string onto the field. After the fix the cavity runs at a flat ~0.45 ms/step and ~3.7 MB RSS. See the Translation-notes callout in CLAUDE.md for the full write-up.

License

GPL-3.0-only (matching the upstream OpenFOAM sources).