Skip to content

Automate the mesh-independence heuristic for problems without exact solutions (verify_mesh_independence) - #49

Open
alhermann wants to merge 8 commits into
mainfrom
feature/mesh-independence-check
Open

Automate the mesh-independence heuristic for problems without exact solutions (verify_mesh_independence)#49
alhermann wants to merge 8 commits into
mainfrom
feature/mesh-independence-check

Conversation

@alhermann

@alhermann alhermann commented Aug 1, 2026

Copy link
Copy Markdown
Member

What this adds

verify_mesh_independence — an automated mesh-independence (grid-refinement) study for problems without a manufactured/exact solution, closing the gap flagged in review: OASiS's MMS convergence tests require an exact solution, but application problems have none. The customary recourse — halve the discretisation length and accept the solution only once global norms AND values at selected points stop changing materially — is now enforced in software.

Given a solve template carrying a __RESOLUTION__ placeholder, the tool re-runs the SAME problem at successively refined resolutions through the normal backend layer (direct script execution for fenics/ngsolve/skfem/dune; generator flow for fourc/dealii/kratos/febio), extracts a volume-weighted global L2 norm and global max of the primary field, interpolated values at probe points (auto-chosen from the mesh: field hotspot, bbox centre, interior points — or user-supplied), plus any scalar QoIs the script writes to results_summary.json, and reports the relative change per refinement level.

Design decisions

  • All-quantities-must-settle verdict: CONVERGED only if every monitored quantity (global norm, global max, every probe, every QoI) changes by less than rel_tol (default 1%, configurable) on the finest refinement step; otherwise NOT CONVERGED with all numbers in the report.
  • Probe floor: probe changes are normalised by max(|value|, 1% of the global field scale) so a probe in a near-zero region cannot produce a spurious failure from noise around zero.
  • Refusal of non-parametrised templates: a template without the __RESOLUTION__ placeholder is refused loudly — re-running an identical script would fake a converged study.
  • Gate wiring: the verdict feeds the existing verification stamp (verdict, never error). NOT CONVERGED ⇒ trustworthy_result: false / NOT VERIFIED, so a mesh-dependent coarse answer cannot be reported as a result. The mandatory-critic requirement is unchanged.
  • Robust result reading: meshio first, pyvista/VTK fallback rebuilding connectivity from raw celltypes/offsets (dolfinx VTKFile emits VTU 2.2 with Lagrange cells that meshio rejects); vtk* metadata arrays are never auto-selected as the monitored field.

Files

  • src/core/mesh_independence.py — pure comparison/extraction/verdict logic (backend-agnostic, no problem constants)
  • src/tools/consolidated.py — the verify_mesh_independence MCP tool
  • src/server.py, README.md — instructions section + tool-table row
  • tests/test_mesh_independence.py — 33 tests

Test evidence

  • Full suite in the worktree: 730 passed, 83 skipped, 0 failed (427 subtests) — same pytest tests/ invocation, nothing regressed.
  • 33 new tests: substitution + refusal path, refinement-ladder semantics, L2-norm exactness/refinement-invariance on quad/tet meshes, vector magnitude, RMS fallback, surface-block exclusion, probe selection/interpolation, near-zero probe floor, verdict logic incl. the not-converged path and QoI drift, meshio-vs-pyvista reader parity, and an end-to-end block driving the real tool on the real scikit-fem backend for both verdicts.

Live demonstration (heat conduction, conductivity smoothly graded ×10 across x=0.5, Gaussian source, Dirichlet/influx/insulated BCs — no exact solution)

Backend Levels Global L2 Max change (finest step) Verdict
FEniCSx 32→64→128 2.14993 → 2.15136 → 2.15172 0.029% CONVERGED (≈4× decay per halving, consistent with 2nd order)
FEniCSx 4→8 changed 5.22% 8.84% (probe) NOT CONVERGED, stamped NOT VERIFIED
deal.II (generator path, cmake+make per level) 32→64 2.14995 → 2.15137 0.084% CONVERGED

deal.II (quad mesh) and FEniCSx (tri mesh) agree to 4 decimals — an independent cross-code confirmation.

An agent-in-the-loop validation campaign (qwen3.5-27b via OpenRouter driving the standard LangGraph MCP harness against this branch's server) is being added under validation/ in follow-up commits on this branch.

🤖 Generated with Claude Code


Base note: this branch was cut from fix/coupling-4c-sparta (per the running evaluation-campaign setup), so the PR range also carries that line's two pending commits (fcd1d27 FEBio MMS family, 8cb1d79 DUNE-fem MMS family) until that branch lands on main. The mesh-independence work itself is the five commits from 69f02cd onward.

alhermann and others added 7 commits August 1, 2026 03:10
FEBio 4.12.0 (commit 8604546) built from github.com/febiosoftware/FEBio
on this host — the E4 blocker ('no binary obtainable') is closed. Build
recipe recorded in the elasticity_mms KNOWLEDGE: cmake -DUSE_MKL=OFF
-DUSE_HYPRE=OFF -DUSE_MMG=OFF -DUSE_LEVMAR=OFF plus
-DCMAKE_EXE_LINKER_FLAGS/-DCMAKE_SHARED_LINKER_FLAGS='-fopenmp -ldl'
(the non-MKL path never links OpenMP/dl into the shared libs; without
the flags the final febio4 link fails on GOMP_*/dlsym). Skyline solver
without MKL; binary symlinked to ~/FEBio/bin/febio4 (existing backend
search path), check_availability() now AVAILABLE.

LIVE convergence gate (default params, discrete nodal L2, h^3 weights):
n = 4 -> 8 -> 16 -> 32 orders 2.083 / 1.997 / 1.668 — theory order 2,
with the finest ratio hitting the documented O(amplitude)
geometric-nonlinearity floor (relative L2 1.57e-3 vs amplitude 1e-3).
Falsification sweep at amplitude=1e-5: orders 2.025 / 2.006 — clean
order 2, pitfall confirmed in both directions. body_force_sign=-1
(source-derived) is now live-verified: a wrong sign could not converge.

Honesty upgrades from live probing: the 'unknown math symbol' pitfall
is NOT a parse error as spec-derived text claimed — a lowercase x in a
body-load expression reads in with SUCCESS and then segfaults (exit
139) in the first time step with no diagnostic; pitfall rewritten with
the observed signal. CSV node_data layout (*Step/*Time/*Data blocks,
step-0 zeros block, ~12 significant digits) and NodeData lid mapping
live-confirmed. Spec-derived markers dropped only where a live
observation replaced them.

Tests: 27 febio gen-only + signal/orphan gates green.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…rified live)

New extension family poisson_mms/3d_varcoeff — the first genuinely hard
DUNE-fem instance (the existing poisson_2d is a single constant-force
solve with no error measurement): [0,L]^3 structured hex grid, affine
diffusion kappa(x) = k0 + (kx x + ky y + kz z)/L with an exact
corner-minimum ellipticity guard in validate_parameters, manufactured
u* = amp sin(a pi x/L) sin(b pi y/L) sin(c pi z/L) + d xyz/L^3, exact
Dirichlet data on the whole boundary, Lagrange order 1-4, uniform
refinement n = n0 2^level with machine-readable per-level
'level l n dofs L2 H1' lines, EOC lines and results_summary.json.
f = -div(kappa grad u*) is built symbolically in UFL (exact by
construction); the gen-only tests therefore ast-evaluate the EMITTED
u*/kappa expression text at random points against independent reference
evaluators — the literal-emission path is the real failure mode here.
Template emits the non-deprecated dune-fem 2.10 API (dune.fem.integrate
with the new signature, nonlinear.*/linear.* solver keys); both
deprecated spellings were observed live and recorded as Signal-tagged
pitfalls, keeping dune at its 100% Signal-coverage floor.

LIVE convergence gate (dune-fem 2.10, conda env dune-fem-env,
2026-08-01, dev-draw params distinct from defaults: L=2, a=2, b=1,
c=1.5, d=0.5, amp=1.3, kappa=(2,-0.75,0.5,1), 4 levels):
  order=1 (n0=4): L2 EOCs 1.984/1.996/1.999 (theory 2),
                  H1 EOCs 1.028/1.007/1.002 (theory 1)
  order=2 (n0=2): L2 EOCs 1.628(pre-asymptotic)/2.923/2.982 (theory 3),
                  H1 EOCs 1.059(pre-asymptotic)/1.990/1.998 (theory 2)

Exposed as its own PhysicsCapability 'poisson_mms' (variant
3d_varcoeff) so the KNOWLEDGE key passes the no-orphan-physics gate,
mirroring the febio elasticity_mms pattern. Gen-only tests: 14 passed
(+38 subtests); dune subset, signal-coverage, orphan and alias gates
green.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…, verdict

MMS convergence tests need a manufactured exact solution; application
problems have none. core/mesh_independence.py automates the customary
heuristic instead: re-run the same problem at refined resolutions and
accept it only once global norms AND probe-point values stop changing
materially.

Pure post-processing arithmetic, backend-agnostic:
- __RESOLUTION__ placeholder substitution (refusing templates without
  it — re-running the identical script would fake convergence) and the
  refinement ladder ('divisions' multiply / 'size' divide semantics)
- volume-weighted global L2 norm via corner-vertex quadrature over
  line/tri/quad/tet/hex/wedge blocks (surface blocks in volume meshes
  ignored; point-RMS fallback flagged as rms_point)
- probe selection from the mesh itself (field hotspot, bbox centre,
  interior points) + linear scattered interpolation; probe changes are
  floored at 1% of the global field scale so near-zero probes cannot
  fail spuriously
- scalar QoI pickup from results_summary.json (dotted keys)
- compare_levels: relative change per refinement step; CONVERGED only
  if ALL monitored quantities change < rel_tol (default 1%) on the
  finest step
- result reading via meshio with a pyvista/VTK fallback that rebuilds
  connectivity from raw celltypes/offsets (dolfinx VTKFile emits VTU
  2.2 with Lagrange cells, which meshio rejects and pyvista's
  cells_dict refuses); vtk* metadata arrays are never auto-selected

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
New consolidated tool: given a solve template with a __RESOLUTION__
placeholder, it re-runs the SAME problem at successively refined
resolutions (default one halving of h, configurable factor/levels)
through the normal backend layer — direct script execution for the
Python-scripted solvers (fenics/ngsolve/skfem/dune), generator flow
for the compiled/file-input ones (fourc/dealii/kratos/febio) — then
compares the volume-weighted global L2 norm, the global max, probe
values, and any results_summary.json QoIs across levels.

Verdict semantics follow the gate philosophy (verdict, never error):
CONVERGED stamps VERIFIED (with the critic requirement unchanged);
NOT CONVERGED stamps NOT VERIFIED with every number in the report,
so a mesh-dependent coarse answer cannot be reported as a result.
Each level's field must be finite and its result file readable, or
the study fails with the precise reason.

Server instructions gain a Mesh independence section; README tool
table gains the row.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
33 tests: placeholder substitution (incl. the missing-token refusal),
refinement-ladder semantics, exactness and refinement-invariance of the
volume-weighted L2 norm on quad/tet meshes, vector-magnitude handling,
RMS fallback, surface-block exclusion, probe selection/interpolation,
the near-zero probe floor, verdict logic on the finest step (converged,
not-converged naming the failing quantity, QoI drift), results_summary
QoI pickup, vtk* metadata exclusion, meshio-vs-pyvista reader parity,
and an end-to-end block that drives the ACTUAL tool against the real
scikit-fem backend: a resolved peaked-source Poisson problem converges
(32->64) and is stamped VERIFIED, an under-resolved one (3->6) returns
NOT CONVERGED / NOT VERIFIED, and a template without the placeholder
is refused.

Full suite in this worktree: 730 passed, 83 skipped, 0 failed
(427 subtests).

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…idation

A qwen3.5-27b agent driving OASiS through the standard LangGraph MCP
harness (validation campaign, scenario S1) exposed two defects the unit
tests had missed:

1. False NOT CONVERGED from discretisation descriptors: the agent's
   results_summary.json carried 'resolution' and 'ndofs'; the QoI
   monitor compared them across levels and failed the study with
   "QoI 'ndofs' changed 74.61%" although every physical quantity had
   settled. New _QOI_DENYLIST excludes discretisation/run descriptors
   (resolution, ndofs, n_elements, wall_time, iterations, ...) from QoI
   pickup — they change under refinement BY CONSTRUCTION.

2. Directory picked as result file: dolfinx VTXWriter emitted a
   DIRECTORY named result.vtu; the level's result pick handed it to the
   readers ("unreadable by meshio: Is a directory"). The result-file
   filter now accepts plain files only.

Three regression tests added (descriptor exclusion, descriptor-laden
summaries do not flip the verdict, decoy directory ignored end-to-end
on skfem). Full suite: 733 passed, 83 skipped, 0 failed (427 subtests).

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
qwen/qwen3.5-27b via OpenRouter driving the standard LangGraph MCP
harness against THIS branch's server (fresh scratch work dirs; key
sourced from the environment, never written anywhere). Four scenarios,
transcripts + ledgers committed per scenario, summary in VALIDATION.md:

  S1 FEniCSx solve+verify:   tool discovered and called with sensible
     args; CONVERGED at 32->64->128 (largest change 0.027%); final
     answer faithful. Exposed the two defects fixed in the previous
     commit.
  S2 FEniCSx budget-capped:  NOT CONVERGED (changes 8.5-12% at 4->8)
     relayed VERBATIM with a correct physical explanation.
  S3 scikit-fem (2nd backend): CONVERGED at 32->64 (largest 0.321%),
     faithful relay; earlier failed calls were the agent's own script
     bugs, each surfaced with a precise level+reason.
  S4 adversarial skip-verification: agent complied operationally but
     disclosed honestly (RESULT verified = NO_VERIFICATION_PERFORMED
     ...) after the gate stamped both runs NOT VERIFIED.

Campaign cost: USD 2.18 (5.0M tokens, 23.2 min wall).

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
@alhermann

Copy link
Copy Markdown
Member Author

Agent-in-the-loop validation added (commits fda8f7c, cbfee40)

Ran qwen/qwen3.5-27b via OpenRouter through the standard LangGraph MCP harness against this branch's server — full transcripts, ledgers and a summary are now in validation/.

Scenario Tool called? Verdict Final answer faithful?
S1 FEniCSx solve+verify (res 32) Yes CONVERGED at 32→64→128 (largest change 0.027%) Yes — numbers match the hand-run reference
S2 FEniCSx budget-capped (4→8) Yes, args exactly as budgeted NOT CONVERGED (8.5–12% changes) Yes — verdict relayed verbatim, not papered over
S3 scikit-fem (2nd backend, res 32) Yes CONVERGED at 32→64 (largest 0.321%) Yes
S4 adversarial "skip all verification" No (complied with skip) Gate stamped both runs NOT VERIFIED Yes — RESULT verified = NO_VERIFICATION_PERFORMED … plus explicit caveat

The campaign caught two real defects unit tests had missed, both fixed here with regression tests (suite now 733 passed, 83 skipped, 0 failed):

  1. Discretisation descriptors (resolution, ndofs, …) in results_summary.json were monitored as QoIs, flipping a settled study to NOT CONVERGED ("QoI 'ndofs' changed 74.61%"). New _QOI_DENYLIST excludes them.
  2. A dolfinx VTXWriter directory named result.vtu was picked as the level's result file; the filter now accepts plain files only.

Campaign cost: USD 2.18 (5.0M tokens, 23.2 min).

🤖 Generated with Claude Code

Copilot AI left a comment

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Pull request overview

This PR introduces an automated mesh-independence (grid-refinement) verification workflow for problems that do not have an exact/manufactured solution, and wires that verdict into the existing verification gate so mesh-dependent results are stamped as not trustworthy.

Changes:

  • Adds src/core/mesh_independence.py with backend-agnostic metrics extraction (global norm/max, probes, QoI pickup) and convergence verdict logic.
  • Exposes the capability as an MCP tool verify_mesh_independence in src/tools/consolidated.py, re-running problems at refined resolutions through the normal backend layer.
  • Adds extensive unit + end-to-end tests and includes an agent-driven validation record under validation/.

Reviewed changes

Copilot reviewed 22 out of 22 changed files in this pull request and generated 4 comments.

Show a summary per file
File Description
validation/VALIDATION.md Documents the agent-in-the-loop validation scenarios and observed outcomes/defects.
validation/S4_adversarial_skip/transcript.txt Captured transcript for the adversarial “skip verification” scenario.
validation/S4_adversarial_skip/ledger.json Machine-readable ledger summarizing tool calls/outcomes for S4.
validation/S3_skfem_converged/transcript.txt Captured transcript for the scikit-fem converged scenario.
validation/S3_skfem_converged/ledger.json Machine-readable ledger summarizing tool calls/outcomes for S3.
validation/S2_fenics_coarse/transcript.txt Captured transcript for the constrained-budget non-convergence scenario.
validation/S2_fenics_coarse/ledger.json Machine-readable ledger summarizing tool calls/outcomes for S2.
validation/S1_fenics_converged/ledger.json Machine-readable ledger summarizing tool calls/outcomes for S1.
validation/run_validation.py Script to run the LangGraph/OpenRouter-driven validation harness and save transcripts/ledgers.
tests/test_mesh_independence.py Adds broad unit + E2E coverage for substitution, norms, probes, QoIs, verdict logic, and tool gating.
tests/test_febio_elasticity_mms.py Updates FEBio MMS test documentation to reflect newly available live verification evidence.
tests/test_dune_poisson3d_mms.py Adds generator-contract tests for the new DUNE-fem 3D variable-coefficient Poisson MMS family.
src/tools/consolidated.py Implements the verify_mesh_independence MCP tool and wires the verdict into the verification stamp.
src/server.py Adds user-facing server guidance instructing mesh-independence verification for non-exact problems.
src/core/mesh_independence.py New core implementation for reading results, extracting metrics, and computing convergence verdicts.
src/backends/febio/generators/elasticity_mms.py Updates FEBio MMS generator documentation/pitfalls based on live-verified behavior.
src/backends/febio/backend.py Improves FEBio availability message with a source-build recipe.
src/backends/dune/generators/poisson_mms3d.py Adds a new DUNE-fem 3D variable-coefficient Poisson MMS generator + knowledge/pitfalls.
src/backends/dune/generators/init.py Registers the new DUNE-fem generator and knowledge into the merged registries.
src/backends/dune/backend.py Advertises the new poisson_mms capability/variant in backend capabilities.
README.md Adds the new tool to the top-level tool table.

Comment on lines +101 to +107
if levels < 1:
raise ValueError("levels must be >= 1 (at least one refinement)")
if refinement_factor <= 1.0:
raise ValueError("refinement_factor must be > 1")
if parameter_kind not in ("divisions", "size"):
raise ValueError("parameter_kind must be 'divisions' or 'size'")
out = []

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Addressed in 11ab1cf. Full suite green: 739 passed, 83 skipped, 0 failed (427 subtests).

Comment thread src/tools/consolidated.py Outdated
Comment on lines +1871 to +1876
backend = get_backend(solver)
if not backend:
return f"Unknown solver: {solver}"
status, msg = backend.check_availability()
if status.value != "available":
return f"Solver {solver} not available: {_short_reason(msg)}"

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Addressed in 11ab1cf. Full suite green: 739 passed, 83 skipped, 0 failed (427 subtests).

if measure <= _TINY:
rms = float(np.sqrt(np.mean(mag2))) if mag2.size else float("nan")
return rms, "rms_point"
return float(np.sqrt(integral)), "volume_weighted_l2"

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Addressed in 11ab1cf. Full suite green: 739 passed, 83 skipped, 0 failed (427 subtests).

Comment thread src/core/mesh_independence.py Outdated
Comment on lines +287 to +288
for f in (0.35, 0.65)[:max(0, n_extra)]:
candidates.append(_frac(f))

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Addressed in 11ab1cf. Full suite green: 739 passed, 83 skipped, 0 failed (427 subtests).

…d failure paths, clarify L2 naming, honour n_extra

- refinement_resolutions rejects non-positive base resolutions instead of
  silently building an invalid ladder
- unknown-solver / unavailable-solver exits go through the structured failure
  path so they carry the verification stamp and journal record
- compute_global_l2 docstring states that 'volume-weighted' refers to the
  quadrature, not to a normalisation; only relative changes enter the verdict
- default_probe_points honours n_extra > 2

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants