PDE Methods — Math Reference

1. At a glance

Partial differential equations (PDEs) describe almost every continuous physical phenomenon: fluid flow, heat conduction, electromagnetic fields, elasticity and stress, quantum mechanics, mathematical finance, biological pattern formation, traffic flow, and image processing. A PDE relates an unknown function u(x, y, z, t, …) of several independent variables to its partial derivatives.

Three core families by character govern most of classical mathematical physics:

• Elliptic — steady-state equilibrium problems governed by the Laplacian; the canonical example is the Poisson equation ∇²u = f. Solutions are globally smooth — information propagates instantaneously across the whole domain.
• Parabolic — diffusion-like processes; the canonical example is the heat equation u_t = α·∇²u. Solutions smooth sharp data over time; information propagates at infinite speed but with exponential damping.
• Hyperbolic — wave-like processes; the canonical example is the wave equation u_tt = c²·∇²u. Solutions preserve singularities and propagate them at finite speed c along characteristics.

Most engineering simulation (CFD, FEA, EM, multiphysics) reduces to discretizing one of these classes — or a mixed system — and solving the resulting large sparse linear or nonlinear algebraic problem.

2. Classification

Order

• First-order: transport, conservation laws (u_t + a·u_x = 0), Hamilton-Jacobi (u_t + H(∇u) = 0), Eikonal (|∇u| = 1/c).
• Second-order: the bulk of classical physics (Laplace, heat, wave, Navier-Stokes momentum, Maxwell after curl-curl manipulation, Schrödinger).
• Higher-order: biharmonic ∇⁴u = f (thin-plate bending, Stokes streamfunction), Kuramoto-Sivashinsky u_t + u·u_x + u_xx + u_xxxx = 0 (flame fronts, thin films).

Linearity

• Linear: coefficients depend only on (x, t), not u — e.g. heat equation with constant α.
• Semilinear: highest-order terms linear, lower-order nonlinear — e.g. reaction-diffusion u_t = D·∇²u + f(u).
• Quasilinear: highest-order terms linear in derivatives but coefficients depend on u and ∇u — e.g. Navier-Stokes (convective (u·∇)u).
• Fully nonlinear: highest-order term itself nonlinear in u — e.g. Monge-Ampère det(D²u) = f, Hamilton-Jacobi.

Discriminant of second-order PDE

For a 2nd-order linear PDE in two variables A·u_xx + B·u_xy + C·u_yy + D·u_x + E·u_y + F·u + G = 0, the discriminant B² − 4AC classifies it:

• B² − 4AC < 0 → elliptic (Laplace u_xx + u_yy = 0: A=1, B=0, C=1, discriminant = −4).
• B² − 4AC = 0 → parabolic (heat u_t = u_xx: treated as u_xx − u_t = 0; A=1, B=0, C=0, discriminant = 0).
• B² − 4AC > 0 → hyperbolic (wave u_tt − c²·u_xx = 0: A=−c², B=0, C=1; discriminant = 4c² > 0).

The character can change spatially (transonic flow: subsonic = elliptic, supersonic = hyperbolic, transition at sonic line).

Boundary and initial conditions

A PDE alone is under-determined — boundary/initial conditions select a unique solution. Standard types:

• Dirichlet: u = g on ∂Ω — temperature held at known value, displacement clamped.
• Neumann: ∂u/∂n = h on ∂Ω — heat flux prescribed, traction-free surface, insulated wall.
• Robin (mixed): α·u + β·∂u/∂n = γ — convective heat-transfer boundary (Newton cooling), absorbing radiation BC.
• Periodic: u(x + L) = u(x) — turbulence DNS in a box, crystallography.
• Initial conditions (parabolic/hyperbolic): u(x, 0) = u₀(x); hyperbolic also needs u_t(x, 0) = v₀(x).

Well-posedness in the Hadamard sense requires existence, uniqueness, and continuous dependence on data. Mismatched BC types (Dirichlet imposed where Neumann is physically correct) are a leading source of garbage numerical solutions.

3. Classical PDEs

Poisson / Laplace

∇²u = f (Poisson, source f) or ∇²u = 0 (Laplace, homogeneous). Models steady-state heat conduction with sources, electrostatics (Gauss’s law in potential form), incompressible inviscid potential flow, gravitational potential outside masses, mesh smoothing (Laplacian smoothing), image inpainting (Perona-Malik), minimal surfaces (linearized). Green’s function G = 1/(4πr) in 3D — solutions integrate sources against this kernel.

Heat / diffusion

u_t = α·∇²u with diffusivity α = k/(ρ·c_p). Models conductive heat transfer, mass diffusion (Fick’s second law), Black-Scholes option pricing after log-transformation, image smoothing (Gaussian blur is exact heat-equation evolution), Brownian motion (Kolmogorov forward equation for the density).

Wave

u_tt = c²·∇²u with wave speed c. Models acoustics (linearized fluid), elastic waves (P and S), EM in vacuum (each component of E and B), seismic waves, vibrating strings/membranes/drums, structural vibration. d’Alembert solution in 1D: u(x, t) = f(x − c·t) + g(x + c·t).

Convection-diffusion

u_t + v·∇u = D·∇²u + S. Combines transport at velocity v with diffusion D. The Péclet number Pe = v·L/D measures the ratio: high Pe → advection-dominated, low Pe → diffusion-dominated. Numerical stability requires upwinding when Pe is large.

Burgers’ equation

Inviscid: u_t + u·u_x = 0 — the simplest nonlinear hyperbolic conservation law; develops shocks in finite time even from smooth data (characteristics cross). Viscous: u_t + u·u_x = ν·u_xx. Used as the toy model for turbulence and shock-capturing scheme development.

Navier-Stokes (incompressible)

Momentum: ρ·(u_t + (u·∇)u) = −∇p + μ·∇²u + ρ·g. Continuity: ∇·u = 0.

Quasilinear, parabolic in viscous terms, hyperbolic in inertia. Covered in [[Engineering/cfd-deep]].

Maxwell’s equations

∇·E = ρ/ε₀, ∇·B = 0, ∇ × E = −∂B/∂t, ∇ × B = μ₀·J + μ₀·ε₀·∂E/∂t.

In source-free regions, each component of E and B satisfies the vector wave equation. FDTD (Yee 1966) is the dominant time-domain method; frequency-domain via FEM (HFSS, COMSOL RF) or MoM.

Schrödinger

i·ℏ·ψ_t = −(ℏ²/2m)·∇²ψ + V(x)·ψ. Linear PDE — but the unknown is complex-valued. Time-independent eigenvalue form H·ψ = E·ψ is the workhorse of quantum chemistry (Hartree-Fock, DFT).

Black-Scholes

V_t + ½·σ²·S²·V_SS + r·S·V_S − r·V = 0. Parabolic, backward in time (terminal payoff condition). Reduces to heat equation by change of variables. Analytical solution for vanilla European options; FD or Monte Carlo for American / exotic.

Reaction-diffusion

u_t = D·∇²u + R(u) — Fisher-KPP, Allen-Cahn, FitzHugh-Nagumo (cardiac), Gray-Scott (Turing patterns), Belousov-Zhabotinsky.

4. Method of Lines (MOL)

The most flexible and widely used modern approach: discretize space only, leaving time continuous. This converts the PDE into a (large) system of ODEs in time:

du/dt = L_h(u)        — semi-discrete form after spatial discretization

Then integrate du/dt = L_h(u) using any standard ODE method (RK4, BDF, Crank-Nicolson, IMEX) covered in [[Math/ode-numerical-methods]]. Most production CFD, FEA-transient, and EM-FDTD codes are MOL under the hood.

Advantages: spatial and temporal accuracy decoupled; easy to swap time integrators; stiff ODE machinery (CVODE, IDA, BDF) applies directly; adaptive time-stepping straightforward.

5. Finite Difference (FD)

Replace partial derivatives with difference quotients evaluated on a regular grid.

Common stencils

• Forward: u_x ≈ (u_{i+1} − u_i)/Δx — O(Δx), first-order.
• Backward: u_x ≈ (u_i − u_{i−1})/Δx — O(Δx).
• Central: u_x ≈ (u_{i+1} − u_{i−1})/(2·Δx) — O(Δx²).
• Second derivative central: u_xx ≈ (u_{i+1} − 2·u_i + u_{i−1})/Δx² — O(Δx²).
• Five-point Laplacian (2D): ∇²u ≈ (u_E + u_W + u_N + u_S − 4·u_C)/Δx².
• Nine-point compact Laplacian — O(Δx⁴) on regular grids.

Upwinding

For convection-dominated problems u_t + a·u_x = 0, central differences are unconditionally unstable when explicit. Upwind (one-sided in the direction information flows) is stable but adds numerical diffusion. High-resolution schemes (flux-corrected transport / FCT, MUSCL, ENO, WENO) maintain second-order or higher accuracy in smooth regions and degrade gracefully near shocks (Harten 1983 TVD, Shu-Osher 1988 ENO/WENO).

Pros and cons

Pros: trivial to implement; high-order stencils on uniform grids easy; very fast on structured grids; matches geophysics/atmospheric models well. Cons: regular grids only — complex geometry requires immersed boundary, cut cells, or coordinate transforms; boundary conditions awkward on curved surfaces; not naturally conservative; weak for unstructured domains.

6. Finite Element (FE)

The Galerkin method projects the PDE onto a finite-dimensional subspace via the weak form. Multiply the PDE by a test function v, integrate by parts to shift one derivative onto v, and require the residual to be orthogonal to a chosen basis.

Weak form (Poisson example)

Strong: −∇²u = f in Ω, u = 0 on ∂Ω. Weak: find u ∈ H¹₀(Ω) such that ∫_Ω ∇u·∇v dx = ∫_Ω f·v dx for all v ∈ H¹₀(Ω).

Approximating u_h = Σ U_j·φ_j(x) with basis functions {φ_j} (typically piecewise polynomials on a mesh of tetrahedra, hexahedra, triangles, or quadrilaterals) yields a sparse linear system K·U = F with stiffness matrix K_{ij} = ∫ ∇φ_i·∇φ_j dx.

Refinement strategies

• h-refinement: subdivide elements (smaller mesh).
• p-refinement: raise polynomial order on existing elements.
• hp-refinement: combine — exponential convergence on smooth regions, locally h-refine near singularities.
• r-refinement: relocate mesh nodes (moving meshes).

Strengths and uses

Excellent for complex geometry, anisotropic materials, mixed BCs, adaptive refinement, multiphysics coupling. Standard in structural mechanics, thermal, electromagnetics (Nédélec edge elements for H(curl)), acoustics. CFD historically favors FV but stabilized FE (SUPG, GLS) and DG are competitive.

Libraries

Open source: deal.II (C++, AMR-heavy, multigrid), FEniCS / dolfinx (Python, automated weak form), MFEM (C++, high-order, GPU), Firedrake (Python, code-generation), Nektar++ (spectral/hp), FreeFEM (DSL), scikit-fem (Python, lightweight). Commercial: Abaqus, ANSYS Mechanical, COMSOL Multiphysics, MSC Nastran, LS-DYNA (explicit dynamics + contact). See [[Engineering/fem-fea]].

7. Finite Volume (FV)

Discretize the integral / conservation form of the PDE over control volumes (cells). For a conservation law u_t + ∇·F(u) = S, integrating over cell V_i and applying the divergence theorem:

d/dt ∫_{V_i} u dV + ∮_{∂V_i} F·n dA = ∫_{V_i} S dV

becomes

dU_i/dt = −(1/|V_i|)·Σ_faces F_face·A_face + S_i

The flux F_face at each interior face is computed by a numerical flux function (Roe 1981, HLL, HLLC, Osher, Godunov, AUSM) using the states on either side. Flux at a shared face is the same for both adjacent cells — so the scheme is conservative by construction.

Strengths and uses

Native conservation makes FV the dominant choice for compressible CFD, reactive flow, multiphase, free-surface, astrophysics (FLASH, Athena), atmospheric (FV3 in NOAA’s GFS / GFDL). Works on arbitrary unstructured meshes (cell-centered or vertex-centered).

Codes

OpenFOAM (open), ANSYS Fluent, STAR-CCM+, NASA FUN3D, SU2, Code_Saturne (EDF), Basilisk (AMR free-surface). See [[Engineering/cfd-deep]].

8. Spectral methods

Global basis expansion: u_N(x) = Σ_{k=0}^{N} a_k·φ_k(x) where {φ_k} are global smooth functions chosen for the geometry:

• Fourier for periodic domains — φ_k = e^{i·k·x}.
• Chebyshev for bounded intervals — T_k(cos θ) = cos(k·θ).
• Legendre, Jacobi — for L² inner products and quadrature.

For smooth (analytic) solutions, spectral methods achieve exponential convergence in N — far faster than any algebraic order. For solutions with limited regularity, convergence drops to algebraic.

Pseudo-spectral collocation

Enforce the PDE pointwise at quadrature nodes (Chebyshev-Gauss-Lobatto, Gauss-Legendre). Compute derivatives via differentiation matrices D so that (D·u)_i ≈ u'(x_i). Or use FFT for Fourier-spectral: O(N log N) per derivative.

Galerkin spectral

Require residual orthogonal to the basis; integrate weak form by parts; yields full but well-conditioned dense matrices (preconditioned).

Applications

Turbulence direct numerical simulation (Kim-Moin-Moser channel flow), climate dynamical cores (spectral T-grid in atmospheric models, before FV3 replaced GFS), quantum chemistry (Gaussian basis is a chemistry-flavored spectral method), seismology (SPECFEM uses spectral elements).

Spectral elements

Patera 1984 combined high-order spectral basis within each element with finite-element-style geometric flexibility — best of both worlds for complex domains. Nektar++ and Deal.II support hp-spectral elements.

9. Discontinuous Galerkin (DG)

Reed-Hill 1973 (neutron transport); Cockburn-Shu 1989 (Runge-Kutta DG for conservation laws). DG uses polynomial bases element-by-element but does not enforce continuity across element boundaries — instead, faces are coupled via numerical fluxes (as in FV).

Advantages: arbitrarily high order, naturally conservative, parallel-friendly (compact stencil), excellent for hyperbolic systems, easy hp-adaptivity. Competitive with FV in CFD (NASA Discontinuous Galerkin codes, FLEXI, Nektar++), dominant in computational electromagnetics (DGTD) and seismology.

Disadvantages: more degrees of freedom per element (no continuity), trickier for elliptic operators (interior penalty, Bassi-Rebay, LDG schemes), stricter CFL for explicit time-stepping (Δt ∝ Δx/(p²·c)).

10. Boundary Element Method (BEM)

For problems with a fundamental solution (Green’s function) — Laplace, Helmholtz, Stokes, linear elasticity — only the boundary is discretized. Reduces dimensionality by one: 3D problem becomes 2D surface mesh.

Strengths: exact treatment of unbounded domains (acoustic scattering into open air, EM radiation, exterior aerodynamics potential flow); no volume mesh needed.

Weaknesses: produces dense matrices (O(N²) memory, O(N³) direct solve); requires fast multipole method (Greengard-Rokhlin 1987) or H-matrices to reach O(N log N); restricted to linear PDEs with known Green’s function.

Codes: BEM++, Bempp-cl (Python), commercial WAON (vibroacoustics), AcuSolve.

11. Meshfree / SPH / MPM

When meshes break down — large deformation, fracture, fragmentation, free-surface — particle methods replace the mesh with a cloud of points:

• Smoothed Particle Hydrodynamics (SPH) — Lucy 1977, Gingold-Monaghan 1977 (astrophysics); now used for free-surface fluids, splashes, dam breaks. Codes: GADGET-2 (astro), DualSPHysics, PySPH, Splishsplash.
• Material Point Method (MPM) — Sulsky 1994; particles carry state, background grid for momentum solve; excellent for snow, sand, soil (Disney’s Frozen snow); large-deformation solids.
• Reproducing Kernel Particle Method (RKPM), Element-Free Galerkin (EFG) — meshless variants of FE.
• Peridynamics (Silling 2000) — nonlocal continuum model — fracture without re-meshing.

12. Time integration for PDE

After spatial discretization (FD/FE/FV/DG/spectral), MOL gives du/dt = L_h(u). Choice of time integrator depends on stiffness, accuracy, and stability needs.

Explicit

Forward Euler, explicit RK4, SSP-RK3 (Shu-Osher, conservation-law friendly). No linear solve. Bounded by CFL conditions:

• Hyperbolic CFL: Δt ≤ C·Δx / max|c| — wave speed dictates step. C ≈ 1 for first-order, smaller for higher-order.
• Parabolic CFL: Δt ≤ C·Δx² / (2·D) — quadratic in Δx. Becomes lethal on fine grids: halving Δx requires 4× more time steps.

Implicit

Backward Euler (1st-order, L-stable), Crank-Nicolson (2nd-order, A-stable but oscillatory for stiff), BDF1-6 (variable-order, stiff-stable), TR-BDF2 (Bank 1985). Each step requires a linear solve (linear PDE) or nonlinear solve (Newton with Jacobian — large Jacobian-free Newton-Krylov via PETSc SNES).

Unconditionally stable — set Δt by accuracy, not stability. Standard in implicit CFD (steady-state, large Δt), reservoir simulation, semiconductor device modeling.

IMEX (Implicit-Explicit)

Treat stiff terms (diffusion, fast linear) implicitly and non-stiff (convection, nonlinear sources) explicitly. ARK3, IMEX-BDF, splitting methods. Used in atmospheric modeling (vertical implicit, horizontal explicit), reaction-diffusion (diffusion implicit, reaction explicit).

Symplectic / structure-preserving

For Hamiltonian PDEs (Schrödinger, wave, Maxwell): symplectic integrators conserve discrete energy; Yee staggered grid for Maxwell preserves ∇·B = 0 exactly. See [[Math/ode-numerical-methods]].

13. Linear solver for PDE-discretized systems

After implicit time-stepping or steady-state FE/FV, you face a large sparse linear system A·x = b with n up to 10⁹.

Direct solvers (sparse LU)

• UMFPACK (Davis) — unsymmetric multifrontal; SuiteSparse.
• Pardiso (Schenk) — Intel MKL ships it; multithreaded.
• MUMPS (Amestoy et al.) — parallel multifrontal; MPI-distributed.
• SuperLU, WSMP.

Pros: robust, predictable, no convergence issues. Cons: memory grows superlinearly with mesh; impractical beyond n ≈ 10⁶ in 3D (fill-in).

Iterative (Krylov) — see [[Math/numerical-linear-algebra]]

• CG — symmetric positive definite (elliptic Poisson, FEM stiffness).
• MINRES — symmetric indefinite (saddle-point Stokes).
• GMRES — general nonsymmetric (convection-dominated, Navier-Stokes).
• BiCGStab, IDR(s), QMR — nonsymmetric alternatives.

All require preconditioning to be competitive.

Multigrid

Optimal O(n) complexity for elliptic problems. Smooth on a fine grid (Jacobi / Gauss-Seidel reduce high-frequency error), restrict residual to coarse grid, solve coarse correction recursively, prolong back. V-cycle, W-cycle, F-cycle, FMG (full multigrid).

• Geometric multigrid (GMG) — needs explicit hierarchy of meshes.
• Algebraic multigrid (AMG) — builds hierarchy from matrix structure (Ruge-Stüben 1987, smoothed aggregation — Vaněk et al. 1996). Black-box for SPD operators. Hypre BoomerAMG, Trilinos ML, PETSc GAMG, NVIDIA AmgX (GPU).

Preconditioners

• Jacobi, block-Jacobi, SOR — cheap, weak.
• Incomplete LU (ILU(k), ILUT) — moderate cost, broadly effective.
• Incomplete Cholesky (IC) — SPD.
• AMG as preconditioner for CG on elliptic systems.
• Block preconditioners for saddle-point systems (Stokes, mixed FE): Schur complement, SIMPLE, PCD, LSC.
• Domain decomposition — Schwarz (additive, multiplicative, RAS), FETI-DP, BDDC — parallel-friendly.

14. Adaptive Mesh Refinement (AMR)

Refine the mesh where the error is large, coarsen where it is small. The loop is estimate → mark → refine → solve → repeat.

Error estimators

• Residual-based (Babuška-Rheinboldt 1978) — evaluate strong-form residual on each element.
• Recovery-based (Zienkiewicz-Zhu 1987) — compare element gradient to smoothed (recovered) gradient.
• Goal-oriented / dual-weighted residual (Becker-Rannacher 2001) — solve adjoint, weight residual by sensitivity to a target functional.

Codes / libraries

• p4est (Burstedde) — parallel forest of octrees; backbone for many AMR codes.
• deal.II AMR — built-in hierarchical refinement.
• MFEM — non-conforming AMR.
• BoxLib / AMReX — block-structured AMR (LBNL).
• Chombo, PARAMESH, FLASH-AMR.

15. Spectral PDE solvers (FFT-based)

For periodic BCs, the FFT diagonalizes the Laplacian: (−∇²)·e^{i·k·x} = |k|²·e^{i·k·x}. Solving Poisson in periodic 3D box becomes:

1. FFT f → f̂_k.
2. Divide by |k|² (set k=0 mode to zero — solvability condition).
3. Inverse FFT → u.

Complexity O(N log N). Used in:

• Turbulence DNS (channel flow Kim-Moin-Moser, isotropic turbulence pseudo-spectral codes by Rogallo, Pope).
• Plasma physics PIC + spectral field solve.
• Climate atmospheric dynamical cores (spectral transform with associated Legendre).
• Phase-field crystal growth, spinodal decomposition (Cahn-Hilliard).

Dealiasing required for nonlinear products (2/3 rule, padding).

16. PINN — Physics-Informed Neural Networks

Raissi-Perdikaris-Karniadakis 2019. Train a neural network u_θ(x, t) to satisfy the PDE residual at sampled collocation points plus the BC/IC data:

Loss = λ_PDE · Σ |L(u_θ)(x_i, t_i) − f|² + λ_BC · Σ |u_θ − g|²_∂Ω + λ_data · Σ |u_θ − u_obs|²

Derivatives computed by autograd (Jax, PyTorch). No mesh, dimension-free in principle.

Useful for: inverse problems (given sparse data, infer PDE coefficients); assimilation with sparse measurements; high-dimensional Hamilton-Jacobi-Bellman; geometry-flexible exploratory work.

Limitations: training cost dwarfs classical solvers on standard forward problems; stiff multiscale problems converge poorly without curriculum / domain decomposition (XPINN, FBPINN); accuracy rarely matches classical FEM/FV; loss-landscape pathologies. Best treated as a complement, not a replacement, for classical methods.

17. Operator learning and neural operators

Learn the solution operator G: a ↦ u mapping PDE coefficients / initial data to solutions — so a single trained model amortizes the cost of many PDE solves with different parameters.

• Fourier Neural Operator (FNO) — Li et al. 2020 — kernel integration via FFT in spectral space; resolution-invariant; demonstrated on Burgers, Darcy, Navier-Stokes.
• DeepONet — Lu et al. 2021 — branch-trunk architecture: branch encodes input function, trunk evaluates at query points.
• Geometric / graph operators — GINO, MeshGraphNets (Pfaff-Battaglia 2021) for unstructured meshes.

Used for: turbulence surrogate models (NVIDIA Modulus, FourCastNet weather), uncertainty quantification (10⁴ Monte Carlo PDE solves replaced by 10⁴ network evaluations), real-time optimal control.

18. Specific applications

• Structural FEM — linear elasticity, plasticity, contact, large deformation. See [[Engineering/fem-fea]], [[Engineering/structural-analysis]].
• CFD — incompressible/compressible Navier-Stokes, turbulence modeling (RANS, LES, DNS). See [[Engineering/cfd-deep]], [[Engineering/fluid-mechanics]].
• Heat transfer — conduction, convection (coupled to flow), radiation. See [[Engineering/heat-transfer]], [[Engineering/Tier3/heat-transfer-correlations]].
• Electromagnetics — FDTD (Yee 1966) for time-domain, MoM for radiation/scattering, FEM for cavities and complex media (HFSS, CST, COMSOL RF, Ansys Maxwell).
• Quantum — Kohn-Sham DFT for materials (VASP, Quantum ESPRESSO, ABINIT, CP2K, CASTEP); time-dependent DFT for excited states.
• Climate — atmospheric + oceanic GCMs (CESM, GFDL CM6, ICON-A, MPAS, NEMO, MOM6). FV3 (finite-volume cubed-sphere) is now standard in NOAA’s UFS.
• Finance — Black-Scholes solver: analytical for European vanillas; FD or Monte Carlo for American, barrier, Asian options; PDE methods scale poorly past ~3 underlyings.

19. Software ecosystem

Python

• FEniCS / dolfinx — automated weak forms via UFL; production FE.
• Firedrake — UFL with code-generation; tightly coupled to PETSc.
• scikit-fem — lightweight pure-Python FE.
• FiPy — FV for PDEs (NIST).
• devito — symbolic FD with autotuned C generation (seismic).
• SfePy — simple FE; PDE prototyping.

C++ / Fortran

• deal.II — large FE library, AMR, multigrid.
• MFEM — high-order FE, GPU, AMR.
• Nektar++ — spectral / hp FE.
• DUNE — modular grid abstraction.
• MOOSE (INL) — multiphysics framework on libMesh.
• libMesh — FE infrastructure.
• Trilinos (Sandia) — solver + discretization toolkit (Belos, Tpetra, Kokkos, MueLu).
• PETSc (ANL) — the de facto foundation: sparse matrices, Krylov solvers, AMG, SNES nonlinear, TS time-steppers, DMPlex unstructured.

Commercial

• ANSYS Mechanical / Fluent / Maxwell / HFSS / CFX.
• Abaqus / Simulia (Dassault).
• COMSOL Multiphysics — single environment for FE + multiphysics.
• Star-CCM+ (Siemens), CD-adapco.
• MSC Nastran / Marc / Dytran, LS-DYNA (Ansys), RADIOSS (Altair).

Open source CFD

• OpenFOAM (ESI / Foundation forks) — FV, very widely used.
• SU2 (Stanford → community) — FV/DG, adjoint-based aerodynamic optimization.
• Code_Saturne (EDF).
• Basilisk (formerly Gerris, Popinet) — adaptive quadtree/octree FV.

GPU

• AmgX (NVIDIA) — algebraic multigrid on GPU.
• Ginkgo — GPU-portable sparse linear algebra.
• Kokkos (Sandia) — performance-portable kernels.
• PETSc + GPU backends — Cuda, HIP, SYCL.

20. Pitfalls

• Forgetting CFL — explicit hyperbolic solver explodes the moment Δt > Δx/c. Cure: use an adaptive Δt tied to wave speed, or switch to implicit.
• Parabolic CFL trap — refining Δx by 2 makes the explicit time step 4× smaller; on fine grids explicit diffusion is murder.
• Mesh too coarse for boundary layer — Reynolds-averaged turbulence needs y⁺ ≈ 1 near walls; coarse y⁺ ≈ 100 mesh gives wrong shear, wrong drag, wrong heat flux.
• Wrong BC type — Dirichlet imposed where Neumann is physical (or vice versa) gives plausible-looking but qualitatively wrong solutions; double-check against energy balance.
• Numerical diffusion in first-order upwind — smears shocks and contact discontinuities; use TVD / MUSCL / WENO when accuracy matters near sharp features.
• Mass non-conservation in non-conservative form — solving u_t + a·u_x = 0 for transport in primitive form vs flux form can yield mass imbalance over long runs; use the conservative form unless you know what you’re doing.
• Ignoring conditioning of stiff systems — high-order FD on stretched grids, ill-conditioned mass matrices, near-singular saddle points: iterative solvers stall without proper preconditioning.
• Spectral pollution — high-frequency aliasing in pseudo-spectral nonlinear terms; needs 2/3 dealiasing.
• Locking in low-order elements — incompressible elasticity / plate bending with Q1 elements locks (zero displacement); use mixed or higher-order elements.
• Aliasing in DG without slope limiters — Gibbs oscillations near shocks; add Krivodonova-style limiter, WENO subcell, or artificial viscosity.

21. Cross-references

• [[Math/ode-numerical-methods]] — RK, BDF, stability for time integrators.
• [[Math/numerical-linear-algebra]] — sparse solvers, CG, GMRES, multigrid, preconditioning.
• [[Math/eigenvalue-problems]] — modal analysis, structural vibration, Schrödinger eigenstates.
• [[Engineering/fem-fea]] — structural FE practice, element libraries, contact.
• [[Engineering/cfd-deep]] — Navier-Stokes discretization, turbulence modeling, CFD codes.
• [[Engineering/heat-transfer]] — conduction / convection / radiation.
• [[Engineering/Tier3/heat-transfer-correlations]] — Nu-Re-Pr correlations.
• [[Engineering/fluid-mechanics]] — continuum, conservation laws, dimensionless groups.
• [[Engineering/structural-analysis]] — beam / plate / shell theory underlying FE.

22. Citations

• Strang, G. and Fix, G. An Analysis of the Finite Element Method. Prentice-Hall, 1973 (revised Wellesley-Cambridge 2008). Foundational FE convergence theory.
• LeVeque, R. Finite Volume Methods for Hyperbolic Problems. Cambridge UP, 2002. Standard FV reference.
• LeVeque, R. Finite Difference Methods for Ordinary and Partial Differential Equations. SIAM, 2007.
• Hughes, T. J. R. The Finite Element Method: Linear Static and Dynamic Finite Element Analysis. Dover, 2000.
• Quarteroni, A. and Valli, A. Numerical Approximation of Partial Differential Equations. Springer, 1994.
• Trefethen, L. N. Spectral Methods in MATLAB. SIAM, 2000.
• Trefethen, L. N. and Embree, M. Spectra and Pseudospectra. Princeton UP, 2005.
• Cockburn, B. and Shu, C.-W. “TVB Runge-Kutta Local Projection Discontinuous Galerkin Finite Element Method for Conservation Laws II.” Mathematics of Computation 1989.
• Reed, W. H. and Hill, T. R. “Triangular Mesh Methods for the Neutron Transport Equation.” LANL report, 1973.
• Patera, A. T. “A Spectral Element Method for Fluid Dynamics: Laminar Flow in a Channel Expansion.” Journal of Computational Physics 54 (1984): 468-488.
• Yee, K. S. “Numerical Solution of Initial Boundary Value Problems Involving Maxwell’s Equations in Isotropic Media.” IEEE Transactions on Antennas and Propagation 14 (1966): 302-307.
• Raissi, M., Perdikaris, P., and Karniadakis, G. E. “Physics-Informed Neural Networks.” Journal of Computational Physics 378 (2019): 686-707.
• Li, Z. et al. “Fourier Neural Operator for Parametric Partial Differential Equations.” ICLR 2021 (arXiv 2010.08895).
• Lu, L., Jin, P., Pang, G., Zhang, Z., Karniadakis, G. E. “Learning Nonlinear Operators via DeepONet.” Nature Machine Intelligence 3 (2021): 218-229.
• Ruge, J. W. and Stüben, K. “Algebraic Multigrid.” In Multigrid Methods, Frontiers in Applied Mathematics, SIAM 1987.
• Greengard, L. and Rokhlin, V. “A Fast Algorithm for Particle Simulations.” Journal of Computational Physics 73 (1987): 325-348.
• Roe, P. L. “Approximate Riemann Solvers, Parameter Vectors, and Difference Schemes.” Journal of Computational Physics 43 (1981): 357-372.
• Harten, A. “High Resolution Schemes for Hyperbolic Conservation Laws.” Journal of Computational Physics 49 (1983): 357-393.
• Sulsky, D., Chen, Z., and Schreyer, H. L. “A Particle Method for History-Dependent Materials.” Computer Methods in Applied Mechanics and Engineering 118 (1994): 179-196.
• Silling, S. A. “Reformulation of Elasticity Theory for Discontinuities and Long-Range Forces.” Journal of the Mechanics and Physics of Solids 48 (2000): 175-209.
• Zienkiewicz, O. C. and Zhu, J. Z. “A Simple Error Estimator and Adaptive Procedure for Practical Engineering Analysis.” International Journal for Numerical Methods in Engineering 24 (1987): 337-357.