diff --git a/tests/integration/analysis/thermal_gauss_oracle.h b/tests/integration/analysis/thermal_gauss_oracle.h
new file mode 100644
index 0000000..e1822e9
--- /dev/null
+++ b/tests/integration/analysis/thermal_gauss_oracle.h
@@ -0,0 +1,70 @@
+// Generated by thermal_oracle_generator.py; exact root/weight certificates.
+#ifndef FESA_TESTS_THERMAL_GAUSS_ORACLE_H_
+#define FESA_TESTS_THERMAL_GAUSS_ORACLE_H_
+#include <array>
+namespace thermal_oracle {
+struct GaussPoint {
+  double x_lo, x_hi, w_lo, w_hi;
+};
+inline constexpr std::array<GaussPoint, 12> kGauss12{{
+    {-0x1.f68f1d8e42e82p-1, -0x1.f68f1d8e42e80p-1, 0x1.8275d9de72d45p-5,
+     0x1.8275d9dedaddep-5},
+    {-0x1.cee874ffb88b5p-1, -0x1.cee874ffb88b3p-1, 0x1.b60602bcd3d56p-4,
+     0x1.b60602bcf85a4p-4},
+    {-0x1.8a30aeed88f37p-1, -0x1.8a30aeed88f35p-1, 0x1.47d7258f1fdcdp-3,
+     0x1.47d7258f25d51p-3},
+    {-0x1.2cb4f05c077fap-1, -0x1.2cb4f05c077f8p-1, 0x1.a0163e6b1a279p-3,
+     0x1.a0163e6b1b46ap-3},
+    {-0x1.78a8d20a8b19ep-2, -0x1.78a8d20a8b19cp-2, 0x1.de3155c256a18p-3,
+     0x1.de3155c256b52p-3},
+    {-0x1.007a5f8f630e5p-3, -0x1.007a5f8f630e3p-3, 0x1.fe40ce6d4f01bp-3,
+     0x1.fe40ce6d4f030p-3},
+    {0x1.007a5f8f630e3p-3, 0x1.007a5f8f630e5p-3, 0x1.fe40ce6d4f01bp-3,
+     0x1.fe40ce6d4f030p-3},
+    {0x1.78a8d20a8b19cp-2, 0x1.78a8d20a8b19ep-2, 0x1.de3155c256a18p-3,
+     0x1.de3155c256b52p-3},
+    {0x1.2cb4f05c077f8p-1, 0x1.2cb4f05c077fap-1, 0x1.a0163e6b1a279p-3,
+     0x1.a0163e6b1b46ap-3},
+    {0x1.8a30aeed88f35p-1, 0x1.8a30aeed88f37p-1, 0x1.47d7258f1fdcdp-3,
+     0x1.47d7258f25d51p-3},
+    {0x1.cee874ffb88b3p-1, 0x1.cee874ffb88b5p-1, 0x1.b60602bcd3d56p-4,
+     0x1.b60602bcf85a4p-4},
+    {0x1.f68f1d8e42e80p-1, 0x1.f68f1d8e42e82p-1, 0x1.8275d9de72d45p-5,
+     0x1.8275d9dedaddep-5},
+}};
+inline constexpr std::array<GaussPoint, 16> kGauss16{{
+    {-0x1.fa92c264d787fp-1, -0x1.fa92c264d787dp-1, 0x1.bcddab41dec5ap-6,
+     0x1.bcddab55197bep-6},
+    {-0x1.e39f56616f9b1p-1, -0x1.e39f56616f9afp-1, 0x1.fdfb1a27efbd3p-5,
+     0x1.fdfb1a30350e0p-5},
+    {-0x1.bb3403514e484p-1, -0x1.bb3403514e482p-1, 0x1.85c4ee78dfd1dp-4,
+     0x1.85c4ee7ab8792p-4},
+    {-0x1.82c45dda4726cp-1, -0x1.82c45dda4726ap-1, 0x1.fe7af2ba91fc7p-4,
+     0x1.fe7af2bb1510fp-4},
+    {-0x1.3c5a466d5e8b9p-1, -0x1.3c5a466d5e8b7p-1, 0x1.325f61bc9dd8dp-3,
+     0x1.325f61bca9bf1p-3},
+    {-0x1.d50259a43a773p-2, -0x1.d50259a43a771p-2, 0x1.5a6ebbb5a6fa2p-3,
+     0x1.5a6ebbb5a7c56p-3},
+    {-0x1.205cae642337dp-2, -0x1.205cae642337bp-2, 0x1.75f8c77e0bf66p-3,
+     0x1.75f8c77e0c0b8p-3},
+    {-0x1.852bd6676a9fap-4, -0x1.852bd6676a9f8p-4, 0x1.83feae80e4df6p-3,
+     0x1.83feae80e4e01p-3},
+    {0x1.852bd6676a9f8p-4, 0x1.852bd6676a9fap-4, 0x1.83feae80e4df6p-3,
+     0x1.83feae80e4e01p-3},
+    {0x1.205cae642337bp-2, 0x1.205cae642337dp-2, 0x1.75f8c77e0bf66p-3,
+     0x1.75f8c77e0c0b8p-3},
+    {0x1.d50259a43a771p-2, 0x1.d50259a43a773p-2, 0x1.5a6ebbb5a6fa2p-3,
+     0x1.5a6ebbb5a7c56p-3},
+    {0x1.3c5a466d5e8b7p-1, 0x1.3c5a466d5e8b9p-1, 0x1.325f61bc9dd8dp-3,
+     0x1.325f61bca9bf1p-3},
+    {0x1.82c45dda4726ap-1, 0x1.82c45dda4726cp-1, 0x1.fe7af2ba91fc7p-4,
+     0x1.fe7af2bb1510fp-4},
+    {0x1.bb3403514e482p-1, 0x1.bb3403514e484p-1, 0x1.85c4ee78dfd1dp-4,
+     0x1.85c4ee7ab8792p-4},
+    {0x1.e39f56616f9afp-1, 0x1.e39f56616f9b1p-1, 0x1.fdfb1a27efbd3p-5,
+     0x1.fdfb1a30350e0p-5},
+    {0x1.fa92c264d787dp-1, 0x1.fa92c264d787fp-1, 0x1.bcddab41dec5ap-6,
+     0x1.bcddab55197bep-6},
+}};
+}  // namespace thermal_oracle
+#endif
diff --git a/tests/integration/analysis/thermal_oracle_generator.py b/tests/integration/analysis/thermal_oracle_generator.py
new file mode 100644
index 0000000..7c643f9
--- /dev/null
+++ b/tests/integration/analysis/thermal_oracle_generator.py
@@ -0,0 +1,244 @@
+"""Generate test-only certified Gauss brackets from exact rational arithmetic.
+
+Newton/Decimal produces candidates only. Exact Fraction signs prove n distinct
+Legendre roots; rational interval recurrence encloses the associated weights.
+No runtime solver dependency, reference modification, or libm ULP premise.
+"""
+from decimal import Decimal, localcontext
+from fractions import Fraction as F
+from pathlib import Path
+import math
+import re
+import json
+
+
+def add(a, b):
+    return a[0] + b[0], a[1] + b[1]
+
+
+def neg(a):
+    return -a[1], -a[0]
+
+
+def mul(a, b):
+    v = [x * y for x in a for y in b]
+    return min(v), max(v)
+
+
+def div(a, b):
+    assert b[0] > 0 or b[1] < 0
+    return mul(a, (1 / b[1], 1 / b[0]))
+
+
+def scale(a, b):
+    return mul(a, (F(b), F(b)))
+
+
+def legendre(n, x):
+    p, q = (F(1), F(1)), x
+    for k in range(2, n + 1):
+        p, q = q, scale(add(scale(mul(x, q), 2 * k - 1), neg(scale(p, k - 1))), F(1, k))
+    derivative = div(scale(add(mul(x, q), neg(p)), n), add(mul(x, x), (-F(1), -F(1))))
+    return q, derivative
+
+
+def bracket(value):
+    return math.nextafter(float(value[0]), -math.inf), math.nextafter(float(value[1]), math.inf)
+
+
+def rule(n):
+    result = []
+    for i in range(n):
+        with localcontext() as context:
+            context.prec = 90
+            z = Decimal(str(math.cos(math.pi * (i + .75) / (n + .5))))
+            for _ in range(30):
+                p, q = Decimal(1), z
+                for k in range(2, n + 1):
+                    p, q = q, ((2 * k - 1) * z * q - (k - 1) * p) / k
+                d = n * (z * q - p) / (z * z - 1)
+                z -= q / d
+        lo, hi = bracket((F(z), F(z)))
+        pl = legendre(n, (F(lo), F(lo)))[0][0]
+        ph = legendre(n, (F(hi), F(hi)))[0][0]
+        assert pl * ph < 0
+        x = F(lo), F(hi)
+        derivative = legendre(n, x)[1]
+        w = div((F(2), F(2)), mul(add((F(1), F(1)), neg(mul(x, x))), mul(derivative, derivative)))
+        result.append((lo, hi, *bracket(w)))
+    result.sort()
+    assert all(a[1] < b[0] for a, b in zip(result, result[1:]))
+    assert sum(F(r[2]) for r in result) < 2 < sum(F(r[3]) for r in result)
+    return result
+
+
+def slab_fixtures():
+    # Approved Appendix A supplies the independently implemented Cholesky /
+    # Jacobi spectral candidates; Appendix E supplies certified transcendentals.
+    document = Path(__file__).resolve().parents[3] / "docs/shell-transient-heat-transfer/numerical-oracle-scratch.md"
+    blocks = re.findall(r"```python\n(.*?)```", document.read_text(encoding="utf-8-sig"), re.S)
+    spectral = {}
+    exec(blocks[0].split('out={"lambda"')[0], spectral)
+    interval_scope = {}
+    exec(blocks[-1].split("\nfor t in [.001")[0], interval_scope)
+    iv, ia, im, isc = (interval_scope[k] for k in ("interval", "add", "mul", "scale"))
+    ex = interval_scope["exneg"]
+    n = 4
+
+    def exact_inverse(matrix):
+        a = [list(row) + [F(i == j) for j in range(n)] for i, row in enumerate(matrix)]
+        for k in range(n):
+            pivot = next(i for i in range(k, n) if a[i][k])
+            a[k], a[pivot] = a[pivot], a[k]
+            divisor = a[k][k]
+            a[k] = [x / divisor for x in a[k]]
+            for i in range(n):
+                if i != k:
+                    factor = a[i][k]
+                    a[i] = [x - factor * y for x, y in zip(a[i], a[k])]
+        return [row[n:] for row in a]
+
+    def exact_mm(a, b):
+        return [[sum((x*y for x, y in zip(row, col)), F(0)) for col in zip(*b)] for row in a]
+
+    def matrix_multiply(a, b):
+        return [[sum_intervals([im(x, y) for x, y in zip(row, col)]) for col in zip(*b)] for row in a]
+
+    def sum_intervals(values):
+        total = iv(0)
+        for v in values:
+            total = ia(total, v)
+        return total
+
+    def apply(a, v):
+        return [sum_intervals([im(x, y) for x, y in zip(row, v)]) for row in a]
+
+    def power(a, count):
+        result = [[iv(i == j) for j in range(n)] for i in range(n)]
+        while count:
+            if count & 1:
+                result = matrix_multiply(result, a)
+            count >>= 1
+            if count:
+                a = matrix_multiply(a, a)
+        return result
+
+    c = [[F(0) for _ in range(5)] for _ in range(5)]
+    k = [[F(0) for _ in range(5)] for _ in range(5)]
+    for offset in (0, 2):
+        for i in range(3):
+            for j in range(3):
+                c[offset+i][offset+j] += F([[4, 2, -1], [2, 16, 2], [-1, 2, 4]][i][j], 60)
+                k[offset+i][offset+j] += F(2, 3)*[[7, -8, 1], [-8, 16, -8], [1, -8, 7]][i][j]
+    cf, kf = [r[1:] for r in c[1:]], [r[1:] for r in k[1:]]
+    inverse = exact_inverse(cf)
+    d = exact_mm(inverse, kf)
+    y = [-1-sum(inverse[i][j]*c[j+1][0] for j in range(n)) for i in range(n)]
+
+    def matrix_exponential(t):
+        # An independent enclosure of exp(-Cff^-1 Kff t) certifies all spectral
+        # eigenpair/coefficient errors together, rather than assuming libm ULPs.
+        # For ||B||inf<=1/2, the 64-term exponential remainder is bounded by
+        # 2*(1/2)^65/65!. Entrywise intervals preserve subsequent squarings.
+        magnitude = max(sum(abs(x*t) for x in row) for row in d)
+        shift = 0
+        while magnitude > F(1, 2):
+            magnitude /= 2
+            shift += 1
+        a = [[iv(-x*t/F(2)**shift) for x in row] for row in d]
+        term = [[iv(i == j) for j in range(n)] for i in range(n)]
+        result = [row[:] for row in term]
+        for order in range(1, 65):
+            term = [[isc(x, F(1, order)) for x in row] for row in matrix_multiply(term, a)]
+            result = [[ia(x, z) for x, z in zip(row, terms)] for row, terms in zip(result, term)]
+        remainder = 2*F(1, 2)**65/math.factorial(65)
+        result = [[ia(x, iv(-remainder, remainder)) for x in row] for row in result]
+        for _ in range(shift):
+            result = matrix_multiply(result, result)
+        return result
+
+    def outputs(error):
+        temperature = [iv(1)] + [ia(iv(1), z) for z in error]
+        derivative = [iv(0)] + apply([[iv(-x) for x in row] for row in d], error)
+        flux = []
+        for offset, s in ((0, -1), (0, 0), (0, 1), (2, -1), (2, 0), (2, 1)):
+            deriv = [F(4)*F(s-F(1, 2)), F(-8*s), F(4)*F(s+F(1, 2))]
+            flux.append(sum_intervals([isc(temperature[offset+i], -deriv[i]) for i in range(3)]))
+        energy = sum_intervals([isc(temperature[j], sum(c[i][j] for i in range(5))) for j in range(5)])
+        reaction = sum_intervals([isc(temperature[j], k[0][j]) for j in range(5)] +
+                                 [isc(derivative[j], c[0][j]) for j in range(5)])
+        return temperature + flux + [energy, reaction]
+
+    def spectral_candidate(t, count):
+        amplitudes = []
+        for eigenvalue, coefficient in zip(spectral["lam"], spectral["coef"]):
+            decay = ex(iv(F(eigenvalue)*t)) if count is None else power_scalar(
+                interval_scope["div"](iv(1), iv(1+F(eigenvalue)*t/count)), count)
+            amplitudes.append(isc(decay, F(coefficient)))
+        return outputs(apply([[iv(F(x)) for x in row] for row in spectral["W"]], amplitudes))
+
+    def power_scalar(x, count):
+        result = iv(1)
+        while count:
+            if count & 1:
+                result = im(result, x)
+            count >>= 1
+            if count:
+                x = im(x, x)
+        return result
+
+    table = [[.141, 10.7, .0304, 10.9], [.0153, 1.28, .00468, .824],
+             [.00598, .253, .00130, .00632], [.00139, .136, .000666, .00953],
+             [.000142, .0284, .0000945, .000526], [.0000414, .00816, .0000159, .0000479]]
+    fixtures, metrics = [], []
+    for time, space in zip((.001, .01, .05, .1, .5, 1.), table):
+        t = F(time)
+        temperatures, flux, energy, _, _ = interval_scope["series_interval"](time)
+        continuum = temperatures + [flux[i] for i in (0, 1, 2, 2, 3, 4)] + [energy, flux[0]]
+        solutions = []
+        spectral_enclosure = 0.
+        for count in (None, 4096, 8192):
+            if count is None:
+                operator = matrix_exponential(t)
+            else:
+                a = [[cf[i][j]+t/count*kf[i][j] for j in range(n)] for i in range(n)]
+                operator = power([[iv(x) for x in row] for row in exact_mm(exact_inverse(a), cf)], count)
+            exact = outputs(apply(operator, [iv(x) for x in y]))
+            candidate = spectral_candidate(t, count)
+            enclosed = [(min(a[0], b[0]), max(a[1], b[1])) for a, b in zip(exact, candidate)]
+            spectral_enclosure = max(spectral_enclosure, max(float(hi-lo) for lo, hi in enclosed))
+            solutions.append(enclosed)
+        fixtures.append((time, space, continuum, *solutions))
+        metrics.append({"time": time, "spectral_and_matrix_function_enclosure": spectral_enclosure})
+    print(json.dumps({"eigenvalues": spectral["lam"], "slab_oracle_metrics": metrics}, indent=2))
+    return fixtures
+
+
+if __name__ == "__main__":
+    output = ["// Generated by thermal_oracle_generator.py; exact root/weight certificates.",
+              "#ifndef FESA_TESTS_THERMAL_GAUSS_ORACLE_H_", "#define FESA_TESTS_THERMAL_GAUSS_ORACLE_H_",
+              "#include <array>", "namespace thermal_oracle {",
+              "struct GaussPoint { double x_lo,x_hi,w_lo,w_hi; };"]
+    for n in (12, 16):
+        rows = rule(n)
+        output.append(f"inline constexpr std::array<GaussPoint,{n}> kGauss{n}{{{{")
+        output.extend("  {" + ",".join(x.hex() for x in row) + "}," for row in rows)
+        output.append("}};")
+    output.extend(["}", "#endif", ""])
+    destination = Path(__file__).with_name("thermal_gauss_oracle.h")
+    destination.write_text("\n".join(output), encoding="utf-8")
+    print("Certified Legendre roots/weights: 12 + 16 distinct roots;", destination.name)
+    fixtures = slab_fixtures()
+    output = ["// Generated by thermal_oracle_generator.py. Outward binary64 brackets.",
+              "#ifndef FESA_TESTS_THERMAL_SLAB_ORACLE_H_", "#define FESA_TESTS_THERMAL_SLAB_ORACLE_H_",
+              "#include <array>", "namespace thermal_oracle {",
+              "struct Bracket { double lo,hi; };",
+              "struct SlabCase { double time; std::array<double,4> space; std::array<Bracket,13> fourier,continuous,be4096,be8192; };",
+              "inline constexpr std::array<SlabCase,6> kSlab{{"]
+    for time, space, *values in fixtures:
+        output.append("{" + time.hex() + ",{{" + ",".join(x.hex() for x in space) + "}},")
+        for rows in values:
+            output.append("{{" + ",".join("{" + ",".join(x.hex() for x in bracket(row)) + "}" for row in rows) + "}},")
+        output.append("},")
+    output.extend(["}};", "}", "#endif", ""])
+    Path(__file__).with_name("thermal_slab_oracle.h").write_text("\n".join(output), encoding="utf-8")
diff --git a/tests/integration/analysis/thermal_physics_test.cpp b/tests/integration/analysis/thermal_physics_test.cpp
new file mode 100644
index 0000000..de76b29
--- /dev/null
+++ b/tests/integration/analysis/thermal_physics_test.cpp
@@ -0,0 +1,1139 @@
+#include <gtest/gtest.h>
+
+#include <algorithm>
+#include <array>
+#include <cmath>
+#include <cstdint>
+#include <functional>
+#include <iomanip>
+#include <iostream>
+#include <limits>
+#include <sstream>
+#include <stdexcept>
+#include <utility>
+#include <vector>
+
+#include "analysis/thermal_transient_algebra.h"
+#include "fesa/analysis/transient_heat_procedure.h"
+#include "fesa/solvers/linear/mkl_pardiso_solver.h"
+#include "thermal_gauss_oracle.h"
+#include "thermal_slab_oracle.h"
+
+namespace {
+namespace ti = fesa::thermal_transient_internal;
+using fesa::Vector;
+using fesa::Vector3;
+
+// Independent test arithmetic: outward endpoints after each scalar operation.
+// No production interval, kernel, transcendental, or tolerance helper is used
+// to construct the analytical expectations.
+double Down(double x) {
+  return std::nextafter(x, -std::numeric_limits<double>::infinity());
+}
+double Up(double x) {
+  return std::nextafter(x, std::numeric_limits<double>::infinity());
+}
+struct I {
+  double lo, hi;
+  I(double x = 0) : lo(x), hi(x) {}
+  I(double lo, double hi) : lo(lo), hi(hi) {}
+  double Mid() const { return lo + (hi - lo) * .5; }
+  double Abs() const { return std::max(std::abs(lo), std::abs(hi)); }
+};
+I operator+(I a, I b) {
+  if (a.lo == 0 && a.hi == 0) return b;
+  if (b.lo == 0 && b.hi == 0) return a;
+  return {Down(a.lo + b.lo), Up(a.hi + b.hi)};
+}
+I operator-(I a, I b) { return a + I{-b.hi, -b.lo}; }
+I operator*(I a, I b) {
+  if ((a.lo == 0 && a.hi == 0) || (b.lo == 0 && b.hi == 0)) return I{0};
+  const std::array<double, 4> v{a.lo * b.lo, a.lo * b.hi, a.hi * b.lo,
+                                a.hi * b.hi};
+  return {Down(*std::min_element(v.begin(), v.end())),
+          Up(*std::max_element(v.begin(), v.end()))};
+}
+I operator/(I a, I b) {
+  if (b.lo <= 0 && b.hi >= 0)
+    throw std::runtime_error{"Oracle division includes zero."};
+  return a * I{Down(1 / b.hi), Up(1 / b.lo)};
+}
+I Sqrt(I a) {
+  if (a.lo < 0)
+    throw std::runtime_error{"Oracle square root includes negative values."};
+  if (a.hi == 0) return I{0};
+  return {a.lo == 0 ? 0 : Down(std::sqrt(a.lo)), Up(std::sqrt(a.hi))};
+}
+I Square(I x) {
+  if (x.lo == 0 && x.hi == 0) return I{0};
+  const double lower =
+      x.lo <= 0 && x.hi >= 0 ? 0 : Down(std::min(x.lo * x.lo, x.hi * x.hi));
+  return {lower, Up(x.Abs() * x.Abs())};
+}
+I Gamma(std::size_t n) {
+  return (I{static_cast<double>(n)} * I{0x1p-53}) /
+         (I{1} - I{static_cast<double>(n)} * I{0x1p-53});
+}
+I Pi() {
+  static const I value = [] {
+    // Machin's formula; alternating arctan remainder after terms 0..63.
+    const auto atan = [](double denominator) {
+      const I q = I{1} / I{denominator}, q2 = q * q;
+      I power = q, sum;
+      for (int j = 0; j < 64; ++j) {
+        sum = sum + power / I{static_cast<double>(2 * j + 1)};
+        power = power * q2 * I{-1};
+      }
+      const double tail = (power / I{129}).Abs();
+      return sum + I{-tail, tail};
+    };
+    return I{16} * atan(5) - I{4} * atan(239);
+  }();
+  return value;
+}
+I Cos(I x) {
+  if (x.hi < 0) x = I{-x.hi, -x.lo};
+  if (x.lo > Pi().hi / 2) {
+    const auto reduced = Cos(Pi() - x);
+    return I{-reduced.hi, -reduced.lo};
+  }
+  // |x|<=pi<4. Thirty-two terms; the omitted absolute term bounds the
+  // alternating tail after decreasing terms begin. All polynomial work is
+  // enclosed independently of a libm accuracy assumption.
+  I term{1}, sum{1};
+  const I x2 = x * x;
+  for (int j = 1; j < 32; ++j) {
+    term = term * x2 * I{-1} / I{static_cast<double>((2 * j - 1) * 2 * j)};
+    sum = sum + term;
+  }
+  const double tail = (term * x2 / I{63. * 64.}).Abs();
+  return sum + I{-tail, tail};
+}
+I ExpNegative(I x) {
+  if (x.lo < 0) throw std::runtime_error{"Negative exponential argument."};
+  int squares = 0;
+  while (x.hi > .5) {
+    x = x / I{2};
+    ++squares;
+  }
+  I term{1}, sum{1};
+  for (int j = 1; j <= 64; ++j) {
+    term = term * x * I{-1} / I{static_cast<double>(j)};
+    sum = sum + term;
+  }
+  const double tail = (term * x / I{65}).Abs();
+  sum = sum + I{-tail, tail};
+  for (int j = 0; j < squares; ++j) sum = sum * sum;
+  return sum;
+}
+I Power(I x, std::uint32_t n) {
+  I result{1};
+  for (; n; n >>= 1) {
+    if (n & 1U) result = result * x;
+    if (n > 1) x = x * x;
+  }
+  return result;
+}
+template <class T>
+T Take(fesa::Result<T> r) {
+  if (!r.HasValue()) {
+    std::ostringstream text;
+    for (const auto& d : r.GetStatus().Diagnostics())
+      text << d.code << ": " << d.message << '\n';
+    throw std::runtime_error{text.str()};
+  }
+  return std::move(r.Value());
+}
+void Checked(const fesa::Status& status) {
+  if (!status.IsOk()) {
+    std::ostringstream text;
+    for (const auto& d : status.Diagnostics())
+      text << d.code << ": " << d.message << '\n';
+    throw std::runtime_error{text.str()};
+  }
+}
+fesa::ModelDefinition Mesh(
+    std::size_t n, double period, std::uint32_t count,
+    const std::function<double(std::size_t, std::size_t)>& initial) {
+  fesa::ModelDefinition m;
+  for (std::size_t j = 0; j <= n; ++j)
+    for (std::size_t i = 0; i <= n; ++i) {
+      const auto index = static_cast<fesa::EntityIndex>(j * (n + 1) + i);
+      m.nodes.push_back(
+          {{"plate", static_cast<std::int64_t>(index + 1),
+            std::to_string(index + 1)},
+           {static_cast<double>(i) / n, static_cast<double>(j) / n, 0},
+           {}});
+      m.thermal_initial_conditions.push_back({index, initial(i, j), {}});
+    }
+  m.thermal_materials = {{{"plate", 1, "1"}, "unit", 1, 1, 1, {}}};
+  m.thermal_sections = {{{"plate", 1, "1"}, "unit", 1, 0, {}}};
+  for (std::size_t j = 0; j < n; ++j)
+    for (std::size_t i = 0; i < n; ++i) {
+      const auto node = static_cast<fesa::EntityIndex>(j * (n + 1) + i);
+      const auto label = static_cast<std::int64_t>(j * n + i + 1);
+      m.thermal_elements.push_back(
+          {{"plate", label, std::to_string(label)},
+           {{node, node + 1, static_cast<fesa::EntityIndex>(node + n + 2),
+             static_cast<fesa::EntityIndex>(node + n + 1)}},
+           0,
+           {}});
+    }
+  m.thermal_step =
+      fesa::ThermalStepDefinition{"physics", period / count, period, count};
+  return m;
+}
+struct Prepared {
+  fesa::Domain domain;
+  fesa::AnalysisModel model;
+  fesa::DofManager dofs;
+  fesa::ThermalAssembly assembly;
+  fesa::ThermalAnalysisState state;
+  Prepared(fesa::ModelDefinition definition, const fesa::ParallelFor& parallel)
+      : domain(Take(fesa::Domain::Create(std::move(definition)))),
+        model(Take(fesa::AnalysisModel::Create(domain))),
+        dofs(Take(fesa::DofManager::CreateThermal(model))),
+        assembly(Take(fesa::AssembleThermal(model, dofs, parallel))),
+        state(Take(fesa::ThermalAnalysisState::Create(model, dofs, assembly))) {
+  }
+};
+struct RunMetrics {
+  double beta = 0, full_beta = 0, capacity_budget = 0;
+  double energy_defect = 0, energy_bound = 0, dissipation_defect = 0,
+         dissipation_bound = 0;
+  std::uint64_t temperature_values = 0, flux_values = 0;
+  std::vector<ti::ThermalStepEvidence> steps;
+  double last_free_residual_l1 = 0, last_free_residual_inf = 0,
+         last_reaction_evaluation_error = 0;
+};
+RunMetrics Advance(Prepared& p, const fesa::ParallelFor& parallel,
+                   double mass_lower, double initial_capacity_budget = 0) {
+  fesa::MklPardisoSolver solver;
+  ti::ThermalFactorization cache;
+  ti::ThermalSystem system{
+      p.assembly.Conductivity(),   p.assembly.Capacity(),
+      p.assembly.HeatInput(),      &p.assembly.ConductivityError(),
+      &p.assembly.CapacityError(), &p.assembly.HeatInputError()};
+  ti::ThermalPartition partition{p.dofs.FreeDofs(), p.dofs.ConstrainedDofs(),
+                                 p.dofs.PrescribedValues()};
+  RunMetrics metrics;
+  I budget{initial_capacity_budget};
+  const std::uint64_t count = p.state.TimeGrid().Count();
+  for (std::uint64_t step = 1; step <= count; ++step) {
+    const auto schedule =
+        Take(p.state.TimeGrid().Frame(static_cast<std::uint32_t>(step)));
+    const auto& old = p.state.CommittedState().temperature;
+    auto next = Take(ti::SolveThermalIncrement(system, partition, old,
+                                               schedule.delta, solver, cache));
+    const auto& e = next.evidence;
+    metrics.steps.push_back(e);
+    EXPECT_LE(e.beta_cert, e.beta_limit);
+    EXPECT_LE(e.full_beta_cert, e.beta_limit);
+    EXPECT_LE(e.energy_defect, e.energy_bound);
+    if (e.dissipation_checked)
+      EXPECT_LE(e.dissipation_defect, e.dissipation_bound);
+    metrics.beta = std::max(metrics.beta, e.beta_cert);
+    metrics.full_beta = std::max(metrics.full_beta, e.full_beta_cert);
+    metrics.energy_defect = std::max(metrics.energy_defect, e.energy_defect);
+    metrics.energy_bound = std::max(metrics.energy_bound, e.energy_bound);
+    metrics.dissipation_defect =
+        std::max(metrics.dissipation_defect, e.dissipation_defect);
+    metrics.dissipation_bound =
+        std::max(metrics.dissipation_bound, e.dissipation_bound);
+    I residual_square;
+    for (auto row : partition.free)
+      residual_square =
+          residual_square + Square(I{next.certified_residual[row]});
+    // Approved contraction for structured homogeneous error, using observed
+    // component residual enclosures including E_C/E_K, never beta_limit.
+    budget = budget +
+             I{schedule.delta} * Sqrt(residual_square) / Sqrt(I{mass_lower});
+    if (step == count) {
+      I free_l1, reaction_error;
+      for (auto row : partition.free) {
+        free_l1 = free_l1 + I{next.certified_residual[row]};
+        metrics.last_free_residual_inf = std::max(
+            metrics.last_free_residual_inf, next.certified_residual[row]);
+      }
+      for (std::size_t c = 0; c < partition.constrained.size(); ++c) {
+        const auto row = partition.constrained[c];
+        I residual = I{-system.q[row]};
+        for (auto entry = system.c.RowOffsets()[row];
+             entry < system.c.RowOffsets()[row + 1]; ++entry) {
+          const auto col = system.c.ColumnIndices()[entry];
+          const I capacity =
+              I{system.c.Values()[entry]} + I{-system.c_error->Values()[entry],
+                                              system.c_error->Values()[entry]};
+          const I conductivity =
+              I{system.k.Values()[entry]} + I{-system.k_error->Values()[entry],
+                                              system.k_error->Values()[entry]};
+          residual = residual +
+                     capacity * ((I{next.temperature[col]} - I{old[col]}) /
+                                 I{schedule.delta}) +
+                     conductivity * I{next.temperature[col]};
+        }
+        residual =
+            residual + I{-(*system.q_error)[row], (*system.q_error)[row]};
+        reaction_error =
+            reaction_error + I{(residual - I{next.reaction[c]}).Abs()};
+      }
+      metrics.last_free_residual_l1 = free_l1.hi;
+      metrics.last_reaction_evaluation_error = reaction_error.hi;
+    }
+    auto flux = Take(p.assembly.RecoverFlux(next.temperature, parallel));
+    auto frame = Take(p.state.MakeCompletedFrame(
+        static_cast<std::uint32_t>(step), std::move(next.temperature),
+        std::move(next.reaction), std::move(flux)));
+    Checked(p.state.CommitAcceptedFrame(frame));
+  }
+  metrics.capacity_budget = budget.hi;
+  EXPECT_TRUE(p.state.IsComplete());
+  EXPECT_EQ(p.state.AcceptedFrames().size(), count + 1);
+  EXPECT_FALSE(p.state.AcceptedFrames()[0].completed.has_value());
+  for (const auto& frame : p.state.AcceptedFrames()) {
+    metrics.temperature_values += frame.temperature.Size();
+    if (frame.completed)
+      metrics.flux_values += frame.completed->global_heat_flux.size() * 72U;
+  }
+  return metrics;
+}
+void Print(const char* name, std::size_t mesh, std::uint32_t count,
+           const RunMetrics& m) {
+  std::cout << std::setprecision(17) << "THERMAL_PHYSICS case=" << name
+            << " n=" << mesh << " steps=" << count << " beta=" << m.beta
+            << " physical_beta=" << m.full_beta
+            << " capacity_budget=" << m.capacity_budget
+            << " energy_defect=" << m.energy_defect
+            << " energy_bound=" << m.energy_bound
+            << " dissipation_defect=" << m.dissipation_defect
+            << " dissipation_bound=" << m.dissipation_bound
+            << " temperature_values=" << m.temperature_values
+            << " flux_values=" << m.flux_values << '\n';
+}
+// Exact tensor mass independently assembled from the rational Q1 and P2
+// element integrals. This checks actual 2-D nodal/station vectors.
+I CapacityNorm(const Vector& values, std::size_t n) {
+  constexpr int q[4][4] = {
+      {4, 2, 1, 2}, {2, 4, 2, 1}, {1, 2, 4, 2}, {2, 1, 2, 4}};
+  constexpr int z[3][3] = {{4, 2, -1}, {2, 16, 2}, {-1, 2, 4}};
+  I sum;
+  for (std::size_t y = 0; y < n; ++y)
+    for (std::size_t x = 0; x < n; ++x) {
+      const std::size_t a = y * (n + 1) + x;
+      const std::array<std::size_t, 4> nodes{a, a + 1, a + n + 2, a + n + 1};
+      for (std::size_t i = 0; i < 4; ++i)
+        for (std::size_t j = 0; j < 4; ++j)
+          for (std::size_t segment = 0; segment < 2; ++segment)
+            for (std::size_t u = 0; u < 3; ++u)
+              for (std::size_t v = 0; v < 3; ++v)
+                sum = sum + I{values[5 * nodes[i] + 2 * segment + u]} *
+                                I{values[5 * nodes[j] + 2 * segment + v]} *
+                                (I{static_cast<double>(q[i][j] * z[u][v])} /
+                                 I{2160. * n * n});
+    }
+  return Sqrt(sum);
+}
+TEST(ThermalPhysics, InsulatedUniformAndZeroRetainAllFramesWithConservation) {
+  fesa::SerialParallelFor parallel;
+  for (double temperature : {0., 2., -2.}) {
+    Prepared p{Mesh(2, .25, 8, [&](auto, auto) { return temperature; }),
+               parallel};
+    const auto metrics = Advance(p, parallel, Down(1. / 8640));
+    Print("uniform", 2, 8, metrics);
+    const double bound = (I{metrics.capacity_budget} / Sqrt(I{1. / 8640})).hi;
+    for (const auto& frame : p.state.AcceptedFrames())
+      for (std::size_t i = 0; i < frame.temperature.Size(); ++i)
+        EXPECT_LE(std::abs(frame.temperature[i] - temperature), bound);
+    if (temperature == 0) {
+      EXPECT_EQ(metrics.capacity_budget, 0);
+      EXPECT_EQ(metrics.energy_defect, 0);
+      EXPECT_EQ(metrics.dissipation_defect, 0);
+    }
+    EXPECT_EQ(metrics.flux_values, 8U * 4U * 72U);
+  }
+}
+void VerifyDiscreteEigenmode(const Prepared& p, const std::array<I, 9>& mode,
+                             I lambda) {
+  constexpr int mass[4][4] = {
+      {4, 2, 1, 2}, {2, 4, 2, 1}, {1, 2, 4, 2}, {2, 1, 2, 4}};
+  constexpr int stiffness[4][4] = {
+      {4, -1, -2, -1}, {-1, 4, -1, -2}, {-2, -1, 4, -1}, {-1, -2, -1, 4}};
+  constexpr int z_mass[3][3] = {{4, 2, -1}, {2, 16, 2}, {-1, 2, 4}};
+  constexpr int z_stiffness[3][3] = {{7, -8, 1}, {-8, 16, -8}, {1, -8, 7}};
+  std::array<I, 405> cv{}, kv{};
+  for (std::size_t y = 0; y < 8; ++y)
+    for (std::size_t x = 0; x < 8; ++x) {
+      const auto node = y * 9 + x;
+      const std::array<std::size_t, 4> nodes{node, node + 1, node + 10,
+                                             node + 9};
+      for (std::size_t i = 0; i < 4; ++i)
+        for (std::size_t j = 0; j < 4; ++j)
+          for (std::size_t segment = 0; segment < 2; ++segment)
+            for (std::size_t a = 0; a < 3; ++a)
+              for (std::size_t b = 0; b < 3; ++b) {
+                const auto row = 5 * nodes[i] + 2 * segment + a;
+                const I m = I{static_cast<double>(mass[i][j])} / I{36 * 64.};
+                const I s = I{static_cast<double>(stiffness[i][j])} / I{6};
+                const I c = I{static_cast<double>(z_mass[a][b])} / I{60};
+                const I k = I{2. * z_stiffness[a][b]} / I{3};
+                cv[row] = cv[row] + m * c * mode[nodes[j] % 9];
+                kv[row] = kv[row] + (s * c + m * k) * mode[nodes[j] % 9];
+              }
+    }
+  const auto& c = p.assembly.Capacity();
+  const auto& k = p.assembly.Conductivity();
+  I independent_norm, actual_norm;
+  double residual_bound = 0;
+  const auto intersects = [](I a, I b) {
+    EXPECT_LE(std::max(a.lo, b.lo), std::min(a.hi, b.hi));
+  };
+  for (std::size_t row = 0; row < 405; ++row) {
+    I actual_c, actual_k;
+    for (auto entry = c.RowOffsets()[row]; entry < c.RowOffsets()[row + 1];
+         ++entry) {
+      const I v = mode[(c.ColumnIndices()[entry] / 5) % 9];
+      actual_c = actual_c + (I{c.Values()[entry]} +
+                             I{-p.assembly.CapacityError().Values()[entry],
+                               p.assembly.CapacityError().Values()[entry]}) *
+                                v;
+      actual_k =
+          actual_k + (I{k.Values()[entry]} +
+                      I{-p.assembly.ConductivityError().Values()[entry],
+                        p.assembly.ConductivityError().Values()[entry]}) *
+                         v;
+    }
+    intersects(actual_c, cv[row]);
+    intersects(actual_k, kv[row]);
+    const I residual = kv[row] - lambda * cv[row];
+    EXPECT_LE(residual.lo, 0);
+    EXPECT_GE(residual.hi, 0);
+    const I assembled_residual = actual_k - lambda * actual_c;
+    EXPECT_LE(assembled_residual.lo, 0);
+    EXPECT_GE(assembled_residual.hi, 0);
+    residual_bound = std::max(residual_bound, assembled_residual.Abs());
+    independent_norm = independent_norm + mode[(row / 5) % 9] * cv[row];
+    actual_norm = actual_norm + mode[(row / 5) % 9] * actual_c;
+  }
+  const I analytical_norm = (I{2} + mode[1]) / I{6};
+  intersects(independent_norm, analytical_norm);
+  intersects(actual_norm, analytical_norm);
+  const I normalized = actual_norm / analytical_norm;
+  EXPECT_LE(normalized.lo, 1);
+  EXPECT_GE(normalized.hi, 1);
+  std::cout << std::setprecision(17)
+            << "THERMAL_EIGENMODE residual_bound=" << residual_bound
+            << " C_norm_squared_lower=" << actual_norm.lo
+            << " C_norm_squared_upper=" << actual_norm.hi
+            << " normalized_lower=" << normalized.lo
+            << " normalized_upper=" << normalized.hi << '\n';
+}
+TEST(ThermalPhysics, DiscreteEigenmodeHasFirstOrderTimeConvergence) {
+  fesa::SerialParallelFor parallel;
+  std::array<I, 9> cosines;
+  double initial_error = 0;
+  for (std::size_t i = 0; i < cosines.size(); ++i) {
+    cosines[i] = Cos(Pi() * I{static_cast<double>(i)} / I{8});
+    initial_error =
+        std::max(initial_error, (cosines[i] - I{cosines[i].Mid()}).Abs());
+  }
+  const I lambda = I{384} * (I{1} - cosines[1]) / (I{2} + cosines[1]);
+  const double period = (I{1} / lambda).Mid();
+  const I exact = ExpNegative(lambda * I{period});
+  double previous_error = 0;
+  for (std::uint32_t count : {16U, 32U, 64U, 128U, 256U}) {
+    Prepared p{Mesh(8, period, count,
+                    [&](std::size_t i, auto) { return cosines[i].Mid(); }),
+               parallel};
+    if (count == 16) VerifyDiscreteEigenmode(p, cosines, lambda);
+    const auto metrics =
+        Advance(p, parallel, Down(1. / (2160 * 64)), initial_error);
+    const auto& actual = p.state.CommittedState().temperature;
+    Vector error{actual.Size()}, mode{actual.Size()}, difference{actual.Size()};
+    const I be = Power(I{1} / (I{1} + lambda * I{period / count}), count);
+    double oracle_width = 0;
+    for (std::size_t node = 0; node < 81; ++node)
+      for (std::size_t a = 0; a < 5; ++a) {
+        const auto i = node % 9;
+        const I v = cosines[i] * exact, expected = cosines[i] * be;
+        error[5 * node + a] = actual[5 * node + a] - v.Mid();
+        mode[5 * node + a] = cosines[i].Mid();
+        difference[5 * node + a] = actual[5 * node + a] - expected.Mid();
+        oracle_width =
+            std::max(oracle_width, (expected - I{expected.Mid()}).Abs());
+      }
+    const I mode_norm = CapacityNorm(mode, 8);
+    const I error_norm = CapacityNorm(error, 8) / mode_norm;
+    const I agreement = CapacityNorm(difference, 8);
+    const I total_budget = I{metrics.capacity_budget} + I{oracle_width};
+    EXPECT_LE(agreement.hi, total_budget.hi);
+    EXPECT_LE((total_budget / mode_norm).hi, error_norm.lo / 1000);
+    if (previous_error > 0) {
+      const double order = std::log2(previous_error / error_norm.Mid());
+      EXPECT_GE(order, .97);
+      EXPECT_LE(order, 1.02);
+      std::cout << std::setprecision(17) << "THERMAL_TIME N=" << count
+                << " error=" << error_norm.Mid() << " order=" << order
+                << " oracle_width=" << oracle_width << '\n';
+    }
+    previous_error = error_norm.Mid();
+    Print("time", 8, count, metrics);
+  }
+}
+
+TEST(ThermalPhysics, SignedOverlappingAndConstrainedPowerBalanceByComponent) {
+  fesa::SerialParallelFor parallel;
+  auto definition = Mesh(1, .5, 8, [](auto, auto) { return 0.; });
+  for (fesa::EntityIndex i = 0; i < 4; ++i) {
+    auto node = definition.nodes[i];
+    node.source_id.instance_name = "other";
+    node.coordinates[0] += 3;
+    definition.nodes.push_back(node);
+    definition.thermal_initial_conditions.push_back({i + 4, 0, {}});
+  }
+  definition.thermal_elements.push_back(
+      {{"other", 1, "1"}, {{4, 5, 7, 6}}, 0, {}});
+  definition.thermal_model_boundaries = {{0, 11, 0, {}}};
+  definition.thermal_step->loads = {{0, 11, 3, {}, 0},
+                                    {1, 15, -3, {}, 1},
+                                    {4, 11, 2, {}, 2},
+                                    {4, 11, 3, {}, 3},
+                                    {6, 15, -5, {}, 4}};
+  auto without_constrained = definition;
+  without_constrained.thermal_step->loads.erase(
+      without_constrained.thermal_step->loads.begin());
+  Prepared p{std::move(definition), parallel},
+      other{std::move(without_constrained), parallel};
+  EXPECT_EQ(p.assembly.HeatInput()[20], 5);
+  EXPECT_EQ(p.assembly.HeatInput()[34], -5);
+  const auto metrics = Advance(p, parallel, Down(1. / 2160));
+  const auto other_metrics = Advance(other, parallel, Down(1. / 2160));
+  Print("signed-component", 1, 8, metrics);
+  Print("without-constrained-power", 1, 8, other_metrics);
+  const std::array<I, 5> weights{I{1} / I{48}, I{1} / I{12}, I{1} / I{24},
+                                 I{1} / I{12}, I{1} / I{48}};
+  for (std::size_t step = 1; step < p.state.AcceptedFrames().size(); ++step) {
+    const auto& frame = p.state.AcceptedFrames()[step];
+    const auto& old = p.state.AcceptedFrames()[step - 1];
+    const auto& other_frame = other.state.AcceptedFrames()[step];
+    std::array<I, 2> energy;
+    for (std::size_t i = 0; i < 40; ++i) {
+      energy[i / 20] =
+          energy[i / 20] +
+          weights[i % 5] * (I{frame.temperature[i]} - I{old.temperature[i]});
+      EXPECT_DOUBLE_EQ(frame.temperature[i], other_frame.temperature[i]);
+    }
+    const double reaction = frame.completed->reaction_heat_rate[0];
+    const I first = energy[0] - I{frame.time_increment} * I{reaction};
+    // Independent exact unit-volume quadrature and component inventory.
+    EXPECT_LE(first.Abs(), metrics.steps[step - 1].energy_bound);
+    EXPECT_LE(energy[1].Abs(), metrics.steps[step - 1].energy_bound);
+    const double other_reaction = other_frame.completed->reaction_heat_rate[0];
+    const I reaction_shift = I{reaction} - I{other_reaction} + I{3};
+    const I subtraction =
+        Gamma(4) * (I{std::abs(reaction)} + I{std::abs(other_reaction)} + I{3});
+    EXPECT_LE(reaction_shift.Abs(), subtraction.hi);
+  }
+}
+
+fesa::SparseMatrix ZeroLike(const fesa::SparseMatrix& a) {
+  std::vector<fesa::CooContribution> coo;
+  for (std::size_t i = 0; i < a.Rows(); ++i)
+    for (std::size_t j = a.RowOffsets()[i]; j < a.RowOffsets()[i + 1]; ++j)
+      coo.push_back({i, a.ColumnIndices()[j], 0, 0, j});
+  return Take(fesa::SparseMatrix::FromCoo(a.Rows(), a.Columns(), std::move(coo),
+                                          {a.RowOffsets(), a.ColumnIndices()}));
+}
+// Independent dense inverse with a posteriori interval residual certificate.
+// Used only for small fixture partitions, never as a production solver.
+double CertifiedInverse(const fesa::SparseMatrix& matrix,
+                        const std::vector<std::size_t>& free) {
+  const auto n = free.size();
+  std::vector<std::vector<double>> a(n, std::vector<double>(n)),
+      z(n, std::vector<double>(n));
+  for (std::size_t i = 0; i < n; ++i) {
+    z[i][i] = 1;
+    for (std::size_t k = matrix.RowOffsets()[free[i]];
+         k < matrix.RowOffsets()[free[i] + 1]; ++k)
+      for (std::size_t j = 0; j < n; ++j)
+        if (matrix.ColumnIndices()[k] == free[j]) a[i][j] = matrix.Values()[k];
+  }
+  const auto original = a;
+  for (std::size_t k = 0; k < n; ++k) {
+    std::size_t pivot = k;
+    for (std::size_t i = k + 1; i < n; ++i)
+      if (std::abs(a[i][k]) > std::abs(a[pivot][k])) pivot = i;
+    if (a[pivot][k] == 0)
+      throw std::runtime_error{"Singular independent inverse."};
+    std::swap(a[k], a[pivot]);
+    std::swap(z[k], z[pivot]);
+    const double scale = a[k][k];
+    for (std::size_t j = 0; j < n; ++j) {
+      a[k][j] /= scale;
+      z[k][j] /= scale;
+    }
+    for (std::size_t i = 0; i < n; ++i)
+      if (i != k) {
+        const double factor = a[i][k];
+        for (std::size_t j = 0; j < n; ++j) {
+          a[i][j] -= factor * a[k][j];
+          z[i][j] -= factor * z[k][j];
+        }
+      }
+  }
+  double theta = 0, znorm = 0;
+  for (std::size_t i = 0; i < n; ++i) {
+    I defect, row;
+    for (std::size_t j = 0; j < n; ++j) {
+      I entry{static_cast<double>(i == j)};
+      for (std::size_t k = 0; k < n; ++k)
+        entry = entry - I{z[i][k]} * I{original[k][j]};
+      defect = defect + I{entry.Abs()};
+      row = row + I{std::abs(z[i][j])};
+    }
+    theta = std::max(theta, defect.hi);
+    znorm = std::max(znorm, row.hi);
+  }
+  if (theta >= 1)
+    throw std::runtime_error{"Independent inverse certificate inconclusive."};
+  return (I{znorm} / (I{1} - I{theta})).hi;
+}
+
+TEST(ThermalPhysics,
+     PrivateFullEndStripRespectsThicknessMaterialJumpAndFoldResistance) {
+  fesa::SerialParallelFor parallel;
+  for (bool folded : {false, true}) {
+    auto m = Mesh(1, 1, 1, [](auto, auto) { return 0.; });
+    m.nodes.push_back(
+        {{"plate", 5, "5"}, {folded ? 1. : 2., 0, folded ? 1. : 0.}, {}});
+    m.nodes.push_back(
+        {{"plate", 6, "6"}, {folded ? 1. : 2., 1, folded ? 1. : 0.}, {}});
+    m.thermal_initial_conditions.push_back({4, 0, {}});
+    m.thermal_initial_conditions.push_back({5, 0, {}});
+    m.thermal_materials.push_back({{"plate", 2, "2"}, "high-k", 10, 1, 1, {}});
+    m.thermal_sections.push_back({{"plate", 2, "2"}, "thin", .5, 1, {}});
+    m.thermal_elements.push_back({{"plate", 2, "2"}, {{1, 4, 5, 3}}, 1, {}});
+    Prepared p{std::move(m), parallel};
+    const auto zero = ZeroLike(p.assembly.Capacity());
+    ti::ThermalPartition partition;
+    std::vector<double> g;
+    for (std::size_t node = 0; node < 6; ++node)
+      for (std::size_t station = 0; station < 5; ++station) {
+        if (node == 1 || node == 3)
+          partition.free.push_back(5 * node + station);
+        else {
+          partition.constrained.push_back(5 * node + station);
+          g.push_back(node >= 4 ? 1. : 0.);
+        }
+      }
+    partition.prescribed = Vector{g.size()};
+    for (std::size_t i = 0; i < g.size(); ++i) partition.prescribed[i] = g[i];
+    ti::ThermalFactorization cache;
+    fesa::MklPardisoSolver solver;
+    const auto result = Take(ti::SolveThermalIncrement(
+        {p.assembly.Conductivity(), zero, p.assembly.HeatInput(),
+         &p.assembly.ConductivityError(), &zero, &p.assembly.HeatInputError()},
+        partition, Vector{30}, 1, solver, cache));
+    const double inverse =
+        CertifiedInverse(p.assembly.Conductivity(), partition.free);
+    double residual = 0, coefficient_error = 0;
+    for (auto i : partition.free) {
+      residual = std::max(residual, result.certified_residual[i]);
+      I row;
+      for (std::size_t k = p.assembly.ConductivityError().RowOffsets()[i];
+           k < p.assembly.ConductivityError().RowOffsets()[i + 1]; ++k)
+        row = row + I{p.assembly.ConductivityError().Values()[k]};
+      coefficient_error = std::max(coefficient_error, row.hi);
+    }
+    const I denominator = I{1} - I{inverse} * I{coefficient_error};
+    ASSERT_GT(denominator.lo, 0);
+    const I bound = I{inverse} / denominator * I{residual};
+    const I mid = I{5} / I{6};  // R=1/(1*1*1)+1/(10*.5*1)=6/5.
+    for (auto i : partition.free)
+      EXPECT_LE((I{result.temperature[i]} - mid).Abs(),
+                bound.hi + (mid - I{mid.Mid()}).Abs());
+    I left, right;
+    for (std::size_t i = 0; i < partition.constrained.size(); ++i)
+      if (partition.constrained[i] / 5 >= 4)
+        right = right + I{result.reaction[i]};
+      else
+        left = left + I{result.reaction[i]};
+    I reaction_budget;
+    for (auto i : partition.constrained) {
+      I row, scale, error;
+      for (std::size_t k = p.assembly.Conductivity().RowOffsets()[i];
+           k < p.assembly.Conductivity().RowOffsets()[i + 1]; ++k) {
+        const auto j = p.assembly.Conductivity().ColumnIndices()[k];
+        row = row + I{std::abs(p.assembly.Conductivity().Values()[k])} * bound;
+        scale = scale + I{std::abs(p.assembly.Conductivity().Values()[k])} *
+                            I{std::abs(result.temperature[j])};
+        error = error + I{p.assembly.ConductivityError().Values()[k]} *
+                            I{std::abs(result.temperature[j])};
+      }
+      const auto terms = p.assembly.Conductivity().RowOffsets()[i + 1] -
+                         p.assembly.Conductivity().RowOffsets()[i];
+      reaction_budget =
+          reaction_budget + row + error + Gamma(2 * terms + 2) * scale;
+    }
+    EXPECT_LE((right - mid).Abs(), reaction_budget.hi);
+    EXPECT_LE((left + mid).Abs(), reaction_budget.hi);
+    EXPECT_LE((left + right).Abs(), result.evidence.energy_bound);
+    auto flux = Take(p.assembly.RecoverFlux(result.temperature, parallel));
+    for (const auto& q : flux[0])
+      EXPECT_LE((I{q.X()} + mid).Abs(), 24 * bound.hi + Gamma(2048).hi * 24);
+    for (const auto& q : flux[1]) {
+      const double component = folded ? q.Z() : q.X();
+      EXPECT_LE((I{component} + I{5} / I{3}).Abs(),
+                400 * bound.hi + Gamma(2048).hi * 400);
+    }
+    std::cout << std::setprecision(17) << "THERMAL_STRIP folded=" << folded
+              << " resistance=1.2 reaction=" << right.Mid()
+              << " forward_bound=" << bound.hi
+              << " beta=" << result.evidence.beta_cert << '\n';
+  }
+}
+
+TEST(ThermalPhysics,
+     PrivateAffineFieldRotationNormalAndUnaveragedMidplaneTraces) {
+  fesa::SerialParallelFor parallel;
+  const auto rotate = [](Vector3 v) { return Vector3{v.Z(), v.X(), v.Y()}; };
+  for (double h : {.5, 1., 2.})
+    for (double k : {1., 10.})
+      for (bool flip : {false, true})
+        for (bool rotated : {false, true}) {
+          SCOPED_TRACE(h);
+          SCOPED_TRACE(k);
+          SCOPED_TRACE(flip);
+          SCOPED_TRACE(rotated);
+          auto m = Mesh(1, 1, 1, [](auto, auto) { return 0.; });
+          m.thermal_sections[0].thickness = h;
+          m.thermal_materials[0].conductivity = k;
+          if (flip) m.thermal_elements[0].node_indices = {{0, 2, 3, 1}};
+          Vector values{20};
+          for (std::size_t node = 0; node < 4; ++node) {
+            const auto& r = m.nodes[node].coordinates;
+            for (std::size_t station = 0; station < 5; ++station) {
+              const double zeta = -1 + .5 * station;
+              values[5 * node + station] =
+                  2 + 3 * r[0] - 2 * r[1] + 4 * (flip ? -1 : 1) * h * zeta / 2;
+            }
+            if (rotated) {
+              const auto r_new = rotate(Vector3{r[0], r[1], r[2]});
+              m.nodes[node].coordinates = {{r_new.X(), r_new.Y(), r_new.Z()}};
+            }
+          }
+          Prepared p{std::move(m), parallel};
+          ti::ThermalPartition partition;
+          for (std::size_t i = 0; i < 20; ++i)
+            partition.constrained.push_back(i);
+          partition.prescribed = values;
+          fesa::MklPardisoSolver solver;
+          ti::ThermalFactorization cache;
+          const auto result = Take(ti::SolveThermalIncrement(
+              {p.assembly.Conductivity(), p.assembly.Capacity(),
+               p.assembly.HeatInput()},
+              partition, values, 1, solver, cache));
+          EXPECT_FALSE(cache.matrix.has_value());
+          const auto flux =
+              Take(p.assembly.RecoverFlux(result.temperature, parallel));
+          Vector3 expected{-3 * k, 2 * k, -4 * k};
+          if (rotated) expected = rotate(expected);
+          double scale = 0;
+          for (std::size_t i = 0; i < 20; ++i)
+            scale = std::max(scale, std::abs(values[i]));
+          const double budget = (Gamma(2048) * I{k * (8 + 16 / h) * scale}).hi;
+          for (std::size_t point = 0; point < 4; ++point) {
+            const auto& r = p.state.FluxLocalFrames()[0][point];
+            for (std::size_t a = 0; a < 3; ++a)
+              for (std::size_t b = 0; b < 3; ++b) {
+                I dot, absolute;
+                for (std::size_t c = 0; c < 3; ++c) {
+                  dot = dot + I{r[a][c]} * I{r[b][c]};
+                  absolute = absolute + I{std::abs(r[a][c] * r[b][c])};
+                }
+                EXPECT_LE((dot - I{static_cast<double>(a == b)}).Abs(),
+                          (Gamma(2048) * absolute).hi);
+              }
+            EXPECT_GT(r[0].Cross(r[1]).Dot(r[2]), 0);
+            for (std::size_t trace = 0; trace < 6; ++trace) {
+              const auto& q = flux[0][6 * point + trace];
+              for (std::size_t c = 0; c < 3; ++c)
+                EXPECT_LE(std::abs(q[c] - expected[c]), budget);
+              for (std::size_t axis = 0; axis < 3; ++axis) {
+                I actual_local, expected_local;
+                for (std::size_t c = 0; c < 3; ++c) {
+                  actual_local = actual_local + I{r[axis][c]} * I{q[c]};
+                  expected_local =
+                      expected_local + I{r[axis][c]} * I{expected[c]};
+                }
+                EXPECT_LE((actual_local - expected_local).Abs(),
+                          budget + Gamma(16).hi * k * (8 + 16 / h) * scale);
+              }
+            }
+          }
+          // The discontinuous derivative at the shared station has two distinct
+          // one-sided traces; prescribed interior stations exist only in this
+          // private algebra fixture, never in the public Domain BC definition.
+          for (std::size_t node = 0; node < 4; ++node) {
+            for (std::size_t station = 0; station < 5; ++station)
+              values[5 * node + station] =
+                  station < 3 ? 0. : static_cast<double>(station - 2);
+          }
+          partition.prescribed = values;
+          const auto trace_result = Take(ti::SolveThermalIncrement(
+              {p.assembly.Conductivity(), p.assembly.Capacity(),
+               p.assembly.HeatInput()},
+              partition, values, 1, solver, cache));
+          const auto traces =
+              Take(p.assembly.RecoverFlux(trace_result.temperature, parallel));
+          for (std::size_t point = 0; point < 4; ++point) {
+            const auto& normal = p.state.FluxLocalFrames()[0][point][2];
+            EXPECT_EQ(traces[0][6 * point + 2].Dot(normal), 0);
+            EXPECT_LE(
+                std::abs(traces[0][6 * point + 3].Dot(normal) + 4 * k / h),
+                (Gamma(2048) * I{2 * k * (8 + 16 / h)}).hi);
+          }
+        }
+}
+
+I ThicknessMass(const std::array<I, 5>& v) {
+  I result;
+  for (std::size_t segment = 0; segment < 2; ++segment) {
+    const I m = (v[2 * segment] + v[2 * segment + 2]) / I{2};
+    const I d = (v[2 * segment + 2] - v[2 * segment]) / I{2};
+    const I middle = v[2 * segment + 1];
+    // Positive square decomposition of the exact two-P2 consistent mass.
+    result = result + (I{6} * Square(m + I{2} / I{3} * middle) +
+                       I{40} / I{3} * Square(middle) + I{10} * Square(d)) /
+                          I{60};
+  }
+  return result;
+}
+I ThicknessGradient(const std::array<I, 5>& v) {
+  I result;
+  for (std::size_t segment = 0; segment < 2; ++segment) {
+    const I left = v[2 * segment] - v[2 * segment + 1],
+            right = v[2 * segment + 2] - v[2 * segment + 1];
+    result = result +
+             I{2} / I{3} *
+                 (I{6} * (Square(left) + Square(right)) + Square(left + right));
+  }
+  return result;
+}
+struct SpatialError {
+  I temperature, gradient;
+};
+template <std::size_t Order>
+SpatialError IntegrateSpatialError(
+    const Vector& values, std::size_t n,
+    const std::array<thermal_oracle::GaussPoint, Order>& rule, I amplitude) {
+  I temperature, gradient;
+  for (std::size_t y = 0; y < n; ++y)
+    for (std::size_t x = 0; x < n; ++x) {
+      const std::size_t node = y * (n + 1) + x;
+      const std::array<std::size_t, 4> nodes{node, node + 1, node + n + 2,
+                                             node + n + 1};
+      for (const auto& gx : rule) {
+        const I u = (I{1} + I{gx.x_lo, gx.x_hi}) / I{2};
+        const I position =
+            (I{static_cast<double>(x)} + u) / I{static_cast<double>(n)};
+        const I phase = Pi() * position;
+        const I exact = Cos(phase) * amplitude;
+        const I exact_gradient =
+            I{-1} * Pi() * Cos(Pi() / I{2} - phase) * amplitude;
+        for (const auto& gy : rule) {
+          const I v = (I{1} + I{gy.x_lo, gy.x_hi}) / I{2};
+          const I weight =
+              I{gx.w_lo, gx.w_hi} * I{gy.w_lo, gy.w_hi} / I{4. * n * n};
+          std::array<I, 5> field, dx, dy, error;
+          for (std::size_t a = 0; a < 5; ++a) {
+            const I a0{values[5 * nodes[0] + a]}, a1{values[5 * nodes[1] + a]},
+                a2{values[5 * nodes[2] + a]}, a3{values[5 * nodes[3] + a]};
+            field[a] = (I{1} - v) * ((I{1} - u) * a0 + u * a1) +
+                       v * ((I{1} - u) * a3 + u * a2);
+            error[a] = field[a] - exact;
+            dx[a] = I{static_cast<double>(n)} *
+                        ((I{1} - v) * (a1 - a0) + v * (a2 - a3)) -
+                    exact_gradient;
+            dy[a] = I{static_cast<double>(n)} *
+                    ((I{1} - u) * (a3 - a0) + u * (a2 - a1));
+          }
+          temperature = temperature + weight * ThicknessMass(error);
+          gradient =
+              gradient + weight * (ThicknessMass(dx) + ThicknessMass(dy) +
+                                   ThicknessGradient(field));
+        }
+      }
+    }
+  return {Sqrt(temperature), Sqrt(gradient)};
+}
+TEST(ThermalPhysics,
+     FullTwoDimensionalMeshesHaveSecondOrderTemperatureAndFirstOrderGradient) {
+  fesa::SerialParallelFor parallel;
+  constexpr double period = .001;
+  constexpr std::uint32_t count = 32;
+  const I continuum_amplitude = ExpNegative(Pi() * Pi() * I{period});
+  double previous_temperature = 0, previous_gradient = 0;
+  for (std::size_t n : {8U, 16U, 32U, 64U}) {
+    SCOPED_TRACE(n);
+    std::vector<I> mode(n + 1);
+    double initial_error = 0;
+    for (std::size_t i = 0; i <= n; ++i) {
+      mode[i] =
+          Cos(Pi() * I{static_cast<double>(i)} / I{static_cast<double>(n)});
+      initial_error =
+          std::max(initial_error, (mode[i] - I{mode[i].Mid()}).Abs());
+    }
+    Prepared p{Mesh(n, period, count,
+                    [&](std::size_t i, auto) { return mode[i].Mid(); }),
+               parallel};
+    const auto metrics =
+        Advance(p, parallel, Down(1. / (2160 * n * n)), initial_error);
+    const auto& actual = p.state.CommittedState().temperature;
+    const auto twelve = IntegrateSpatialError(
+        actual, n, thermal_oracle::kGauss12, continuum_amplitude);
+    const auto sixteen = IntegrateSpatialError(
+        actual, n, thermal_oracle::kGauss16, continuum_amplitude);
+    EXPECT_LE((twelve.temperature - sixteen.temperature).Abs(),
+              sixteen.temperature.lo * 1e-8);
+    EXPECT_LE((twelve.gradient - sixteen.gradient).Abs(),
+              sixteen.gradient.lo * 1e-8);
+    const I lambda = I{6. * n * n} * (I{1} - mode[1]) / (I{2} + mode[1]);
+    const I semidiscrete = ExpNegative(lambda * I{period});
+    const I be = Power(I{1} / (I{1} + lambda * I{period / count}), count);
+    Vector spatial{actual.Size()}, difference{actual.Size()},
+        mode_values{actual.Size()};
+    double oracle_width = 0;
+    for (std::size_t node = 0; node < (n + 1) * (n + 1); ++node)
+      for (std::size_t a = 0; a < 5; ++a) {
+        const auto i = node % (n + 1);
+        const I expected = mode[i] * be;
+        spatial[5 * node + a] = (mode[i] * semidiscrete).Mid();
+        mode_values[5 * node + a] = mode[i].Mid();
+        difference[5 * node + a] = actual[5 * node + a] - expected.Mid();
+        oracle_width =
+            std::max(oracle_width, (expected - I{expected.Mid()}).Abs());
+      }
+    const I agreement = CapacityNorm(difference, n);
+    const I forward = I{metrics.capacity_budget} + I{oracle_width};
+    EXPECT_LE(agreement.hi, forward.hi);
+    EXPECT_LE(forward.hi, sixteen.temperature.lo / 1000);
+    EXPECT_LE((Sqrt(I{24. * n * n + 1280}) * forward).hi,
+              sixteen.gradient.lo / 1000);
+    const auto spatial_error = IntegrateSpatialError(
+        spatial, n, thermal_oracle::kGauss12, continuum_amplitude);
+    const I temporal =
+        I{(be - semidiscrete).Abs()} * CapacityNorm(mode_values, n);
+    EXPECT_LE(temporal.hi, spatial_error.temperature.lo / 32);
+    double temperature_order = 0, gradient_order = 0;
+    if (previous_temperature > 0) {
+      temperature_order =
+          std::log2(previous_temperature / sixteen.temperature.Mid());
+      gradient_order = std::log2(previous_gradient / sixteen.gradient.Mid());
+      EXPECT_GE(temperature_order, 1.95);
+      EXPECT_LE(temperature_order, 2.05);
+      EXPECT_GE(gradient_order, .98);
+      EXPECT_LE(gradient_order, 1.03);
+    }
+    previous_temperature = sixteen.temperature.Mid();
+    previous_gradient = sixteen.gradient.Mid();
+    Print("space", n, count, metrics);
+    std::cout << std::setprecision(17) << "THERMAL_SPACE n=" << n
+              << " T_error=" << sixteen.temperature.Mid()
+              << " gradient_error=" << sixteen.gradient.Mid()
+              << " T_order=" << temperature_order
+              << " gradient_order=" << gradient_order << " time_spatial_ratio="
+              << (temporal / spatial_error.temperature).hi
+              << " quadrature_T_gap="
+              << (twelve.temperature - sixteen.temperature).Abs()
+              << " quadrature_gradient_gap="
+              << (twelve.gradient - sixteen.gradient).Abs()
+              << " modal_agreement=" << agreement.hi
+              << " roundoff_budget=" << forward.hi << '\n';
+    if (n == 64)
+      EXPECT_EQ(8 * (metrics.temperature_values + metrics.flux_values),
+                81074472U);
+  }
+}
+I Bracket(thermal_oracle::Bracket x) { return {x.lo, x.hi}; }
+class ThermalSlab
+    : public testing::TestWithParam<std::pair<int, std::uint32_t>> {};
+TEST_P(ThermalSlab, ActualMeshesMatchCertifiedP2AndSeparateFourierError) {
+  const auto [case_index, count] = GetParam();
+  const auto& oracle = thermal_oracle::kSlab[case_index];
+  const auto& be = count == 4096 ? oracle.be4096 : oracle.be8192;
+  fesa::SerialParallelFor parallel;
+  std::array<I, 13> first_mesh;
+  std::array<double, 4> first_budget{};
+  for (std::size_t n : {1U, 2U}) {
+    auto definition =
+        Mesh(n, oracle.time, count, [](auto, auto) { return 0.; });
+    for (std::size_t node = 0; node < definition.nodes.size(); ++node)
+      definition.thermal_model_boundaries.push_back(
+          {static_cast<fesa::EntityIndex>(node), 11, 1, {}});
+    Prepared p{std::move(definition), parallel};
+    const double mass_lower = Down(1. / (2160 * n * n));
+    for (std::size_t row = 0; row < p.state.CommittedState().temperature.Size();
+         ++row)
+      EXPECT_EQ(p.state.CommittedState().temperature[row], 0);
+    const auto metrics = Advance(p, parallel, mass_lower);
+    const auto& final = p.state.CommittedState();
+    // Independently certify the small last-step inverse as well as the
+    // structured C-norm propagation used for the full 4096/8192-step history.
+    const auto& c = p.assembly.Capacity();
+    const auto& k = p.assembly.Conductivity();
+    const double delta = final.time_increment;
+    std::vector<fesa::CooContribution> coo;
+    double operator_error = 0;
+    for (std::size_t row = 0; row < c.Rows(); ++row) {
+      I error;
+      for (auto entry = c.RowOffsets()[row]; entry < c.RowOffsets()[row + 1];
+           ++entry) {
+        const double center = c.Values()[entry] / delta + k.Values()[entry];
+        coo.push_back({row, c.ColumnIndices()[entry], center, 0, entry});
+        const I exact = (I{c.Values()[entry]} +
+                         I{-p.assembly.CapacityError().Values()[entry],
+                           p.assembly.CapacityError().Values()[entry]}) /
+                            I{delta} +
+                        I{k.Values()[entry]} +
+                        I{-p.assembly.ConductivityError().Values()[entry],
+                          p.assembly.ConductivityError().Values()[entry]};
+        error = error + I{(exact - I{center}).Abs()};
+      }
+      operator_error = std::max(operator_error, error.hi);
+    }
+    const auto a =
+        Take(fesa::SparseMatrix::FromCoo(c.Rows(), c.Columns(), std::move(coo),
+                                         {c.RowOffsets(), c.ColumnIndices()}));
+    const double inverse = CertifiedInverse(a, p.dofs.FreeDofs());
+    const I inverse_denominator = I{1} - I{inverse} * I{operator_error};
+    ASSERT_GT(inverse_denominator.lo, 0);
+    const I local_forward =
+        I{inverse} / inverse_denominator * I{metrics.last_free_residual_inf};
+    const I temperature_budget =
+        I{metrics.capacity_budget} / Sqrt(I{mass_lower});
+    double temperature_scale = 0;
+    for (std::size_t row = 0; row < final.temperature.Size(); ++row)
+      temperature_scale =
+          std::max(temperature_scale, std::abs(final.temperature[row]));
+    const I flux_budget = I{8. * n + 16} * (temperature_budget +
+                                            Gamma(2048) * I{temperature_scale});
+    I energy, reaction;
+    constexpr int thickness_weight[5] = {1, 4, 2, 4, 1};
+    for (std::size_t j = 0; j <= n; ++j)
+      for (std::size_t i = 0; i <= n; ++i) {
+        const int plane_weight =
+            (i == 0 || i == n ? 1 : 2) * (j == 0 || j == n ? 1 : 2);
+        for (std::size_t a = 0; a < 5; ++a)
+          energy =
+              energy +
+              I{final.temperature[5 * (j * (n + 1) + i) + a]} *
+                  I{static_cast<double>(plane_weight * thickness_weight[a])} /
+                  I{48. * n * n};
+      }
+    for (std::size_t row = 0; row < final.completed->reaction_heat_rate.Size();
+         ++row)
+      reaction = reaction + I{final.completed->reaction_heat_rate[row]};
+    // Exact P2 Schur reaction row L = (-400/51,56/17,-80/51,28/51),
+    // L Cff^-1 L^T = 6447104/14739. Plane integration is a contraction
+    // in the product C norm; Ccf Cff^-1 has infinity norm 3/17.
+    const I reaction_budget =
+        Sqrt(I{6447104} / I{14739}) * I{metrics.capacity_budget} +
+        (I{3} / I{17}) * I{metrics.last_free_residual_l1} +
+        I{metrics.last_reaction_evaluation_error} +
+        I{(reaction - I{reaction.Mid()}).Abs()};
+    const std::array<double, 4> solver_budget{
+        temperature_budget.hi, flux_budget.hi,
+        (I{metrics.capacity_budget} + I{(energy - I{energy.Mid()}).Abs()}).hi,
+        reaction_budget.hi};
+    std::array<double, 4> agreement{}, physical{}, oracle_budget{},
+        time_error{}, space_error{};
+    const auto observe = [&](std::size_t quantity, std::size_t group, I value) {
+      agreement[group] =
+          std::max(agreement[group], (value - Bracket(be[quantity])).Abs());
+      physical[group] = std::max(
+          physical[group], (value - Bracket(oracle.fourier[quantity])).Abs());
+    };
+    for (std::size_t quantity = 0; quantity < 13; ++quantity) {
+      const std::size_t group = quantity < 5    ? 0
+                                : quantity < 11 ? 1
+                                                : quantity - 9;
+      oracle_budget[group] = std::max(
+          oracle_budget[group],
+          (I{be[quantity].hi} - I{be[quantity].lo} +
+           I{oracle.continuous[quantity].hi} -
+           I{oracle.continuous[quantity].lo} + I{oracle.fourier[quantity].hi} -
+           I{oracle.fourier[quantity].lo})
+              .hi);
+      time_error[group] = std::max(
+          time_error[group],
+          (Bracket(be[quantity]) - Bracket(oracle.continuous[quantity])).Abs());
+      space_error[group] =
+          std::max(space_error[group], (Bracket(oracle.continuous[quantity]) -
+                                        Bracket(oracle.fourier[quantity]))
+                                           .Abs());
+    }
+    for (std::size_t row = 0; row < final.temperature.Size(); ++row)
+      observe(row % 5, 0, I{final.temperature[row]});
+    for (const auto& element : final.completed->global_heat_flux)
+      for (std::size_t sample = 0; sample < 24; ++sample) {
+        EXPECT_LE(std::abs(element[sample][0]), flux_budget.hi);
+        EXPECT_LE(std::abs(element[sample][1]), flux_budget.hi);
+        observe(5 + sample % 6, 1, I{element[sample][2]});
+      }
+    observe(11, 2, energy);
+    observe(12, 3, reaction);
+    std::array<I, 13> samples;
+    for (std::size_t a = 0; a < 5; ++a) samples[a] = I{final.temperature[a]};
+    for (std::size_t trace = 0; trace < 6; ++trace)
+      samples[5 + trace] = I{final.completed->global_heat_flux[0][trace][2]};
+    samples[11] = energy;
+    samples[12] = reaction;
+    if (n == 1) {
+      first_mesh = samples;
+      first_budget = solver_budget;
+    } else {
+      for (std::size_t quantity = 0; quantity < 13; ++quantity) {
+        const std::size_t group = quantity < 5    ? 0
+                                  : quantity < 11 ? 1
+                                                  : quantity - 9;
+        EXPECT_LE((samples[quantity] - first_mesh[quantity]).Abs(),
+                  (I{first_budget[group]} + I{solver_budget[group]}).hi);
+      }
+    }
+    Print("slab", n, count, metrics);
+    for (std::size_t group = 0; group < 4; ++group) {
+      const I roundoff = I{solver_budget[group]} + I{oracle_budget[group]};
+      EXPECT_LE(roundoff.hi, oracle.space[group] / 1000);
+      EXPECT_LE(agreement[group], roundoff.hi);
+      EXPECT_LE(physical[group],
+                (I{oracle.space[group]} + I{time_error[group]} + roundoff).hi);
+      EXPECT_LE((I{physical[group]} - I{space_error[group]}).Abs(),
+                (I{time_error[group]} + roundoff).hi);
+      const bool plateau =
+          (I{time_error[group]} + roundoff).hi <= space_error[group] / 10;
+      if (case_index < 4 && count == 8192) EXPECT_TRUE(plateau);
+      std::cout << std::setprecision(17) << "THERMAL_SLAB time=" << oracle.time
+                << " n=" << n << " N=" << count << " quantity=" << group
+                << " finite_BE_gap=" << agreement[group]
+                << " solver_bound=" << solver_budget[group]
+                << " oracle_bound=" << oracle_budget[group]
+                << " physical_error=" << physical[group]
+                << " space_error=" << space_error[group]
+                << " D_space=" << oracle.space[group]
+                << " D_time=" << time_error[group] << " plateau=" << plateau
+                << " last_step_inverse_bound=" << inverse
+                << " last_step_forward_bound=" << local_forward.hi
+                << " penetration=" << std::sqrt(oracle.time)
+                << " station_spacing=0.25 segment_length=0.5" << '\n';
+    }
+  }
+}
+INSTANTIATE_TEST_SUITE_P(
+    FixedStations, ThermalSlab,
+    testing::Values(std::make_pair(0, 4096U), std::make_pair(0, 8192U),
+                    std::make_pair(1, 4096U), std::make_pair(1, 8192U),
+                    std::make_pair(2, 4096U), std::make_pair(2, 8192U),
+                    std::make_pair(3, 4096U), std::make_pair(3, 8192U),
+                    std::make_pair(4, 4096U), std::make_pair(4, 8192U),
+                    std::make_pair(5, 4096U), std::make_pair(5, 8192U)));
+}  // namespace
diff --git a/tests/integration/analysis/thermal_slab_oracle.h b/tests/integration/analysis/thermal_slab_oracle.h
new file mode 100644
index 0000000..01577f9
--- /dev/null
+++ b/tests/integration/analysis/thermal_slab_oracle.h
@@ -0,0 +1,359 @@
+// Generated by thermal_oracle_generator.py. Outward binary64 brackets.
+#ifndef FESA_TESTS_THERMAL_SLAB_ORACLE_H_
+#define FESA_TESTS_THERMAL_SLAB_ORACLE_H_
+#include <array>
+namespace thermal_oracle {
+struct Bracket {
+  double lo, hi;
+};
+struct SlabCase {
+  double time;
+  std::array<double, 4> space;
+  std::array<Bracket, 13> fourier, continuous, be4096, be8192;
+};
+inline constexpr std::array<SlabCase, 6> kSlab{{
+    {
+        0x1.0624dd2f1a9fcp-10,
+        {{0x1.20c49ba5e353fp-3, 0x1.5666666666666p+3, 0x1.f212d77318fc5p-6,
+          0x1.5cccccccccccdp+3}},
+        {{{0x1.fffffffffffffp-1, 0x1.0000000000001p+0},
+          {0x1.85b89453c09acp-26, 0x1.85b89453c09aep-26},
+          {0x1.fe4543738e03ap-95, 0x1.0ff4c2f026ce0p-94},
+          {-0x1.64845561fdc38p-101, 0x1.c12311737a124p-99},
+          {-0x1.68eca14e83537p-101, 0x1.c008fe7858ae5p-99},
+          {0x1.1d75b94accf81p+4, 0x1.1d75b94accf83p+4},
+          {0x1.88168e094c2d6p-19, 0x1.88168e094c2d8p-19},
+          {0x1.e433884766565p-87, 0x1.02b7b8a2bf31ap-86},
+          {0x1.e433884766565p-87, 0x1.02b7b8a2bf31ap-86},
+          {-0x1.3a94ff682132fp-91, 0x1.b25320f2bf337p-92},
+          {-0x1.09df47f0c0666p-91, 0x1.09df47f0c0666p-91},
+          {0x1.244f96d3c38d6p-5, 0x1.244f96d3c38d8p-5},
+          {0x1.1d75b94accf81p+4, 0x1.1d75b94accf83p+4}}},
+        {{{0x1.fffffffffffffp-1, 0x1.0000000000001p+0},
+          {-0x1.d37c7f414f866p-4, -0x1.d37c7f414f859p-4},
+          {0x1.1eb81bc756156p-3, 0x1.1eb81bc75615ap-3},
+          {-0x1.17e5fa341d20bp-6, -0x1.17e5fa341d110p-6},
+          {0x1.2d98feedf9757p-5, 0x1.2d98feedf9788p-5},
+          {0x1.cc5b11a49f520p+2, 0x1.cc5b11a49f523p+2},
+          {0x1.b851f90e2a7a9p+0, 0x1.b851f90e2a7abp+0},
+          {-0x1.e0642a3b1429cp+1, -0x1.e0642a3b14296p+1},
+          {0x1.0ce0640ae3c9bp+0, 0x1.0ce0640ae3cbbp+0},
+          {0x1.a6a3b817af6eap-3, 0x1.a6a3b817af708p-3},
+          {-0x1.466eec09efdf1p-1, -0x1.466eec09efdc2p-1},
+          {0x1.0e51bf4794d1ep-4, 0x1.0e51bf4794d3ap-4},
+          {0x1.be34bece7dbcbp+2, 0x1.be34bece7dbd0p+2}}},
+        {{{0x1.fffffffffffffp-1, 0x1.0000000000001p+0},
+          {-0x1.d37cd7bcda5f0p-4, -0x1.d37cd7bcda5e4p-4},
+          {0x1.1eb8818ac9073p-3, 0x1.1eb8818ac9077p-3},
+          {-0x1.17e77bedcde83p-6, -0x1.17e77bedcdd88p-6},
+          {0x1.2d9a6a607059ep-5, 0x1.2d9a6a60705d0p-5},
+          {0x1.cc5b231047dc3p+2, 0x1.cc5b231047dc6p+2},
+          {0x1.b851df9d4dbe2p+0, 0x1.b851df9d4dbe4p+0},
+          {-0x1.e064668341fa9p+1, -0x1.e064668341fa4p+1},
+          {0x1.0ce0f74bd7864p+0, 0x1.0ce0f74bd7884p+0},
+          {0x1.a6a3cde559e00p-3, 0x1.a6a3cde559e1dp-3},
+          {-0x1.467007a5021f8p-1, -0x1.467007a5021c9p-1},
+          {0x1.0e51b2b4d8688p-4, 0x1.0e51b2b4d86a3p-4},
+          {0x1.be34d80c5223ap+2, 0x1.be34d80c5223fp+2}}},
+        {{{0x1.fffffffffffffp-1, 0x1.0000000000001p+0},
+          {-0x1.d37cab7f333f2p-4, -0x1.d37cab7f333e5p-4},
+          {0x1.1eb84ea93c6aep-3, 0x1.1eb84ea93c6b2p-3},
+          {-0x1.17e6bb11c7d2ep-6, -0x1.17e6bb11c7c33p-6},
+          {0x1.2d99b4a80421ep-5, 0x1.2d99b4a80424fp-5},
+          {0x1.cc5b1a5a7a2e7p+2, 0x1.cc5b1a5a7a2eap+2},
+          {0x1.b851ec55b0e53p+0, 0x1.b851ec55b0e55p+0},
+          {-0x1.e064485f43780p+1, -0x1.e064485f4377ap+1},
+          {0x1.0ce0adaba68aep+0, 0x1.0ce0adaba68cdp+0},
+          {0x1.a6a3c2fe76c36p-3, 0x1.a6a3c2fe76c53p-3},
+          {-0x1.466f79d811b70p-1, -0x1.466f79d811b41p-1},
+          {0x1.0e51b8fe3293cp-4, 0x1.0e51b8fe32958p-4},
+          {0x1.be34cb6d72726p+2, 0x1.be34cb6d7272bp+2}}},
+    },
+    {
+        0x1.47ae147ae147bp-7,
+        {{0x1.f559b3d07c84bp-7, 0x1.47ae147ae147bp+0, 0x1.32b55ef1fddecp-8,
+          0x1.a5e353f7ced91p-1}},
+        {{{0x1.fffffffffffffp-1, 0x1.0000000000001p+0},
+          {0x1.3bcd133aa0ffbp-4, 0x1.3bcd133aa0ffdp-4},
+          {0x1.aab859b20ac9ep-12, 0x1.aab859b20aca0p-12},
+          {0x1.e87470e505fc1p-24, 0x1.e87470e505fc3p-24},
+          {0x1.b0c1a759f773cp-39, 0x1.b0c1a759f773ep-39},
+          {0x1.6914d24534247p+2, 0x1.6914d24534249p+2},
+          {0x1.2ebf3dcc9f22ep+0, 0x1.2ebf3dcc9f230p+0},
+          {0x1.64e3dcd3af4bap-7, 0x1.64e3dcd3af4bcp-7},
+          {0x1.64e3dcd3af4bap-7, 0x1.64e3dcd3af4bcp-7},
+          {0x1.27c272dde935ep-18, 0x1.27c272dde9360p-18},
+          {-0x1.64c3500000001p-300, 0x1.64c3500000001p-300},
+          {0x1.ce2f21a042be1p-4, 0x1.ce2f21a042be3p-4},
+          {0x1.6914d24534247p+2, 0x1.6914d24534249p+2}}},
+        {{{0x1.fffffffffffffp-1, 0x1.0000000000001p+0},
+          {0x1.7a435a2984e32p-4, 0x1.7a435a2984e3ap-4},
+          {0x1.9dc400057caf6p-12, 0x1.9dc400057d765p-12},
+          {0x1.ad9feeec8c48dp-7, 0x1.ad9feeec8c649p-7},
+          {-0x1.986d06319bd56p-7, -0x1.986d06319bcd8p-7},
+          {0x1.50c482dacf8f7p+2, 0x1.50c482dacf8f9p+2},
+          {0x1.ffcc477fff504p+0, 0x1.ffcc477fff508p+0},
+          {-0x1.43797c6b3f9d5p+0, -0x1.43797c6b3f9d2p+0},
+          {-0x1.05044c3c6933ap-3, -0x1.05044c3c69290p-3},
+          {0x1.a55b2631c7b93p-6, 0x1.a55b2631c7baep-6},
+          {0x1.6e5b15c8db17bp-3, 0x1.6e5b15c8db21fp-3},
+          {0x1.e15645f1d96e9p-4, 0x1.e15645f1d9703p-4},
+          {0x1.345d88b145b6ap+2, 0x1.345d88b145b6ep+2}}},
+        {{{0x1.fffffffffffffp-1, 0x1.0000000000001p+0},
+          {0x1.7a374829b1404p-4, 0x1.7a374829b140cp-4},
+          {0x1.af3900931af0ap-12, 0x1.af3900931bb78p-12},
+          {0x1.ad81aefddf88dp-7, 0x1.ad81aefddfa49p-7},
+          {-0x1.983e6ad119e83p-7, -0x1.983e6ad119e05p-7},
+          {0x1.50c690c2ce70cp+2, 0x1.50c690c2ce70ep+2},
+          {0x1.ffca18dfed9c8p+0, 0x1.ffca18dfed9cbp+0},
+          {-0x1.4386114b5e8a2p+0, -0x1.4386114b5e89ep+0},
+          {-0x1.04baf9d759bb3p-3, -0x1.04baf9d759b09p-3},
+          {0x1.a5b832d5b2be1p-6, 0x1.a5b832d5b2bfcp-6},
+          {0x1.6e29068cc6608p-3, 0x1.6e29068cc66acp-3},
+          {0x1.e154626664a52p-4, 0x1.e154626664a6cp-4},
+          {0x1.34600de90e320p+2, 0x1.34600de90e323p+2}}},
+        {{{0x1.fffffffffffffp-1, 0x1.0000000000001p+0},
+          {0x1.7a3d511a6a009p-4, 0x1.7a3d511a6a012p-4},
+          {0x1.a67e9a3b3c161p-12, 0x1.a67e9a3b3cdd0p-12},
+          {0x1.ad90cec9c0f6fp-7, 0x1.ad90cec9c112bp-7},
+          {-0x1.9855b85b10049p-7, -0x1.9855b85b0ffcbp-7},
+          {0x1.50c589d18499bp+2, 0x1.50c589d18499ep+2},
+          {0x1.ffcb302cb8985p+0, 0x1.ffcb302cb8989p+0},
+          {-0x1.437fc6eca1366p+0, -0x1.437fc6eca1362p+0},
+          {-0x1.04dfa2a190d26p-3, -0x1.04dfa2a190c7cp-3},
+          {0x1.a589ad2ce9e39p-6, 0x1.a589ad2ce9e54p-6},
+          {0x1.6e420deccb411p-3, 0x1.6e420deccb4b5p-3},
+          {0x1.e1555429f7bf1p-4, 0x1.e1555429f7c0bp-4},
+          {0x1.345ecb5088dc7p+2, 0x1.345ecb5088dcbp+2}}},
+    },
+    {
+        0x1.999999999999ap-5,
+        {{0x1.87e7c06e19b91p-8, 0x1.03126e978d4fep-2, 0x1.54c985f06f694p-10,
+          0x1.9e30014f8b589p-8}},
+        {{{0x1.fffffffffffffp-1, 0x1.0000000000001p+0},
+          {0x1.b77efb25b8dd2p-2, 0x1.b77efb25b8dd4p-2},
+          {0x1.d252b304ed681p-4, 0x1.d252b304ed683p-4},
+          {0x1.235c8909400c5p-6, 0x1.235c8909400c7p-6},
+          {0x1.9a5c5e03484dep-9, 0x1.9a5c5e03484e0p-9},
+          {0x1.42f6019270650p+1, 0x1.42f6019270652p+1},
+          {0x1.d89101a46e842p+0, 0x1.d89101a46e844p+0},
+          {0x1.721a47b716779p-1, 0x1.721a47b71677bp-1},
+          {0x1.721a47b716779p-1, 0x1.721a47b71677bp-1},
+          {0x1.343beb6102cc6p-3, 0x1.343beb6102cc8p-3},
+          {-0x1.3f24600000001p-301, 0x1.3f24600000001p-301},
+          {0x1.025e67b9cdad7p-2, 0x1.025e67b9cdad9p-2},
+          {0x1.42f6019270650p+1, 0x1.42f6019270652p+1}}},
+        {{{0x1.fffffffffffffp-1, 0x1.0000000000001p+0},
+          {0x1.b8566cd59dd36p-2, 0x1.b8566cd59dd3ap-2},
+          {0x1.b9d5c5e3870abp-4, 0x1.b9d5c5e3870c1p-4},
+          {0x1.1587ab02cfe81p-6, 0x1.1587ab02cff38p-6},
+          {0x1.1ad4cd8131fbap-7, 0x1.1ad4cd8131fd0p-7},
+          {0x1.6346ef889a9d2p+1, 0x1.6346ef889a9d4p+1},
+          {0x1.c8c547438f1e7p+0, 0x1.c8c547438f1ecp+0},
+          {0x1.95f95eebd2052p-1, 0x1.95f95eebd205ep-1},
+          {0x1.0ed51015fadbfp-1, 0x1.0ed51015fade0p-1},
+          {0x1.967b2c3360cb1p-3, 0x1.967b2c3360cc9p-3},
+          {-0x1.0e5de7f129e18p-3, -0x1.0e5de7f129d6ep-3},
+          {0x1.010a6df46a07bp-2, 0x1.010a6df46a082p-2},
+          {0x1.43c4ee9a8ce2bp+1, 0x1.43c4ee9a8ce2fp+1}}},
+        {{{0x1.fffffffffffffp-1, 0x1.0000000000001p+0},
+          {0x1.b84e9bd4911f6p-2, 0x1.b84e9bd4911fap-2},
+          {0x1.b9ccf55bd5095p-4, 0x1.b9ccf55bd50abp-4},
+          {0x1.15b1a7159d30ap-6, 0x1.15b1a7159d3c1p-6},
+          {0x1.1b178ee6bb5d6p-7, 0x1.1b178ee6bb5ecp-7},
+          {0x1.634e33812c311p+1, 0x1.634e33812c313p+1},
+          {0x1.c8c66154855eap+0, 0x1.c8c66154855eep+0},
+          {0x1.95e0b74d64b63p-1, 0x1.95e0b74d64b6ep-1},
+          {0x1.0ec60ab6ae53dp-1, 0x1.0ec60ab6ae55ep-1},
+          {0x1.966a037efd9d7p-3, 0x1.966a037efd9f0p-3},
+          {-0x1.0e4423dcbe1c4p-3, -0x1.0e4423dcbe11ap-3},
+          {0x1.0108815a6658cp-2, 0x1.0108815a66593p-2},
+          {0x1.43cb51041b9e3p+1, 0x1.43cb51041b9e7p+1}}},
+        {{{0x1.fffffffffffffp-1, 0x1.0000000000001p+0},
+          {0x1.b85284545e3fep-2, 0x1.b85284545e402p-2},
+          {0x1.b9d15d4ee5b37p-4, 0x1.b9d15d4ee5b4dp-4},
+          {0x1.159ca92438c35p-6, 0x1.159ca92438cebp-6},
+          {0x1.1af6311dbcc58p-7, 0x1.1af6311dbcc6ep-7},
+          {0x1.634a9180901b3p+1, 0x1.634a9180901b5p+1},
+          {0x1.c8c5d45623495p+0, 0x1.c8c5d4562349ap+0},
+          {0x1.95ed0b564cb88p-1, 0x1.95ed0b564cb94p-1},
+          {0x1.0ecd8d3b0bfa0p-1, 0x1.0ecd8d3b0bfc1p-1},
+          {0x1.9672972b2e1aap-3, 0x1.9672972b2e1c2p-3},
+          {-0x1.0e510695d3bacp-3, -0x1.0e510695d3b02p-3},
+          {0x1.010977a63dfebp-2, 0x1.010977a63dff1p-2},
+          {0x1.43c81fc9071c0p+1, 0x1.43c81fc9071c4p+1}}},
+    },
+    {
+        0x1.999999999999ap-4,
+        {{0x1.6c61522a6f3f5p-10, 0x1.16872b020c49cp-3, 0x1.5d2d01c0ca601p-11,
+          0x1.38476f2a5a46ap-7}},
+        {{{0x1.fffffffffffffp-1, 0x1.0000000000001p+0},
+          {0x1.270906ee5a6cfp-1, 0x1.270906ee5a6d1p-1},
+          {0x1.0eb16bef11388p-2, 0x1.0eb16bef1138ap-2},
+          {0x1.945c98b83c928p-4, 0x1.945c98b83c92ap-4},
+          {0x1.9f4a5c2b49fb7p-5, 0x1.9f4a5c2b49fb9p-5},
+          {0x1.c8b1bdc9d3f27p+0, 0x1.c8b1bdc9d3f29p+0},
+          {0x1.8673042a10047p+0, 0x1.8673042a10049p+0},
+          {0x1.e5a6c4f8749dfp-1, 0x1.e5a6c4f8749e1p-1},
+          {0x1.e5a6c4f8749dfp-1, 0x1.e5a6c4f8749e1p-1},
+          {0x1.9af613adf14c0p-2, 0x1.9af613adf14c2p-2},
+          {-0x1.c357400000001p-302, 0x1.c357400000001p-302},
+          {0x1.6d631d0d902e4p-2, 0x1.6d631d0d902e6p-2},
+          {0x1.c8b1bdc9d3f27p+0, 0x1.c8b1bdc9d3f29p+0}}},
+        {{{0x1.fffffffffffffp-1, 0x1.0000000000001p+0},
+          {0x1.26f65a0fd4a73p-1, 0x1.26f65a0fd4a75p-1},
+          {0x1.0e6d0af83d0acp-2, 0x1.0e6d0af83d0b2p-2},
+          {0x1.916bd5fe070d9p-4, 0x1.916bd5fe07101p-4},
+          {0x1.93ecedc1f0d7fp-5, 0x1.93ecedc1f0d89p-5},
+          {0x1.eb5d1d3ccbe86p+0, 0x1.eb5d1d3ccbe89p+0},
+          {0x1.78c97a83e17a7p+0, 0x1.78c97a83e17aap+0},
+          {0x1.0635d7caf70c6p+0, 0x1.0635d7caf70cdp+0},
+          {0x1.cc58e8a2ee2c2p-1, 0x1.cc58e8a2ee2dfp-1},
+          {0x1.b7deda7ffddf7p-2, 0x1.b7deda7ffde03p-2},
+          {-0x1.47a0e22f04e5ep-5, -0x1.47a0e22f04c0dp-5},
+          {0x1.6cb487831db4ap-2, 0x1.6cb487831db50p-2},
+          {0x1.cb21e9c5f325fp+0, 0x1.cb21e9c5f3263p+0}}},
+        {{{0x1.fffffffffffffp-1, 0x1.0000000000001p+0},
+          {0x1.26f24a61a564bp-1, 0x1.26f24a61a564dp-1},
+          {0x1.0e667df8508a0p-2, 0x1.0e667df8508a6p-2},
+          {0x1.9170e359fb821p-4, 0x1.9170e359fb849p-4},
+          {0x1.9417f976522f8p-5, 0x1.9417f97652303p-5},
+          {0x1.eb6a157592b20p+0, 0x1.eb6a157592b23p+0},
+          {0x1.78ccc103d7badp+0, 0x1.78ccc103d7bb1p+0},
+          {0x1.062f6c921cc39p+0, 0x1.062f6c921cc3fp+0},
+          {0x1.cc4595bdc0604p-1, 0x1.cc4595bdc0621p-1},
+          {0x1.b7c6fd930c880p-2, 0x1.b7c6fd930c88cp-2},
+          {-0x1.47e982ab3d9f7p-5, -0x1.47e982ab3d7a5p-5},
+          {0x1.6cb19975e1802p-2, 0x1.6cb19975e1808p-2},
+          {0x1.cb2cb028081f8p+0, 0x1.cb2cb028081fcp+0}}},
+        {{{0x1.fffffffffffffp-1, 0x1.0000000000001p+0},
+          {0x1.26f45238988fdp-1, 0x1.26f45238988ffp-1},
+          {0x1.0e69c46af3672p-2, 0x1.0e69c46af3678p-2},
+          {0x1.916e5ca776061p-4, 0x1.916e5ca776088p-4},
+          {0x1.940274204c9cdp-5, 0x1.940274204c9d8p-5},
+          {0x1.eb63995317741p+0, 0x1.eb63995317743p+0},
+          {0x1.78cb1dca864c4p+0, 0x1.78cb1dca864c8p+0},
+          {0x1.0632a241f5246p+0, 0x1.0632a241f524cp+0},
+          {0x1.cc4f3f1d6dc15p-1, 0x1.cc4f3f1d6dc32p-1},
+          {0x1.b7d2ebcdd3a6fp-2, 0x1.b7d2ebcdd3a7bp-2},
+          {-0x1.47c534f9a1c17p-5, -0x1.47c534f9a19c6p-5},
+          {0x1.6cb3107b2e3dep-2, 0x1.6cb3107b2e3e4p-2},
+          {0x1.cb274cef49caep+0, 0x1.cb274cef49cb1p+0}}},
+    },
+    {
+        0x1.0000000000000p-1,
+        {{0x1.29cbab649d389p-13, 0x1.d14e3bcd35a86p-6, 0x1.8c5c9a34ca0c3p-14,
+          0x1.13c68661ae70cp-11}},
+        {{{0x1.fffffffffffffp-1, 0x1.0000000000001p+0},
+          {0x1.b7590ca5d56abp-1, 0x1.b7590ca5d56adp-1},
+          {0x1.79c27556059fcp-1, 0x1.79c27556059fep-1},
+          {0x1.509c59ccb1964p-1, 0x1.509c59ccb1966p-1},
+          {0x1.422975f1d50c3p-1, 0x1.422975f1d50c5p-1},
+          {0x1.2a37abf469800p-1, 0x1.2a37abf469802p-1},
+          {0x1.1382380277aa0p-1, 0x1.1382380277aa2p-1},
+          {0x1.a5b316788e428p-2, 0x1.a5b316788e42ap-2},
+          {0x1.a5b316788e428p-2, 0x1.a5b316788e42ap-2},
+          {0x1.c86932802133fp-3, 0x1.c869328021341p-3},
+          {-0x1.26ee800000001p-303, 0x1.26ee800000001p-303},
+          {0x1.87247f6ca05d8p-1, 0x1.87247f6ca05dap-1},
+          {0x1.2a37abf469800p-1, 0x1.2a37abf469802p-1}}},
+        {{{0x1.fffffffffffffp-1, 0x1.0000000000001p+0},
+          {0x1.b7553624e5316p-1, 0x1.b7553624e5318p-1},
+          {0x1.79b5b3c71720bp-1, 0x1.79b5b3c71720dp-1},
+          {0x1.509288faf814fp-1, 0x1.509288faf8152p-1},
+          {0x1.4216f7b692db3p-1, 0x1.4216f7b692db6p-1},
+          {0x1.38c1b66704b5dp-1, 0x1.38c1b66704b63p-1},
+          {0x1.0c949871d1be7p-1, 0x1.0c949871d1be9p-1},
+          {0x1.c0cef4f93d8ddp-2, 0x1.c0cef4f93d8e5p-2},
+          {0x1.b3b7bc7fdfa53p-2, 0x1.b3b7bc7fdfa69p-2},
+          {0x1.bcf5e084222b4p-3, 0x1.bcf5e084222bdp-3},
+          {0x1.27c480885099dp-7, 0x1.27c4808850cf9p-7},
+          {0x1.87181cfb04848p-1, 0x1.87181cfb0484ap-1},
+          {0x1.2a7c94e9a3a3dp-1, 0x1.2a7c94e9a3a41p-1}}},
+        {{{0x1.fffffffffffffp-1, 0x1.0000000000001p+0},
+          {0x1.b751be282259dp-1, 0x1.b751be28225a0p-1},
+          {0x1.79af4d5fe0284p-1, 0x1.79af4d5fe0286p-1},
+          {0x1.508a2ff978465p-1, 0x1.508a2ff978468p-1},
+          {0x1.420defb8a1c5ap-1, 0x1.420defb8a1c5dp-1},
+          {0x1.38d0a97ead811p-1, 0x1.38d0a97ead816p-1},
+          {0x1.0ca165403faf5p-1, 0x1.0ca165403faf8p-1},
+          {0x1.c0e44203a3bafp-2, 0x1.c0e44203a3bb7p-2},
+          {0x1.b3cc5fc984930p-2, 0x1.b3cc5fc984946p-2},
+          {0x1.bd0aed39f3147p-3, 0x1.bd0aed39f3150p-3},
+          {0x1.27d1ae0dd005dp-7, 0x1.27d1ae0dd03b9p-7},
+          {0x1.87125aea910c9p-1, 0x1.87125aea910cbp-1},
+          {0x1.2a8ad7437b2aep-1, 0x1.2a8ad7437b2b1p-1}}},
+        {{{0x1.fffffffffffffp-1, 0x1.0000000000001p+0},
+          {0x1.b7537a2158170p-1, 0x1.b7537a2158173p-1},
+          {0x1.79b28088e2c65p-1, 0x1.79b28088e2c67p-1},
+          {0x1.508e5c6b105ecp-1, 0x1.508e5c6b105efp-1},
+          {0x1.421273a69f88ep-1, 0x1.421273a69f891p-1},
+          {0x1.38c9300704d3bp-1, 0x1.38c9300704d40p-1},
+          {0x1.0c9afeee3a733p-1, 0x1.0c9afeee3a736p-1},
+          {0x1.c0d99baae0253p-2, 0x1.c0d99baae025bp-2},
+          {0x1.b3c20e5419820p-2, 0x1.b3c20e5419836p-2},
+          {0x1.bd00671219eacp-3, 0x1.bd00671219eb4p-3},
+          {0x1.27cb17c00cee7p-7, 0x1.27cb17c00d243p-7},
+          {0x1.87153be8d5e91p-1, 0x1.87153be8d5e93p-1},
+          {0x1.2a83b62aa8393p-1, 0x1.2a83b62aa8396p-1}}},
+    },
+    {
+        0x1.0000000000000p+0,
+        {{0x1.5b49d2b1e91bbp-15, 0x1.0b630a91537a0p-7, 0x1.0ac1fae1b30dep-16,
+          0x1.91d0777c385fcp-15}},
+        {{{0x1.fffffffffffffp-1, 0x1.0000000000001p+0},
+          {0x1.ead7f87070237p-1, 0x1.ead7f87070239p-1},
+          {0x1.d8e87b7ae7a62p-1, 0x1.d8e87b7ae7a64p-1},
+          {0x1.ccec8c4ec5ef1p-1, 0x1.ccec8c4ec5ef3p-1},
+          {0x1.c8b73b9ac9a1dp-1, 0x1.c8b73b9ac9a1fp-1},
+          {0x1.5b5c75846aaf9p-3, 0x1.5b5c75846aafbp-3},
+          {0x1.40eb7cceb0482p-3, 0x1.40eb7cceb0484p-3},
+          {0x1.eb3e2ce365fd5p-4, 0x1.eb3e2ce365fd7p-4},
+          {0x1.eb3e2ce365fd5p-4, 0x1.eb3e2ce365fd7p-4},
+          {0x1.09dbd573f51e4p-4, 0x1.09dbd573f51e6p-4},
+          {-0x1.586e000000001p-305, 0x1.586e000000001p-305},
+          {0x1.dcce118e2abfap-1, 0x1.dcce118e2abfcp-1},
+          {0x1.5b5c75846aaf9p-3, 0x1.5b5c75846aafbp-3}}},
+        {{{0x1.fffffffffffffp-1, 0x1.0000000000001p+0},
+          {0x1.eada37c43dbcdp-1, 0x1.eada37c43dbd0p-1},
+          {0x1.d8eb0b5ab5aebp-1, 0x1.d8eb0b5ab5aeep-1},
+          {0x1.ccf1f94448ecfp-1, 0x1.ccf1f94448ed1p-1},
+          {0x1.c8badb25ce1b3p-1, 0x1.c8badb25ce1b6p-1},
+          {0x1.6c11624df5d8dp-3, 0x1.6c11624df5d96p-3},
+          {0x1.38a7a52a52899p-3, 0x1.38a7a52a5289dp-3},
+          {0x1.053de806af3a2p-3, 0x1.053de806af3a6p-3},
+          {0x1.fb41824cb7384p-4, 0x1.fb41824cb739ap-4},
+          {0x1.0303034e7937cp-4, 0x1.0303034e79384p-4},
+          {0x1.58908a0766c2ap-9, 0x1.58908a0767016p-9},
+          {0x1.dcd0247fc702ap-1, 0x1.dcd0247fc702dp-1},
+          {0x1.5b75854f419a5p-3, 0x1.5b75854f419aap-3}}},
+        {{{0x1.fffffffffffffp-1, 0x1.0000000000001p+0},
+          {0x1.ead630c226440p-1, 0x1.ead630c226443p-1},
+          {0x1.d8e399ff30e90p-1, 0x1.d8e399ff30e93p-1},
+          {0x1.cce84031702f6p-1, 0x1.cce84031702f8p-1},
+          {0x1.c8b05490c3638p-1, 0x1.c8b05490c363bp-1},
+          {0x1.6c56b7b4bec58p-3, 0x1.6c56b7b4bec60p-3},
+          {0x1.38e3300678b71p-3, 0x1.38e3300678b76p-3},
+          {0x1.056fa85832a89p-3, 0x1.056fa85832a8dp-3},
+          {0x1.fba21c895611fp-4, 0x1.fba21c8956135p-4},
+          {0x1.033456e6d857dp-4, 0x1.033456e6d8585p-4},
+          {0x1.58d2288b53912p-9, 0x1.58d2288b53cfep-9},
+          {0x1.dcc97107ca959p-1, 0x1.dcc97107ca95cp-1},
+          {0x1.5bb7b0fadf43fp-3, 0x1.5bb7b0fadf445p-3}}},
+        {{{0x1.fffffffffffffp-1, 0x1.0000000000001p+0},
+          {0x1.ead8344139715p-1, 0x1.ead8344139717p-1},
+          {0x1.d8e752a94ecadp-1, 0x1.d8e752a94ecafp-1},
+          {0x1.cced1cb61a4edp-1, 0x1.cced1cb61a4f0p-1},
+          {0x1.c8b597d621e34p-1, 0x1.c8b597d621e37p-1},
+          {0x1.6c340d23482b2p-3, 0x1.6c340d23482bap-3},
+          {0x1.38c56ab589a8cp-3, 0x1.38c56ab589a90p-3},
+          {0x1.0556c847cb264p-3, 0x1.0556c847cb268p-3},
+          {0x1.fb71cf9a50858p-4, 0x1.fb71cf9a5086ep-4},
+          {0x1.031bad32ce78bp-4, 0x1.031bad32ce792p-4},
+          {0x1.58b159698d54ep-9, 0x1.58b159698d939p-9},
+          {0x1.dccccac08134cp-1, 0x1.dccccac08134fp-1},
+          {0x1.5b969b4572366p-3, 0x1.5b969b457236cp-3}}},
+    },
+}};
+}  // namespace thermal_oracle
+#endif
diff --git a/tests/reference/thermal_reference_cases_test.cpp b/tests/reference/thermal_reference_cases_test.cpp
index 77c2aba..1484203 100644
--- a/tests/reference/thermal_reference_cases_test.cpp
+++ b/tests/reference/thermal_reference_cases_test.cpp
@@ -13,7 +13,7 @@ std::string Bytes(const std::filesystem::path& path) {
   std::ifstream in{path, std::ios::binary};
   return {std::istreambuf_iterator<char>{in}, std::istreambuf_iterator<char>{}};
 }
-TEST(ThermalReference, OriginalAbaqusNt11AtExactTime37) {
+TEST(ThermalReference, OriginalAbaqusNt11RemainsAnExplicitFailingDiagnostic) {
   const auto folder =
       std::filesystem::path{FESA_TEST_SOURCE_DIR} / "reference/shell-thermal";
   const auto input = folder / "shell-thermal.inp";
@@ -52,9 +52,17 @@ TEST(ThermalReference, OriginalAbaqusNt11AtExactTime37) {
             << " B=" << r.family.near_zero_band
             << " RMS/S=" << r.family.relative_rms
             << " max_abs=" << r.family.maximum_absolute_error
-            << " failed_rows=" << failed << " passed=" << r.passed << '\n';
+            << " failed_rows=" << failed << " numerical_passed=" << r.passed
+            << " classification=diagnostic blocking=" << r.family.blocking
+            << '\n';
   EXPECT_EQ(r.family.reference_scale, 1000);
   EXPECT_EQ(r.family.near_zero_band, 10);
-  EXPECT_TRUE(r.passed);
+  EXPECT_FALSE(r.family.blocking);
+  EXPECT_FALSE(r.passed);
+  EXPECT_FALSE(r.family.rms_passed);
+  EXPECT_EQ(failed, 87U);
+  EXPECT_GT(r.family.relative_rms,
+            fesa::test::ReferenceTolerancePolicy::kRelativeRmsTolerance);
+  EXPECT_GT(r.family.maximum_absolute_error, r.family.near_zero_band);
 }
 }  // namespace
diff --git a/tests/reference/thermal_reference_comparison.cpp b/tests/reference/thermal_reference_comparison.cpp
index d57b15e..6056972 100644
--- a/tests/reference/thermal_reference_comparison.cpp
+++ b/tests/reference/thermal_reference_comparison.cpp
@@ -154,14 +154,14 @@ Result<ThermalReferenceReport> CompareThermalSnapshots(
     ReferenceToleranceFamily family{
         "shell-thermal/Step-1/time=37/NT11/temperature/user-consistent",
         {"NT11"},
-        true,
+        false,
         {}};
     for (auto label : labels) family.rows.push_back({a.at(label), r.at(label)});
     const auto evaluation = ReferenceTolerancePolicy::Evaluate({family});
     Require(evaluation.valid && evaluation.families.size() == 1,
             "Invalid tolerance evaluation.");
     ThermalReferenceReport report{
-        evaluation.passed, {}, evaluation.families[0]};
+        evaluation.families[0].passed, {}, evaluation.families[0]};
     std::size_t index = 0;
     for (auto label : labels)
       report.rows.push_back({instance, label, 11, a.at(label), r.at(label),
diff --git a/tests/reference/thermal_reference_comparison.h b/tests/reference/thermal_reference_comparison.h
index 58ab065..0d02343 100644
--- a/tests/reference/thermal_reference_comparison.h
+++ b/tests/reference/thermal_reference_comparison.h
@@ -32,6 +32,8 @@ struct ThermalReferenceRowDecision {
   ReferenceToleranceRowDecision tolerance;
 };
 struct ThermalReferenceReport {
+  // Numerical family verdict, retained even when the original artifact is
+  // classified as nonblocking diagnostic by the approved option A decision.
   bool passed = false;
   std::vector<ThermalReferenceRowDecision> rows;
   ReferenceToleranceFamilyDecision family;
diff --git a/tests/reference/thermal_reference_comparison_test.cpp b/tests/reference/thermal_reference_comparison_test.cpp
index a4845c4..2d469fe 100644
--- a/tests/reference/thermal_reference_comparison_test.cpp
+++ b/tests/reference/thermal_reference_comparison_test.cpp
@@ -217,7 +217,7 @@ TEST(ThermalReferenceComparison,
             fesa::test::ReferenceToleranceBranch::kRelative);
   EXPECT_EQ(result.Value().family.reference_scale, 1000);
   EXPECT_EQ(result.Value().family.near_zero_band, 10);
-  EXPECT_TRUE(result.Value().family.blocking);
+  EXPECT_FALSE(result.Value().family.blocking);
 }
 TEST(ThermalReferenceComparison,
      RejectsMissingExtraDuplicateAndNonfiniteBeforeTolerance) {