Brittle Numerics Root Cause

Brittle Numerics Root Cause

Core idea

A calculation whose result changes under compiler flags, hardware, or optimization level is not "flaky" the way a race condition is. It's almost always exposing an existing ill-conditioned spot in the algorithm — a near-singular linear solve, a slowly- (algebraically, not exponentially-) converging series truncated at an arbitrary cutoff, or a subtraction of two large near-equal numbers — where the compiler's FP rounding choices tip the balance. Genuinely fixing it means finding and repairing that ill-conditioned spot, not suppressing where it happens to surface.

The near-singular-solve case also surfaces through solver internals, not just compiler/hardware choices: PETSc can select a different default preconditioner depending on run configuration (e.g. ILU in serial vs. block-Jacobi-with-ILU-on-the-blocks once the run spans more than one rank), and if the linear system is ill-conditioned enough, different PCs that each satisfy the same convergence tolerance can still land on distinguishably different solutions — because a converged residual only bounds the solution error by roughly condition number x solver tolerance. -ksp_view (or -snes_view) on the passing and failing runs confirms whether a PC/solver-type change is what's different before chasing per-iteration instrumentation.

Anti-patterns (don't reach for these first)

• Pinning or loosening compiler flags to make the test pass again.
• Regolding the test to whatever the new answer happens to be, without understanding why it moved.
• Loosening the test's tolerance until the diff disappears.

Any of these can be a legitimate final step, but only after you understand why the value is unstable — otherwise you're just relocating where the instability bites next time.

Procedure

1. Reproduce deterministically first. Confirm the failure is deterministic given fixed inputs/flags (same build + same run always gives the same wrong answer). This rules out actual nondeterminism (uninitialized memory, races, order-dependent floating sums) and confirms you're chasing numerical conditioning, not a different bug class.
2. Instrument, don't guess. Add temporary per-iteration/per-mode debug prints of the internal quantities actually accumulating error (partial sums, determinants, convergence-criterion values), or step through with a debugger/conditional breakpoints when the quantity of interest is easier to inspect live than to log — either is "instrumenting," pick whichever gets you real numbers fastest. If logging, write one small reusable parser for the log format rather than a fresh one-off regex per question — you'll ask several follow-up questions of the same data.
3. Find the mechanism from the data, not from intuition. Let the instrumented values tell you the shape of the problem instead of assuming it. For example: measure the decay rate of a truncated series' terms empirically rather than assuming it's exponential; check whether a "near-singular" determinant is really near zero or just the result of large-magnitude cancellation. The specific check depends on the algorithm — the point is to derive it from what the numbers actually do, not from what's typical. For a linear or nonlinear solve where different builds or PCs converge to different answers, check whether the discrepancy is within roughly condition number x solver tolerance — if so, the fix is tightening a tolerance, prioritized by problem type. For a nonlinear problem, the Executioner's nl_rel_tol/ nl_abs_tol is the tolerance that ultimately governs the final answer's robustness: because Newton re-linearizes and self-corrects every iteration, a sufficiently tight nonlinear tolerance drives the accepted answer to that residual regardless of which PC/linear-solver path got there, as long as the linear solves converge well enough for Newton to keep making progress. l_tol/ l_abs_tol mainly governs whether Newton converges at all, since a too-loose linear tolerance can stall it before it reaches a tight nl_rel_tol. Tightening l_tol directly, or replacing the linear solve with a direct factorization (e.g. -pc_type lu -pc_factor_mat_solver_type mumps), are equally effective, fast ways to remove that risk — the direct factorization just avoids having to search for a tolerance value that's tight enough. If the divergence persists under a tight nl_rel_tol and a direct linear factorization, that's a sign the matrix is severely enough conditioned that roundoff itself is the exposed mechanism — the main procedure's case. Only fall back to the tolerance-widening path below if the near-singularity itself is physically expected and tightening the solve isn't warranted.
4. Look for precedent in sibling code. If a related algorithm/branch in the same codebase already handles this class of problem (e.g., a fallback for non-convergence), extend that pattern rather than inventing a new one.
5. Validate the fix against independent data — and treat "brute force" skeptically too. Cross-check a proposed closed-form correction (tail sum, limit, etc.) against a much less-truncated run of the same algorithm (e.g., rerun once, offline, with an absurdly large iteration cap). That run is not automatically ground truth: it's still inside the same slow-converging/cancellation-prone regime, just further along it, so it can carry its own residual error. Treat it as a second, independently-derived estimate to check trend and order-of-magnitude agreement against — not an authoritative reference — and corroborate with another angle (e.g., does the correction's sign and shape match physical expectation?) before trusting either number alone.
6. Test whether "more correct" complexity actually pays for itself. An enhancement (e.g. an early-exit heuristic) can look principled but not help the actual pathological case if it's gated on a coupled quantity that converges slower than the one you're optimizing for. Measure it on the failing case specifically before keeping it.
7. Anticipate blast radius before regolding. A fix to shared numerical code changes results for every consumer of that code, not just the test that was failing. Before mass-regolding: (a) enumerate every test exercising the changed code path, (b) sample across parameter regimes that stress the fix differently (e.g. a much smaller iteration cap, a much faster-converging case) to confirm it degrades gracefully rather than blowing up outside the regime it was tuned against, (c) only then regold in bulk, checking the diffs are smooth/deterministic/physically sensible rather than erratic.
8. Separate "already-known artifact in the gold" from "new bug." If something ugly turns up mid-investigation (e.g. a value that should be physically non-negative but isn't), diff against old logs/gold to check whether it predates your change before treating it as something you caused.

When targeted tolerance widening is legitimate, not an anti-pattern

Not every diff that varies with build/flags/hardware/thread count is an algorithmic conditioning bug. A quantity computed as the residual of two much larger, comparable-magnitude terms can be genuinely near a physical zero for reasons unrelated to any conditioning problem in the production code — e.g. a test deliberately suppresses the process that quantity measures, a physics residual that's genuinely near zero at a converged/equilibrium state, or it's a release/growth quantity still ramping up from zero in the first few steps. There, ordinary floating-point rounding noise — from any source: compiler codegen, FMA/vectorization, thread or MPI reduction order — dominates the relative error simply because the true value is tiny, while every other reported quantity in the file stays tight. That's a legitimate case for widening only that column's tolerance (e.g. MOOSE's override_columns/override_rel_err for CSVDiff, or custom_cmp for Exodiff — see below, and run-tests), with a comment recording the physical reason and the measured noise floor — not for a blanket rel_err bump on the whole test. Confirm before doing this that (a) the widened floor stays far enough below any physically meaningful violation that the test would still catch one — even when the near-zero column is itself what the test is validating (e.g. a mass-balance residual asserted to be ~0) — and (b) the column is near-zero for a documented physical/test-setup reason rather than because an upstream calculation is ill-conditioned — otherwise you're back in this skill's main case.

A related but distinct sub-case is a solver-telemetry column — a linear/nonlinear iteration count (and any cumulative postprocessor built on one), or a Picard iteration count — rather than a physical quantity. A non-smooth solve (e.g. frictional/complementarity contact) can take one more or fewer Newton or linear iterations on a different platform/build without the converged solution changing at all. Before reaching for a tolerance, though, ask whether the column belongs in the comparison at all: unless the test's own purpose is to validate solver performance (e.g. a dedicated iteration-count regression test), the telemetry postprocessor usually shouldn't be gold-compared in the first place, and dropping it from the comparison — rather than widening a tolerance to tolerate a quantity the test was never meant to pin down — is the simpler, more robust fix. Reserve the rest of this sub-case for tests that do intend to track solver performance. There, this isn't the "near-zero" case above — the count itself is O(10-100), not tiny — so a floor doesn't fit; widen that column's relative tolerance instead. Before doing so, confirm the actual solution variables (nodal variables, and any other global postprocessors) show zero or negligible diff on the same comparison — that's what tells you the extra iteration changed only the solver's path, not the answer.

A blanket file-wide rel_err bump introduced alongside an unrelated fix is a red flag worth git-archaeology (git log/git show on the test spec) even when it isn't currently failing: it often means this exact situation was mishandled by loosening everything instead of the one or two affected columns.

Exodiff: per-variable floor/relative error via custom_cmp

MOOSE's Exodiff tester has no override_columns equivalent; instead point custom_cmp at a command file passed to exodiff -f. Build it, don't hand-write it:

• exodiff -summary gold/<file>.e prints every variable with its file-wide peak magnitude, already formatted as NODAL VARIABLES/ELEMENT VARIABLES blocks — use this as the starting template (it also tells you each column's peak, which the sizing rule below needs).
• Give every variable-type block (all) (GLOBAL VARIABLES, NODAL VARIABLES, ELEMENT VARIABLES alike — exodiff's Parse_Variables in ED_SystemInterface.C parses the (all) flag and the indented per-variable override list independently of each other and of block type) with no explicit relative/floor on the header line. With no tolerance of its own, (all) falls through to the -F <abs_zero> -t <rel_err> the tester already passes on the command line from the test spec — so the rest of the file's comparison stays exactly as strict as before, and it stays correct if the spec's abs_zero/rel_err change later. Because the two are independent, this holds even for a block containing some affected variables: list only the affected variable's override line(s) under (all) and leave every unaffected variable out entirely, rather than enumerating them by hand just because they share a block with an affected one. (all) also keeps the file resilient to new output variables being added later.
• Only add an explicit floor (near-zero-residual case) or relative (solver-telemetry case, and only when the test's purpose is to track solver performance — otherwise omit the telemetry postprocessor from the compared variables instead) on the specific affected variable line(s), with a comment stating the physical reason. For a floor, keep -summary's auto-generated trailing # min: ... max: ... comment on that line (and on every other line, for consistency) rather than dropping it — it's tool-sourced evidence backing the floor, not hand-typed, so it corroborates the hand-written physical-reason comment instead of just asserting it. For a solver-telemetry relative override, look for an existing precedent value in a sibling .cmp before picking a number.
• Size that floor as roughly 1e-8 × that variable's own peak magnitude (from the retained -summary comment) — a principled, reproducible "many orders of magnitude below the values that matter" rule, rather than hand-tuning to the exact noise observed on one build. Then sanity-check it against the actually-measured noise from step 2's instrumentation: it should clear the worst observed noise by at least ~10x, and dumping the raw variable's full time/node distribution (e.g. via scipy.io.netcdf_file, since Exodus is classic netCDF) to confirm there's no real, physically-meaningful value of that same variable sitting just below the chosen floor.
• Look for a sibling .cmp file already in the test suite before inventing the format — e.g. a neighboring test's custom_cmp file that already floors a different near-zero column is the precedent to extend (per the Procedure's "look for precedent in sibling code" step above), not a new convention to invent.