diff --git a/.github/workflows/publish-test-badges.yml b/.github/workflows/publish-test-badges.yml index 73d54f20..90558ab5 100644 --- a/.github/workflows/publish-test-badges.yml +++ b/.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(): diff --git a/.gitignore b/.gitignore index 0dd9f39d..c48844bc 100644 --- a/.gitignore +++ b/.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] diff --git a/Makefile b/Makefile index 5a3730e4..35c003b7 100644 --- a/Makefile +++ b/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} diff --git a/docs/Explanations/testing.md b/docs/Explanations/testing.md index df7cb457..da0c44c2 100644 --- a/docs/Explanations/testing.md +++ b/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). diff --git a/docs/Validation/bc.md b/docs/Validation/bc.md index 71325039..db4f2e83 100644 --- a/docs/Validation/bc.md +++ b/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. diff --git a/docs/Validation/energy.md b/docs/Validation/energy.md index 24bd71b2..ae718c8e 100644 --- a/docs/Validation/energy.md +++ b/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. diff --git a/docs/Validation/eos.md b/docs/Validation/eos.md index 660a39ec..33bbc8a2 100644 --- a/docs/Validation/eos.md +++ b/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. diff --git a/docs/Validation/ic.md b/docs/Validation/ic.md index 355f8f0c..8cebd42f 100644 --- a/docs/Validation/ic.md +++ b/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. diff --git a/ic.c b/ic.c index e615cd34..50550d4a 100644 --- a/ic.c +++ b/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; in_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]; diff --git a/pyproject.toml b/pyproject.toml index c896ef20..f1c5dd47 100644 --- a/pyproject.toml +++ b/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 diff --git a/tests/c/test_dimensionalisablefield.c b/tests/c/test_dimensionalisablefield.c new file mode 100644 index 00000000..437e5530 --- /dev/null +++ b/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 + +#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; +} diff --git a/tests/conftest.py b/tests/conftest.py index f19d1b0f..dc882cea 100644 --- a/tests/conftest.py +++ b/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 diff --git a/tests/opts/reaction_library.opts b/tests/opts/reaction_library.opts new file mode 100644 index 00000000..443e9e85 --- /dev/null +++ b/tests/opts/reaction_library.opts @@ -0,0 +1,247 @@ +# Exercises the full named-reaction library: the IVTANTHERMO water and +# carbon dioxide pairs, the IVTANTHERMO methane reaction, and the ammonia +# reaction, over seven volatiles. Initial partial pressures are solved +# from the per-volatile total abundances so the initial state satisfies +# every reaction's equilibrium constraint. + +# time stepping + +# number of macro steps +-nstepsmacro 10 + +# macro step size in years +-dtmacro 1000 + +# tolerances +-ts_sundials_atol 1.0e-10 +-ts_sundials_rtol 1.0e-10 + +# mesh +# number of basic nodes +#-n 500 +# reduced mesh size for simple testing +-n 100 + +-MASS_COORDINATES + +-activate_poststep +-activate_rollback +-tsurf_poststep_change 1050 + +# all units are SI unless indicated otherwise + +# primary scaling parameters +# - these are used to non-dimensionalise (scale) the numerical problem +# and do not represent the physical quantities that are actually +# used (physical parameters are listed below) +-radius0 63710000.0 # Earth radius +-entropy0 2993.025100070677 # overturn of liquidus in Stixrude et al. (2009) +-time0 1.0E5 +-pressure0 10.0E5 +-volatile0 1.0E-6 + +# energy fluxes and sources to include +# [0] is off +# [1] is on +-CONDUCTION 1 # conduction +-CONVECTION 1 # convection +-MIXING 1 # mixing (latent heat transport) +-SEPARATION 1 # gravitational separation of solid and melt phase +-HTIDAL 0 # tidal heating (not currently implemented) + +# output directory +-outputDirectory output + +-phase_names melt,solid + +-melt_TYPE 1 +-melt_alpha_filename_rel_to_src lookup_data/1TPa-dK09-elec-free/thermal_exp_melt.dat +-melt_cp_filename_rel_to_src lookup_data/1TPa-dK09-elec-free/heat_capacity_melt.dat +-melt_dTdPs_filename_rel_to_src lookup_data/1TPa-dK09-elec-free/adiabat_temp_grad_melt.dat +-melt_rho_filename_rel_to_src lookup_data/1TPa-dK09-elec-free/density_melt.dat +-melt_temp_filename_rel_to_src lookup_data/1TPa-dK09-elec-free/temperature_melt.dat +-melt_log10visc 2.0 +-melt_cond 4.0 # conductivity of melt +-melt_phase_boundary_filename_rel_to_src lookup_data/1TPa-dK09-elec-free/liquidus_A11_H13.dat + +-solid_TYPE 1 +-solid_alpha_filename_rel_to_src lookup_data/1TPa-dK09-elec-free/thermal_exp_solid.dat +-solid_cp_filename_rel_to_src lookup_data/1TPa-dK09-elec-free/heat_capacity_solid.dat +-solid_dTdPs_filename_rel_to_src lookup_data/1TPa-dK09-elec-free/adiabat_temp_grad_solid.dat +-solid_rho_filename_rel_to_src lookup_data/1TPa-dK09-elec-free/density_solid.dat +-solid_temp_filename_rel_to_src lookup_data/1TPa-dK09-elec-free/temperature_solid.dat +-solid_log10visc 22.0 +-solid_cond 4.0 # conductivity of solid +-solid_phase_boundary_filename_rel_to_src lookup_data/1TPa-dK09-elec-free/solidus_A11_H13.dat + +# initial condition +# [1] is set dS/dr to constant everywhere +# [2] is restart from file +-IC_INTERIOR 1 +#-ic_interior_filename output/dSdr_b_aug_0.m +-ic_adiabat_entropy 2600.0 +-ic_dsdr -1.0E-5 +# set surface and core entropy (-1 is a flag to ignore, positive values are used) +-ic_surface_entropy -1 +-ic_core_entropy -1 + +# mixing length +# [1] is conventional, i.e. distance from the nearest boundary +# which is either the surface or core-mantle boundary +# [2] is constant, i.e. 1/4 x mantle depth, which is the mean +# of conventional theory. This can help to prevent a rigid +# lid from forming at the surface +# [3] is conventional as [1] above, but additionally allows you to +# specify a mid-mantle interface + +-mixing_length 2 +#-mixing_length_layer_radius 0.0 # non-dimensional radius + +# physical parameters + +# radius of the planet +-radius 6371000.0 # Earth radius + +# static pressure profile derived from Adams-Williamson equation of state +# these parameters are from fitting PREM in the lower mantle (for Earth) +-adams_williamson_rhos 4078.95095544 # surface density +-adams_williamson_beta 1.1115348931000002e-07 # beta parameter +-gravity -10.0 # gravity + +# eddy diffusivity +# if negative, this value is adopted (units m^2/s) +# if positive, this value is used to scale the internally calculated eddy diffusivity +-eddy_diffusivity_thermal 1.0 +-eddy_diffusivity_chemical 1.0 + +# material properties +-grain 1.0E-1 # grain size (m) + +# smoothing of material properties across liquidus and solidus +# units of melt fraction (non-dimensional) +-matprop_smooth_width 1.0E-2 + +# viscosity +-phi_critical 0.4 # transition melt fraction (non-dimensional) +-phi_width 0.15 # transition width (non-dimensional) + +# core-mantle boundary condition +# [1] core-cooling +# [2] heat flux (prescribe value using core_bc_value) +# [3] entropy +-CORE_BC 1 +-core_bc_value 0.0 + +# core parameters for CORE_BC=1 +-coresize 0.55 # fractional radius of core-mantle boundary +-rho_core 10738.332568062382 # core density +-cp_core 880.0 # core heat capacity + +# surface boundary condition +# [1] grey-body, i.e. sigma*emiss*(T^4-Teqm^4), with constant emissivity +# [2] steam atmosphere parameterisation (Zahnle et al., 1988) +# [3] self-consistent volatile evolution (Bower et al., 2019, Bower et al., 2021) +# [4] heat flux (prescribe value using surface_bc_value) +# [5] entropy +-SURFACE_BC 3 +-surface_bc_value 0.0 + +# if SURFACE_BC=1 or 2, the following constant emissivity is used +-emissivity0 1.0 +# equilibrium temperature of the planet +-teqm 273.0 +# parameterise the upper thermal boundary layer, where +# dT = param_utbl_const * [Surface temperature]**3.0 +# this is useful for high melt fraction dynamics where the ultra-thin thermal boundary layer +# (often only a few cms thick) cannot be resolved by the choice of mesh +-PARAM_UTBL 0 # flag to turn on parameterisation +-param_utbl_const 1.0E-7 # value of parameterisation + +-IC_ATMOSPHERE 1 + +# atmosphere reference pressure (Pa) +-P0 101325.0 # Pa (= 1 atm) + +# Define the volatiles +-volatile_names H2O,H2,CO2,CO,CH4,NH3,N2 + +-CO2_SOLUBILITY 3 +-CO2_kdist 0.0 +-CO2_kabs 1E-4 +-CO2_henry 0.0 +-CO2_henry_pow 1.0 +-CO2_molar_mass 0.04401 # kg/mol + +-H2O_SOLUBILITY 1 +-H2O_kdist 0.0 +-H2O_kabs 0.01 +-H2O_henry 1.6886562705299146 +-H2O_henry_pow 2.0 +-H2O_molar_mass 0.01801528 # kg/mol + +-CO_SOLUBILITY 1 +-CO_kdist 0.0 +-CO_kabs 1.0E-5 +-CO_henry 0.0 +-CO_henry_pow 1.0 +-CO_molar_mass 0.02801 # kg/mol + +-H2_SOLUBILITY 1 +-H2_kdist 0.0 +-H2_kabs 5.0E-5 # absorption (m^2/kg) +-H2_henry 0.0 +-H2_henry_pow 1.0 +-H2_molar_mass 0.00201588 # kg/mol + +# IW buffer +-OXYGEN_FUGACITY 7 + +# note that an offset of 0.5 is also the default, following +# Sossi et al., 2020 +-OXYGEN_FUGACITY_offset 0.5 + +-reaction_water_IVTANTHERMO +-reaction_carbondioxide_IVTANTHERMO +-reaction_methane_IVTANTHERMO +-reaction_ammonia1 + +# Initial volatile abundances (ppmw of the mantle) for the IC solve +-H2O_initial_total_abundance 100.0 +-H2_initial_total_abundance 20.0 +-CO2_initial_total_abundance 50.0 +-CO_initial_total_abundance 10.0 +-CH4_initial_total_abundance 1.0 +-NH3_initial_total_abundance 0.5 +-N2_initial_total_abundance 2.0 + +-CH4_SOLUBILITY 1 +-CH4_kdist 0.0 +-CH4_kabs 1.0E-5 +-CH4_henry 0.0 +-CH4_henry_pow 1.0 +-CH4_molar_mass 0.01604 # kg/mol + +-NH3_SOLUBILITY 1 +-NH3_kdist 0.0 +-NH3_kabs 1.0E-5 +-NH3_henry 0.0 +-NH3_henry_pow 1.0 +-NH3_molar_mass 0.017031 # kg/mol + +-N2_SOLUBILITY 1 +-N2_kdist 0.0 +-N2_kabs 1.0E-5 +-N2_henry 0.0 +-N2_henry_pow 1.0 +-N2_molar_mass 0.028014 # kg/mol + +# Initial guesses for the partial-pressure solve (Pa), near the +# equilibrium of the abundances above so the Newton iteration converges +# from the same starting point on every platform. +-H2O_initial_atmos_pressure 2.78e4 +-H2_initial_atmos_pressure 2.61e4 +-CO2_initial_atmos_pressure 3.27e5 +-CO_initial_atmos_pressure 3.33e6 +-CH4_initial_atmos_pressure 5.1e-4 +-NH3_initial_atmos_pressure 1.2 +-N2_initial_atmos_pressure 2.07e5 diff --git a/tests/test_bc.py b/tests/test_bc.py index 4694b994..451a1a11 100644 --- a/tests/test_bc.py +++ b/tests/test_bc.py @@ -99,3 +99,94 @@ def test_mantle_cools_from_the_top_between_surface_and_cmb(blackbody_short): # Scale guard on the deep end: blackbody50 reaches ~1.38e11 Pa at # the CMB; a nondimensional leak would sit near unity. assert 1e10 < pressure_b[-1] < 1e12 + + +@pytest.fixture(scope='module') +def core_flux_run(cached_spider_run): + """Two-macro-step blackbody50 run with a prescribed core heat flux.""" + return cached_spider_run( + overrides=('-nstepsmacro', '2', '-CORE_BC', '2', '-core_bc_value', '1.0e5'), + name='core_flux', + ) + + +@pytest.mark.physics_invariant +def test_prescribed_core_flux_slows_cmb_cooling(core_flux_run, blackbody_short): + """CORE_BC 2 injects heat at the CMB and slows the deep cooling. + + Against the default core-cooling run at identical initial + conditions and output times, a prescribed 1e5 W/m^2 inflow keeps + the deepest staggered node hotter after 200 years (observed 2515 + vs 2511 J/kg/K) and roughly doubles the total heat flux carried at + the core-mantle boundary node. + """ + flux_0 = data_si(read_output(core_flux_run, 0), 'S_s') # J/kg/K + cool_0 = data_si(read_output(blackbody_short, 0), 'S_s') # J/kg/K + # Edge case, identical IC: the boundary condition acts only during + # integration, so t = 0 must agree between the two runs exactly. + np.testing.assert_allclose(flux_0, cool_0, rtol=1e-12) + + flux_200 = data_si(read_output(core_flux_run, 200), 'S_s') + cool_200 = data_si(read_output(blackbody_short, 200), 'S_s') + # Direction discrimination: heating from below leaves the CMB node + # hotter than core cooling does; the observed gap (3.9 J/kg/K) is + # far above the CVODE tolerance floor. + assert flux_200[-1] > cool_200[-1] + assert flux_200[-1] - cool_200[-1] > 1.0 # J/kg/K + + # The prescribed inflow roughly doubles the CMB heat flux relative + # to core cooling (observed 1.9e5 vs 9.2e4 W/m^2). + jtot_flux = data_si(read_output(core_flux_run, 200), 'Jtot_b') # W/m^2 + jtot_cool = data_si(read_output(blackbody_short, 200), 'Jtot_b') # W/m^2 + assert jtot_flux[-1] > 1.5 * jtot_cool[-1] + # Positivity: entropy stays physical under the modified BC. + assert np.all(flux_200 > 0) + + +@pytest.fixture(scope='module') +def entropy_bc_run(cached_spider_run): + """One-step run with constant-entropy boundaries and a steady IC. + + The prescribed boundary entropies (2550 and 2650 J/kg/K) bracket + the 2600 J/kg/K adiabat; the steady-state energy solve keeps the + perturbed initial condition integrable within one macro step. + """ + return cached_spider_run( + overrides=( + '-nstepsmacro', '1', + '-dtmacro', '10', + '-ic_surface_entropy', '2550', + '-ic_core_entropy', '2650', + '-ic_steady_state_energy', + ), + name='entropy_bc_steady', + ) + + +@pytest.mark.physics_invariant +@pytest.mark.reference_pinned +def test_constant_entropy_bc_pins_the_initial_surface(entropy_bc_run, blackbody_short): + """ic_surface_entropy fixes the surface basic node exactly. + + Anchor: analytical identity. The constant-entropy boundary writes + the option value (2550 J/kg/K) onto the surface basic node, and the + subsequent steady-state solve preserves it, so the initial output + carries the prescribed value exactly. The default run's surface + (2599.7 J/kg/K, from the 2600 adiabat and the ic_dsdr gradient) + discriminates the option not being applied. + """ + s_b = data_si(read_output(entropy_bc_run, 0), 'S_b') # J/kg/K + + # rel=1e-10: the value is assigned, not solved for. + assert s_b[0] == pytest.approx(2550.0, rel=1e-10) + + # Discrimination guard: the default IC surface value differs by + # ~50 J/kg/K, five orders above the pin tolerance. + s_b_default = data_si(read_output(blackbody_short, 0), 'S_b') + assert abs(s_b_default[0] - 2550.0) > 40.0 # J/kg/K + + # The core-side prescription seeds the deep profile before the + # steady-state solve rebalances it; the result must stay physical + # and within the entropy range of the loaded tables. + assert np.all(s_b > 0) + assert 2400.0 < s_b[-1] < 2700.0 # J/kg/K diff --git a/tests/test_dimensionalisablefield.py b/tests/test_dimensionalisablefield.py new file mode 100644 index 00000000..920e46b2 --- /dev/null +++ b/tests/test_dimensionalisablefield.py @@ -0,0 +1,65 @@ +"""Tests for dimensionalisablefield.c (the scaled-field wrapper). + +Driven through the C test executable tests/c/test_dimensionalisablefield, +which builds a five-node field with scaling 2.5, fills it with 1..5, and +exercises scale, unscale, duplicate, the local-vector factory, and the +scaling query. Contract clauses exercised: scaling multiplies by the +stored factor, unscaling inverts it exactly, the duplicate carries the +scaling and domain count, and the serial local vector matches the mesh +size. See docs/How-to/build_tests.md for the tier system. +""" + +from __future__ import annotations + +import pytest + +pytestmark = [pytest.mark.unit, pytest.mark.timeout(30)] + +# Values hard-coded in the C executable. +SCALING = 2.5 +FIRST_VALUE = 1.0 +MESH_SIZE = 5 + + +def test_scale_unscale_round_trip_is_exact(c_test): + """Scaling multiplies by the stored factor and unscaling inverts it. + + The field starts unscaled with first value 1.0; after Scale the + stored entry is 2.5, and after Unscale it returns to 1.0 exactly + (a multiply and divide by the same factor, no accumulation). The + unit scaling factor is the edge the wrapper must not special-case: + with 2.5 the scaled and raw values differ by 150 percent, so a + no-op Scale cannot pass. + """ + out = c_test('test_dimensionalisablefield') + + # rel=1e-15: one multiplication, roundoff only. + assert out['scaled0'] == pytest.approx(SCALING * FIRST_VALUE, rel=1e-15) + # State-change guard: scaling must actually change the stored values. + assert abs(out['scaled0'] - out['original0']) > 1.0 + # Error contract: scaling an already-scaled field is the documented + # no-op guard; it warns and leaves every value untouched, so a + # second Scale call must not multiply again (which would give 6.25). + assert out['rescaled0'] == pytest.approx(out['scaled0'], rel=1e-15) + assert abs(out['rescaled0'] - SCALING**2 * FIRST_VALUE) > 1.0 + # Round trip: unscale exactly inverts scale. + assert out['roundtrip0'] == pytest.approx(out['original0'], rel=1e-15) + assert out['original0'] == pytest.approx(FIRST_VALUE, rel=1e-15) + + +def test_duplicate_carries_scaling_and_layout(c_test): + """A duplicated field reports the parent's scaling and domain count. + + The duplicate is created over the same serial DMDA, so it has one + domain, the parent's 2.5 scaling, and a local vector of the mesh + size (the limit input: a single-domain, single-rank layout with no + ghost points). + """ + out = c_test('test_dimensionalisablefield') + + assert out['dup_scaling0'] == pytest.approx(SCALING, rel=1e-15) + # Scaling discrimination: the duplicate must carry 2.5, not the + # unit scaling a fresh field would default to. + assert abs(out['dup_scaling0'] - 1.0) > 1.0 + assert int(out['num_domains']) == 1 + assert int(out['nlocal']) == MESH_SIZE diff --git a/tests/test_energy.py b/tests/test_energy.py index 28d8e445..6772a369 100644 --- a/tests/test_energy.py +++ b/tests/test_energy.py @@ -154,3 +154,210 @@ def test_surface_node_flux_is_the_grey_body_atmospheric_flux(blackbody_short): # Scale guard: magma-ocean surface fluxes are 1e2 to 1e9 W/m^2; a # nondimensional value leaking through unscaled would sit near unity. assert 1e2 < fatm < 1e9 # W/m^2 + + +# Two-isotope internal heating configuration. The half-lives (100 and +# 200 years) are deliberately short so the decay resolves within the +# 200-year run; the production rates and concentrations are of the +# order of al26 and k40 in early solar system material. +HEATING_OVERRIDES = ( + '-nstepsmacro', '2', + # dtmacro is pinned here because the decay pin below requires the + # second output at exactly t = 200 years (two al26 half-lives, one + # k40 half-life); the test must not inherit it silently from the + # shared options file. + '-dtmacro', '100', + '-radionuclide_names', 'al26,k40', + '-al26_t0', '0.0', + '-al26_abundance', '1.0', + '-al26_concentration', '10.0', + '-al26_heat_production', '0.3568', + '-al26_half_life', '100.0', + '-k40_t0', '0.0', + '-k40_abundance', '1.0', + '-k40_concentration', '300.0', + '-k40_heat_production', '2.92e-5', + '-k40_half_life', '200.0', + '-HTIDAL', '1', + '-htidal_value', '1.0e-7', +) +# Specific heating of each isotope at t0: concentration (mass fraction) +# times isotopic abundance times heat production. al26 dominates. +H_AL26 = 10.0e-6 * 1.0 * 0.3568 # 3.568e-6 W/kg +H_K40 = 300.0e-6 * 1.0 * 2.92e-5 # 8.76e-9 W/kg +HTIDAL_VALUE = 1.0e-7 # W/kg, prescribed constant tidal heating + + +@pytest.fixture(scope='module') +def heating_run(cached_spider_run): + """Two-macro-step blackbody50 run with radiogenic and tidal heating.""" + return cached_spider_run(overrides=HEATING_OVERRIDES, name='heating') + + +@pytest.mark.physics_invariant +@pytest.mark.reference_pinned +def test_radiogenic_heating_pins_two_isotope_decay(heating_run): + """Radiogenic heating sums both isotopes and follows exponential decay. + + Anchor: analytical limit. Each radionuclide contributes + concentration * abundance * heat_production * 2^(-(t - t0)/half_life) + of specific heating, uniform across the mantle. At t = 0 the sum is + 3.57676e-6 W/kg, and after 200 years (two al26 half-lives, one k40 + half-life) the total falls to 0.2506 of its initial value. + """ + h_0 = data_si(read_output(heating_run, 0), 'Hradio_s') # W/kg + h_200 = data_si(read_output(heating_run, 200), 'Hradio_s') # W/kg + + # Radiogenic heating is spatially uniform: the same isotopes and + # concentration heat every shell (rel 1e-12 is roundoff headroom). + np.testing.assert_allclose(h_0, h_0[0], rtol=1e-12) + + # rel=1e-6: the sum is evaluated in closed form from the option + # values, so only scaling round-trip error enters. + assert h_0[0] == pytest.approx(H_AL26 + H_K40, rel=1e-6) + # Isotope-sum guard: dropping k40 shifts the total by 2.4e-3 + # relative, three orders above the pin tolerance. + assert abs(h_0[0] - H_AL26) > 1e-3 * h_0[0] + + # Decay pin: the two-isotope weighted decay factor after 200 years. + expected_ratio = (H_AL26 * 0.25 + H_K40 * 0.5) / (H_AL26 + H_K40) + ratio = h_200[0] / h_0[0] + # rel=1e-4: the run writes output exactly at t = 200 years, so the + # closed-form factor applies with only output rounding on top. + assert ratio == pytest.approx(expected_ratio, rel=1e-4) + # Decay guards: no decay gives 1.0, a single-isotope slip gives + # exactly 0.25 (al26 only) or 0.5 (k40 only); the two-isotope value + # 0.25061 is resolved from all three far beyond the tolerance. + assert abs(ratio - 1.0) > 0.7 + assert abs(ratio - 0.25) > 5e-4 + assert abs(ratio - 0.5) > 0.2 + # Sign guard: heat production is strictly positive. + assert np.all(h_0 > 0) and np.all(h_200 > 0) + + +@pytest.mark.physics_invariant +def test_prescribed_tidal_heating_is_uniform_and_steady(heating_run): + """HTIDAL 1 applies a constant specific heating that does not decay. + + The prescribed scalar must appear unchanged at every node and at + every output time, in contrast to the radiogenic contribution in + the same run, which decays by a factor of four over the window. + """ + ht_0 = data_si(read_output(heating_run, 0), 'Htidal_s') # W/kg + ht_200 = data_si(read_output(heating_run, 200), 'Htidal_s') # W/kg + + # rel=1e-12: the option value is stored and re-emitted, not + # recomputed. + np.testing.assert_allclose(ht_0, HTIDAL_VALUE, rtol=1e-12) + np.testing.assert_allclose(ht_200, HTIDAL_VALUE, rtol=1e-12) + + # Scale guard: a nondimensional leak (the value divided by the + # heat-generation scaling) would land far outside this bracket. + assert 5e-8 < ht_0[0] < 2e-7 # W/kg + + # Contrast with the decaying radiogenic term in the same run: the + # tidal column is time independent while Hradio drops fourfold. + h_0 = data_si(read_output(heating_run, 0), 'Hradio_s') + h_200 = data_si(read_output(heating_run, 200), 'Hradio_s') + assert h_200[0] < 0.3 * h_0[0] + assert ht_200[0] == pytest.approx(ht_0[0], rel=1e-12) + + +@pytest.mark.physics_invariant +def test_tidal_heating_profile_from_file_is_recovered( + blackbody_short, run_spider, tmp_path +): + """HTIDAL 2 reads a per-node heating profile and applies it verbatim. + + A linear ramp from 1e-8 to 5e-8 W/kg across the staggered nodes is + written in the Interp1d file format (header with row count, scaling + line, then pressure and heating columns) and must be recovered + node for node. A file with the wrong number of rows must be + refused, exercising the documented length check. + """ + pressure_s = data_si(read_output(blackbody_short, 0), 'pressure_s') # Pa + n = len(pressure_s) + ramp = np.linspace(1.0e-8, 5.0e-8, n) # W/kg + + lines = [ + f'# 5 {n}', + '# Pressure, Htidal', + '# column * scaling factor should be SI units', + '# scaling factors (constant) for each column given on line below', + '# 1.0 1.0', + ] + lines += [f'{p:.18e} {h:.18e}' for p, h in zip(pressure_s, ramp)] + profile = tmp_path / 'htidal_ramp.dat' + profile.write_text('\n'.join(lines) + '\n') + + outdir = run_spider( + overrides=('-nstepsmacro', '1', '-HTIDAL', '2', '-htidal_filename', str(profile)), + name='htidal_file', + ) + ht = data_si(read_output(outdir, 0), 'Htidal_s') # W/kg + + # rel=1e-10: the file values are read by node index with unit + # scalings, so recovery is exact to I/O rounding. + np.testing.assert_allclose(ht, ramp, rtol=1e-10) + # Orientation guard: a reversed node indexing would flip the ramp; + # the endpoints differ by a factor of five. + assert ht[0] == pytest.approx(1.0e-8, rel=1e-10) + assert ht[-1] == pytest.approx(5.0e-8, rel=1e-10) + assert np.all(np.diff(ht) > 0) + + # Error contract: the reader checks the header's row count against + # the staggered mesh and refuses a profile of the wrong length. + short_lines = [f'# 5 {5}'] + lines[1:5] + short_lines += [f'{pressure_s[0]:.18e} 1.0e-8'] * 5 + short = tmp_path / 'htidal_short.dat' + short.write_text('\n'.join(short_lines) + '\n') + with pytest.raises(RuntimeError): + run_spider( + overrides=('-nstepsmacro', '1', '-HTIDAL', '2', '-htidal_filename', str(short)), + name='htidal_short', + ) + + +@pytest.fixture(scope='module') +def steady_state_run(cached_spider_run): + """One-step blackbody50 run with the steady-state initial condition.""" + return cached_spider_run( + overrides=('-nstepsmacro', '1', '-dtmacro', '10', '-ic_steady_state_energy'), + name='steady_ic', + ) + + +@pytest.mark.physics_invariant +@pytest.mark.reference_pinned +def test_steady_state_ic_balances_interior_and_surface(steady_state_run): + """The steady-state IC equalises the energy flow through the mantle. + + Anchor: analytical limit. With -ic_steady_state_energy the surface + radiation balance is solved together with a steady interior, so at + t = 0 the energy flow Etot is the same through every basic node + (observed relative spread 2e-11) and the surface flux satisfies the + grey-body law sigma * (T_surf^4 - teqm^4) with emissivity 1. + """ + doc = read_output(steady_state_run, 0) + etot = data_si(doc, 'Etot_b') # W + fatm = atmosphere_si(doc, 'Fatm') # W/m^2 + t_surf = atmosphere_si(doc, 'temperature_surface') # K + + # Steady state: constant energy flow through the interior. The + # 1e-8 bound sits three orders above the observed 2e-11 spread yet + # far below the order-one spread of the default (non-steady) IC. + assert (etot.max() - etot.min()) / abs(etot.mean()) < 1e-8 + assert np.all(etot > 0) + + # Grey-body pin with the sigma hard-coded in constants.c. + sigma = 5.670367e-8 # W m^-2 K^-4 + expected = sigma * (t_surf**4 - 273.0**4) + # rel=1e-5: the radiation balance is solved to snes_rtol 1e-9; the + # margin absorbs output rounding of T_surf. + assert fatm == pytest.approx(expected, rel=1e-5) + # Exponent guard: a T^3 slip misses by roughly a factor of T_surf + # (~2000 at the ~1960 K of this state), far beyond any tolerance. + assert abs(fatm - sigma * (t_surf**3 - 273.0**3)) > 0.5 * fatm + # Sign and scale guards. + assert fatm > 0 + assert 1e2 < fatm < 1e9 # W/m^2 diff --git a/tests/test_eos.py b/tests/test_eos.py index 8247934d..a55abe7e 100644 --- a/tests/test_eos.py +++ b/tests/test_eos.py @@ -84,3 +84,129 @@ def test_phase_boundary_matches_independent_table_read(c_test): # positive and of order a few thousand J/kg/K. assert out['phase_boundary_S'][1] > 0 assert 1.5e3 < out['phase_boundary_S'][1] < 4.0e3 + + +def _spaargaren_prefactor(mg_si): + """log10 viscosity prefactor of Spaargaren et al. (2020), as in eos.c.""" + if mg_si <= 1.0: + return 0.5185 * (1 - mg_si) / 0.3 + if mg_si <= 1.25: + return -1.4815 * (mg_si - 1) / 0.25 + if mg_si <= 1.5: + return -2 + 0.5185 * (1.5 - mg_si) / 0.25 + return -2.0 + + +@pytest.mark.physics_invariant +@pytest.mark.reference_pinned +def test_compositional_viscosity_matches_spaargaren(c_test): + """The Mg/Si viscosity prefactor reproduces Spaargaren et al. (2020). + + Anchor: the piecewise log10 prefactor of Spaargaren et al. (2020). + Two evaluations cover all four Mg/Si branches: 1.4 against a 1.6 + reference (upper-intermediate and Fp-rich branches) and 1.2 against + a 0.9 reference (lower-intermediate and Mg-rich branches). The + shift adds to the configured melt log10visc of 2.0. + """ + probe = ((1.0e10,), (2600.0,)) + high = _eval_with( + c_test, ('-melt_visc_comp', '1.4', '-melt_visc_ref_comp', '1.6'), *probe + ) + low = _eval_with( + c_test, ('-melt_visc_comp', '1.2', '-melt_visc_ref_comp', '0.9'), *probe + ) + + shift_high = _spaargaren_prefactor(1.4) - _spaargaren_prefactor(1.6) + shift_low = _spaargaren_prefactor(1.2) - _spaargaren_prefactor(0.9) + # Transcription pins: the helper re-encodes the published piecewise + # coefficients independently of eos.c; these literals (hand-derived + # from Spaargaren et al. 2020) guard the helper itself, so a shared + # transcription error cannot silently cancel. + assert shift_high == pytest.approx(0.2074, abs=1e-4) + assert shift_low == pytest.approx(-1.35803, abs=1e-4) + # abs=1e-9: the prefactor is a closed-form expression of the two + # option values, so only roundoff enters. + assert high['log10visc'][0] == pytest.approx(2.0 + shift_high, abs=1e-9) + assert low['log10visc'][0] == pytest.approx(2.0 + shift_low, abs=1e-9) + + # Branch discrimination: the two configurations shift in opposite + # directions and differ by more than 1.5 decades, so a collapsed + # piecewise (any single branch applied throughout) fails. + assert shift_high > 0 > shift_low + assert abs(high['log10visc'][0] - low['log10visc'][0]) > 1.5 + + # Edge case: with the compositional term disabled (no visc_comp + # option) the melt viscosity is the bare configured constant. + base = _eval(c_test, 'melt', *probe) + assert base['log10visc'][0] == pytest.approx(2.0, abs=1e-10) + + +def _eval_with(c_test, extra, p_list, s_list): + """Evaluate the melt EOS with additional viscosity options.""" + return c_test( + 'test_eos', + ( + *OPTS, + *extra, + '-eos_prefix', 'melt', + '-P_si', ','.join(str(p) for p in p_list), + '-S_si', ','.join(str(s) for s in s_list), + ), + ) + + +@pytest.mark.physics_invariant +def test_activation_terms_shape_the_viscosity_profile(c_test): + """Activation volume and energy bend the viscosity with P and T. + + The unpinned activation-volume term adds V * P * exp(-P/Ps) / T to + the log viscosity: it vanishes identically at the P = 0 surface + limit (edge case), stiffens the deep mantle, and is damped by the + pressure scale Ps. The pinned activation-energy term with a + reference temperature below the actual temperature softens the + melt, discriminating the sign convention of dT. + """ + probes_p = (0.0, 6.55e10) # Pa; surface limit and deep mantle + probes_s = (2600.0, 2600.0) # J/kg/K + + vol = _eval_with( + c_test, ('-melt_activation_volume', '1.0e-6'), probes_p, probes_s + ) + # Edge limit: no pressure, no activation-volume contribution; the + # bare constant 2.0 returns exactly. + assert vol['log10visc'][0] == pytest.approx(2.0, abs=1e-10) + # Deep mantle stiffens: positive V and P give a positive shift. + assert vol['log10visc'][1] > 2.0 + 1e-6 + + # The pressure scale damps the deep contribution: with Ps = 10 GPa + # the 65.5 GPa probe carries exp(-6.55) = 1.4e-3 of the undamped + # term, so the shift shrinks by orders of magnitude but stays + # positive. + damped = _eval_with( + c_test, + ( + '-melt_activation_volume', '1.0e-6', + '-melt_activation_volume_pressure_scale', '1.0e10', + ), + probes_p, + probes_s, + ) + undamped_shift = vol['log10visc'][1] - 2.0 + damped_shift = damped['log10visc'][1] - 2.0 + assert 0 < damped_shift < 0.1 * undamped_shift + + # Pinned activation energy with a reference temperature below the + # melt temperature: dT < 0, so the melt is softer than the bare + # constant, and the shift grows with temperature contrast (the + # hotter 3000 J/kg/K probe sits further from the reference). + energy = _eval_with( + c_test, + ('-melt_activation_energy', '1.0e5', '-melt_visc_ref_temp', '1500.0'), + (1.0e10, 1.0e10), + (2500.0, 3000.0), + ) + assert energy['log10visc'][0] < 2.0 + assert energy['log10visc'][1] < 2.0 + # Temperature dependence: the two probes shift by different + # amounts, so a T-independent (constant) offset fails. + assert abs(energy['log10visc'][0] - energy['log10visc'][1]) > 1e-4 diff --git a/tests/test_ic.py b/tests/test_ic.py index 0d5cc0b8..e7aa40fc 100644 --- a/tests/test_ic.py +++ b/tests/test_ic.py @@ -15,7 +15,7 @@ import numpy as np import pytest -from tests._json_utils import data_si, read_output, solution_entries +from tests._json_utils import data_si, field_si, read_output, solution_entries pytestmark = [pytest.mark.smoke, pytest.mark.timeout(120)] @@ -98,3 +98,214 @@ def test_initial_adiabat_is_nearly_flat_until_cooling_imprints_structure(blackbo s_200 = data_si(read_output(blackbody_short, 200), 'S_s') assert s_200[0] < s_0[0] assert s_0[0] - s_200[0] > 10.0 # J/kg/K + + +# Molar masses (kg/mol) of the reaction.opts volatiles, matching the +# per-volatile molar_mass options in the configuration. +MOLAR_MASS = { + 'H2O': 0.01801528, + 'H2': 0.00201588, + 'CO2': 0.04401, + 'CO': 0.02801, +} +# Elemental composition (atoms per molecule) for the mole bookkeeping. +STOICHIOMETRY = { + 'H2O': {'H': 2, 'O': 1}, + 'H2': {'H': 2}, + 'CO2': {'C': 1, 'O': 2}, + 'CO': {'C': 1, 'O': 1}, +} +# Requested inventories (ppm by mass of the mantle) for the abundance IC. +# The values sit 20 to 40 percent off the three-ocean equilibrium implied +# by the configured guess pressures (984, 0.45, 5.0, 22.6 ppm), so the +# Newton solve stays within its convergence basin on every platform while +# the reactions still have to redistribute a resolvable amount of mass. +ABUNDANCE_PPM = {'H2O': 800.0, 'H2': 0.6, 'CO2': 7.0, 'CO': 18.0} +# Requested inventories (Earth oceans) for the ocean-moles IC, matched to +# the abundance targets above; the mole count of one ocean is the +# OCEAN_MOLES constant in constants.c. +OCEAN_MOLES = 7.68894973907177e22 # mol per Earth ocean of H2O (or H2) +MOLES_OCEANS = {'H2O': 2.43, 'H2': 0.0163, 'CO2': 0.0087, 'CO': 0.0352} + + +def _volatile_reservoirs_kg(doc, volatile): + """SI masses of the liquid, solid, and atmosphere reservoirs.""" + block = doc['atmosphere'][volatile] + return { + key: float(field_si(block[key])[0]) + for key in ('initial_kg', 'liquid_kg', 'solid_kg', 'atmosphere_kg', 'physical_kg') + } + + +def _element_moles(masses_kg, element): + """Total moles of one element across the four volatile inventories.""" + return sum( + masses_kg[v] / MOLAR_MASS[v] * STOICHIOMETRY[v].get(element, 0) + for v in MOLAR_MASS + ) + + +@pytest.fixture(scope='module') +def abundance_ic_run(cached_spider_run): + """Zero-step reaction run with the abundance-based atmosphere IC. + + IC_ATMOSPHERE 1 solves the initial partial pressures and reaction + masses from the per-volatile total abundances, subject to the + water and carbon dioxide equilibrium constraints. Zero macro steps: + the tests read only the initial condition, so the run skips time + integration entirely. + """ + return cached_spider_run( + opts_file='reaction.opts', + overrides=( + '-IC_ATMOSPHERE', '1', + '-nstepsmacro', '0', + '-n', '50', + '-H2O_initial_total_abundance', '800.0', + '-H2_initial_total_abundance', '0.6', + '-CO2_initial_total_abundance', '7.0', + '-CO_initial_total_abundance', '18.0', + ), + name='ic_abundance', + ) + + +@pytest.mark.physics_invariant +@pytest.mark.reference_pinned +def test_abundance_ic_realises_the_requested_inventory(abundance_ic_run): + """The abundance IC conserves the requested elemental inventories. + + Anchor: mass-balance identity through the initial-pressure solve. + Each volatile's initial_kg must equal the requested ppm of the + mantle mass exactly, and because the water and carbon dioxide + reactions conserve hydrogen and carbon, the realised reservoirs + carry the same H and C mole totals as the request. Oxygen is NOT + conserved among the volatiles: the reactions exchange it with the + melt's oxygen-fugacity buffer (observed 2e-3 relative for these + targets), which discriminates an inert fO2 pathway. + """ + doc = read_output(abundance_ic_run, 0) + mantle_kg = float(field_si(doc['atmosphere']['mass_mantle'])[0]) + + requested = {} + realised = {} + for volatile, ppm in ABUNDANCE_PPM.items(): + res = _volatile_reservoirs_kg(doc, volatile) + requested[volatile] = ppm * 1e-6 * mantle_kg + realised[volatile] = res['liquid_kg'] + res['solid_kg'] + res['atmosphere_kg'] + # rel=1e-9: initial_kg is the requested abundance times the + # mantle mass, stored rather than solved. + assert res['initial_kg'] == pytest.approx(requested[volatile], rel=1e-9) + # physical_kg is defined as the reservoir sum. + assert res['physical_kg'] == pytest.approx(realised[volatile], rel=1e-10) + + # Elemental conservation through the equilibrium solve: H and C + # totals match the request (observed to 3e-15 relative; rel=1e-9 + # leaves platform headroom). This is the discriminating check that + # the solve redistributed mass without creating or destroying it. + for element in ('H', 'C'): + assert _element_moles(realised, element) == pytest.approx( + _element_moles(requested, element), rel=1e-9 + ) + + # Oxygen exchange: the realised volatile O inventory departs from + # the requested one (observed 2.2e-3 relative; the threshold sits + # twentyfold below) because the fO2 buffer participates in both + # reactions. An inert fO2 pathway would conserve O to the same + # 1e-9 the H and C totals meet. + o_req = _element_moles(requested, 'O') + o_real = _element_moles(realised, 'O') + assert abs(o_real - o_req) > 1e-4 * o_req + + # Edge case: the reactions moved water mass at the IC (observed + # 3.6e-3 relative), so the realised H2O reservoirs sit resolvably + # off the bare request. + assert abs(realised['H2O'] - requested['H2O']) > 1e-4 * requested['H2O'] + + +@pytest.fixture(scope='module') +def ocean_moles_ic_run(cached_spider_run): + """Zero-step reaction run with the ocean-moles atmosphere IC. + + The mole counts are chosen so the implied abundances match the + abundance-IC configuration, keeping the equilibrium solve within + its convergence basin. Zero macro steps: the tests read only the + initial condition. + """ + return cached_spider_run( + opts_file='reaction.opts', + overrides=( + '-IC_ATMOSPHERE', '4', + '-nstepsmacro', '0', + '-n', '50', + '-H2O_initial_ocean_moles', '2.43', + '-H2_initial_ocean_moles', '0.0163', + '-CO2_initial_ocean_moles', '0.0087', + '-CO_initial_ocean_moles', '0.0352', + ), + name='ic_ocean_moles', + ) + + +@pytest.mark.physics_invariant +@pytest.mark.reference_pinned +def test_ocean_moles_ic_converts_moles_to_mass(ocean_moles_ic_run): + """IC_ATMOSPHERE 4 converts Earth-ocean mole counts to inventory mass. + + Anchor: the OCEAN_MOLES constant in constants.c (7.68894974e22 mol + per Earth ocean). Each volatile's initial_kg must equal + moles * OCEAN_MOLES * molar_mass, so the pin verifies both the + constant and the per-volatile molar mass wiring. + """ + doc = read_output(ocean_moles_ic_run, 0) + + for volatile, oceans in MOLES_OCEANS.items(): + res = _volatile_reservoirs_kg(doc, volatile) + expected = oceans * OCEAN_MOLES * MOLAR_MASS[volatile] # kg + # rel=1e-9: a closed-form product of three stored constants. + assert res['initial_kg'] == pytest.approx(expected, rel=1e-9) + + # Molar-mass discrimination: reading the CO2 inventory with the CO + # molar mass shifts the expected mass by the mass ratio (36%), far + # beyond the pin tolerance. + co2 = _volatile_reservoirs_kg(doc, 'CO2') + wrong = MOLES_OCEANS['CO2'] * OCEAN_MOLES * MOLAR_MASS['CO'] + assert abs(co2['initial_kg'] - wrong) > 0.3 * co2['initial_kg'] + + # Edge case: the smallest inventory (CO2, 0.0087 oceans) still + # produces a positive, finite mass at the expected scale. + assert 0 < co2['initial_kg'] < 1e21 # kg + assert np.isfinite(co2['initial_kg']) + + +@pytest.mark.physics_invariant +def test_abundance_and_ocean_mole_solves_agree(abundance_ic_run, ocean_moles_ic_run): + """The two atmosphere-IC routes converge to the same equilibrium. + + The ocean-mole counts are chosen to imply the abundance targets to + about 0.1 percent, so both configurations pose the same equilibrium + problem through different entry points and their SOLVED partial + pressures must agree. This validates the converged state itself: + a solve that stopped early, hit a floor, or landed on a different + root would separate the two routes far beyond the input rounding. + """ + doc_a = read_output(abundance_ic_run, 0) + doc_m = read_output(ocean_moles_ic_run, 0) + + for volatile in MOLAR_MASS: + p_a = float(field_si(doc_a['atmosphere'][volatile]['atmosphere_bar'])[0]) + p_m = float(field_si(doc_m['atmosphere'][volatile]['atmosphere_bar'])[0]) + # Positivity of the solved state itself (the in-code guard + # rejects negative solutions; this pins it from the outside). + assert p_a > 0 and p_m > 0, volatile + # rel=5e-3: the mole counts imply the abundances to about 7e-4 + # (observed agreement 3.5e-4 to 6.7e-4); the margin is sevenfold + # while a wrong-branch or non-converged solve differs at the + # tens-of-percent level. + assert p_m == pytest.approx(p_a, rel=5e-3), volatile + + # Scale guard: the equilibrium sits at bar-scale pressures for this + # inventory (observed 1.7 to 14 bar); a nondimensional leak or a + # Pa/bar slip would leave this bracket by orders of magnitude. + p_h2o = float(field_si(doc_a['atmosphere']['H2O']['atmosphere_bar'])[0]) + assert 0.1 < p_h2o < 100.0 # bar diff --git a/tests/test_poststep.py b/tests/test_poststep.py new file mode 100644 index 00000000..e399d79e --- /dev/null +++ b/tests/test_poststep.py @@ -0,0 +1,98 @@ +"""Tests for poststep.c (post-step event checks and rollback). + +Driven through short blackbody50 runs of the spider binary with the +post-step machinery active. Contract clauses exercised: a surface +temperature drop beyond tsurf_poststep_change rolls the offending step +back and stops the run early with the rolled-back state as the final +output, and activating the post-step check without rollback support is +refused. See docs/How-to/build_tests.md for the tier system. +""" + +from __future__ import annotations + +import pytest + +from tests._json_utils import atmosphere_si, read_output + +pytestmark = [pytest.mark.smoke, pytest.mark.timeout(120)] + +# blackbody50 cools by hundreds of kelvin in the first 100-year macro +# step, so a 1 K ceiling guarantees the event fires on step one. +TSURF_CEILING = 1.0 # K +NOMINAL_STEPS = 6 +NOMINAL_END = 600 # years, 6 steps of the default 100-year dtmacro + + +@pytest.fixture(scope='module') +def rollback_run(cached_spider_run): + """Six-step run whose first step exceeds the tsurf-change ceiling.""" + return cached_spider_run( + overrides=( + '-nstepsmacro', str(NOMINAL_STEPS), + '-activate_poststep', + '-activate_rollback', + '-tsurf_poststep_change', str(TSURF_CEILING), + ), + name='rollback_event', + ) + + +def test_tsurf_event_rolls_back_and_stops_early(rollback_run): + """A tsurf drop beyond the ceiling ends the run at the rolled-back state. + + The run must terminate before its nominal six macro steps: the + output directory holds the initial condition plus exactly one + further state, written at the rolled-back time inside the first + macro step, and the surface cooling across that window exceeds the + 1 K ceiling that triggered the event. + """ + outputs = sorted( + int(p.stem) for p in rollback_run.glob('*.json') if p.stem.isdigit() + ) + + # Early stop: the initial condition plus the rolled-back state. + assert len(outputs) == 2 + assert outputs[0] == 0 + # The final output falls inside the first nominal macro step, far + # short of the 600-year end of an uninterrupted run. + assert 0 < outputs[-1] < NOMINAL_END + assert outputs[-1] < 100 # years, within macro step one + + # Rollback semantics: the offending step is discarded, so the + # written final state is from BEFORE the ceiling was crossed. The + # cooling across the written window is therefore nonzero but + # bounded by the ceiling (observed 0.99 K against the 1 K limit); + # a run that wrote the over-limit state would exceed it by + # hundreds of kelvin. + t_start = atmosphere_si(read_output(rollback_run, 0), 'temperature_surface') + t_end = atmosphere_si(read_output(rollback_run, outputs[-1]), 'temperature_surface') + drop = t_start - t_end # K + assert 0 < drop <= TSURF_CEILING + # The integrator ran close to the ceiling before the event fired. + assert drop > 0.1 * TSURF_CEILING + # Positivity: the rolled-back state is still physical. + assert t_end > 0 # K + + +def test_poststep_without_rollback_is_refused(run_spider): + """The post-step check requires rollback support and errors without it. + + poststep.c raises a hard error when the first post-step evaluation + finds rollback inactive, so the run must exit nonzero rather than + silently continue without the safety net. + """ + with pytest.raises(RuntimeError) as excinfo: + run_spider( + overrides=( + '-nstepsmacro', '2', + '-activate_poststep', + '-tsurf_poststep_change', str(TSURF_CEILING), + ), + name='poststep_no_rollback', + ) + message = str(excinfo.value) + # The documented refusal fired, not an arbitrary crash: the PETSc + # error text names the missing option. (The exit code itself + # differs by platform, so the message is the stable contract.) + assert 'You must run with -activate_rollback' in message + assert 'spider exited with code' in message diff --git a/tests/test_reaction.py b/tests/test_reaction.py index 3c12164f..bf4d11a1 100644 --- a/tests/test_reaction.py +++ b/tests/test_reaction.py @@ -147,3 +147,131 @@ def test_elements_conserved_across_each_reaction_pair(reaction_run): # inventory; both bounds discriminate unit slips. assert 1e16 < abs(h2o['reaction_kg']) < 1e20 assert abs(h2o['reaction_kg']) < 0.01 * h2o['initial_kg'] + + +# The full named-reaction library configuration: seven volatiles and +# the IVTANTHERMO water, IVTANTHERMO carbon dioxide, IVTANTHERMO +# methane, and ammonia reactions (tests/opts/reaction_library.opts). +LIBRARY_VOLATILES = ('H2O', 'H2', 'CO2', 'CO', 'CH4', 'NH3', 'N2') +LIBRARY_MOLAR_MASS = { + 'H2O': 0.01801528, + 'H2': 0.00201588, + 'CO2': 0.04401, + 'CO': 0.02801, + 'CH4': 0.01604, + 'NH3': 0.017031, + 'N2': 0.028014, +} +LIBRARY_STOICHIOMETRY = { + 'H2O': {'H': 2, 'O': 1}, + 'H2': {'H': 2}, + 'CO2': {'C': 1, 'O': 2}, + 'CO': {'C': 1, 'O': 1}, + 'CH4': {'C': 1, 'H': 4}, + 'NH3': {'N': 1, 'H': 3}, + 'N2': {'N': 2}, +} + + +@pytest.fixture(scope='module') +def library_run(cached_spider_run): + """One macro step of the full named-reaction library configuration.""" + return cached_spider_run( + opts_file='reaction_library.opts', + overrides=('-nstepsmacro', '1', '-n', '50'), + name='reaction_library', + ) + + +def _library_element_moles(doc, element): + """Moles of one element summed over all seven volatile inventories.""" + total = 0.0 + for volatile in LIBRARY_VOLATILES: + block = doc['atmosphere'][volatile] + kg = sum( + float(field_si(block[key])[0]) + for key in ('liquid_kg', 'solid_kg', 'atmosphere_kg') + ) + moles = kg / LIBRARY_MOLAR_MASS[volatile] + total += moles * LIBRARY_STOICHIOMETRY[volatile].get(element, 0) + return total + + +@pytest.mark.physics_invariant +def test_reaction_library_conserves_h_c_n(library_run): + """The four-reaction network conserves hydrogen, carbon, and nitrogen. + + The methane reaction exchanges C and H, the ammonia reaction N and + H, and the water and carbon dioxide reactions H and C, but none of + them creates or destroys those elements, so the molar totals over + all seven volatiles are the same before and after the macro step. + """ + doc_0 = read_output(library_run, 0) + doc_1 = read_output(library_run, 1000) + + # rel=1e-5: the coupled volatile solve conserves the totals to its + # own tolerance (observed drifts at or below 1.5e-6). + for element in ('H', 'C', 'N'): + total_0 = _library_element_moles(doc_0, element) + total_1 = _library_element_moles(doc_1, element) + assert total_1 == pytest.approx(total_0, rel=1e-5), element + assert total_0 > 0 + + # Equilibrium partitioning discrimination: at these conditions the + # ammonia reaction pushes nearly all nitrogen into N2, so the + # realised NH3 share of the N inventory is tiny even though 0.5 ppm + # of NH3 was requested. An inert ammonia reaction would leave the + # requested 20 percent molar share in place. + n2_block = doc_0['atmosphere']['N2'] + n2_kg = sum( + float(field_si(n2_block[key])[0]) + for key in ('liquid_kg', 'solid_kg', 'atmosphere_kg') + ) + n_total = _library_element_moles(doc_0, 'N') + assert 2.0 * n2_kg / LIBRARY_MOLAR_MASS['N2'] > 0.99 * n_total + + # Edge case: the trace species (NH3) still carries a positive, + # finite inventory through the equilibrium. + nh3_kg = sum( + float(field_si(doc_0['atmosphere']['NH3'][key])[0]) + for key in ('liquid_kg', 'solid_kg', 'atmosphere_kg') + ) + assert 0 < nh3_kg < 1e21 # kg + + +@pytest.mark.physics_invariant +def test_oxygen_flows_through_the_melt_buffer(library_run): + """Oxygen is exchanged with the melt fO2 buffer, not conserved. + + The water, carbon dioxide, and methane reactions carry an + oxygen-fugacity stoichiometry (the ammonia reaction does not), so + the volatile O inventory drifts over a step (observed 5e-3 + relative in 1000 years) while H stays conserved in the same + window. A conserved O total would mean the fO2 coupling is inert, + which this test is designed to expose. + """ + doc_0 = read_output(library_run, 0) + doc_1 = read_output(library_run, 1000) + + o_0 = _library_element_moles(doc_0, 'O') + o_1 = _library_element_moles(doc_1, 'O') + assert o_0 > 0 + # The fO2 exchange must move a resolvable amount of oxygen. + assert abs(o_1 - o_0) / o_0 > 1e-4 + + # Contrast: hydrogen is conserved in the very same step, so the + # oxygen drift is the fO2 pathway and not a global mass leak. + h_0 = _library_element_moles(doc_0, 'H') + h_1 = _library_element_moles(doc_1, 'H') + assert h_1 == pytest.approx(h_0, rel=1e-5) + + # Positivity to solver tolerance: the coupled volatile solve has no + # hard positivity clamp, so a trace reservoir can undershoot zero + # by roundoff (observed -7e9 kg of CH4 against its 4.2e18 kg + # inventory, i.e. -2e-9 relative). Bound the undershoot at 1e-6 of + # the species inventory instead of asserting a hard zero. + for volatile in LIBRARY_VOLATILES: + block = doc_1['atmosphere'][volatile] + floor = -1e-6 * float(field_si(block['initial_kg'])[0]) + for key in ('liquid_kg', 'solid_kg', 'atmosphere_kg'): + assert float(field_si(block[key])[0]) >= floor, (volatile, key) diff --git a/tools/generate_test_badges.py b/tools/generate_test_badges.py index b6121f65..5b2007ed 100644 --- a/tools/generate_test_badges.py +++ b/tools/generate_test_badges.py @@ -6,11 +6,18 @@ {"schemaVersion": 1, "label": "