Benchmarks and verification

MILET's verification story has three layers: fast smoke tests that pin the pointwise physics to analytic results, solver gates that require the iterative and reduced formulations to reproduce the direct/reference solution to tight tolerances, and community benchmarks (SolCx, Blankenbach, Glerum, Busse) with published reference values. Everything below lives in test/ and runs against the same source files as production (test/runtests.jl includes the identical module graph as main.jl).

# smoke + solver gates (minutes)
julia --project=. test/runtests.jl

# additionally run the short physics benchmarks
RUN_BENCHMARKS=1 julia --project=. test/runtests.jl

# convergence studies (longer)
julia --project=. test/validation/solcx_convergence.jl
julia --project=. test/validation/blankenbach_convergence.jl

# laptop-budget long bench (≈ 30–45 min)
julia --project=. test/long_bench.jl

SolCx — variable-viscosity Stokes convergence

test/validation/solcx_convergence.jl runs the SolCx benchmark (Duretz et al. 2011) — a sharp 10⁶ viscosity jump at x = 0.5 with the analytic body force $\rho g_y=\sin (\pi y)\cos (\pi x)$ and free-slip walls — at nx = 33, 65, 129 and estimates the convergence order of $|v{|}_{max}$ by Richardson extrapolation between successive refinements (h halves → an error ratio of 4 implies second order). The reference peak velocity is $|v{|}_{max}\approx 6.9\times {10}^{-9}$. The case definitions are cases/solcx.toml and the fixed-resolution series cases/solcx_fixed_n24…n128.toml; the short gated variant (test/benchmarks/run_solcx.jl, RUN_BENCHMARKS=1) asserts finiteness and brackets vmax. Because the viscosity jump is element-aligned for even cell counts, this benchmark exercises exactly the discontinuous-coefficient regime the Q2/Q1 Taylor–Hood pair (see Discretization) must handle.

Blankenbach 1A — thermal convection references

test/validation/blankenbach_convergence.jl runs Case 1A of Blankenbach et al. (1989) — isoviscous Ra = 10⁴ convection, nondimensional unit box — at nx = 41/61/81 to quasi-steady state and compares against the canonical steady-state references:

$$\mathrm{Nu}=4.884,V_{\mathrm{r}\mathrm{m}\mathrm{s}}=42.86.$$

Both diagnostics are FE-integral-accurate (nusselt and vrms in src/io/Diagnostics.jl): Nu is computed as $\mathrm{Nu}=⟨\partial T/\partial y{|}_{y=0}⟩H/\mathrm{\Delta }T$ (positive in MILET's y-down convention) over the surface boundary measure, and $V_{\mathrm{r}\mathrm{m}\mathrm{s}}=\sqrt{\int u\cdot ud\mathrm{\Omega }/|\mathrm{\Omega }|}$. The script prints the percent error per resolution. The benchmark validates the entire PIC heat-transport loop — marker_T_mode = "pic_increment" (markers carry and advect temperature, the FE step supplies the diffusion increment) plus the buoyancy_from_T Stokes coupling. The cases/blankenbach_sg_n41/65/97 series repeats the study with Gerya & Yuen (2003) subgrid diffusion enabled (do_subgrid_diffusion = true), the configuration used in the methods paper. The 3D analogue is cases/busse_1a.toml (Busse et al. 1993: Nu ≈ 3.539, V_rms ≈ 40.999).

Geometric multigrid gates

test/smoke/test_gmg.jl contains three gates for the GMG solver stack (src/solvers/GMG.jl, src/solvers/BlockMG.jl; see Geometric multigrid):

1. 2D V-cycle on the real Q2 velocity space with a 4-decade viscosity contrast (stiff block in a soft matrix): the Galerkin hierarchy ($A_{\ell +1}=P_{\ell }^T{A}_{\ell }{P}_{\ell }$) must contract the residual by ≥ 10⁴ within 8 V-cycles with a mean per-cycle rate < 0.5, and gmg_update! (new matrix values, same pattern — the Picard-iteration path) must keep converging.

2. 3D hierarchy: the same machinery on a 16³ Q2 vector space with a viscosity contrast — relative residual < 10⁻³ after 8 cycles. This validates transfers, Galerkin products, and Chebyshev smoothing in 3D.

3. End-to-end Stokes: GMG-preconditioned FGMRES must reproduce the MUMPS direct solution of a buoyancy-driven variable-viscosity flow to a relative velocity difference < 2×10⁻³ (measured agreement ≈ 4×10⁻⁴ — the gate reflects the backward-error limit of a fully dimensional system with η ~ 10²⁰ Pa·s, not discretization error). The test runs n_picard = 2 deliberately so the cached in-place reassembly (assemble_stokes_cached!), the MUMPS numeric-only refactorization, and gmg_update! — the paths production runs live on — are all exercised.

The runnable demos for the iterative path are cases/plume3d_cube_gmg.toml and cases/plume3d_slab_gmg.toml.

Two-phase flow — the φ = 0 gate

test/smoke/test_twophase.jl enforces the degeneracy requirement of the McKenzie three-field formulation (src/physics/TwoPhase.jl; see Two-phase flow): as $\varphi \to 0$ the permeability $k_{\varphi }\to$ tiny and $1/{\zeta }_{\varphi }\to$ tiny, so the fluid-pressure equation degenerates to $-\mathrm{\nabla }\cdot u=0$ and $p_c\to 0$ — the system must reproduce single-phase Stokes. The test solves the same falling-block problem with picard_stokes! and with solve_twophase! at φ = 0 and asserts

• relative velocity difference < 5×10⁻³, and

• relative difference of the total pressure $p_f+p_c$ (mean-free, since the split between the two multipliers is regularization-dependent at φ = 0) against the single-phase pressure < 5×10⁻³.

A companion testset checks the constitutive closures pointwise (permeability power law, compaction viscosity, melt weakening $\eta e^{-{\alpha }_{\varphi }\varphi }$, finite compaction length $\delta =\sqrt{k_{\varphi }{\zeta }_{\varphi }/\mu }$).

Grain-size evolution — analytic tests

test/smoke/test_grainsize.jl pins the grain-size module (src/materials/GrainSize.jl, design in doc/design/grain_size.md; see Grain-size evolution) to closed-form results:

• Growth-only analytic: with ${\mathrm{\Psi }}_{\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{l}}=0$ the ODE integrates exactly to $d(t{)}^p=d_0^p+G_{0}t{e}^{-E_g/RT}$; update_grain_size must match within 1 % (it does so in a single substep, because the growth sub-flow is integrated in closed form).

• Wattmeter steady state: integrating from 10× above and 0.1× below the Austin & Evans (2007) fixed point $d_{ss}=(c\gamma K/(p\lambda \mathrm{\Psi }){)}^{1/(p+1)}$ must converge to the same value within 5 % of the analytic $d_{ss}$ (and to each other within 10⁻⁶ — a unique fixed point).

• Composite-rheology identity: with diffusion creep enabled, the pre-clamp harmonic identity $1/{\eta }_{\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{p}}=1/{\eta }_{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}+1/{\eta }_{\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{l}}$ must hold to 10⁻¹⁰; diff_A = 0 must regress exactly to the single-mechanism path.

• Mechanism-partition consistency (viscous_mechanism_partition): stress continuity requires ${\dot{\epsilon }}_{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}+{\dot{\epsilon }}_{\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{l}}={\dot{\epsilon }}_{II}$ and ${\mathrm{\Psi }}_{\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{l}}=2{\eta }_v{\dot{\epsilon }}_{II}{\dot{\epsilon }}_{\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{l}}$, both to 10⁻¹⁰.

Softening and healing tests

test/smoke/test_softening.jl covers the strain-softening additions (design in doc/design/softening_healing.md; see Softening & healing):

• Healing decay: ${\epsilon }_p(t)={\epsilon }_p(0)e^{-t/\tau }$ halves after $t=\tau \ln 2$ at $T=T_{\mathrm{r}\mathrm{e}\mathrm{f}}$; the Arrhenius τ(T) must shrink for hotter and grow for colder material.

• Viscous strain-softening ramp: visc_soft_factor = 0.1 must leave fresh material at η₀, reach exactly 0.1 η₀ beyond strain_thresh_weak, and be monotonic in between; the default factor 1.0 must be a no-op.

• Configurable yield cap: at 3 GPa confining pressure the implied stress $2\eta {\dot{\epsilon }}_{II}$ must respect yield_stress_max = 1e8, and disabling the cap (≤ 0) must let the Drucker–Prager yield grow with depth.

• Yielded flag: compute_VEP_viscosity must report plasticity governing (the flag that feeds the per-cell yield fraction gating strain accumulation) for a low-cohesion/high-η rock, and never for a purely viscous one.

Other smoke tests

• test/smoke/test_rheology.jl — pointwise rheology invariants: constant Newtonian rock returns exactly η₀; Arrhenius ∂η/∂T < 0; Maxwell VEP gives Z ∈ (0, 1], χ + Z = 1, η_VEP ≤ η_viscous; strain weakening is monotonic. Plus density: ρ(T_ref) = ρ₀, thermal contraction, additive phase Δρ.

• test/smoke/test_mesh.jl, test/smoke/test_fespaces.jl — 2D/3D mesh construction, wall tagging, and Taylor–Hood space sanity.

• test/smoke/test_phase_melt.jl — depth-triggered phase change and the Katz-lite solidus/liquidus melt fraction.

Gated short benchmarks and the long bench

With RUN_BENCHMARKS=1, test/runtests.jl additionally runs short end-to-end models (test/benchmarks/): SolCx (analytic ρ/η), the falling block (50 steps), a 20-step plume, and a short Blankenbach 1A that logs Nu and V_rms against 4.884 / 42.86. test/long_bench.jl is the curated laptop-budget suite (SolCx at two resolutions, falling block, plume, Blankenbach with ~200 steps) used to track wall-clock and physics regressions together.

Related pages

• The case files under cases/ — the runnable TOML definitions behind each benchmark.

• Testing — how the test harness is organised.

• Linear solvers, Geometric multigrid — the solver architecture the GMG gates protect.