cfd-fem-solver

How it is built

Architecture overview

An engineering simulation solver turns a physical problem (air over a wing, stress in a bracket, heat in a block) into a large sparse system of equations and solves it. The geometry is cut into a mesh of cells or elements; a discretisation (finite volume in most CFD codes, finite elements in most structural and multiphysics codes) turns the partial differential equations into algebraic ones per cell; the contributions are assembled into a sparse matrix; a solver stack (Krylov iteration, preconditioner, multigrid or a direct factorisation) solves it; and a time loop and a nonlinear loop repeat this until the answer stops changing. At scale the mesh is split across processes, and each iteration exchanges a thin layer of boundary values with neighbours. The one idea that makes it work is that the problem is local (a cell only touches its neighbours), so the matrix is sparse, the work is mostly sparse matrix-vector products, and the parallel cost is a surface-to-volume ratio. In almost every real code, the linear solve takes most of the time, and the pressure equation of an incompressible flow is the hardest part of it.

Parts:

• Mesh and geometry — builds and stores the mesh — delaunay-triangulation (no page yet), exact-geometric-predicates (no page yet), half-edge-mesh (no page yet), unstructured-mesh-representation (no page yet), space-filling-curves (no page yet)
• Finite-volume discretisation — conservative per-cell balances of fluxes — finite-volume-method (no page yet), flux-limiters-upwinding (no page yet), Finite differences
• Finite-element discretisation — weak form, element integrals, assembly — Finite elements, Gaussian quadrature
• Time and nonlinear loops — marches in time, iterates on nonlinearity — Explicit Euler, Runge-Kutta 4, Stiff solvers (implicit BDF), Newton root finding
• Pressure-velocity coupling — enforces incompressibility — pressure-projection-method (no page yet), Multigrid, Conjugate gradient
• Linear algebra — stores and solves the sparse systems — Sparse matrices (CSR), Conjugate gradient, gmres (no page yet), Incomplete-factorization preconditioners, LU and Cholesky, Sparse direct solver, Fill-reducing ordering, thomas-tridiagonal-solver (no page yet)
• Parallel layer — splits the mesh, exchanges boundaries, reduces globally — Graph partitioning, Halo exchange, Collective communication, eager-rendezvous-protocol (no page yet), rdma (no page yet), Zero-copy I/O
• Supporting numerics — priority queues and random numbers inside the above — d-ary heap, PRNG (PCG, xoshiro)

 geometry ─► mesh generation ─► partition ─► per-rank mesh + ghost cells
                                                   │
         ┌─────────────── time loop ───────────────▼───────────────────────────┐
         │  nonlinear loop (Picard / Newton / SIMPLE-PISO sweeps)              │
         │    assemble: cell/element loop ─► CSR matrix A, vector b            │
         │    solve:  Krylov (CG / GMRES) + preconditioner (ILU / multigrid)   │
         │              each iteration: SpMV ──halo exchange──► neighbours     │
         │                             dot products ──allreduce──► all ranks   │
         │    update fields, check residual                                    │
         └──────────────────────────────────────────────────────────────────────┘
                                                   │ write fields (parallel I/O) ─► post-processing

How the parts fit together

Mesh and geometry

Everything downstream depends on mesh quality: skewed or tiny cells degrade accuracy and the conditioning of the matrix. Triangle and tetrahedral meshes come from delaunay-triangulation (no page yet), which maximises the minimum angle; it is only reliable with exact-geometric-predicates (no page yet), because a wrong orientation or in-circle test on nearly degenerate input corrupts the mesh (Shewchuk's adaptive-precision predicates are the standard fix). Editable meshes (refinement, smoothing) use a half-edge-mesh (no page yet). The solver does not use that structure: it uses unstructured-mesh-representation (no page yet), cell-face-node connectivity in compressed arrays, because the solver loops over faces and cells and wants contiguous memory. OpenFOAM's mesh is polyhedral and built from points, faces and cells, with no restriction on the number of faces per cell (source: OpenFOAM user guide, Mesh description); a flux is computed once per face and added to the two cells that share it, so assembly is one pass over faces. Cells are renumbered along space-filling-curves (no page yet) or with a bandwidth-reducing order so neighbouring cells sit near each other in memory; this raises cache hits in every sparse operation and also gives a cheap first partition.

Finite-volume discretisation

finite-volume-method (no page yet) writes each conservation law as a balance on a cell: change of the stored quantity equals the net flux through its faces. Conservation is exact by construction because the flux leaving one cell is the flux entering its neighbour, which is why CFD uses it. Each flux needs a face value from the surrounding cells. Central interpolation is accurate but oscillates in convection-dominated flow; pure upwind is stable but smears. flux-limiters-upwinding (no page yet) blend the two so that the scheme is upwind near sharp gradients and higher order in smooth regions (Sweby's flux-limiter analysis is the reference). Gradients at cell centres use stencils built on the same ideas as Finite differences. The output of this step per equation is one sparse row per cell: a diagonal and one off-diagonal entry per face neighbour.

Finite-element discretisation

Finite elements start from the weak form: multiply the equation by a test function and integrate by parts. The domain integral becomes a sum of element integrals, each computed by Gaussian quadrature on a reference element and mapped to the real one. Each element yields a small dense matrix and vector, which assembly adds into the global CSR matrix using the local-to-global dof map. Higher-order elements give more accuracy per unknown at the cost of denser rows; matching the quadrature order to the integrand matters, because too low an order loses accuracy and too high wastes time. FEniCS expresses the weak form in a Unified Form Language and generates the element kernels from it, so the assembler is a generic loop over cells that calls generated code.

Time and nonlinear loops

Transient problems step in time. Explicit Euler and Runge-Kutta 4 need no linear solve but are limited by a stability bound tied to the cell size (the CFL condition), so fine meshes force tiny steps. Stiff problems (fast chemistry, diffusion on fine meshes) use implicit schemes such as BDF, which take larger steps at the cost of a nonlinear solve per step. That solve is Newton: linearise, solve J dx = -F, repeat, with each Newton step ending in a sparse linear solve that can be made inexact (stop the inner solve early) as long as the outer residual still falls (source: Dembo, Eisenstat and Steihaug 1982 on inexact Newton methods). For steady problems the same loop runs as pseudo-time or under-relaxed iterations.

Pressure-velocity coupling

Incompressible flow has no equation for pressure; it is whatever field makes the velocity divergence-free. The pressure-projection-method (no page yet) (Chorin 1968) advances velocity without the pressure term, then solves a Poisson equation for the pressure that projects the result onto a divergence-free field, then corrects the velocity. SIMPLE (Patankar and Spalding 1972) applies the same split to steady problems with under-relaxation; PISO and PIMPLE add correction loops for transient runs, and OpenFOAM's solution-control file configures SIMPLE, PISO and PIMPLE loops and GAMG (geometric-algebraic multigrid) solvers (source: OpenFOAM user guide, fvSolution). The pressure equation is symmetric and elliptic, so it is solved with Conjugate gradient preconditioned by Multigrid, whose cost grows about linearly with the number of cells. A code that needs to solve many pressure equations per time step spends most of its time here, so this is the part to optimise first.

Linear algebra

Matrices live in CSR (or a blocked variant). A sparse matrix-vector product reads each nonzero once and does 2 flops on it; with an 8-byte value and a 4-byte index that is about 0.17 flop per byte, far below what a modern CPU or GPU can sustain, so SpMV is limited by memory bandwidth and not arithmetic (derived from the format; see the roofline block for the model). The solver menu:

• Conjugate gradient for symmetric positive-definite systems (pressure, elasticity, diffusion).
• gmres (no page yet) for nonsymmetric systems (convection-dominated momentum, Newton Jacobians); it needs memory for a growing basis, so it is restarted (Saad and Schultz 1986).
• Incomplete-factorization preconditioners: ILU or IC approximate the inverse cheaply and cut the iteration count; they are local, so they parallelise as block-Jacobi or additive Schwarz.
• Multigrid uses a hierarchy of coarser problems to remove smooth error that a smoother cannot; algebraic variants build the hierarchy from the matrix alone, so they work on unstructured meshes.
• LU and Cholesky through a Sparse direct solver: robust for hard or small problems; the permutation from Fill-reducing ordering decides how much fill the factors get and therefore memory and time. Direct solvers scale worse than multigrid for large 3-D problems, which is why iterative solvers dominate big CFD runs.
• thomas-tridiagonal-solver (no page yet) solves tridiagonal systems in O(n); it is the kernel of line-implicit smoothers and alternating-direction schemes on structured blocks.

Parallel layer

Graph partitioning cuts the cell graph into one piece per process, balancing cells per part and minimising cut faces, since cut faces are communication (Karypis and Kumar's multilevel scheme in METIS is the standard). A d-ary heap holds the gain priorities in refinement passes, and randomised tie-breaking from a seeded PRNG keeps runs reproducible. Each process owns its cells and keeps a ghost layer of copies of neighbours' boundary cells. Before every SpMV, Halo exchange fills the ghost values with point-to-point messages. Small messages are sent eagerly and large ones by a rendezvous handshake (eager-rendezvous-protocol, no page yet), and on high-speed fabrics the transfer is RDMA (rdma, no page yet) so the NIC writes straight into the destination buffer (Zero-copy I/O). Every dot product in CG or GMRES needs a global sum, an allreduce from Collective communication; it is latency-bound, which is why at high core counts Krylov methods stop scaling before the stencil work does.

Data and control flow

Flow 1 — one SIMPLE-type iteration of a steady flow, per process

1. Assemble the momentum equation: loop over faces, compute fluxes with the limited scheme, add to the CSR matrix and source; under-relax. Memory-bound, one pass over the mesh.
2. Solve momentum approximately (a few smoother or Krylov iterations; accuracy is not needed yet). Each iteration costs one SpMV plus a halo exchange.
3. Assemble the pressure equation from the momentum diagonal and the face fluxes.
4. Solve pressure to a loose tolerance with multigrid-preconditioned CG. Typically the largest share of time.
5. Correct face fluxes and cell velocities with the pressure gradient; update turbulence and scalar equations the same way.
6. Compute residuals (an allreduce), stop if they fall below tolerance, otherwise repeat.

Flow 2 — finite-element solve, from weak form to field

1. Read or generate the mesh, number the degrees of freedom, build the sparsity pattern (the CSR structure) once.
2. Loop over elements: gather the nodal data, evaluate basis functions at quadrature points, form the element matrix and vector, scatter-add into the global system.
3. Apply boundary conditions (modify rows and the right-hand side).
4. If nonlinear, form the Jacobian and residual, solve for the Newton update, line search, repeat until the residual norm is small.
5. Solve each linear system with the chosen Krylov method and preconditioner, or a direct factorisation, then write the field for post-processing.

Flow 3 — start of a parallel run

1. Partition the mesh graph; renumber the cells within each part for locality.
2. Distribute cells, build ghost layers and the send/receive lists used by the halo exchange.
3. Run the solve loops; time goes to SpMV bandwidth, the halo exchange latency and allreduces, in a ratio that shifts toward communication as parts get smaller.

Scaling

In the order the limits appear:

1. Memory bandwidth per node. SpMV and assembly are bandwidth-bound; renumber for locality, use blocked or lower-precision storage where accuracy allows, and put more nodes on the problem.
2. Iteration count. Unpreconditioned CG needs more iterations as the mesh refines; multigrid keeps it nearly flat, so a good preconditioner is worth more than faster hardware.
3. Communication surface. Surface-to-volume ratio worsens as parts shrink; use better partitions, overlap halo exchange with interior computation, and avoid over-decomposing.
4. Global reductions. Each dot product is an allreduce whose latency grows with the process count; use variants that fuse or defer reductions, and reduce Krylov iterations with a stronger preconditioner.
5. Memory of direct solvers and GMRES. Switch to multigrid or restarted methods for the largest runs.
6. I/O and meshing. Writing fields and generating huge meshes serially become the serial fraction; use parallel I/O and partition the mesh during generation.

Failure modes

Real implementations

Sources

• Weller, Tabor, Jasak, Fureby, A tensorial approach to computational continuum mechanics using object-oriented techniques, Computers in Physics, 1998. paper
• Chorin, Numerical solution of the Navier-Stokes equations, Mathematics of Computation, 1968. paper
• Patankar, Spalding, A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows, Int. J. Heat Mass Transfer, 1972. paper
• Sweby, High resolution schemes using flux limiters for hyperbolic conservation laws, SIAM J. Numer. Anal., 1984. paper
• Saad, Schultz, GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems, SIAM J. Sci. Stat. Comput., 1986. paper
• Karypis, Kumar, A fast and high quality multilevel scheme for partitioning irregular graphs, SIAM J. Sci. Comput., 1998. paper
• Shewchuk, Adaptive precision floating-point arithmetic and fast robust geometric predicates, Discrete and Computational Geometry, 1997. paper
• Dembo, Eisenstat, Steihaug, Inexact Newton methods, SIAM J. Numer. Anal., 1982. paper
• OpenFOAM Foundation, OpenFOAM User Guide (v11), Mesh description and fvSolution. docs · docs
• PETSc manual. docs
• DOLFINx documentation. docs · Dokken, The FEniCSx tutorial. tutorial