Raise C line coverage past the 90 percent ecosystem target

.github/workflows/publish-test-badges.yml

@@ -64,7 +64,14 @@ jobs:
           from pathlib import Path
 
           payload = Path(sys.argv[1])
-          expected = ('tests-total.json', 'tests-fast.json', 'tests-nightly.json')
+          expected = (
+              'tests-total.json',
+              'tests-unit.json',
+              'tests-smoke.json',
+              'tests-integration.json',
+              'tests-fast.json',
+              'tests-nightly.json',
+          )
           for name in expected:
               path = payload / name
               if not path.is_file():

.gitignore

@@ -16,6 +16,12 @@ spider
 # dependency files
 *.d
 
+# coverage instrumentation and counters (from --coverage builds)
+*.gcno
+*.gcda
+coverage.xml
+coverage-*.json
+
 # Vim swap files
 .*.sw[omp]
 

Makefile

@@ -88,6 +88,7 @@ TEST_C_SRC = \
         tests/c/test_interp.c \
         tests/c/test_eos.c \
         tests/c/test_eos_composite.c \
+        tests/c/test_dimensionalisablefield.c \
 
 TEST_C_EXE = ${TEST_C_SRC:%.c=%}
 TEST_C_O = ${TEST_C_SRC:%.c=%.o}

docs/Explanations/testing.md

@@ -2,8 +2,9 @@
 
 [![Tests](https://img.shields.io/github/actions/workflow/status/FormingWorlds/SPIDER/ci.yml?branch=main&label=Tests)](https://github.com/FormingWorlds/SPIDER/actions/workflows/ci.yml)
 [![tests](https://img.shields.io/endpoint?url=https://raw.githubusercontent.com/FormingWorlds/SPIDER/badges/tests-total.json)](https://proteus-framework.org/testing)
-[![fast tests](https://img.shields.io/endpoint?url=https://raw.githubusercontent.com/FormingWorlds/SPIDER/badges/tests-fast.json)](https://github.com/FormingWorlds/SPIDER/actions/workflows/ci.yml)
-[![nightly tests](https://img.shields.io/endpoint?url=https://raw.githubusercontent.com/FormingWorlds/SPIDER/badges/tests-nightly.json)](https://github.com/FormingWorlds/SPIDER/actions/workflows/nightly.yml)
+[![unit tests](https://img.shields.io/endpoint?url=https://raw.githubusercontent.com/FormingWorlds/SPIDER/badges/tests-unit.json)](https://github.com/FormingWorlds/SPIDER/actions/workflows/ci.yml)
+[![smoke tests](https://img.shields.io/endpoint?url=https://raw.githubusercontent.com/FormingWorlds/SPIDER/badges/tests-smoke.json)](https://github.com/FormingWorlds/SPIDER/actions/workflows/ci.yml)
+[![integration tests](https://img.shields.io/endpoint?url=https://raw.githubusercontent.com/FormingWorlds/SPIDER/badges/tests-integration.json)](https://github.com/FormingWorlds/SPIDER/actions/workflows/nightly.yml)
 
 SPIDER's tests run under pytest in four tiers. The unit tier drives small C test executables that evaluate pure functions (interpolation, equation-of-state lookups, phase blending) at probe points and checks the results in Python; the smoke tier runs the real `spider` binary for a few macro steps and asserts physical invariants on the JSON output (mass closure, positivity, monotonicity, boundary-condition identities); the integration tier compares full runs against frozen reference output; and the slow tier is reserved for long validation runs. Every physics source file has a companion test file, and each is pinned against a published benchmark, an analytical limit, or an independent cross-check, inventoried under [Validation](../Validation/index.md).
 

docs/Validation/bc.md

@@ -11,3 +11,11 @@
 **Discrimination guards**: The grey-body flux at the same surface temperature differs from the prescribed value by far more than tolerance (dispatch guard); temperature ordering across the mantle (CMB hotter than surface) and the near-zero surface pressure are asserted on the default configuration.
 
 The boundary-condition dispatch is pinned through the prescribed-flux pathway PROTEUS uses in coupled runs: the value handed in through the option must reappear identically in the output, with the emissivity consistently inverted from it.
+
+**Reference-pinned test**: `tests/test_bc.py::test_constant_entropy_bc_pins_the_initial_surface`
+
+**Anchor**: Analytical identity: the constant-entropy boundary writes ic_surface_entropy (2550 J/kg/K) onto the surface basic node, preserved through the steady-state energy solve into the initial output.
+
+**Tolerance**: rel 1e-10 (assigned, not solved).
+
+**Discrimination guards**: The default configuration's surface value (2599.7 J/kg/K from the 2600 adiabat) sits about 50 J/kg/K away; positivity and table-range bounds on the core-side value after the steady rebalance. The prescribed-core-flux comparison against the core-cooling baseline (same file) discriminates the CORE_BC dispatch by the sign of the CMB entropy difference and the doubled CMB heat flux.

docs/Validation/energy.md

@@ -11,3 +11,19 @@
 **Discrimination guards**: The complementary mechanism stays nonzero when its sibling is disabled (wrong-path guard); the total flux is asserted to be the sum of the four components at every node in the full-physics configuration; scale bounds on the magma-ocean flux magnitudes.
 
 The flux assembly is validated through its superposition structure: each transport term can be switched off through its runtime option, and the total must respond exactly. This pins the dispatch wiring between the options, the individual flux routines, and the total used by the time integration.
+
+**Reference-pinned test**: `tests/test_energy.py::test_radiogenic_heating_pins_two_isotope_decay`
+
+**Anchor**: Analytical limit: radiogenic heating is the sum over isotopes of concentration times abundance times heat production, decaying as 2^(-(t - t0)/half_life). A two-isotope configuration with 100-year and 200-year half-lives pins the sum at t = 0 and the weighted decay factor 0.2506 after 200 years.
+
+**Tolerance**: rel 1e-6 on the initial sum, rel 1e-4 on the decay factor.
+
+**Discrimination guards**: A dropped isotope shifts the sum by 2.4e-3 relative; the two-isotope decay factor is resolved from the single-isotope values 0.25 and 0.5 and from the no-decay value 1.0.
+
+**Reference-pinned test**: `tests/test_energy.py::test_steady_state_ic_balances_interior_and_surface`
+
+**Anchor**: Analytical limit: the steady-state initial condition equalises the energy flow through every basic node and closes the surface radiation balance sigma (T_surf^4 - teqm^4) with the Stefan-Boltzmann constant of constants.c.
+
+**Tolerance**: Relative energy-flow spread below 1e-8 (observed 2e-11); grey-body pin at rel 1e-5.
+
+**Discrimination guards**: A T^3 exponent slip misses the flux by roughly a factor of the surface temperature (about 2000 here); sign and scale bounds on the flux.

docs/Validation/eos.md

@@ -11,3 +11,11 @@
 **Discrimination guards**: The solidus value at 10 GPa (about 1524 J/kg/K) sits far from the pinned liquidus value (about 2128 J/kg/K), so a swapped-boundary regression fails; sign and scale guards bound the entropies; a companion test pins the per-phase viscosity and conductivity dispatch with a 19-decade melt/solid contrast.
 
 eos.c owns the option parsing shared by all EOS implementations and loads the phase boundaries. The boundary interpolation is compared against numpy.interp on the same file, and the per-phase transport constants configured in the options file are recovered exactly from the evaluation.
+
+**Reference-pinned test**: `tests/test_eos.py::test_compositional_viscosity_matches_spaargaren`
+
+**Anchor**: The piecewise Mg/Si log10 viscosity prefactor of Spaargaren et al. (2020). The test re-encodes the published coefficients independently of eos.c and additionally pins the two branch-pair shifts as hand-derived literals (0.2074 and -1.35803), so a transcription error in either copy of the formula surfaces.
+
+**Tolerance**: abs 1e-9 on the evaluated shifts (closed-form option arithmetic); abs 1e-4 on the literal transcription pins.
+
+**Discrimination guards**: The two configurations shift in opposite directions across more than 1.5 decades, covering all four piecewise branches; a collapsed piecewise or a disabled compositional term fails, and the bare-constant edge case pins the term switched off. The activation-term test (same file) pins the zero-pressure limit exactly and the pressure-scale damping direction.

docs/Validation/ic.md

@@ -11,3 +11,19 @@
 **Discrimination guards**: A run with ic_adiabat_entropy = 2400 shifts the profile by exactly 200 J/kg/K (wrong-value discrimination); cooling at later outputs breaks the initial flatness, confirming the IC is not a fixed point of the output pipeline.
 
 The initial condition writes the prescribed adiabat onto the mesh. The top-node entropy, the integrated gradient ramp, and the response to a different prescribed value are pinned from the t = 0 output.
+
+**Reference-pinned test**: `tests/test_ic.py::test_abundance_ic_realises_the_requested_inventory`
+
+**Anchor**: Mass-balance identity through the initial partial-pressure solve: initial_kg equals the requested ppm of the mantle mass per volatile, and the realised reservoirs carry the requested hydrogen and carbon mole totals (conserved by the water and carbon dioxide reactions).
+
+**Tolerance**: rel 1e-9 on the stored inventories and the H and C totals (observed agreement 3e-15).
+
+**Discrimination guards**: The volatile oxygen total drifts by about 2e-3 relative against the request because the reactions exchange O with the melt's fO2 buffer, and the realised H2O reservoirs sit about 4e-3 relative off the bare request; the test thresholds sit at 1e-4, twenty-fold below the observed signals, while an inert reaction network would match the request to the 1e-9 the H and C totals meet.
+
+**Reference-pinned test**: `tests/test_ic.py::test_ocean_moles_ic_converts_moles_to_mass`
+
+**Anchor**: The OCEAN_MOLES constant of constants.c (7.68894973907177e22 mol per Earth ocean): initial_kg equals moles times OCEAN_MOLES times the volatile molar mass.
+
+**Tolerance**: rel 1e-9 (closed-form product of stored constants).
+
+**Discrimination guards**: Reading the CO2 inventory with the CO molar mass shifts the expected mass by 36 percent; positivity and scale bounds on the smallest inventory.

ic.c

@@ -932,9 +932,12 @@ static PetscErrorCode solve_for_initial_partial_pressure( Ctx *E )
        is way off.  The newtontr seems to help, and particularly increasing the size
        of delta0 to 10.0.  But I'm sure further optimisations are possible */
     ierr = PetscOptionsSetValue(NULL,"-atmosic_snes_type","newtontr");CHKERRQ(ierr);
-    /* Inform the nonlinear solver to generate a finite-difference approximation
-       to the Jacobian */
-    ierr = PetscOptionsSetValue(NULL,"-atmosic_snes_mf",NULL);CHKERRQ(ierr);
+    /* Build an explicit finite-difference Jacobian. The system is tiny
+       (one dof per volatile plus one per reaction), so the dense FD
+       Jacobian is cheap, and the trust-region solver needs a matrix
+       that supports transpose products, which the matrix-free
+       approximation does not provide. */
+    ierr = PetscOptionsSetValue(NULL,"-atmosic_snes_fd",NULL);CHKERRQ(ierr);
     /* Turn off convergence based on step size */
     ierr = PetscOptionsSetValue(NULL,"-atmosic_snes_stol","0");CHKERRQ(ierr);
     /* Turn off convergenced based on trust region tolerance */
@@ -973,10 +976,12 @@ static PetscErrorCode solve_for_initial_partial_pressure( Ctx *E )
 
     ierr = VecGetArray(x,&xx);CHKERRQ(ierr);
     for (i=0; i<Ap->n_volatiles; ++i) {
-        if( A->volatiles[i].p < 0.0 ){
-            /* Sanity check on solution (since it's non-unique) */
+        /* Sanity check on the solved value (the solution is non-unique):
+           test the converged unknown itself, not the stale entry left in
+           the atmosphere struct by the last residual evaluation. */
+        if( xx[i] < 0.0 ){
             SETERRQ2(PetscObjectComm((PetscObject)snes),PETSC_ERR_CONV_FAILED,
-                "Unphysical initial volatile partial pressure: volatile %d, x: %g",i,A->volatiles[i].p);
+                "Unphysical initial volatile partial pressure: volatile %d, x: %g",i,xx[i]);
         }
         else{
             A->volatiles[i].p = xx[i];

pyproject.toml

@@ -18,19 +18,17 @@ markers = [
 ]
 
 # C line coverage gates, measured with gcovr on a --coverage build.
-# Both gates ratchet upward via tools/update_coverage_threshold.py
-# (capped at the 90.0 ecosystem ceiling) and are never manually
-# decreased. The fast gate covers the unit + smoke tiers on every PR;
-# the full gate covers all tiers in the nightly workflow.
-# Initial floors sit a few points below the locally measured values
-# (79.6 fast, 81.2 full with llvm instrumentation) to absorb the
-# difference between compiler instrumentation variants; the ratchet
-# raises them from CI-measured coverage.
+# The fast gate covers the unit + smoke tiers on every PR; the full
+# gate covers all tiers in the nightly workflow. Both sit at the 90.0
+# ecosystem ceiling: the Linux PR runners measure 91 percent on the
+# fast tier and the full suite measures above 92, so each gate holds
+# a point or two of margin against instrumentation variants. Gates
+# are never decreased (the CI ratchet guard rejects any lowering).
 [tool.spider.coverage_fast]
-fail_under = 75.0
+fail_under = 90.0
 
 [tool.spider.coverage_full]
-fail_under = 77.0
+fail_under = 90.0
 
 [tool.ruff]
 line-length = 96

tests/c/test_dimensionalisablefield.c

@@ -0,0 +1,92 @@
+static const char help[] =
+    "Exercises the DimensionalisableField wrapper on a small DMDA:\n"
+    "create, fill, scale, unscale, duplicate, local vector, and the\n"
+    "scaling query. Prints one JSON object to stdout; the pytest\n"
+    "wrapper owns all assertions.\n";
+
+#include <petsc.h>
+
+#include "dimensionalisablefield.h"
+
+#define TEST_N 5
+#define TEST_SCALING 2.5
+
+int main(int argc, char **argv)
+{
+  PetscErrorCode         ierr;
+  DM                     dm;
+  DimensionalisableField f, fdup;
+  Vec                    v, vlocal;
+  PetscScalar            scalings[1] = {TEST_SCALING};
+  PetscScalar            dupScaling[1];
+  const PetscScalar      *arr;
+  PetscScalar            original0, scaled0, rescaled0, roundtrip0;
+  PetscInt               i, numDomains, nlocal;
+
+  ierr = PetscInitialize(&argc, &argv, NULL, help);
+  if (ierr) return (int)ierr;
+
+  ierr = DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, TEST_N, 1, 1, NULL, &dm);CHKERRQ(ierr);
+  ierr = DMSetUp(dm);CHKERRQ(ierr);
+
+  ierr = DimensionalisableFieldCreate(&f, dm, scalings, PETSC_FALSE);CHKERRQ(ierr);
+  ierr = DimensionalisableFieldSetName(f, "test field");CHKERRQ(ierr);
+  ierr = DimensionalisableFieldSetUnits(f, "test units");CHKERRQ(ierr);
+
+  /* fill the global vector with 1..N so index slips are visible */
+  ierr = DimensionalisableFieldGetGlobalVec(f, &v);CHKERRQ(ierr);
+  for (i = 0; i < TEST_N; ++i) {
+    ierr = VecSetValue(v, i, (PetscScalar)(i + 1), INSERT_VALUES);CHKERRQ(ierr);
+  }
+  ierr = VecAssemblyBegin(v);CHKERRQ(ierr);
+  ierr = VecAssemblyEnd(v);CHKERRQ(ierr);
+
+  ierr = VecGetArrayRead(v, &arr);CHKERRQ(ierr);
+  original0 = arr[0];
+  ierr = VecRestoreArrayRead(v, &arr);CHKERRQ(ierr);
+
+  /* scale multiplies by the stored scaling, unscale inverts it */
+  ierr = DimensionalisableFieldScale(f);CHKERRQ(ierr);
+  ierr = VecGetArrayRead(v, &arr);CHKERRQ(ierr);
+  scaled0 = arr[0];
+  ierr = VecRestoreArrayRead(v, &arr);CHKERRQ(ierr);
+
+  /* scaling an already-scaled field is the documented no-op guard:
+     it warns and must leave the values untouched */
+  ierr = DimensionalisableFieldScale(f);CHKERRQ(ierr);
+  ierr = VecGetArrayRead(v, &arr);CHKERRQ(ierr);
+  rescaled0 = arr[0];
+  ierr = VecRestoreArrayRead(v, &arr);CHKERRQ(ierr);
+
+  ierr = DimensionalisableFieldUnscale(f);CHKERRQ(ierr);
+  ierr = VecGetArrayRead(v, &arr);CHKERRQ(ierr);
+  roundtrip0 = arr[0];
+  ierr = VecRestoreArrayRead(v, &arr);CHKERRQ(ierr);
+
+  /* the duplicate carries the scaling and domain layout */
+  ierr = DimensionalisableFieldDuplicate(f, &fdup);CHKERRQ(ierr);
+  numDomains = 1;
+  ierr = DimensionalisableFieldGetScaling(fdup, &numDomains, dupScaling);CHKERRQ(ierr);
+
+  /* the local vector matches the serial DMDA layout */
+  ierr = DimensionalisableFieldCreateLocalVec(f, &vlocal);CHKERRQ(ierr);
+  ierr = VecGetLocalSize(vlocal, &nlocal);CHKERRQ(ierr);
+
+  ierr = PetscPrintf(PETSC_COMM_WORLD, "{\n");CHKERRQ(ierr);
+  ierr = PetscPrintf(PETSC_COMM_WORLD, "\"original0\": %.17g,\n", (double)original0);CHKERRQ(ierr);
+  ierr = PetscPrintf(PETSC_COMM_WORLD, "\"scaled0\": %.17g,\n", (double)scaled0);CHKERRQ(ierr);
+  ierr = PetscPrintf(PETSC_COMM_WORLD, "\"rescaled0\": %.17g,\n", (double)rescaled0);CHKERRQ(ierr);
+  ierr = PetscPrintf(PETSC_COMM_WORLD, "\"roundtrip0\": %.17g,\n", (double)roundtrip0);CHKERRQ(ierr);
+  ierr = PetscPrintf(PETSC_COMM_WORLD, "\"scaling\": %.17g,\n", (double)TEST_SCALING);CHKERRQ(ierr);
+  ierr = PetscPrintf(PETSC_COMM_WORLD, "\"dup_scaling0\": %.17g,\n", (double)dupScaling[0]);CHKERRQ(ierr);
+  ierr = PetscPrintf(PETSC_COMM_WORLD, "\"num_domains\": %d,\n", (int)numDomains);CHKERRQ(ierr);
+  ierr = PetscPrintf(PETSC_COMM_WORLD, "\"nlocal\": %d\n", (int)nlocal);CHKERRQ(ierr);
+  ierr = PetscPrintf(PETSC_COMM_WORLD, "}\n");CHKERRQ(ierr);
+
+  ierr = VecDestroy(&vlocal);CHKERRQ(ierr);
+  ierr = DimensionalisableFieldDestroy(&fdup);CHKERRQ(ierr);
+  ierr = DimensionalisableFieldDestroy(&f);CHKERRQ(ierr);
+  ierr = DMDestroy(&dm);CHKERRQ(ierr);
+  ierr = PetscFinalize();
+  return 0;
+}

tests/conftest.py

@@ -73,8 +73,17 @@ def _run(
         (outdir / 'spider_stdout.log').write_text(proc.stdout)
         (outdir / 'spider_stderr.log').write_text(proc.stderr)
         if proc.returncode != 0:
-            tail = (proc.stderr or proc.stdout)[-2000:]
-            raise RuntimeError(f'spider exited with code {proc.returncode} for {cmd}:\n{tail}')
+            # The PETSc error block sits near the top of a long stderr
+            # stream while the tail is MPI abort boilerplate; keep both
+            # so tests can assert on the actual error message.
+            stream = proc.stderr or proc.stdout
+            petsc = '\n'.join(
+                line for line in stream.splitlines() if 'PETSC ERROR' in line
+            )[:2000]
+            tail = stream[-2000:]
+            raise RuntimeError(
+                f'spider exited with code {proc.returncode} for {cmd}:\n{petsc}\n{tail}'
+            )
         return outdir
 
     return _run