Verification and validation

Nothing in this library is checked only against itself. Four independent references are used, and each caught defects the others did not.

Against MFEM’s own integrators

Every block generated by automatic differentiation is compared against the hand-written MFEM integrator for the same bilinear form, on the same integration rule (set explicitly with SetIntRule, or the comparison measures the rule difference instead).

form	MFEM integrator	agreement
`inner(grad u, grad p)`	`DiffusionIntegrator`	2e-16
`u p`	`MassIntegrator`	2e-16
elasticity	`ElasticityIntegrator`	5.5e-16
`dot(u, p)`, H(curl)/H(div)	`VectorFEMassIntegrator`	1.3e-15
`curl u . curl p`	`CurlCurlIntegrator`	1.0e-15
`div u div p`	`DivDivIntegrator`	1.4e-15
boundary `u p`	`MassIntegrator` as a boundary integrator	2.6e-16
boundary `dot(grad u, n) p`	`BoundaryNormalLFIntegrator`	2.5e-16

on triangles, quadrilaterals, tetrahedra and hexahedra, P1 to P3, scalar and vector, including mixed-geometry meshes and the H(curl) case that needs a DofTransformation.

Blocks taken between two families are checked the same way: the RT mass and the L2 x RT divergence of a mixed Darcy density against VectorFEMassIntegrator and VectorFEDivergenceIntegrator (worst 1e-12), with the two off-diagonal blocks required to be each other’s negative transpose, and the diagonal blocks of a poroelasticity density against ElasticityIntegrator and a mass-plus-diffusion pair.

Against hIPPYlibx on a shared discrete problem

validation/run_validation.sh solves the same discrete problem in hIPPyMFEM and in hIPPYlibx (hIPPYlib on FEniCSx): the same mesh arrays, the same quadrature degree, the same observation points, and the same noise realization read from a file written by whichever driver runs first. The two run under different interpreters because they need different environments.

Twenty-seven quantities agree. The decisive one is the exact dense generalized spectrum of the misfit Hessian, evaluated at a fixed analytic point so that neither optimizer enters the comparison: 3.7e-13. Costs, gradients, Hessian quadratic forms, the MAP point and the MAP fields agree to 1e-13 to 1e-15.

Internal consistency

adjoint identities, and symmetry of every second-order block;
ReducedHessian against FDHessian;
modelVerify slopes of exactly 1 for the gradient and the Hessian;
BFGS reproducing the Newton-CG minimizer to 1e-5, from the cost and gradient alone;
randomized eigensolvers against a dense eigh (1e-15);
\(\sqrt{M}\sqrt{M}^{\top} = M\) exactly;
MCMC against an exact Gaussian posterior, and the autocorrelation time against an AR(1) process whose value is known analytically;
the full Hessian positive definite at the MAP point over several directions, though not everywhere, because away from a minimum it is legitimately indefinite.

Implementation equivalences

Where a faster path replaced a slower one, the slower one is kept as the reference and the two are required to agree exactly:

direct-CSR assembly against MFEM’s per-element callback route: exactly zero difference for every block, every geometry, every diagonal policy, on 1, 2 and 4 ranks;
GPU against CPU element kernels and assembled operators: 4.3e-16 and 3.6e-16;
the deterministic device scatter against the atomics-based one: 1e-16;
exact parallel solves against iterative ones on the same inverse problem.

Partition independence

Random vectors are bit-identical on 1, 2 and 4 ranks, because the random stream is indexed by global position rather than by rank. Prior samples, synthetic data and the MAP point pass through linear solves, so they agree across rank counts to the tolerance of those solves: to round-off with the exact solvers, and to 1e-14 for a prior sample whose Krylov solve is run to that tolerance. This is the check that catches ghost-dof and true-dof mistakes, which symmetry tests do not.

The test suites

./run_tests.sh            # every suite on one rank, a few minutes
./run_tests.sh 1 2 4      # and on two and four ranks
python -m pytest          # the same suites through pytest, on one and two ranks

suite	what it checks
`test_vectors`	vectors, operators, MultiVector, the parallel RNG; bit-identical random vectors on 1 to 4 ranks
`test_kernels`	every AD-generated block against MFEM’s own integrators on triangles, quadrilaterals, tetrahedra and hexahedra, P1 and P2, scalar and vector; finite-difference consistency of second and third derivatives
`test_solves`	linear and nonlinear forward solves, adjoint and incremental systems, the adjoint gradient against finite differences
`test_modeling`	prior sample covariance against \(R^{-1}\), the vector and mollified priors against MFEM’s forms, observation operators, the multiplicative-noise and multi-state misfits against their formulas, `modelVerify` slopes, `ReducedHessian` against a finite-difference Hessian
`test_optimization`	trust-region CG, randomized eigensolvers against dense `eigh`, BFGS reproducing the Newton-CG minimizer, steepest descent against a dense solve, the Laplace approximation
`test_timedependent`	time-dependent inversion: first-order gradient, an exact Hessian for the quadratic cost
`test_uq`	pCN, gpCN, MALA and importance sampling against an exact Gaussian posterior, autocorrelation times, QoI derivatives, Taylor moments, variance-reduced Monte Carlo
`test_assembly`	the direct-CSR assembly against MFEM’s callback route: every block, every geometry, exactly zero difference
`test_solvers`	exact parallel solves against dense ones, the same answer on 1, 2 and 4 ranks, the lumped-mass and transpose solvers, the bookkeeping of hypre’s device memory pool, and the PETSc bridge
`test_boundary`	boundary (`ds`) integrals against MFEM’s boundary integrators, a Robin problem, parameter-dependent boundary blocks
`test_facets`	interior-penalty DG against MFEM’s `DGDiffusionIntegrator`, faces shared between ranks, a DG inverse problem’s gradient and Hessian
`test_vectorfe`	H(curl) and H(div) blocks against MFEM, mixed Darcy and Biot poroelasticity blocks, an H(curl) inverse problem
`test_nb`	the plotting helpers, exact at every sample point on every rank count
`test_gpu`, `test_device`	the GPU path against the CPU path at round-off, and MFEM and hypre on the device with a CUDA or HIP build of PyMFEM