From fbd38535d753d4b75715f03af9c04417103819c8 Mon Sep 17 00:00:00 2001 From: timlichtenberg Date: Sun, 12 Jul 2026 19:44:24 +0200 Subject: [PATCH 1/6] Fix the initial atmosphere solve for abundance and ocean-moles starts The initial partial-pressure solver forced a matrix-free Jacobian onto its trust-region Newton method, but PETSc 3.19's trust-region step needs a matrix that supports transpose products, which the matrix-free approximation does not provide. Every run using IC_ATMOSPHERE 1 (abundances) or 4 (ocean moles) therefore aborted with "Matrix type mffd does not have a multiply transpose defined" before taking a step. The system is tiny, one unknown per volatile plus one per reaction, so build an explicit finite-difference Jacobian instead; it supports the transpose product and costs a handful of residual evaluations. --- ic.c | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/ic.c b/ic.c index e615cd34..d5eb791b 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 */ From 4f691d6ba6f3efca86dfc5e89ba44219bd2f1641 Mon Sep 17 00:00:00 2001 From: timlichtenberg Date: Sun, 12 Jul 2026 19:45:17 +0200 Subject: [PATCH 2/6] Cover the optional physics paths with invariant-pinned smoke tests The suite exercised only the default blackbody and reaction configurations, leaving the option-dispatched physics dark: radiogenic and tidal heating, the steady-state initial condition, constant-entropy and prescribed-core-flux boundaries, the abundance and ocean-moles atmosphere starts, the IVTANTHERMO and ammonia reaction constructors, and the post-step rollback event. Line coverage of the C sources rises from 81 to 91 percent. New cached smoke runs pin each path against an analytical anchor: two-isotope radiogenic heating against the exponential decay law, prescribed and file-based tidal heating against the option value and a per-node ramp, the steady-state start against a node-independent energy flow and the grey-body surface balance, the entropy boundary against the exact option value, the core-flux boundary against the core-cooling baseline, the abundance start against elemental hydrogen and carbon conservation through the equilibrium solve (with oxygen exchanged through the melt fO2 buffer, which the tests require rather than forbid), the ocean-moles start against the OCEAN_MOLES constant, and a four-reaction, seven-volatile configuration against elemental conservation. The rollback test verifies that a surface-temperature event ends the run at the rolled-back state, and that the post-step check refuses to run without rollback support. The new anchors are recorded on the validation pages for energy.c, bc.c, and ic.c. Coverage counters and reports from instrumented builds are ignored. --- .gitignore | 6 + docs/Validation/bc.md | 8 ++ docs/Validation/energy.md | 16 +++ docs/Validation/ic.md | 16 +++ tests/opts/reaction_library.opts | 240 +++++++++++++++++++++++++++++++ tests/test_bc.py | 91 ++++++++++++ tests/test_energy.py | 201 ++++++++++++++++++++++++++ tests/test_ic.py | 170 +++++++++++++++++++++- tests/test_poststep.py | 98 +++++++++++++ tests/test_reaction.py | 127 ++++++++++++++++ 10 files changed, 972 insertions(+), 1 deletion(-) create mode 100644 tests/opts/reaction_library.opts create mode 100644 tests/test_poststep.py 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/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..44be00f3 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 a factor of about 300; sign and scale bounds on the flux. diff --git a/docs/Validation/ic.md b/docs/Validation/ic.md index 355f8f0c..2bf68ba4 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 moves by a factor of about two against the request because the reactions exchange O with the melt's fO2 buffer; the realised H2O reservoirs differ from the bare request by more than half, so an inert reaction network fails both checks. + +**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/tests/opts/reaction_library.opts b/tests/opts/reaction_library.opts new file mode 100644 index 00000000..eccdbf36 --- /dev/null +++ b/tests/opts/reaction_library.opts @@ -0,0 +1,240 @@ +# 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_initial_atmos_pressure 1.8412864823803783E5 +-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_initial_atmos_pressure 3.3703436080366105E5 +-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_initial_atmos_pressure 15.555831058257253E5 +-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_initial_atmos_pressure 4.27087847230405E5 +-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 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_energy.py b/tests/test_energy.py index 28d8e445..5c98b951 100644 --- a/tests/test_energy.py +++ b/tests/test_energy.py @@ -154,3 +154,204 @@ 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', + '-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 at ~1960 K misses by a factor of ~300. + 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_ic.py b/tests/test_ic.py index 0d5cc0b8..98a4088c 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,171 @@ 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. +ABUNDANCE_PPM = {'H2O': 100.0, 'H2': 20.0, 'CO2': 50.0, 'CO': 10.0} +# Requested inventories (Earth oceans) for the ocean-moles IC; 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': 0.29, 'H2': 0.52, 'CO2': 0.0595, 'CO': 0.0187} + + +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): + """One-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. + """ + return cached_spider_run( + opts_file='reaction.opts', + overrides=( + '-IC_ATMOSPHERE', '1', + '-nstepsmacro', '1', + '-n', '50', + '-H2O_initial_total_abundance', '100.0', + '-H2_initial_total_abundance', '20.0', + '-CO2_initial_total_abundance', '50.0', + '-CO_initial_total_abundance', '10.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, roughly doubling the volatile O + inventory, 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 sits far from + # the requested one (observed factor 2.1) because the fO2 buffer + # participates in both reactions. + o_req = _element_moles(requested, 'O') + o_real = _element_moles(realised, 'O') + assert abs(o_real - o_req) > 0.5 * o_req + + # Edge case: the reactions moved water mass at the IC, so the + # realised H2O reservoirs differ strongly from the bare request. + assert abs(realised['H2O'] - requested['H2O']) > 0.5 * requested['H2O'] + + +@pytest.fixture(scope='module') +def ocean_moles_ic_run(cached_spider_run): + """One-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. + """ + return cached_spider_run( + opts_file='reaction.opts', + overrides=( + '-IC_ATMOSPHERE', '4', + '-nstepsmacro', '1', + '-n', '50', + '-H2O_initial_ocean_moles', '0.29', + '-H2_initial_ocean_moles', '0.52', + '-CO2_initial_ocean_moles', '0.0595', + '-CO_initial_ocean_moles', '0.0187', + ), + 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 (CO, 0.0187 oceans) still + # produces a positive, finite mass at the expected scale. + co = _volatile_reservoirs_kg(doc, 'CO') + assert 0 < co['initial_kg'] < 1e21 # kg + assert np.isfinite(co['initial_kg']) diff --git a/tests/test_poststep.py b/tests/test_poststep.py new file mode 100644 index 00000000..e87e2670 --- /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 refusal is the documented unsupported-configuration error + # (PETSC_ERR_SUP, exit code 56), not a solver convergence crash + # (which aborts with code 76) or a clean exit. + assert 'spider exited with code 56' in message + assert '-activate_poststep' in message # the offending option, echoed in the command diff --git a/tests/test_reaction.py b/tests/test_reaction.py index 3c12164f..3f31f0c0 100644 --- a/tests/test_reaction.py +++ b/tests/test_reaction.py @@ -147,3 +147,130 @@ 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. + + Every reaction in the library carries an oxygen-fugacity + stoichiometry, 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) From 39060a2d03ddee37949297e27f7217c533a8733d Mon Sep 17 00:00:00 2001 From: timlichtenberg Date: Sun, 12 Jul 2026 19:57:39 +0200 Subject: [PATCH 3/6] Cover the field wrapper and viscosity laws; publish per-tier badge counts A C test executable exercises the DimensionalisableField wrapper (scale, unscale, duplicate, local vector, scaling query) with exact round-trip pins, and the equation-of-state tests pin the compositional viscosity prefactor against all four Mg/Si branches of Spaargaren et al. (2020) plus the activation-volume and activation-energy behavior, including the zero-pressure surface limit. Fast-tier line coverage reaches 91 percent locally and the full suite 92 percent, so the fast gate rises to 88 and the full gate to the 90 percent ecosystem ceiling. The badge generator now also publishes per-tier counts (tests-unit, tests-smoke, tests-integration) alongside the total and the CI-split fast and nightly files, matching the badge scheme of the other ecosystem modules; the publish workflow validates all six files and the testing explainer shows the per-tier badges. --- .github/workflows/publish-test-badges.yml | 9 +- Makefile | 1 + docs/Explanations/testing.md | 5 +- pyproject.toml | 4 +- tests/c/test_dimensionalisablefield.c | 84 +++++++++++++++ tests/test_dimensionalisablefield.py | 60 +++++++++++ tests/test_eos.py | 120 ++++++++++++++++++++++ tools/generate_test_badges.py | 18 +++- 8 files changed, 292 insertions(+), 9 deletions(-) create mode 100644 tests/c/test_dimensionalisablefield.c create mode 100644 tests/test_dimensionalisablefield.py 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/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/pyproject.toml b/pyproject.toml index c896ef20..b96ab854 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -27,10 +27,10 @@ markers = [ # difference between compiler instrumentation variants; the ratchet # raises them from CI-measured coverage. [tool.spider.coverage_fast] -fail_under = 75.0 +fail_under = 88.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..ce278e54 --- /dev/null +++ b/tests/c/test_dimensionalisablefield.c @@ -0,0 +1,84 @@ +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, 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); + + 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, "\"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/test_dimensionalisablefield.py b/tests/test_dimensionalisablefield.py new file mode 100644 index 00000000..ecd57559 --- /dev/null +++ b/tests/test_dimensionalisablefield.py @@ -0,0 +1,60 @@ +"""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) + # No-op guard: scaling must actually change the stored values. + assert abs(out['scaled0'] - out['original0']) > 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_eos.py b/tests/test_eos.py index 8247934d..143f9ebb 100644 --- a/tests/test_eos.py +++ b/tests/test_eos.py @@ -84,3 +84,123 @@ 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) # +0.2074 + shift_low = _spaargaren_prefactor(1.2) - _spaargaren_prefactor(0.9) # -1.3580 + # 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/tools/generate_test_badges.py b/tools/generate_test_badges.py index b6121f65..1efddc97 100644 --- a/tools/generate_test_badges.py +++ b/tools/generate_test_badges.py @@ -6,11 +6,18 @@ {"schemaVersion": 1, "label": "