Tier 3 — COMSOL differential tests (hours, nightly)

Tier 3 — COMSOL differential tests (hours, nightly)

The oracle tier, available because the author retains a COMSOL licence.

Scope (SPECIFICATION.md §8.2.4 D6). The comparison that gates v0.7 is against the paper's published current–voltage relationships and in-pore averages (VAL-16, VAL-17), which need no licence. The field comparison below is kept, but not required. What a published-results comparison cannot do is localise: see the ablation paragraph below, and knowledge/04-clya-geometry-and-charge/06-validation-targets.md#65-correction-ablation-targets-t-sectransportcomparison_corrections, for the part of localisation it can recover.

For a frozen set of cases spanning the envelope, compare field-by-field against exported COMSOL solutions — φ, c_i, u, p sampled on a common probe grid — plus every scalar QoI.

Acceptance: < 1 % relative L² on fields, < 0.5 % on integrated QoIs.

Golden files stored as compressed arrays with the generating COMSOL model archived alongside.

Treat the oracle as a depreciating asset. Generate and archive the full reference set early, in the solver-core phase, so the project stays verifiable if licence access ever lapses.

Differential-test design. Because the corrections are independently switchable, the most informative comparisons are ablations: classical PNP-NS, then each correction family enabled alone, then all. A discrepancy that appears only when one correction is on localises the bug immediately. This is the strongest argument for making PNP-vs-ePNP a first-class configuration rather than a code branch.

Measured: the attribution ladder on the VER-11 benchmark pore [tested: tests/tier3/test_comsol_comparison.py::test_val01_full_ladder_against_a_self_golden]

Four configurations of one model against a self-golden generated by the last of them, on the 708-element VER-11 pore (2 nm radius, 13 nm membrane, 0.5 M, +50 mV, −0.05 C/m², validated ePNP-NS with the flow on). E_k is the r-weighted relative L² on a 157-point probe grid; Δ_resid = E₃ is zero by construction because rung 3 generated the golden. Four full ladder solves take 75 s on this mesh.

field	Δ_total (none, P2/P1)	Δ_transport (→ supg)	Δ_flow (→ reference)	Δ_pair (→ P1/P1)
potential	4.01e-2	−3.55e-2	−3.40e-3	−1.23e-3
c_Na+	2.35e-3	−1.38e-3	−9.20e-4	−4.25e-5
c_Cl-	7.27e-4	+1.27e-4	−8.12e-4	−4.21e-5
pressure	1.89e-1	+3.97e-2	−1.98e-1	−3.09e-2
velocity_r	8.53e-1	+2.71e-3	−7.69e-1	−8.71e-2
velocity_z	5.77e-2	−8.53e-3	−2.08e-2	−2.84e-2

Two things this says, and neither is visible in a three-rung ladder.

The transport stabilisation dominates the potential and the flow GLS dominates the velocity. Δ_transport carries 88 % of the potential's distance and 0.3 % of velocity_r's; Δ_flow carries 90 % of velocity_r's. Collapsing the two into one Δ_stab, as a three-rung ladder does, would report a single number that is 88 % transport for one field and 90 % flow operator for another — and attribute neither.

Δ_pair is not negligible for the flow. The element pair alone moves velocity_r by 8.7 % and the pressure by 3.1 %, against 0.004 % for c_Na+. A comparison against a P1/P1 reference that ran at Taylor-Hood would carry that as an unexplained residual.

These are machinery measurements on a coarse benchmark, not statements about the reference implementation: the golden is one of our own solutions and the report says golden_source: self.

The transposed table is silent, and non-square patches are the cheap gate

A COMSOL %Data block is one row per z sample and one value per r sample, which is RadialGrid's own [i_z, i_r] layout. The reader checks that every row has the same width, which catches a ragged file and not a cleanly transposed square one: the shape check passes and the field is wrong everywhere except on the diagonal r = z, where it is exactly right. Requiring n_r ≠ n_z on every probe patch turns that case into a shape error at ingest and costs nothing — no physically motivated patch needs equal sample counts.

The r-weighted norm cannot see the axis [tested: tests/tier1/test_comparison_norms.py::test_val01_axis_defect_is_invisible_to_the_weighted_norm] [tested: tests/tier1/test_comparison_norms.py::test_val01_the_same_defect_at_the_wall_is_visible_to_both]

On the checked-in ClyA probe grid, a unit perturbation confined to the axis_line patch — 1 610 of 6 205 points, all within 0.045 nm of the symmetry axis — gives rel_l2 = 0.51 and rel_L2_r = 4.0e-4. Three orders. The weight is w_p · r_p, and r_p is 0.005 nm there against 60 nm at the far edge of the reservoir patch. The axis is where the 1/r forms of §6.2 and the NUM-06 natural condition are fragile, so a Tier-3 comparison reporting the weighted norm alone would be least sensitive exactly where the implementation is most likely to be wrong. The same perturbation moved to the reservoir patch gives rel_L2_r > 0.1, which is what shows the contrast is the axis and not the patch's point count.

The telescoping identity catches a slip, not a stale golden

Δ_total + Δ_transport + Δ_flow + Δ_pair = E₀ + (E₁−E₀) + (E₂−E₁) + (E₃−E₂) = E₃ = Δ_resid holds for any four numbers. It is worth asserting — it fails if the deltas are formed or assigned wrongly — but it cannot detect a rung compared against a re-exported golden or an edited probe grid, which telescope exactly as cleanly. Those need the hashes compared directly, and the retained point count per field compared between rungs. The WP13 plan's decision row claimed the identity did that job; it does not.

Its round-off is not absolute. An E_k is a relative error and exceeds 10 wherever the golden is small, and the round-off of the telescoped sum grows with it: over 100k random draws an absolute 1e-14 failed 0 % for E_k ∈ [0, 2], 2.7 % for [10, 100] and 16.1 % for [100, 1000], aborting the report on correct input. The tolerance is therefore taken relative to max(1, |E_k|) (codebase review of 2026-09-28, CR-7). [tested: tests/tier1/test_attribution.py::test_val01_ladder_identity_holds_on_synthetic_errors] [tested: tests/tier1/test_knowledge_arithmetic.py::test_telescoping_identity_round_off_scales_with_error_magnitude]

A norm over no points is zero, which reads as agreement [tested: tests/tier1/test_comparison_norms.py::test_val01_a_norm_over_no_retained_point_aborts]

field_error over an all-false mask returned rel_L2_r = rel_l2 = max_abs_rel = 0.0, and a mask gate comparing two empty masks passes. A self-golden made by a probe grid whose masks came back empty — which _inside's catch-all made one NGSolve API change away — would pass the Tier-2 round trip on every field while comparing nothing. An empty mask is now refused (CR-5).

A %Grid table round-trips float64 exactly at %.17g [tested: tests/tier1/test_comsol_golden.py::test_val03_grid_table_round_trips_through_npz] [tested: tests/tier2/test_self_golden_round_trip.py::test_val01_self_golden_reproduces_its_own_solution]

Seventeen significant decimal digits round-trip an IEEE-754 double by definition, so a solution written to a %Grid table and read back through density.grid.read_grid reproduces its array bitwise, NaN included (f"{nan:.17g}" is "nan", which np.fromstring parses back). A self-golden round trip therefore asserts exactly 0, not a tolerance — any difference at all is a defect in the transport format rather than a rounding budget being spent.

---