diff --git a/AGENTS.md b/AGENTS.md index 0a8e034e..34a40d7c 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -137,11 +137,13 @@ numbers in both runs; keep it that way and do not reseed from the clock. | `hfls` | latent heat flux, positive away from the surface | W/m2 | | `hfss` | sensible heat flux, positive away from the surface | W/m2 | | `hfg` | ground heat flux at the actual soil surface, positive into the ground; a flux taken below the surface must be corrected for heat storage above that depth | W/m2 | +| `rlus` | total upward longwave radiation at the surface: emission plus reflected downward longwave, positive away from the surface; row i's value is at row i's `time`, the same instant as row i's `rlds` | W/m2 | | `mrso` | soil water storage | mm | | `snw` | snow water equivalent | mm | | `canopy` | canopy interception storage | mm | | `gw` | groundwater storage below the soil column | mm | | `channel` | water generated as runoff but not yet released by the model's routing | mm | +| `ts` | surface (skin) temperature at row i's `time`, the same instant as row i's `rlds`; a diagnostic, declared under `emits.diagnostics`, neither integrated nor differenced by any budget | K | For `hfg`, an adapter mapping a plate-depth or deeper-boundary flux must use `G_surface = G_depth + (E_above_end - E_above_start) / dt`, with downward @@ -200,6 +202,19 @@ cannot be put to the model, so it counts neither way. Fabricating an `evspsbl` column to avoid `INCOMPLETE` produces `VIOLATION` instead, which is worse and is also dishonest. +`needs_forcing` and `needs_static` declare inputs the adapter cannot run +without; a missing input makes the case `N/A (INCOMPATIBLE)`. +`uses_forcing` and `uses_static` declare optional inputs: the adapter must +consume them whenever supplied, but can run without them using a documented +fallback. A probe's `requires.forcing` and `requires.static` accept either +declaration. For example, `energy/radiation-consistency` requires consumption +of `rlds` and `eps`; a model that does not declare it consumes both is +`N/A (INCOMPATIBLE)` because it may be computing its own sky or emissivity. + +Declare diagnostic outputs under the optional key `diagnostics: [ts]`. +Criteria read diagnostics, but budgets never integrate or difference them; +a surface temperature must not be summed into water storage. + ## Verify before you submit ```bash @@ -211,7 +226,7 @@ ht run --model # the actual evaluation use, and checks the shape of what came back. Get that green before looking at any residual. Without `--probe` it checks the closure probe, or the first probe the model can consume when it cannot consume that one: a step it does -not declare, a forcing the probe does not generate, or a window that drops a +not declare, a forcing or static input the probe does not generate, or a window that drops a stretch the probe scores makes a probe N/A (INCOMPATIBLE) for the model. When that is true of the probe named with `--probe`, or of every probe, the adapter is not run and the command exits 1, as `ht run` does for a model no diff --git a/CITATION.cff b/CITATION.cff index d03f9d2f..a4e6fb2c 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -39,6 +39,13 @@ authors: Laboratory of Catchment Hydrology and Geomorphology, École Polytechnique Fédérale de Lausanne (EPFL), 1951 Sion, Switzerland orcid: "https://orcid.org/0009-0001-9658-4375" + - given-names: Xin + family-names: Lan + affiliation: >- + Center for Systems Integration and Sustainability, Department of + Fisheries and Wildlife, Michigan State University, East Lansing, MI + 48823, USA + orcid: "https://orcid.org/0000-0002-0607-2270" repository-code: "https://github.com/Flood-Lab/HydroTuring" url: "https://flood-lab.github.io/HydroTuring/" license: MIT diff --git a/CONTRIBUTORS.md b/CONTRIBUTORS.md index 79afeef0..51677d37 100644 --- a/CONTRIBUTORS.md +++ b/CONTRIBUTORS.md @@ -35,6 +35,7 @@ long. | `energy/latent-heat-et-consistency` | Changming Li (SCUT) | | `energy/evaporative-partition` | Changming Li (SCUT) | | `energy/surface-energy-closure` | Han Wang (The Hong Kong University of Science and Technology) | +| `energy/radiation-consistency` | Xin Lan (Michigan State University) | | `momentum/routing-conservation` | Zhi Li (CU Boulder) | ## Models diff --git a/README.md b/README.md index 1ddcb7e7..1824903f 100644 --- a/README.md +++ b/README.md @@ -68,15 +68,17 @@ not. ## The probes -Twenty-one: sixteen under mass, four under energy and one under momentum. Each was -merged only after the acceptance gate saw it pass four physical models, a -bucket that conserves water exactly, two hand-written FLEX models and the -NWS's SAC-SMA with Snow-17, and fail a purpose-built broken one on the -named criterion. A probe that fails a +Twenty-two: sixteen under mass, five under energy and one under momentum. +Each was merged only after the acceptance gate saw it pass its declared +exact reference and fail a purpose-built broken one on the named criterion. +Four physical models, a bucket that conserves water exactly, two +hand-written FLEX models and the NWS's SAC-SMA with Snow-17, must pass +every probe that can ask them anything; four of the five energy probes need +outputs they do not report and are not scored for them. A probe that fails a physical model is examined before the model is; that is the first thing done -with any probe pull request. Eleven of the twenty-one can be scored on a model +with any probe pull request. Eleven of the twenty-two can be scored on a model that reports runoff and nothing else. `ht list` prints them; -[ROADMAP.md](ROADMAP.md#probes-we-want) has the eleven more we want, all +[ROADMAP.md](ROADMAP.md#probes-we-want) has the ten more we want, all unclaimed. | Probe | Law | What it asks | The broken model it catches | @@ -101,6 +103,7 @@ unclaimed. | [`energy/latent-heat-et-consistency`](probes/energy/latent-heat-et-consistency) | energy | The evaporation a model reports as water and the evaporation implied by the latent heat it reports: are they the same evaporation? | `reference_two_head`, `reference_constant_lambda`, `reference_sublimation_blind`, `reference_energy_leak` | | [`energy/evaporative-partition`](probes/energy/evaporative-partition) | energy | One summer without rain under net radiation that did not change: the latent heat a drying surface gives up has to warm the air. | `reference_two_head`, `reference_ground_dodge` | | [`energy/surface-energy-closure`](probes/energy/surface-energy-closure) | energy | Does each day and night close its hourly surface energy budget, without opposite errors cancelling? | `reference_diurnal_bias` | +| [`energy/radiation-consistency`](probes/energy/radiation-consistency) | energy | The surface temperature a model reports and the upward longwave it reports: do they describe one surface, hour by hour, at the emissivity it was given? | `reference_air_emitter`, `reference_no_reflection` | | [`momentum/routing-conservation`](probes/momentum/routing-conservation) | momentum | The channel store is never negative and never holds more than its hydrograph can. | `reference_stuck_router` | ## Models @@ -113,11 +116,12 @@ hand-written conceptual models from [chrimerss/HydrologicModels](https://github.com/chrimerss/HydrologicModels) and the NWS's SAC-SMA with Snow-17 must pass every probe that can ask them anything, so a probe that fails one is wrong until shown otherwise. All four -report water and no energy, so all four are N/A, with reason INCOMPLETE, on the three probes -that need the latent, sensible and ground heat fluxes, -`energy/latent-heat-et-consistency`, `energy/evaporative-partition` and -`energy/surface-energy-closure`, -rather than passing them: they have not +report water and no energy, so all four are N/A, with reason INCOMPLETE, on +the four probes that need an energy output: the three that need the latent, +sensible and ground heat fluxes, `energy/latent-heat-et-consistency`, +`energy/evaporative-partition` and `energy/surface-energy-closure`, and +`energy/radiation-consistency`, which needs a surface temperature and its +upward longwave. Not scored rather than passing: they have not violated conservation of energy, they have declined to be falsifiable about it, exactly as `reference_streamflow_only` does on the budget probes. A broken model is broken in one specific way, so that no criterion goes untested. @@ -134,8 +138,8 @@ between models. | [`wflow_sbm`](models/wflow_sbm) | submitted | Deltares' Wflow.jl SBM: a soil column with unsaturated and saturated stores, interception, snow and kinematic-wave routing, adapted in Julia and run on one representative cell at the resolution of Wflow's Moselle model. | **FAIL (VIOLATION)**, 15 of 18 probes passed. Its water budget closes to 2e-4 of the rain and it passes the area, routing, step, calendar, causality, extreme-rain, phase, precipitation-counterfactual, warming and human-abstraction probes, the last by taking the prescribed withdrawal through Wflow's own water demand and allocation. That 2e-4 is water its river kinematic wave creates when open-water evaporation exceeds the rain on the river and the wave floors its discharge rather than drying the channel; on `mass/extreme-event-closure` it is 0.0016 mm on one 0.006 mm drizzle day on one seed, past that probe's 0.001 mm floor, the only one of the probe's events to fail. Its other two failures both come from its soil column: a storm that does not fill the column makes almost no runoff, and under the shipped mapping a wet month's extra water has already been evaporated down to the rooting depth when the storm arrives, so the wetter catchment runs off barely more than the dry one, about a thousandth of the storm where the probe asks for two hundredths; and once the root zone dries in a rainless spell, drainage from the unsaturated store raises a water table that sits below the roots, and lateral flow rises with no rain. Both move with how the stated soil capacity is mapped onto soil thickness and roots. | | [`summa`](models/summa) | submitted | SUMMA 4.0.0, a land model that solves the water and energy balances of canopy, snow, soil and aquifer with one implicit solver, run as one lumped HRU from its shipped test-case setup, with the radiation and humidity it needs mocked from each forcing row. | **FAIL (VIOLATION)**, 8 of 21 probes passed. Its water budget closes to rounding and its surface energy budget to 0.2% of its own net radiation; what fails is a latent heat of vaporisation held at its 0 C value, a net radiation of its own that is not the probe's `rn` (every hourly phase, and a drought that warms the surface), a canopy that, as the shipped setup configures SUMMA, intercepts all rain and keeps what freezes near 0 C, holding up to 27 mm of ice against a 2 mm capacity on the probes that score it (also the only failure on the precipitation counterfactual), leaf area that keeps its seasons under constant weather, runoff that follows the melt rather than the rain on one seed, a 5.2% runoff response to turning snow into rain on another, Green-Ampt runoff that appears only at the hourly step, and no human water use, so the prescribed withdrawal on the human-abstraction probe is never taken | | [`cwatm`](models/cwatm) | submitted | CWatM 1.11, IIASA's Community Water Model and an ISIMIP global hydrological model, run on one grid cell with every store it carries reported. | **FAIL (VIOLATION)**, 14 of 18 probes passed. Its budget closes to 0.03%, and its own water-demand module pumps a prescribed withdrawal out of groundwater to within 0.004%. That 0.03% is water its capillary rise creates, a few thousandths of a millimetre a day, and on `mass/extreme-event-closure` it fails 3 to 8 one-day drizzle events of under 0.07 mm per seed, where 0.001 to 0.005 mm more leaves or is stored than fell, against that probe's allowance of 5% of the rain or 0.001 mm, whichever is larger. It also fails on a groundwater reservoir with no dt that drains 24 times too fast at an hourly step (52% of the rain); on a 0.29 mm/day dip after an added storm, where preferential flow turns surface runoff into interflow that leaves through a slower runoff-concentration lag; and on evaporation on wet soil at 0.70 of demand, near the cap its crop coefficients set. The last two are packaging choices as much as results: this package follows the CWatM-Earth-30min template its parameters come from, and with `preferentialFlow = False`, the setting of the pinned model repository's own 30′ templates, both pass | -| [`lisflood`](models/lisflood) | submitted | LISFLOOD 5.0.0, the EC Joint Research Centre's distributed model behind EFAS and GloFAS, stepped through its own Python framework on one representative 5 km cell and reporting every store its own water balance module counts. | **FAIL (ERROR)**, 15 of 18 probes passed. Its budget closes to 1e-13 mm per step; it reports no heat fluxes, so the three energy-flux probes cannot ask it anything and are N/A. It fails the step probe because its potential infiltration is a pore-space storage multiplied by the step length, so rain falling within hours runs off at the hourly step (13.0% of the rain between PT1H and PT1D). The two ten-year probes take about 90 s per run under amd64 emulation against a 60 s budget, so their rows are ERROR from the host's speed, and those two ERROR rows alone make the verdict FAIL (ERROR). Run outside the limit, `mass/precipitation-counterfactual` passes every criterion. On `mass/human-abstraction`, LISFLOOD's own water-use rule takes the groundwater share in full and the rest only from channel water above an environmental-flow reserve, recording what the channel cannot give as shortage. On this one-cell water region it withdraws 17 to 33% of the prescription, from the sourced reserve to none, and leaves 67 to 83% where the probe allows 5% | -| `reference_bucket` | exact | conserves water exactly by construction | must pass every probe that can ask it anything; INCOMPLETE on the three energy-flux probes, which need fluxes it does not report | +| [`lisflood`](models/lisflood) | submitted | LISFLOOD 5.0.0, the EC Joint Research Centre's distributed model behind EFAS and GloFAS, stepped through its own Python framework on one representative 5 km cell and reporting every store its own water balance module counts. | **FAIL (ERROR)**, 15 of 18 probes passed. Its budget closes to 1e-13 mm per step; it reports no heat fluxes and no surface temperature, so the four energy probes that need an energy output cannot ask it anything and are N/A. It fails the step probe because its potential infiltration is a pore-space storage multiplied by the step length, so rain falling within hours runs off at the hourly step (13.0% of the rain between PT1H and PT1D). The two ten-year probes take about 90 s per run under amd64 emulation against a 60 s budget, so their rows are ERROR from the host's speed, and those two ERROR rows alone make the verdict FAIL (ERROR). Run outside the limit, `mass/precipitation-counterfactual` passes every criterion. On `mass/human-abstraction`, LISFLOOD's own water-use rule takes the groundwater share in full and the rest only from channel water above an environmental-flow reserve, recording what the channel cannot give as shortage. On this one-cell water region it withdraws 17 to 33% of the prescription, from the sourced reserve to none, and leaves 67 to 83% where the probe allows 5% | +| `reference_bucket` | exact | conserves water exactly by construction | must pass every probe that can ask it anything; N/A on the four energy probes that need outputs it does not report | | [`flex_lumped`](models/flex_lumped) | physical | lumped FLEX/HBV: interception, beta-partitioned unsaturated store, fast and slow reservoirs, triangular lag | must pass every probe that can ask it anything; **PASS**, 18 of 18 | | [`flex_topo`](models/flex_topo) | physical | FLEX-Topo: plateau, hillslope and wetland units on real Wark fractions sharing one groundwater store | must pass every probe that can ask it anything; **PASS**, 18 of 18 | | [`sacsma_snow17`](models/sacsma_snow17) | physical | the NWS's SAC-SMA with Snow-17 and a gamma unit hydrograph, ported from the legacy Fortran and checked against it | must pass every probe that can ask it anything; **PASS**, 18 of 18 | @@ -146,6 +150,9 @@ between models. | `reference_sublimation_blind` | broken | coherent for liquid water, but converts snow sublimation at the latent heat of vaporisation instead of sublimation | caught by `flux_identity` | | `reference_ground_dodge` | broken | coherent, both budgets close, and the sensible flux never reads the soil: when the soil dries the ground flux absorbs the whole shift | caught by `partition_shift` | | `reference_diurnal_bias` | broken | shifts sensible heat so the energy residual is +20 W/m2 by day and -20 W/m2 by night; the full-record residual cancels | caught by `energy_closure_by_phase` | +| `reference_radiative` | exact | the coupled reference with a skin: temperature diagnosed from the sensible heat flux through a fixed conductance, upward longwave that skin's emission plus the reflected sky at the emissivity it was given | must pass every criterion of `energy/radiation-consistency` | +| `reference_air_emitter` | broken | reports the skin's temperature but emits at the air's, reflected sky unchanged | caught by `radiative_identity` | +| `reference_no_reflection` | broken | reports emission alone as the total upward longwave, the reflected sky left out: 3 to 23 W/m2 that a 5% tolerance would not reliably see | caught by `radiative_identity` | | `reference_energy_leak` | broken | discards 15% of net radiation; the energy counterpart of `reference_leaky` | caught by `energy_closure` | | `reference_leaky` | broken | hides a silent 15% sink | caught by `closure` | | `reference_cheater` | broken | solves for storage as whatever balances the budget | caught by `state_bounds` | @@ -196,7 +203,8 @@ a pass nor a fail. That happens two ways: - `INCOMPLETE` the model never reported enough to be checked. Every streamflow-only model lands here on the budget probes. It has not violated conservation; it has declined to be falsifiable. -- `INCOMPATIBLE` the model and probe disagree on timestep, required forcing or +- `INCOMPATIBLE` the model and probe disagree on timestep, on a forcing or + static input one needs and the other does not supply or declare, or on paired-perturbation support, so running them would not be meaningful. A model passes when at least one probe could be put to it and every probe @@ -331,7 +339,7 @@ where that conversation happens, before and alongside the issues. ## Status -Suite `0.1.0`, pre-release. Twenty-one probes, sixteen mass, four energy, one momentum, synthetic track only. More +Suite `0.1.0`, pre-release. Twenty-two probes, sixteen mass, five energy, one momentum, synthetic track only. More energy and momentum probes, and the real-data track, are next. The harness runs paired cases and scores labelled regimes, so the generalisation probes on the roadmap — extrapolation in space and time, counterfactual response, invariance — are diff --git a/ROADMAP.md b/ROADMAP.md index af9ba9d9..46653236 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -195,11 +195,16 @@ must not occur while the pack is below freezing. *Discriminates:* models that melt snow on a warm day regardless of whether the pack has the energy to melt. -### `energy/radiation-consistency` · standard · **unclaimed** -Outgoing longwave must be consistent with the reported surface temperature -through Stefan-Boltzmann, given emissivity. -*Discriminates:* models that predict surface temperature and radiation with -separate heads that never have to agree. +### `energy/radiation-consistency` · **merged** +The surface temperature a model reports and the upward longwave it reports +must describe one gray surface at the emissivity it was given, hour by hour: +`rlus = eps sigma ts^4 + (1 - eps) rlds` within 0.5 percent of the reported +flux, with a 0.5 W/m2 floor. A same-instant identity, so nothing cancels +across hours; net radiation stays a prescribed forcing. +*Discriminates:* a model whose temperature and radiation heads never have to +agree, through `reference_air_emitter`, which emits at the air temperature, +and `reference_no_reflection`, which drops the reflected sky. +Contributed by Xin Lan (Michigan State University). --- diff --git a/docs/adapting-a-model.md b/docs/adapting-a-model.md index e2ed602b..76a71323 100644 --- a/docs/adapting-a-model.md +++ b/docs/adapting-a-model.md @@ -88,7 +88,7 @@ from a malformed table tells you nothing. It runs on the closure probe, or on the first probe your model can consume when it cannot consume that one: a model driven by net radiation is checked on an energy probe, because the closure probe generates none. A probe the -model cannot consume, at its step, with its forcing or on its window, is N/A +model cannot consume, at its step, with its forcing or static inputs, or on its window, is N/A (INCOMPATIBLE) for it. When that is true of the probe named with `--probe`, or of every probe, the adapter is not run and the command exits 1, as `ht run` does for a model no probe could score. Exit 0 is a contract that diff --git a/docs/writing-a-probe.md b/docs/writing-a-probe.md index 74a7e136..abfe3064 100644 --- a/docs/writing-a-probe.md +++ b/docs/writing-a-probe.md @@ -113,6 +113,7 @@ Every one is binary. | `phase_invariance` | the same water as rain instead of snow leaves the integrated volumes within a share of the rain | paired runs | | `demand_consistency` | evaporation reaches demand when the model's own soil is wettest, never exceeds it, and falls when driest | one run | | `energy_closure_by_phase` | the mean absolute surface-energy residual in each contiguous day or night stays within the larger of the relative radiation tolerance and the absolute flux floor | one run, contiguous day/night blocks | +| `radiative_identity` | upward longwave equals what the reported surface temperature emits plus the reflected downward longwave, at every step, within the larger of a relative tolerance and an absolute floor; emissivity comes from `static.json` | one run, instantaneous values | | `routing_conservation` | the channel store is non-negative and never exceeds `max_lag_days` of the largest recent runoff | one run | Picking a denominator for `closure` and `regime_transfer`: @@ -267,6 +268,9 @@ The reference models available today: | `reference_energy_leak` | discards 15% of net radiation | `energy_closure` | | `reference_ground_dodge` | coherent and closed, but its sensible flux never reads the soil, so the ground flux absorbs a drydown's whole shift | `partition_shift` | | `reference_diurnal_bias` | shifts sensible heat to leave opposite day and night energy residuals that cancel over the full record | `energy_closure_by_phase` | +| `reference_radiative` | the coupled reference with a skin: temperature from the sensible flux through a fixed conductance, upward longwave from that skin at the given emissivity | nothing, it must pass the radiation probe | +| `reference_air_emitter` | reports the skin's temperature but emits at the air's, reflected sky unchanged | `radiative_identity` | +| `reference_no_reflection` | reports emission alone as the total upward longwave, the reflected sky left out | `radiative_identity` | If your probe needs a broken model that does not exist yet, add it under `models/` alongside the probe. A criterion with nothing that trips it is diff --git a/models/_template/model.yaml b/models/_template/model.yaml index 94618969..4c3c04b1 100644 --- a/models/_template/model.yaml +++ b/models/_template/model.yaml @@ -18,8 +18,16 @@ timestep: PT1D emits: fluxes: [pr, evspsbl, mrro] states: [mrso, snw, canopy] + # Outputs no budget integrates, such as a surface temperature. + # diagnostics: [ts] +# Inputs the adapter cannot run without. needs_forcing: [pr, tas, pet] +# needs_static: [soil_capacity_mm] +# Optional inputs consumed whenever supplied; document the fallback if absent. +# Either needs_* or uses_* satisfies a probe's required input declaration. +# uses_forcing: [rlds] +# uses_static: [eps] supports: perturbation: true resources: {cpu: 2, memory_gb: 4, gpu: false} diff --git a/models/cwatm/README.md b/models/cwatm/README.md index a97cd758..a4bc48b8 100644 --- a/models/cwatm/README.md +++ b/models/cwatm/README.md @@ -52,7 +52,7 @@ as VIOLATION, and they are not alike. | Probe | Result | Mechanism | | --- | --- | --- | -| `energy/evaporative-partition`, `energy/latent-heat-et-consistency`, `energy/surface-energy-closure` | N/A (INCOMPLETE) | CWatM reports no latent, sensible or ground heat flux (and the last probe is hourly), so these probes cannot ask it anything | +| `energy/evaporative-partition`, `energy/latent-heat-et-consistency`, `energy/surface-energy-closure`, `energy/radiation-consistency` | N/A (INCOMPLETE) | CWatM reports no latent, sensible or ground heat flux, and no surface temperature or upward longwave (and the last two probes are hourly), so these probes cannot ask it anything | | `energy/pet-consistency` | VIOLATION: evaporation on the wettest fifth of soil days is 0.695 of demand on the worst seed (0.695–0.705; at least 0.7) | a land cover's transpiration, bare-soil and interception evaporation together cannot exceed `crop_correct × cropKC × ETRef`, and the fraction-weighted crop coefficient is 0.689; snow evaporation is added on top of that cap, which is why the ratio sits just above 0.689 (Sensitivity). **A packaging choice**: with `preferentialFlow = False` it passes at 0.708–0.712, and higher crop coefficients or more forest pass too | | `mass/resolution-invariance` | VIOLATION: runoff differs by 52.1 % of the rain between PT1H and PT1D (38.6–52.1 %), evaporation by 0.6 % | CWatM has no dt (below); its groundwater reservoir releases `recessionCoeff × storage` per step, so at PT1H it drains 24 times too fast | | `mass/response-nonnegativity` | VIOLATION: runoff 0.29 mm/day below the control on 2002-02-12, eight days after 120 mm was added (seed 1713476937; the other two never dip) | on wetter soil preferential flow takes a larger share of a later storm, and its interflow part replaces surface runoff; runoff concentration releases interflow through a slower kernel than surface runoff, so the next day carries less (below). **A packaging choice**: with `preferentialFlow = False` it passes (largest dip 0.048 mm/day, within tolerance), and the verdict also moves with land cover and with the slope that sets the lag | diff --git a/models/lisflood/README.md b/models/lisflood/README.md index 26c762a4..07aad535 100644 --- a/models/lisflood/README.md +++ b/models/lisflood/README.md @@ -394,7 +394,7 @@ of their own. ## Result -**FAIL (ERROR), 15 of 18 probes passed, 3 N/A (INCOMPLETE).** These are the +**FAIL (ERROR), 15 of 18 probes passed, 4 N/A (INCOMPLETE).** These are the rows of the full gate-seed run of `5.0.0-onecell.5`, made on the emulated host described under "Native re-run". The verdict is ERROR because two probes ran out of time on that host. Any ERROR among the scored probes makes the verdict @@ -409,9 +409,10 @@ budget for each case, and every wet event closes to 1e-13 mm. are to be replaced by a native x86-64 evaluation; "Native re-run" gives the commands and the row replacement. -- **N/A (INCOMPLETE), 3, not scored:** `energy/evaporative-partition`, - `energy/latent-heat-et-consistency` and `energy/surface-energy-closure`. - LISFLOOD reports no heat fluxes, so these probes cannot ask it anything. +- **N/A (INCOMPLETE), 4, not scored:** `energy/evaporative-partition`, + `energy/latent-heat-et-consistency`, `energy/surface-energy-closure` and + `energy/radiation-consistency`. LISFLOOD reports no heat fluxes and no + surface temperature, so these probes cannot ask it anything. They are neither a pass nor a fail, and do not decide the verdict. - **ERROR, 2:** `mass/precipitation-counterfactual` and `mass/human-abstraction`. The container exceeded the 60 s budget, because a @@ -426,7 +427,7 @@ commands and the row replacement. `closure` and `state_bounds` pass. - On a host fast enough for the budget, the first should PASS and the second be VIOLATION. The model's verdict would then be FAIL (VIOLATION), with 16 - of 18 probes passed and 3 N/A. + of 18 probes passed and 4 N/A. - **VIOLATION, 1:** `mass/resolution-invariance`. Rain that falls within an hour runs off, so `mrro` differs by 13.0% of `pr` between PT1H and PT1D, against a 10% limit. @@ -451,11 +452,12 @@ Against the `.2` rows: - `mass/precipitation-counterfactual` is ERROR on this host in both. - `mass/human-abstraction` is new since `.2`. - The three energy-flux probes were FAIL (INCOMPLETE) under the earlier - roll-up and are N/A (INCOMPLETE) under main's. + roll-up and are N/A (INCOMPLETE) under main's; `energy/radiation-consistency`, + new since `.2`, is N/A too. - No other probe's verdict moved. The standing counts passes out of the 18 probes that could score LISFLOOD; the -three N/A energy-flux probes are in neither number. +four N/A energy probes are in neither number. Against the `.3` rows, `.4` changes the environmental-flow reserve and the channel's bottom width, bankfull depth and gradient to the headwater values. diff --git a/models/reference_air_emitter/ht_adapter.py b/models/reference_air_emitter/ht_adapter.py new file mode 100644 index 00000000..af82faec --- /dev/null +++ b/models/reference_air_emitter/ht_adapter.py @@ -0,0 +1,193 @@ +#!/usr/bin/env python3 +"""Negative control: report skin temperature but emit longwave at air temperature.""" + +from __future__ import annotations + +import argparse +import csv +import json +import sys +from pathlib import Path + +MODE = "air_emitter" +MODEL = {"name": "reference_air_emitter", "version": "1.0.0"} + +COLUMNS = [ + "time", "pr", "evspsbl", "mrro", "sbl", "hfls", "hfss", "hfg", "rlus", + "mrso", "snw", "canopy", "channel", "ts", +] + +EVAP_SHAPE = 0.5 # soil moisture at which evaporation reaches its potential rate +SUBL_SHARE = 0.35 # share of the remaining demand a snowpack can meet +GROUND_SHARE = 0.10 # share of net radiation conducted into the ground + +# Latent heat of vaporisation as A + B*T (T in degC) and of fusion, J kg-1. +LAMBDA_A, LAMBDA_B = 2.501e6, -2361.0 +LAMBDA_F = 3.337e5 +SECONDS_PER_DAY = 86400.0 + +# Stefan-Boltzmann constant (CODATA 2018) and bulk heat-transfer coefficient. +STEFAN_BOLTZMANN = 5.670374419e-8 +CONDUCTANCE = 20.0 # W m-2 K-1 +KELVIN = 273.15 + +TIMESTEP_DAYS = { + "PT1D": 1.0, + "PT1H": 1.0 / 24.0, + "PT15M": 1.0 / 96.0, + "PT5M": 1.0 / 288.0, + "PT1M": 1.0 / 1440.0, +} + + +def partition_energy(rn, tas, liquid_mm, sublimated_mm, dt_days): + """Keep reference_coupled's phase-specific latent heat and closed surface budget.""" + seconds = dt_days * SECONDS_PER_DAY + lam_v = LAMBDA_A + LAMBDA_B * tas + lam_s = LAMBDA_A + LAMBDA_F + ground = GROUND_SHARE * rn + latent = (lam_v * liquid_mm + lam_s * sublimated_mm) / seconds + return latent, rn - ground - latent, ground + + +def radiate(tas, sensible, rlds, eps): + """Diagnose Ts from H = g_H * (Ts - Ta), then total upward longwave. + + The interval-mean sensible flux is treated as holding at the sampling + instant. The sign of H sets whether skin is warmer or cooler than air. + Net radiation remains prescribed independently of longwave. + MODE changes only the emitting temperature or the reflected-sky term.""" + air_k = tas + KELVIN + skin_k = air_k + sensible / CONDUCTANCE + emitter_k = air_k if MODE == "air_emitter" else skin_k + upward = eps * STEFAN_BOLTZMANN * emitter_k ** 4 + if MODE != "no_reflection": + upward += (1.0 - eps) * rlds + return skin_k, upward + + +def simulate(forcing, static, dt_days=1.0): + """Run reference_coupled's water and energy core, adding Ts and rlus. + + Snow sublimation remains part of evaporation and closes the water budget.""" + soil_cap = static["soil_capacity_mm"] + canopy_cap = static["canopy_capacity_mm"] + ddf = static["degree_day_factor_mm_per_C_day"] + k_base = static["baseflow_coefficient"] + t_snow = static["snow_threshold_degC"] + # Use the supplied emissivity; otherwise assume a black surface. + eps = static.get("eps", 1.0) + + soil = 0.5 * soil_cap + swe = 0.0 + canopy = 0.0 + rows = [] + + for step in forcing: + pr_rate, tas, pet_rate = step["pr"], step["tas"], step["pet"] + rn = step.get("rn", 0.0) + rlds = step.get("rlds", 0.0) + pr = pr_rate * dt_days + pet = pet_rate * dt_days + + snowfall = pr if tas < t_snow else 0.0 + rain = 0.0 if tas < t_snow else pr + + swe += snowfall + melt = min(swe, ddf * max(tas - t_snow, 0.0) * dt_days) + swe -= melt + + water_in = rain + melt + intercepted = min(canopy_cap - canopy, water_in) + canopy += intercepted + throughfall = water_in - intercepted + + canopy_evap = min(canopy, pet) + canopy -= canopy_evap + pet_left = pet - canopy_evap + + sublimation = min(swe, SUBL_SHARE * pet_left) if swe > 0.0 else 0.0 + swe -= sublimation + pet_left -= sublimation + + soil += throughfall + surface = max(0.0, soil - soil_cap) + soil -= surface + baseflow = k_base * soil * dt_days + soil -= baseflow + soil_evap = min(soil, pet_left * min(1.0, soil / (EVAP_SHAPE * soil_cap))) + soil -= soil_evap + + liquid = canopy_evap + soil_evap + latent, sensible, ground = partition_energy(rn, tas, liquid, sublimation, dt_days) + skin_k, upward = radiate(tas, sensible, rlds, eps) + + rows.append({ + "time": step["time"], + "pr": pr_rate, + "evspsbl": (liquid + sublimation) / dt_days, + # The sublimating share of the evaporation above, so a + # criterion never has to infer which kilograms left as ice. + "sbl": sublimation / dt_days, + "mrro": (surface + baseflow) / dt_days, + "hfls": latent, + "hfss": sensible, + "hfg": ground, + "rlus": upward, + "mrso": soil, + "snw": swe, + "canopy": canopy, + # Runoff leaves in the same step; no water remains in routing. + "channel": 0.0, + "ts": skin_k, + }) + return rows + + +def read_request(path: Path) -> tuple[dict, Path]: + request = json.loads(path.read_text()) + return request, path.parent + + +def read_forcing(path: Path) -> list[dict]: + with open(path, newline="") as fh: + rows = list(csv.DictReader(fh)) + for row in rows: + for key in ("pr", "tas", "pet", "rn", "rlds"): + if key in row: + row[key] = float(row[key]) + return rows + + +def write_result(path: Path, rows: list[dict]) -> None: + path.parent.mkdir(parents=True, exist_ok=True) + with open(path, "w", newline="") as fh: + writer = csv.DictWriter(fh, fieldnames=COLUMNS) + writer.writeheader() + writer.writerows(rows) + + +def main() -> int: + parser = argparse.ArgumentParser() + parser.add_argument("--request", required=True) + args = parser.parse_args() + + request_path = Path(args.request).resolve() + request, io_dir = read_request(request_path) + forcing = read_forcing(io_dir / request["input"]["forcing"]) + static = json.loads((io_dir / request["input"]["static"]).read_text()) + + timestep = request.get("timestep", "PT1D") + if timestep not in TIMESTEP_DAYS: + raise SystemExit(f"unsupported timestep {timestep!r}") + rows = simulate(forcing, static, TIMESTEP_DAYS[timestep]) + + write_result(io_dir / request["output"]["table"], rows) + (io_dir / request["output"]["run"]).write_text( + json.dumps({"status": "ok", "model": MODEL, "n_steps": len(rows)}, indent=2) + ) + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/models/reference_air_emitter/model.yaml b/models/reference_air_emitter/model.yaml new file mode 100644 index 00000000..053cf797 --- /dev/null +++ b/models/reference_air_emitter/model.yaml @@ -0,0 +1,21 @@ +name: reference_air_emitter +version: "1.0.0" +description: > + Negative control for energy/radiation-consistency. Uses reference_radiative's + skin temperature but emits longwave at air temperature, with reflected sky + unchanged. Only upward longwave differs from the positive control. +authors: ["Xin Lan"] +license: PolyForm-Noncommercial-1.0.0 +runner: subprocess +entrypoint: ["python3", "ht_adapter.py"] +timestep: [PT1H] +emits: + fluxes: [pr, evspsbl, mrro, sbl, hfls, hfss, hfg, rlus] + states: [mrso, snw, canopy, channel] + diagnostics: [ts] +needs_forcing: [pr, tas, pet, rn] +uses_forcing: [rlds] +uses_static: [eps] +supports: + perturbation: true +resources: {cpu: 1, memory_gb: 1, gpu: false} diff --git a/models/reference_no_reflection/ht_adapter.py b/models/reference_no_reflection/ht_adapter.py new file mode 100644 index 00000000..adaff718 --- /dev/null +++ b/models/reference_no_reflection/ht_adapter.py @@ -0,0 +1,193 @@ +#!/usr/bin/env python3 +"""Negative control: omit reflected sky from total upward longwave.""" + +from __future__ import annotations + +import argparse +import csv +import json +import sys +from pathlib import Path + +MODE = "no_reflection" +MODEL = {"name": "reference_no_reflection", "version": "1.0.0"} + +COLUMNS = [ + "time", "pr", "evspsbl", "mrro", "sbl", "hfls", "hfss", "hfg", "rlus", + "mrso", "snw", "canopy", "channel", "ts", +] + +EVAP_SHAPE = 0.5 # soil moisture at which evaporation reaches its potential rate +SUBL_SHARE = 0.35 # share of the remaining demand a snowpack can meet +GROUND_SHARE = 0.10 # share of net radiation conducted into the ground + +# Latent heat of vaporisation as A + B*T (T in degC) and of fusion, J kg-1. +LAMBDA_A, LAMBDA_B = 2.501e6, -2361.0 +LAMBDA_F = 3.337e5 +SECONDS_PER_DAY = 86400.0 + +# Stefan-Boltzmann constant (CODATA 2018) and bulk heat-transfer coefficient. +STEFAN_BOLTZMANN = 5.670374419e-8 +CONDUCTANCE = 20.0 # W m-2 K-1 +KELVIN = 273.15 + +TIMESTEP_DAYS = { + "PT1D": 1.0, + "PT1H": 1.0 / 24.0, + "PT15M": 1.0 / 96.0, + "PT5M": 1.0 / 288.0, + "PT1M": 1.0 / 1440.0, +} + + +def partition_energy(rn, tas, liquid_mm, sublimated_mm, dt_days): + """Keep reference_coupled's phase-specific latent heat and closed surface budget.""" + seconds = dt_days * SECONDS_PER_DAY + lam_v = LAMBDA_A + LAMBDA_B * tas + lam_s = LAMBDA_A + LAMBDA_F + ground = GROUND_SHARE * rn + latent = (lam_v * liquid_mm + lam_s * sublimated_mm) / seconds + return latent, rn - ground - latent, ground + + +def radiate(tas, sensible, rlds, eps): + """Diagnose Ts from H = g_H * (Ts - Ta), then total upward longwave. + + The interval-mean sensible flux is treated as holding at the sampling + instant. The sign of H sets whether skin is warmer or cooler than air. + Net radiation remains prescribed independently of longwave. + MODE changes only the emitting temperature or the reflected-sky term.""" + air_k = tas + KELVIN + skin_k = air_k + sensible / CONDUCTANCE + emitter_k = air_k if MODE == "air_emitter" else skin_k + upward = eps * STEFAN_BOLTZMANN * emitter_k ** 4 + if MODE != "no_reflection": + upward += (1.0 - eps) * rlds + return skin_k, upward + + +def simulate(forcing, static, dt_days=1.0): + """Run reference_coupled's water and energy core, adding Ts and rlus. + + Snow sublimation remains part of evaporation and closes the water budget.""" + soil_cap = static["soil_capacity_mm"] + canopy_cap = static["canopy_capacity_mm"] + ddf = static["degree_day_factor_mm_per_C_day"] + k_base = static["baseflow_coefficient"] + t_snow = static["snow_threshold_degC"] + # Use the supplied emissivity; otherwise assume a black surface. + eps = static.get("eps", 1.0) + + soil = 0.5 * soil_cap + swe = 0.0 + canopy = 0.0 + rows = [] + + for step in forcing: + pr_rate, tas, pet_rate = step["pr"], step["tas"], step["pet"] + rn = step.get("rn", 0.0) + rlds = step.get("rlds", 0.0) + pr = pr_rate * dt_days + pet = pet_rate * dt_days + + snowfall = pr if tas < t_snow else 0.0 + rain = 0.0 if tas < t_snow else pr + + swe += snowfall + melt = min(swe, ddf * max(tas - t_snow, 0.0) * dt_days) + swe -= melt + + water_in = rain + melt + intercepted = min(canopy_cap - canopy, water_in) + canopy += intercepted + throughfall = water_in - intercepted + + canopy_evap = min(canopy, pet) + canopy -= canopy_evap + pet_left = pet - canopy_evap + + sublimation = min(swe, SUBL_SHARE * pet_left) if swe > 0.0 else 0.0 + swe -= sublimation + pet_left -= sublimation + + soil += throughfall + surface = max(0.0, soil - soil_cap) + soil -= surface + baseflow = k_base * soil * dt_days + soil -= baseflow + soil_evap = min(soil, pet_left * min(1.0, soil / (EVAP_SHAPE * soil_cap))) + soil -= soil_evap + + liquid = canopy_evap + soil_evap + latent, sensible, ground = partition_energy(rn, tas, liquid, sublimation, dt_days) + skin_k, upward = radiate(tas, sensible, rlds, eps) + + rows.append({ + "time": step["time"], + "pr": pr_rate, + "evspsbl": (liquid + sublimation) / dt_days, + # The sublimating share of the evaporation above, so a + # criterion never has to infer which kilograms left as ice. + "sbl": sublimation / dt_days, + "mrro": (surface + baseflow) / dt_days, + "hfls": latent, + "hfss": sensible, + "hfg": ground, + "rlus": upward, + "mrso": soil, + "snw": swe, + "canopy": canopy, + # Runoff leaves in the same step; no water remains in routing. + "channel": 0.0, + "ts": skin_k, + }) + return rows + + +def read_request(path: Path) -> tuple[dict, Path]: + request = json.loads(path.read_text()) + return request, path.parent + + +def read_forcing(path: Path) -> list[dict]: + with open(path, newline="") as fh: + rows = list(csv.DictReader(fh)) + for row in rows: + for key in ("pr", "tas", "pet", "rn", "rlds"): + if key in row: + row[key] = float(row[key]) + return rows + + +def write_result(path: Path, rows: list[dict]) -> None: + path.parent.mkdir(parents=True, exist_ok=True) + with open(path, "w", newline="") as fh: + writer = csv.DictWriter(fh, fieldnames=COLUMNS) + writer.writeheader() + writer.writerows(rows) + + +def main() -> int: + parser = argparse.ArgumentParser() + parser.add_argument("--request", required=True) + args = parser.parse_args() + + request_path = Path(args.request).resolve() + request, io_dir = read_request(request_path) + forcing = read_forcing(io_dir / request["input"]["forcing"]) + static = json.loads((io_dir / request["input"]["static"]).read_text()) + + timestep = request.get("timestep", "PT1D") + if timestep not in TIMESTEP_DAYS: + raise SystemExit(f"unsupported timestep {timestep!r}") + rows = simulate(forcing, static, TIMESTEP_DAYS[timestep]) + + write_result(io_dir / request["output"]["table"], rows) + (io_dir / request["output"]["run"]).write_text( + json.dumps({"status": "ok", "model": MODEL, "n_steps": len(rows)}, indent=2) + ) + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/models/reference_no_reflection/model.yaml b/models/reference_no_reflection/model.yaml new file mode 100644 index 00000000..d15c3aac --- /dev/null +++ b/models/reference_no_reflection/model.yaml @@ -0,0 +1,21 @@ +name: reference_no_reflection +version: "1.0.0" +description: > + Negative control for energy/radiation-consistency. Reports surface emission + alone as total upward longwave, omitting (1 - eps) times the downward + longwave. Only upward longwave differs from reference_radiative. +authors: ["Xin Lan"] +license: PolyForm-Noncommercial-1.0.0 +runner: subprocess +entrypoint: ["python3", "ht_adapter.py"] +timestep: [PT1H] +emits: + fluxes: [pr, evspsbl, mrro, sbl, hfls, hfss, hfg, rlus] + states: [mrso, snw, canopy, channel] + diagnostics: [ts] +needs_forcing: [pr, tas, pet, rn] +uses_forcing: [rlds] +uses_static: [eps] +supports: + perturbation: true +resources: {cpu: 1, memory_gb: 1, gpu: false} diff --git a/models/reference_radiative/ht_adapter.py b/models/reference_radiative/ht_adapter.py new file mode 100644 index 00000000..696102d4 --- /dev/null +++ b/models/reference_radiative/ht_adapter.py @@ -0,0 +1,193 @@ +#!/usr/bin/env python3 +"""Standard-library adapter with a skin satisfying the radiation identity.""" + +from __future__ import annotations + +import argparse +import csv +import json +import sys +from pathlib import Path + +MODE = "radiative" +MODEL = {"name": "reference_radiative", "version": "1.0.0"} + +COLUMNS = [ + "time", "pr", "evspsbl", "mrro", "sbl", "hfls", "hfss", "hfg", "rlus", + "mrso", "snw", "canopy", "channel", "ts", +] + +EVAP_SHAPE = 0.5 # soil moisture at which evaporation reaches its potential rate +SUBL_SHARE = 0.35 # share of the remaining demand a snowpack can meet +GROUND_SHARE = 0.10 # share of net radiation conducted into the ground + +# Latent heat of vaporisation as A + B*T (T in degC) and of fusion, J kg-1. +LAMBDA_A, LAMBDA_B = 2.501e6, -2361.0 +LAMBDA_F = 3.337e5 +SECONDS_PER_DAY = 86400.0 + +# Stefan-Boltzmann constant (CODATA 2018) and bulk heat-transfer coefficient. +STEFAN_BOLTZMANN = 5.670374419e-8 +CONDUCTANCE = 20.0 # W m-2 K-1 +KELVIN = 273.15 + +TIMESTEP_DAYS = { + "PT1D": 1.0, + "PT1H": 1.0 / 24.0, + "PT15M": 1.0 / 96.0, + "PT5M": 1.0 / 288.0, + "PT1M": 1.0 / 1440.0, +} + + +def partition_energy(rn, tas, liquid_mm, sublimated_mm, dt_days): + """Keep reference_coupled's phase-specific latent heat and closed surface budget.""" + seconds = dt_days * SECONDS_PER_DAY + lam_v = LAMBDA_A + LAMBDA_B * tas + lam_s = LAMBDA_A + LAMBDA_F + ground = GROUND_SHARE * rn + latent = (lam_v * liquid_mm + lam_s * sublimated_mm) / seconds + return latent, rn - ground - latent, ground + + +def radiate(tas, sensible, rlds, eps): + """Diagnose Ts from H = g_H * (Ts - Ta), then total upward longwave. + + The interval-mean sensible flux is treated as holding at the sampling + instant. The sign of H sets whether skin is warmer or cooler than air. + Net radiation remains prescribed independently of longwave. + MODE changes only the emitting temperature or the reflected-sky term.""" + air_k = tas + KELVIN + skin_k = air_k + sensible / CONDUCTANCE + emitter_k = air_k if MODE == "air_emitter" else skin_k + upward = eps * STEFAN_BOLTZMANN * emitter_k ** 4 + if MODE != "no_reflection": + upward += (1.0 - eps) * rlds + return skin_k, upward + + +def simulate(forcing, static, dt_days=1.0): + """Run reference_coupled's water and energy core, adding Ts and rlus. + + Snow sublimation remains part of evaporation and closes the water budget.""" + soil_cap = static["soil_capacity_mm"] + canopy_cap = static["canopy_capacity_mm"] + ddf = static["degree_day_factor_mm_per_C_day"] + k_base = static["baseflow_coefficient"] + t_snow = static["snow_threshold_degC"] + # Use the supplied emissivity; otherwise assume a black surface. + eps = static.get("eps", 1.0) + + soil = 0.5 * soil_cap + swe = 0.0 + canopy = 0.0 + rows = [] + + for step in forcing: + pr_rate, tas, pet_rate = step["pr"], step["tas"], step["pet"] + rn = step.get("rn", 0.0) + rlds = step.get("rlds", 0.0) + pr = pr_rate * dt_days + pet = pet_rate * dt_days + + snowfall = pr if tas < t_snow else 0.0 + rain = 0.0 if tas < t_snow else pr + + swe += snowfall + melt = min(swe, ddf * max(tas - t_snow, 0.0) * dt_days) + swe -= melt + + water_in = rain + melt + intercepted = min(canopy_cap - canopy, water_in) + canopy += intercepted + throughfall = water_in - intercepted + + canopy_evap = min(canopy, pet) + canopy -= canopy_evap + pet_left = pet - canopy_evap + + sublimation = min(swe, SUBL_SHARE * pet_left) if swe > 0.0 else 0.0 + swe -= sublimation + pet_left -= sublimation + + soil += throughfall + surface = max(0.0, soil - soil_cap) + soil -= surface + baseflow = k_base * soil * dt_days + soil -= baseflow + soil_evap = min(soil, pet_left * min(1.0, soil / (EVAP_SHAPE * soil_cap))) + soil -= soil_evap + + liquid = canopy_evap + soil_evap + latent, sensible, ground = partition_energy(rn, tas, liquid, sublimation, dt_days) + skin_k, upward = radiate(tas, sensible, rlds, eps) + + rows.append({ + "time": step["time"], + "pr": pr_rate, + "evspsbl": (liquid + sublimation) / dt_days, + # The sublimating share of the evaporation above, so a + # criterion never has to infer which kilograms left as ice. + "sbl": sublimation / dt_days, + "mrro": (surface + baseflow) / dt_days, + "hfls": latent, + "hfss": sensible, + "hfg": ground, + "rlus": upward, + "mrso": soil, + "snw": swe, + "canopy": canopy, + # Runoff leaves in the same step; no water remains in routing. + "channel": 0.0, + "ts": skin_k, + }) + return rows + + +def read_request(path: Path) -> tuple[dict, Path]: + request = json.loads(path.read_text()) + return request, path.parent + + +def read_forcing(path: Path) -> list[dict]: + with open(path, newline="") as fh: + rows = list(csv.DictReader(fh)) + for row in rows: + for key in ("pr", "tas", "pet", "rn", "rlds"): + if key in row: + row[key] = float(row[key]) + return rows + + +def write_result(path: Path, rows: list[dict]) -> None: + path.parent.mkdir(parents=True, exist_ok=True) + with open(path, "w", newline="") as fh: + writer = csv.DictWriter(fh, fieldnames=COLUMNS) + writer.writeheader() + writer.writerows(rows) + + +def main() -> int: + parser = argparse.ArgumentParser() + parser.add_argument("--request", required=True) + args = parser.parse_args() + + request_path = Path(args.request).resolve() + request, io_dir = read_request(request_path) + forcing = read_forcing(io_dir / request["input"]["forcing"]) + static = json.loads((io_dir / request["input"]["static"]).read_text()) + + timestep = request.get("timestep", "PT1D") + if timestep not in TIMESTEP_DAYS: + raise SystemExit(f"unsupported timestep {timestep!r}") + rows = simulate(forcing, static, TIMESTEP_DAYS[timestep]) + + write_result(io_dir / request["output"]["table"], rows) + (io_dir / request["output"]["run"]).write_text( + json.dumps({"status": "ok", "model": MODEL, "n_steps": len(rows)}, indent=2) + ) + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/models/reference_radiative/model.yaml b/models/reference_radiative/model.yaml new file mode 100644 index 00000000..b86ed5f0 --- /dev/null +++ b/models/reference_radiative/model.yaml @@ -0,0 +1,23 @@ +name: reference_radiative +version: "1.0.0" +description: > + Positive control for energy/radiation-consistency. Retains the water and + energy core of reference_coupled, diagnoses skin temperature from sensible + heat through a fixed conductance, and reports emission plus reflected sky + at the supplied emissivity. The radiation identity holds by construction; + the diagnosed temperature is an approximation. +authors: ["Xin Lan"] +license: PolyForm-Noncommercial-1.0.0 +runner: subprocess +entrypoint: ["python3", "ht_adapter.py"] +timestep: [PT1H] +emits: + fluxes: [pr, evspsbl, mrro, sbl, hfls, hfss, hfg, rlus] + states: [mrso, snw, canopy, channel] + diagnostics: [ts] +needs_forcing: [pr, tas, pet, rn] +uses_forcing: [rlds] +uses_static: [eps] +supports: + perturbation: true +resources: {cpu: 1, memory_gb: 1, gpu: false} diff --git a/models/result.csv b/models/result.csv index e1a223fc..f9183247 100644 --- a/models/result.csv +++ b/models/result.csv @@ -299,3 +299,12 @@ run_date,model,version,suite_version,runner,probe,verdict,reason,window,seeds,de 2026-09-12,summa,4.0.0-f787fa5.4,0.1.0,docker,mass/extreme-event-closure,FAIL,VIOLATION,full record,492502526 999454640 1506406754 2013358868 372827335,"state_bounds: canopy leaves [0, 2] on 68 steps (range 0.0092 to 25)" 2026-09-12,cwatm,1.11-5baaadd.4,0.1.0,docker,mass/extreme-event-closure,FAIL,VIOLATION,full record,492502526 999454640 1506406754 2013358868 372827335,"event_water_closure: worst event residual 0.00491841 mm / 0.001 mm allowed (114.5948% of event precipitation; 4/1257 events fail; allowance=max(5.0% of rain, 0.001 mm))" 2026-09-12,lisflood,5.0.0-onecell.5,0.1.0,docker,mass/extreme-event-closure,PASS,OK,7300-day flood event,492502526 999454640 1506406754 2013358868 372827335,"event_water_closure: worst event residual 1.35003e-13 mm / 0.001 mm allowed (0.0000% of event precipitation; 0/1299 events fail; allowance=max(5.0% of rain, 0.001 mm))" +2026-09-12,cwatm,1.11-5baaadd.4,0.1.0,docker,energy/radiation-consistency,N/A,INCOMPLETE,not run,,"does not report rlus, ts; model timestep PT1D does not cover the probe's PT1H; model does not declare that it consumes forcing rlds; model does not declare that it consumes static eps" +2026-09-12,dhbv2,0.5.4-hbv2ep100.3,0.1.0,docker,energy/radiation-consistency,N/A,INCOMPLETE,not run,,"does not report rlus, ts; model timestep PT1D does not cover the probe's PT1H; model does not declare that it consumes forcing rlds; model does not declare that it consumes static eps" +2026-09-12,flex_lumped,1.0.0,0.1.0,subprocess,energy/radiation-consistency,N/A,INCOMPLETE,not run,,"does not report rlus, ts; model does not declare that it consumes forcing rlds; model does not declare that it consumes static eps" +2026-09-12,flex_topo,1.0.0,0.1.0,subprocess,energy/radiation-consistency,N/A,INCOMPLETE,not run,,"does not report rlus, ts; model does not declare that it consumes forcing rlds; model does not declare that it consumes static eps" +2026-09-12,google_flood_forecast,0.1.0-828dfc5,0.1.0,docker,energy/radiation-consistency,N/A,INCOMPLETE,not run,,"does not report rlus, ts; model timestep PT1D does not cover the probe's PT1H; model does not declare that it consumes forcing rlds; model does not declare that it consumes static eps" +2026-09-12,sacsma_snow17,1.0.0,0.1.0,subprocess,energy/radiation-consistency,N/A,INCOMPLETE,not run,,"does not report rlus, ts; model does not declare that it consumes forcing rlds; model does not declare that it consumes static eps" +2026-09-12,summa,4.0.0-f787fa5.4,0.1.0,docker,energy/radiation-consistency,N/A,INCOMPLETE,not run,,"does not report rlus, ts; model does not declare that it consumes forcing rlds; model does not declare that it consumes static eps" +2026-09-12,wflow_sbm,1.0.4-ht.4,0.1.0,docker,energy/radiation-consistency,N/A,INCOMPLETE,not run,,"does not report rlus, ts; model does not declare that it consumes forcing rlds; model does not declare that it consumes static eps" +2026-09-12,lisflood,5.0.0-onecell.5,0.1.0,docker,energy/radiation-consistency,N/A,INCOMPLETE,not run,,"does not report rlus, ts; model does not declare that it consumes forcing rlds; model does not declare that it consumes static eps" diff --git a/models/summa/README.md b/models/summa/README.md index 3dc3c953..cb73d7d7 100644 --- a/models/summa/README.md +++ b/models/summa/README.md @@ -12,7 +12,7 @@ vegetation canopy, a layered snowpack, a layered soil column and an aquifer with one implicit solver, and it reports its latent, sensible and ground heat fluxes. It is the first submission that can be asked about both budgets and the identity between them, rather than being N/A (INCOMPLETE) on the -energy probes. +energy probes that need the heat fluxes. v4.0.1 was released on the day this was packaged. It adds brackets around `iden_ice/iden_water` in four files (a last-bit change in floating point) and @@ -261,7 +261,7 @@ SUMMA wants shortwave and longwave down, wind, pressure and specific humidity. They are mocked from each forcing row alone: nothing reads the calendar, the clock or another row, and `run.json` labels every one. -**Net radiation of the row.** Where the probe supplies `rn` (the three energy +**Net radiation of the row.** Where the probe supplies `rn` (the four energy probes that need it), that is used. Otherwise the row's `pet` is multiplied by the Priestley-Taylor (alpha 1.26) conversion from potential evaporation to net radiation. The conversion is evaluated at one fixed reference temperature, @@ -522,6 +522,8 @@ It is out of 21 since `mass/extreme-event-closure` merged. SUMMA fails that probe on the canopy's `state_bounds` alone, holding up to 25 mm of ice against a 2 mm capacity on 68 steps, while every wet event's water budget closes to within 3e-4 of its allowance. +`energy/radiation-consistency`, merged since, is N/A (INCOMPLETE): the adapter +reports no `rlus` or `ts`, so the count stays out of 21. In `4.0.0-f787fa5.3` the threshold translation flipped no probe verdict against `4.0.0-f787fa5.2`. It moved one failure inside a probe and changed the diff --git a/models/wflow_sbm/README.md b/models/wflow_sbm/README.md index 2be74a4a..63797234 100644 --- a/models/wflow_sbm/README.md +++ b/models/wflow_sbm/README.md @@ -270,6 +270,7 @@ case on the full record (`window_days: full`). | `energy/evaporative-partition` | N/A | INCOMPLETE | does not report `hfls`, `hfss`, `hfg` | | `energy/latent-heat-et-consistency` | N/A | INCOMPLETE | does not report `hfls`, `hfss`, `hfg` | | `energy/pet-consistency` | PASS | OK | evaporation 0.96 of demand when the soil is wettest, 0.15 when driest | +| `energy/radiation-consistency` | N/A | INCOMPLETE | does not report `rlus`, `ts` | | `energy/surface-energy-closure` | N/A | INCOMPLETE | does not report `hfls`, `hfss`, `hfg` | | `mass/antecedent-monotonicity` | FAIL | VIOLATION | a wet month before the storm adds almost no runoff (0.0013 and 0.0005 of the storm on 2 of 3 seeds, where 0.02 is asked) | | `mass/area-invariance` | PASS | OK | identical to floating point at ten times the area | @@ -289,8 +290,9 @@ case on the full record (`window_days: full`). | `mass/warming-response` | PASS | OK | runoff falls by 0.27 to 0.31 per unit of added demand | | `momentum/routing-conservation` | PASS | OK | the channel holds at most 0.17 of what a 15-day hydrograph of recent runoff allows | -The three energy-flux probes are N/A (INCOMPLETE) because wflow_sbm computes no latent, -sensible or ground heat flux; that is the model declining to be asked, not a failure. +The four energy probes that need an energy output are N/A (INCOMPLETE) because wflow_sbm +computes no latent, sensible or ground heat flux and no surface temperature; that is the +model declining to be asked, not a failure. `mass/human-abstraction` passes because the adapter now takes the prescribed withdrawal through Wflow's own water demand and allocation (see [A prescribed withdrawal](#a-prescribed-withdrawal)). diff --git a/probes/energy/radiation-consistency/README.md b/probes/energy/radiation-consistency/README.md new file mode 100644 index 00000000..8462c35c --- /dev/null +++ b/probes/energy/radiation-consistency/README.md @@ -0,0 +1,135 @@ +# energy/radiation-consistency + +## What it checks + +For a uniform, opaque, snow-free gray surface with fixed emissivity, total +upward longwave must equal surface emission plus reflected downward longwave: + +``` +r_i = rlus_i - [eps * sigma * ts_i^4 + (1 - eps) * rlds_i] +|r_i| <= max(0.005 * |rlus_i|, 0.5 W m-2) +``` + +`ts` is skin temperature in K; `rlus` and `rlds` are upward and downward +longwave in W m-2. `sigma` is the Stefan-Boltzmann constant. The criterion +reads `eps` from the supplied `static.json`, never from model output. + +Every scored step must pass. Opposite errors cannot cancel, and interval +means are out of scope: mean(T^4) is not mean(T)^4. The relative tolerance +uses the **reported** upward flux. `floor_steps` counts steps where the +absolute floor governs; none occur in the reference cases. + +The 0.5 percent tolerance and 0.5 W m-2 floor follow +[proposal #22](https://github.com/Flood-Lab/HydroTuring/issues/22). +They are consistency thresholds, not observational uncertainty. For example, +at 290 K, `eps = 0.98` and `rlds = 300 W m-2`, reflected longwave is 6 W m-2, +about 1.5 percent of total upward flux. Omitting it passes a 5 percent bound +but fails this one. + +## Case and controls + +The generator supplies 768 hourly rows starting at local-solar 06:00: +two spinup days and 30 scored days. `min_window_days: 30` preserves all +scored day/night cycles for submitted models, including dry hours outside +the default seven-day flood window. + +`tas`, `rlds`, `ts` and `rlus` are instantaneous values at each row's `time`. +`pr`, `pet` and `rn` are means over the following hour and drive the reference +water/energy budgets. +Net radiation is prescribed and is **not reconciled with the longwave +components**. + +Air temperature spans 16-30 degC with a daily weather offset and a 10 K +diurnal range. Hourly cloud transmission (0.45-1.0) both dims the sun and +raises daily clear-sky emissivity (0.74-0.82) towards one. Downward longwave +follows from that sky emissivity and air temperature; the scan below spans +297-456 W m-2. Surface emissivity is drawn once per seed in [0.95, 0.99], +rounded to four decimals, and held fixed. + +| Model | Construction | Expected result | +| --- | --- | --- | +| `reference_radiative` | `Ts = Ta + H / g_H`, with `g_H = 20 W m-2 K-1`; longwave uses this Ts, supplied eps and reflected sky | PASS | +| `reference_air_emitter` | Reports the same Ts but emits at air temperature; reflection unchanged | FAIL on `radiative_identity` | +| `reference_no_reflection` | Reports surface emission alone as total upward longwave | FAIL on `radiative_identity` | + +The positive control retains `reference_coupled`'s water and energy columns +exactly. It treats interval-mean sensible heat as holding at the sampling +instant: an explicit approximation, not a solved instantaneous energy +balance. Its radiation identity holds by construction; the diagnosed +temperature is approximate. Each negative control changes the radiation +mode; tests verify that only `rlus` differs from the positive control. + +All three references consume supplied `rlds` and `eps` through `uses_forcing` +and `uses_static`. When absent, they default to zero downward longwave and +a black surface (`eps = 1`), so they can still run on other compatible +probes. This probe always supplies both inputs; the defaults do not apply. + +## Limits + +A wrong Ts paired with longwave computed from that Ts can pass. This checks +output consistency, not temperature accuracy, energy balance or thermal +inertia. Adapters must expose native model outputs, not manufacture `rlus` +from `ts`. Missing either output yields INCOMPLETE, as for every archived +physical and submitted model on this probe. + +A model using its own land-cover emissivity cannot be scored here unless +its adapter consumes the case's `eps`, as well as its `rlds`. + +Around 290 K under a 300 W m-2 sky, temperature differences below about +0.4 K can pass; resolution is coarser at warmer temperatures. Unit tests +also document a linearised Stefan-Boltzmann law: around 288 K, excursions +of 5 K pass while 10 K fails, with a neglected second-order term near +2.8 W m-2. That optional control is not part of the gate. + +In a sweep of the five gate seeds plus seeds 0-199, a blackbody approximation +(`rlus = sigma * ts^4`) passes 11 of 205 cases: every case whose drawn `eps` +is 0.9881 or more, and none at 0.9864 or less, so it passes above about +`eps` 0.987. The gate seeds draw `eps` up to 0.9844, so it fails all five of +them today, but it passes all 205 cases with `eps` forced to 0.99. + +Emission linearized around the previous step's Ts passes 8 of those 205 +cases. Its maximum residual-to-tolerance ratio ranges from 0.88 to 1.80 +(1.25-1.59 on the gate seeds), depending on hourly skin-temperature jumps, +which reach 12.6 K in this sweep. It remains outside the gate. + +## Reproduce + +Run from the repository root after installing HydroTuring: + +```bash +ht validate +ht gate --probe energy/radiation-consistency +ht run --model reference_radiative --probe energy/radiation-consistency --seed 0 +ht run --model reference_no_reflection --probe energy/radiation-consistency --seed 0 +pytest -q tests/test_radiation_consistency.py tests/test_radiation_probe.py tests/test_radiation_references.py +python3 scripts/radiation_margins.py +``` + +The first model must pass; the second must fail on `radiative_identity`. +The gate makes 15 adapter invocations across five seeds. The scan adds +seeds 0-19, seed 36 (eps 0.9891, highest among seeds 0-49), and the five +gate cases forced to each emissivity endpoint: 36 cases in total. + +Ratios below are absolute residual divided by tolerance. The positive +column bounds the largest ratio in the group: the value is rounding, of +order 1e-13 or smaller, varying with the numpy and pandas installed. +Negative columns give the **minimum across cases of each case's maximum +ratio**; the design target is at least 1.3 for both negative controls, +with the criterion threshold unchanged at 1. + +| Group | Cases | Positive, bound | Air emitter | No reflection | +| --- | --- | --- | --- | --- | +| gate | 5 | < 1e-12 | 45.7 | 3.10 | +| additional | 20 | < 1e-12 | 44.5 | 3.12 | +| drawn eps 0.9891 | 1 | < 1e-12 | 49.9 | 2.15 | +| forced eps 0.99 | 5 | < 1e-12 | 46.1 | 1.95 | +| forced eps 0.95 | 5 | < 1e-12 | 44.7 | 10.2 | + +The air emitter violates 679-690 of 720 scored steps per case: every night +step and all daytime steps except some 06:00 and 17:00 rows near Ts = Ta. +No reflection violates all 720 steps. At eps 0.99 its ratio is +`2 * rlds / rlus`: per-case maxima are 1.95-1.98, while the minimum over +**all individual steps** is 1.20. Skin temperature stays within 287-320 K +and never freezes. The script also checks unchanged water/energy columns +and that negative controls differ only in `rlus`; an optional output path +writes the per-case results as JSON. diff --git a/probes/energy/radiation-consistency/generate.py b/probes/energy/radiation-consistency/generate.py new file mode 100644 index 00000000..dfd1274b --- /dev/null +++ b/probes/energy/radiation-consistency/generate.py @@ -0,0 +1,95 @@ +"""Seeded hourly forcing for a warm, snow-free gray surface. + +Timestamps use local solar time. tas and rlds are instantaneous, as are the +requested ts and rlus. Water/energy inputs pr, pet and rn are means over +[time, time + 1 h). Net radiation is prescribed independently of longwave.""" + +from __future__ import annotations + +import numpy as np +import pandas as pd + +PERIOD_DAYS = 30 +SPINUP_DAYS = 2 +STEPS_PER_DAY = 24 +N_STEPS = (PERIOD_DAYS + SPINUP_DAYS) * STEPS_PER_DAY + +# Stefan-Boltzmann constant, W m-2 K-4 (CODATA 2018). +STEFAN_BOLTZMANN = 5.670374419e-8 +# Shared bounds for shortwave transmission and derived cloud fraction. +TRANSMISSION = (0.45, 1.0) + +STATIC = { + "area_km2": 250.0, + "soil_capacity_mm": 180.0, + "canopy_capacity_mm": 0.0, + "degree_day_factor_mm_per_C_day": 3.2, + "baseflow_coefficient": 0.006, + "snow_threshold_degC": 0.0, + "latitude_deg": 0.0, +} + + +def generate(seed: int) -> tuple[pd.DataFrame, dict]: + rng = np.random.default_rng(seed) + time = pd.date_range("2000-03-20T06:00:00", periods=N_STEPS, freq="h") + hour = time.hour.to_numpy() + daytime = (hour >= 6) & (hour < 18) + solar_day = np.arange(N_STEPS) // STEPS_PER_DAY + n_days = PERIOD_DAYS + SPINUP_DAYS + + # Exact hourly mean of the positive sine between sunrise at 06:00 and + # sunset at 18:00, as in energy/surface-energy-closure. Cloud attenuation + # is drawn at every hour, night included, because it also sets how much + # the sky radiates downward. + omega = np.pi / 12.0 + daylight = np.where( + daytime, + (np.cos(omega * (hour - 6)) - np.cos(omega * (hour - 5))) / omega, + 0.0, + ) + peak_shortwave = rng.uniform(420.0, 540.0, n_days)[solar_day] + transmission = rng.uniform(*TRANSMISSION, N_STEPS) + longwave_cooling = rng.uniform(25.0, 40.0, n_days)[solar_day] + rn = peak_shortwave * transmission * daylight - longwave_cooling + + # Sample warm air at the radiation timestamp, not as an hourly mean. + weather = rng.uniform(-2.0, 2.0, n_days)[solar_day] + tas = 23.0 + weather + 5.0 * np.cos(2.0 * np.pi * (hour - 14.0) / 24.0) + + # The same cloud draw dims the sun and raises sky emissivity. Overcast + # conditions close three quarters of the gap from clear sky to a black body. + cloud = (TRANSMISSION[1] - transmission) / (TRANSMISSION[1] - TRANSMISSION[0]) + eps_clear = rng.uniform(0.74, 0.82, n_days)[solar_day] + eps_sky = eps_clear + (1.0 - eps_clear) * 0.75 * cloud + rlds = eps_sky * STEFAN_BOLTZMANN * (tas + 273.15) ** 4 + + # Water-side inputs for the reference model, in mm/day at every hour. + pet = 0.12 + 3.2 * daylight * transmission + wet = rng.random(N_STEPS) < 0.025 + pr = np.where(wet, rng.uniform(12.0, 48.0, N_STEPS), 0.0) + + # Supply one fixed surface emissivity for the entire case. + static = dict(STATIC) + static["eps"] = round(float(rng.uniform(0.95, 0.99)), 4) + + forcing = pd.DataFrame( + { + "time": time.strftime("%Y-%m-%dT%H:%M:%S"), + "pr": np.round(pr, 6), + "tas": np.round(tas, 6), + "pet": np.round(pet, 6), + "rn": np.round(rn, 6), + "rlds": np.round(rlds, 6), + } + ) + return forcing, static + + +if __name__ == "__main__": + frame, static = generate(20260911) + print(f"{len(frame)} hourly rows: {SPINUP_DAYS} spinup + {PERIOD_DAYS} scored days") + print(f"Air temperature: {frame['tas'].min():.1f} to {frame['tas'].max():.1f} degC") + print(f"Downward longwave: {frame['rlds'].min():.0f} to {frame['rlds'].max():.0f} W/m2") + print(f"Net radiation: {frame['rn'].min():.0f} to {frame['rn'].max():.0f} W/m2") + print(f"Surface emissivity: {static['eps']}") diff --git a/probes/energy/radiation-consistency/probe.yaml b/probes/energy/radiation-consistency/probe.yaml new file mode 100644 index 00000000..dcd26499 --- /dev/null +++ b/probes/energy/radiation-consistency/probe.yaml @@ -0,0 +1,76 @@ +id: energy/radiation-consistency +title: Surface temperature and upward longwave must describe one surface +law: energy +track: synthetic +version: 1 + +authors: + - name: Xin Lan + affiliation: >- + Center for Systems Integration and Sustainability, Department of + Fisheries and Wildlife, Michigan State University, East Lansing, MI + 48823, USA + orcid: "0000-0002-0607-2270" + github: xinlan-technology + +description: > + For a uniform, opaque, snow-free gray surface, total upward longwave must + equal emission at the reported skin temperature plus reflected downward + longwave, using emissivity supplied in static.json. Every instantaneous + hourly sample must agree within 0.5 percent of reported upward flux or + 0.5 W m-2, whichever is larger. This checks consistency, not temperature + accuracy. Net radiation is prescribed independently of these components. + +requires: + fluxes: [rlus] + states: [] + diagnostics: [ts] + # The verdict rests on the case's sky and emissivity. A model that does + # not declare it consumes them may be using its own, and is INCOMPATIBLE + # rather than judged against values it never read. + forcing: [rlds] + static: [eps] + +case: + generator: generate.py + n_seeds: 5 + timestep: PT1H + period_days: 30 + spinup_days: 2 + # Keep all scored day/night cycles, including dry hours outside the + # default seven-day window around the largest flood event. + min_window_days: 30 + max_output_mb: 2.0 + max_runtime_s: 60 + +criteria: + - radiative_identity: + upward: rlus + temperature: ts + downward: rlds + # Emissivity belongs to the case, not the model output. + emissivity: eps + # 0.5 percent of the model's own reported upward flux, not the suite's + # 5 percent: two outputs of one model are held to one equation. At + # eps 0.98, a 290 K surface under a 300 W m-2 sky reflects 6 W m-2, + # about 1.5 percent of its upward flux, so a model that drops the + # reflected sky would pass at 5 percent. The floor stops the bound + # shrinking towards zero at weak fluxes; at this warm surface the + # relative bound always governs. + rel_tol: 0.005 + abs_floor: 0.5 + +baselines: + must_pass: [reference_radiative] + must_fail: + reference_air_emitter: radiative_identity + reference_no_reflection: radiative_identity + +provenance: > + Synthetic, generated at run time from the recorded seed. Hourly cloud + attenuation sets both the shortwave and the sky emissivity behind the + downward longwave; daily shortwave amplitude, clear-sky emissivity, + longwave cooling and a weather offset vary by seed, as does the surface + emissivity drawn in [0.95, 0.99] and held fixed for the run. Timestamps + are local-solar instants. Proposal and agreed scope: + https://github.com/Flood-Lab/HydroTuring/issues/22. diff --git a/schemas/model.schema.json b/schemas/model.schema.json index b65bbc76..3a5eb53b 100644 --- a/schemas/model.schema.json +++ b/schemas/model.schema.json @@ -39,10 +39,29 @@ "states": { "type": "array", "items": {"type": "string"}, "uniqueItems": true, "description": "Absolute storage values, never tendencies. The harness differences them." + }, + "diagnostics": { + "type": "array", "items": {"type": "string"}, "uniqueItems": true, + "description": "Diagnostic outputs that are neither fluxes nor stores, such as a surface temperature. Read by criteria, never integrated or differenced." } } }, - "needs_forcing": {"type": "array", "items": {"type": "string"}}, + "needs_forcing": { + "type": "array", "items": {"type": "string"}, + "description": "Forcing columns the adapter cannot run without. A missing column makes the case INCOMPATIBLE." + }, + "needs_static": { + "type": "array", "items": {"type": "string"}, "uniqueItems": true, + "description": "static.json keys the adapter cannot run without. A missing key makes the case INCOMPATIBLE." + }, + "uses_forcing": { + "type": "array", "items": {"type": "string"}, "uniqueItems": true, + "description": "Optional forcing columns the adapter consumes whenever supplied. Their absence does not prevent running. These and needs_forcing satisfy a probe's requires.forcing." + }, + "uses_static": { + "type": "array", "items": {"type": "string"}, "uniqueItems": true, + "description": "Optional static.json keys the adapter consumes whenever supplied. Their absence does not prevent running. These and needs_static satisfy a probe's requires.static." + }, "supports": { "type": "object", "additionalProperties": false, diff --git a/schemas/probe.schema.json b/schemas/probe.schema.json index 8628b825..46298b74 100644 --- a/schemas/probe.schema.json +++ b/schemas/probe.schema.json @@ -94,6 +94,30 @@ "type": "string" }, "uniqueItems": true + }, + "diagnostics": { + "type": "array", + "items": { + "type": "string" + }, + "uniqueItems": true, + "description": "Diagnostic outputs a criterion reads but no budget integrates or differences, such as a surface temperature." + }, + "forcing": { + "type": "array", + "items": { + "type": "string" + }, + "uniqueItems": true, + "description": "Forcing columns the verdict rests on. A model must list them in needs_forcing or uses_forcing, or it is INCOMPATIBLE: it may be using its own estimate instead." + }, + "static": { + "type": "array", + "items": { + "type": "string" + }, + "uniqueItems": true, + "description": "static.json keys the verdict rests on. A model must list them in needs_static or uses_static, or it is INCOMPATIBLE." } } }, diff --git a/scripts/radiation_margins.py b/scripts/radiation_margins.py new file mode 100644 index 00000000..4e67fb5d --- /dev/null +++ b/scripts/radiation_margins.py @@ -0,0 +1,157 @@ +#!/usr/bin/env python3 +"""Reproduce the radiation probe's margin table through real adapters. + +Run five gate seeds, twenty additional seeds, seed 36 near the emissivity +ceiling, and gate seeds at both emissivity bounds. Check output parity and +report per-case maximum/minimum residual-to-tolerance ratios and day/night +violation counts. Each negative control's maximum ratio must reach 1.3; +the criterion's threshold remains 1. This extended sweep is separate from +the fast regression tests. + + python3 scripts/radiation_margins.py [output.json] +""" + +from __future__ import annotations + +import json +import sys +import tempfile +from dataclasses import replace +from pathlib import Path + +import numpy as np +import pandas as pd + +from hydroturing import registry +from hydroturing.criteria import get +from hydroturing.criteria.radiation import STEFAN_BOLTZMANN +from hydroturing.harness import build_case +from hydroturing.runner import get_runner +from hydroturing.seeds import gate_seeds + +PROBE = registry.find_probe("energy/radiation-consistency") +POSITIVE = "reference_radiative" +NEGATIVES = ("reference_air_emitter", "reference_no_reflection") +CRITERION = get("radiative_identity") +PARAMS = dict(PROBE.criteria[0].params) +KELVIN = 273.15 +MARGIN = 1.3 + + +def run(name, case): + # Read back what the model declares, as verify-adapter does, so that the + # coupled reference can be run on a probe that asks for more than it has. + model = registry.find_model(name) + asks = replace( + PROBE, requires_fluxes=model.emits_fluxes, requires_states=model.emits_states, + requires_diagnostics=model.emits_diagnostics, + ) + with tempfile.TemporaryDirectory(prefix="ht-margins-") as workdir: + return get_runner(model).run(model, asks, case, Path(workdir)) + + +def score(result): + r = CRITERION(result, PROBE, PARAMS) + case = result.case + forcing = case.forcing.iloc[case.spinup_steps:].reset_index(drop=True) + hours = pd.to_datetime(forcing["time"]).dt.hour + day = ((hours >= 6) & (hours < 18)).to_numpy() + # The criterion reports the worst step; the split needs every step, so + # the residual is recomputed here with the criterion's own constants. + scored = result.table.iloc[case.spinup_steps:].reset_index(drop=True) + eps = case.static[PARAMS["emissivity"]] + sigma = PARAMS.get("sigma", STEFAN_BOLTZMANN) + expected = eps * sigma * scored["ts"] ** 4 + (1.0 - eps) * forcing["rlds"] + slack = ((scored["rlus"] - expected).abs() + / np.maximum(PARAMS["rel_tol"] * scored["rlus"].abs(), PARAMS["abs_floor"])).to_numpy() + violating = slack > 1.0 + return { + "status": r.status, + "worst_slack": r.diagnostics.get("worst_slack", float("nan")), + # Distinguish the weakest single-step margin from the worst step. + "min_slack": float(slack.min()), + "violating": int(violating.sum()), + "violating_day": int((violating & day).sum()), + "violating_night": int((violating & ~day).sum()), + } + + +def sweep(label, seeds, eps=None): + rows = [] + for seed in seeds: + case = build_case(PROBE, seed) + if eps is not None: + case.static["eps"] = eps + tables = {} + row = {"label": label, "seed": seed, "eps": case.static["eps"]} + for name in (POSITIVE, *NEGATIVES): + result = run(name, case) + tables[name] = result.table + row[name] = score(result) + coupled = run("reference_coupled", case).table + + positive = tables[POSITIVE] + row["positive_matches_coupled"] = all( + np.array_equal(positive[c].to_numpy(), coupled[c].to_numpy()) for c in coupled.columns + ) + for name in NEGATIVES: + same = [ + c for c in positive.columns + if np.array_equal(positive[c].to_numpy(), tables[name][c].to_numpy()) + ] + row[f"{name}_differs_only_in_rlus"] = set(positive.columns) - set(same) == {"rlus"} + ts = positive["ts"].to_numpy() + row["ts_min_k"], row["ts_max_k"] = float(ts.min()), float(ts.max()) + row["ts_frozen_steps"] = int((ts <= KELVIN).sum()) + row["ok"] = ( + row["positive_matches_coupled"] + and all(row[f"{n}_differs_only_in_rlus"] for n in NEGATIVES) + and row["ts_frozen_steps"] == 0 + and row[POSITIVE]["status"] == "pass" + and all(row[n]["status"] == "fail" and row[n]["worst_slack"] >= MARGIN for n in NEGATIVES) + ) + rows.append(row) + return rows + + +def main() -> int: + gate = gate_seeds(PROBE.id, PROBE.n_seeds) + groups = [ + ("gate", gate, None), + ("additional", list(range(20)), None), + # Seed 36 draws the highest emissivity among seeds 0-49, 0.9891: a + # drawn case near the ceiling. Forced 0.99 also checks the endpoint. + ("hi-eps", [36], None), + ("eps=0.99", gate, 0.99), + ("eps=0.95", gate, 0.95), + ] + rows = [row for label, seeds, eps in groups for row in sweep(label, seeds, eps)] + if len(sys.argv) > 1: + Path(sys.argv[1]).write_text(json.dumps(rows, indent=1)) + + print(f"{'group':10} {'seed':>10} {'eps':>6} | {'positive':>9} | " + f"{'air_emitter':>11} {'viol':>4} {'day/night':>9} | " + f"{'no_reflection':>13} {'viol':>4} {'day/night':>9} |") + for r in rows: + p, a, n = r[POSITIVE], r[NEGATIVES[0]], r[NEGATIVES[1]] + print(f"{r['label']:10} {r['seed']:>10} {r['eps']:>6.4f} | {p['worst_slack']:>9.1e} | " + f"{a['worst_slack']:>11.2f} {a['violating']:>4} {a['violating_day']:>4}/{a['violating_night']:<4} | " + f"{n['worst_slack']:>13.2f} {n['violating']:>4} {n['violating_day']:>4}/{n['violating_night']:<4} | " + f"{'ok' if r['ok'] else 'PROBLEM'}") + + for label, _, _ in groups: + sub = [r for r in rows if r["label"] == label] + print(f"\n{label}: {len(sub)} seed(s); positive worst step " + f"{max(r[POSITIVE]['worst_slack'] for r in sub):.1e}; " + f"min over seeds of the worst step: air_emitter " + f"{min(r[NEGATIVES[0]]['worst_slack'] for r in sub):.2f}, no_reflection " + f"{min(r[NEGATIVES[1]]['worst_slack'] for r in sub):.2f}; " + f"smallest single step: air_emitter " + f"{min(r[NEGATIVES[0]]['min_slack'] for r in sub):.2f}, no_reflection " + f"{min(r[NEGATIVES[1]]['min_slack'] for r in sub):.2f}; " + f"skin {min(r['ts_min_k'] for r in sub):.1f} to {max(r['ts_max_k'] for r in sub):.1f} K") + return 0 if all(r["ok"] for r in rows) else 1 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/site/index.html b/site/index.html index b8755268..a3c83de9 100644 --- a/site/index.html +++ b/site/index.html @@ -318,8 +318,10 @@

How it works

partition day / night - - probes wanted + + radiation + + probes wanted MOMENTUM routing conserves @@ -635,6 +637,7 @@

Credit

"infra.probe.latent":"latent heat", "infra.probe.partition":"partition", "infra.probe.diurnal":"day / night", + "infra.probe.radiation":"radiation", "infra.probe.routing":"routing", "infra.flow":"seed → case → run → one bit", "infra.s1a":"seeded", @@ -714,6 +717,7 @@

Credit

energy/latent-heat-et-consistencyenergyThe evaporation a model reports as water and the one implied by its latent heat: are they the same evaporation?merged energy/evaporative-partitionenergyOne summer without rain under radiation that did not change: the latent heat a drying surface gives up has to warm the air.merged energy/surface-energy-closureenergyDoes each day and night close its hourly surface energy budget, without opposite errors cancelling?merged +energy/radiation-consistencyenergyDo the surface temperature and the upward longwave a model reports describe one surface, hour by hour, at the emissivity it was given?merged momentum/routing-conservationmomentumThe channel store is never negative and never holds more than its hydrograph can.merged mass/extreme-event-closuremassEvent water budgets under overlapping storms targeting synthetic 100-year rainfall depths.merged coupled/snowmelt-energy-watercoupledMelt in the water budget must match the energy spent melting.wanted @@ -722,7 +726,6 @@

Credit

mass/routing-network-closuremassA branching network: mass must close reach by reach.wanted mass/ungauged-basin-closuremassClosure on catchments outside the range models are fitted to.wanted energy/snowpack-cold-contentenergyThe full snowpack energy budget, cold content and phase change included.wanted -energy/radiation-consistencyenergyOutgoing longwave must be consistent with the reported surface temperature.wanted momentum/channel-routing-massmomentumInflow minus outflow minus the change in channel storage, per reach.wanted momentum/stage-discharge-monotonicmomentumA steady-flow rating must be monotonic.wanted momentum/wave-celerity-boundsmomentumKinematic wave celerity must be positive and near the Manning expectation.wanted @@ -839,6 +842,7 @@

Credit

"infra.probe.latent":"calor latente", "infra.probe.partition":"partición", "infra.probe.diurnal":"día / noche", + "infra.probe.radiation":"radiación", "infra.probe.routing":"tránsito", "infra.flow":"semilla → caso → ejecución → un bit", "infra.s1a":"generador", @@ -918,6 +922,7 @@

Credit

energy/latent-heat-et-consistencyenergíaLa evaporación que un modelo reporta como agua y la que implica su calor latente: ¿son la misma evaporación?integrada energy/evaporative-partitionenergíaUn verano sin lluvia bajo una radiación que no cambió: el calor latente que cede una superficie que se seca tiene que calentar el aire.integrada energy/surface-energy-closureenergía¿Cierra el balance energético horario en cada día y cada noche, sin que se cancelen errores opuestos?integrada +energy/radiation-consistencyenergía¿Describen la temperatura superficial y la onda larga ascendente que reporta un modelo una misma superficie, hora a hora, con la emisividad que se le dio?integrada momentum/routing-conservationmomentoEl almacén del cauce nunca es negativo ni retiene más de lo que cabe en su hidrograma.integrada mass/extreme-event-closuremasaBalances de agua por evento bajo tormentas superpuestas con objetivos sintéticos de lluvia de 100 años.integrada coupled/snowmelt-energy-wateracopladaEl deshielo del balance hídrico debe coincidir con la energía gastada en fundir.se busca @@ -926,7 +931,6 @@

Credit

mass/routing-network-closuremasaUna red ramificada: la masa debe cerrar tramo a tramo.se busca mass/ungauged-basin-closuremasaCierre en cuencas fuera del rango al que se ajustan los modelos.se busca energy/snowpack-cold-contentenergíaEl balance energético completo del manto de nieve, con contenido frío y cambio de fase.se busca -energy/radiation-consistencyenergíaLa onda larga saliente debe ser coherente con la temperatura superficial reportada.se busca momentum/channel-routing-massmomentoEntrada menos salida menos el cambio de almacenamiento en el cauce, por tramo.se busca momentum/stage-discharge-monotonicmomentoUna curva de gasto en flujo estacionario debe ser monótona.se busca momentum/wave-celerity-boundsmomentoLa celeridad de la onda cinemática debe ser positiva y cercana a la esperada por Manning.se busca @@ -1043,6 +1047,7 @@

Credit

"infra.probe.latent":"潜热一致", "infra.probe.partition":"能量分配", "infra.probe.diurnal":"昼夜闭合", + "infra.probe.radiation":"辐射一致", "infra.probe.routing":"汇流演算", "infra.flow":"种子 → 案例 → 运行 → 一位判定", "infra.s1a":"带种子的", @@ -1122,6 +1127,7 @@

Credit

energy/latent-heat-et-consistency能量模型作为水量报告的蒸发与其潜热通量所隐含的蒸发：是否是同一个蒸发？已合并 energy/evaporative-partition能量净辐射不变的情况下一个无雨的夏天：干燥地表让出的潜热必须去加热空气。已合并 energy/surface-energy-closure能量每个白天和夜晚的逐小时地表能量预算是否闭合，而没有依赖正负误差抵消？已合并 +energy/radiation-consistency能量模型报告的地表温度与上行长波辐射，在给定发射率下是否逐小时描述同一个表面？已合并 momentum/routing-conservation动量河道蓄水永不为负，也永不超过其单位线所能容纳的量。已合并 mass/extreme-event-closure质量叠加降雨以超过合成气候的百年一遇雨量阈值，逐事件检查水量闭合。已合入 coupled/snowmelt-energy-water耦合水量收支中的融雪量必须与融化耗能一致。征集中 @@ -1130,7 +1136,6 @@

Credit

mass/routing-network-closure质量分叉河网：质量须逐河段闭合。征集中 mass/ungauged-basin-closure质量在模型拟合范围之外的流域上保持闭合。征集中 energy/snowpack-cold-content能量完整的积雪能量收支，含冷储量与相变。征集中 -energy/radiation-consistency能量出射长波辐射必须与所报告的地表温度一致。征集中 momentum/channel-routing-mass动量逐河段：入流减出流减河道蓄量变化。征集中 momentum/stage-discharge-monotonic动量恒定流的水位流量关系必须单调。征集中 momentum/wave-celerity-bounds动量运动波波速必须为正，且接近曼宁公式的预期。征集中 diff --git a/src/hydroturing/criteria/__init__.py b/src/hydroturing/criteria/__init__.py index 503e7647..d30b6f92 100644 --- a/src/hydroturing/criteria/__init__.py +++ b/src/hydroturing/criteria/__init__.py @@ -20,6 +20,7 @@ event_water_closure, human, limits, + radiation, regime, response, stress, diff --git a/src/hydroturing/criteria/radiation.py b/src/hydroturing/criteria/radiation.py new file mode 100644 index 00000000..aabe5fb6 --- /dev/null +++ b/src/hydroturing/criteria/radiation.py @@ -0,0 +1,169 @@ +"""Check instantaneous surface temperature against total upward longwave. + +For a uniform, opaque, snow-free gray surface with a fixed emissivity, the +upward longwave radiation at an instant is fixed by the surface temperature +and the sky: + + rlus = eps * sigma * ts**4 + (1 - eps) * rlds + +Emission and reflected sky must agree with the reported flux at every step; +interval means and cancellation between steps are not allowed. The default +tolerance is max(0.005 * abs(rlus), 0.5 W m-2). Emissivity comes only from +static.json. This checks output consistency, not temperature accuracy. +""" + +from __future__ import annotations + +import numpy as np + +from hydroturing.criteria.base import FAIL, PASS, CriterionResult, criterion, make_window +from hydroturing.protocol import RunResult +from hydroturing.spec import ProbeSpec + +# Stefan-Boltzmann constant, W m-2 K-4 (CODATA 2018). +STEFAN_BOLTZMANN = 5.670374419e-8 + + +@criterion("radiative_identity") +def radiative_identity(run: RunResult, probe: ProbeSpec, params: dict) -> CriterionResult: + """Upward longwave must equal what the reported surface temperature emits + plus the reflected share of the downward longwave, at every step.""" + upward = str(params.get("upward", "rlus")) + temperature = str(params.get("temperature", "ts")) + downward = str(params.get("downward", "rlds")) + emissivity_key = str(params.get("emissivity", "eps")) + sigma = float(params.get("sigma", STEFAN_BOLTZMANN)) + rel_tol = float(params.get("rel_tol", 0.005)) + abs_floor = float(params.get("abs_floor", 0.5)) + # NaN comparisons and infinite allowances can silently accept bad output. + if not (np.isfinite(sigma) and sigma > 0.0): + raise ValueError(f"radiative_identity needs a finite sigma > 0, not {sigma!r}") + if not (np.isfinite(rel_tol) and rel_tol >= 0.0): + raise ValueError(f"radiative_identity needs a finite rel_tol >= 0, not {rel_tol!r}") + if not (np.isfinite(abs_floor) and abs_floor > 0.0): + raise ValueError(f"radiative_identity needs a finite abs_floor > 0, not {abs_floor!r}") + + w = make_window(run, probe) + for var in (upward, temperature): + if var not in w.table.columns: + raise ValueError(f"radiative_identity needs '{var}' in the model result") + if downward not in w.forcing.columns: + raise ValueError( + f"radiative_identity needs forcing column '{downward}'; this probe's " + "generator does not produce it" + ) + + # Output-derived emissivity would let a model force its own consistency. + static = run.case.static + if emissivity_key not in static: + raise ValueError( + f"radiative_identity needs '{emissivity_key}' in the case's static " + "attributes; this probe's generator does not supply it" + ) + eps = float(static[emissivity_key]) + if not 0.0 < eps <= 1.0: + raise ValueError(f"emissivity '{emissivity_key}' is {eps!r}, outside (0, 1]") + + rlds = w.forcing[downward].to_numpy(dtype=float) + if not np.isfinite(rlds).all(): + raise ValueError(f"forcing column '{downward}' has non-finite values") + + rlus = w.table[upward].to_numpy(dtype=float) + ts = w.table[temperature].to_numpy(dtype=float) + finite = np.isfinite(rlus) & np.isfinite(ts) + if not finite.all(): + n_bad = int((~finite).sum()) + return CriterionResult( + name="radiative_identity", + status=FAIL, + message=f"non-finite surface temperature or upward longwave on {n_bad} scored steps", + diagnostics={"non_finite_steps": n_bad}, + ) + # The fourth power hides the sign: -290 K would otherwise emit as +290 K. + below_zero = ts <= 0.0 + if below_zero.any(): + n_bad = int(below_zero.sum()) + return CriterionResult( + name="radiative_identity", + status=FAIL, + message=( + f"'{temperature}' is at or below 0 K on {n_bad} scored steps; the surface " + "temperature must be reported in kelvin" + ), + diagnostics={"non_positive_kelvin_steps": n_bad}, + ) + + # Finite inputs can overflow. Infinite allowances or NaN ratios can + # silently pass; rejecting non-finite ratios also keeps diagnostics finite. + with np.errstate(over="ignore", invalid="ignore"): + expected = eps * sigma * ts**4 + (1.0 - eps) * rlds + residual = rlus - expected + allowance = np.maximum(rel_tol * np.abs(rlus), abs_floor) + slack = np.abs(residual) / allowance + if not np.isfinite(allowance).all(): + raise ValueError("radiative_identity needs a finite allowance; rel_tol overflows") + if not np.isfinite(slack).all(): + return CriterionResult( + name="radiative_identity", + status=FAIL, + message="surface radiation calculation overflowed on scored steps", + diagnostics={"non_finite_calculation_steps": int((~np.isfinite(slack)).sum())}, + ) + violating = slack > 1.0 + + n_bad = int(violating.sum()) + worst = int(slack.argmax()) + times = w.forcing["time"].astype(str).to_numpy() + max_residual = float(np.abs(residual).max()) + # Scale before averaging so finite residuals do not overflow in the sum. + residual_scale = max(1.0, max_residual) + + detail = ( + f"worst step {slack[worst]:.2f} of tolerance at {times[worst]} " + f"(residual {residual[worst]:+.3g} W m-2 against {allowance[worst]:.3g} allowed; " + f"eps {eps:.4g})" + ) + if n_bad == 0: + message = ( + f"upward longwave agrees with the reported surface temperature on all " + f"{len(rlus)} steps: {detail}" + ) + else: + first = int(violating.argmax()) + message = ( + f"upward longwave contradicts the reported surface temperature on {n_bad} of " + f"{len(rlus)} steps, first at index {first} ({times[first]}): " + f"{rlus[first]:.3f} W m-2 reported, {expected[first]:.3f} expected from " + f"{temperature} {ts[first]:.2f} K and {downward} {rlds[first]:.1f} W m-2; {detail}" + ) + + return CriterionResult( + name="radiative_identity", + status=PASS if n_bad == 0 else FAIL, + value=float(slack[worst]), + threshold=1.0, + message=message, + diagnostics={ + "violating_steps": n_bad, + "scored_steps": len(rlus), + "worst_slack": float(slack[worst]), + "worst_step": { + "index": worst, + "time": times[worst], + "upward_w_m2": float(rlus[worst]), + "expected_w_m2": float(expected[worst]), + "residual_w_m2": float(residual[worst]), + "allowance_w_m2": float(allowance[worst]), + "temperature_k": float(ts[worst]), + "downward_w_m2": float(rlds[worst]), + }, + "max_abs_residual_w_m2": max_residual, + "mean_abs_residual_w_m2": float((np.abs(residual) / residual_scale).mean() * residual_scale), + # Identify weak fluxes where the absolute floor sets the bound. + "floor_steps": int((rel_tol * np.abs(rlus) < abs_floor).sum()), + "emissivity": eps, + "sigma_w_m2_k4": sigma, + "rel_tol": rel_tol, + "abs_floor_w_m2": abs_floor, + }, + ) diff --git a/src/hydroturing/harness.py b/src/hydroturing/harness.py index c0703bb9..3da292e4 100644 --- a/src/hydroturing/harness.py +++ b/src/hydroturing/harness.py @@ -307,6 +307,7 @@ def compatibility_issues( case: Case | None = None, *, check_perturbation: bool = True, + check_declared_inputs: bool = True, ) -> list[str]: """Explain why a model cannot be meaningfully run on a probe. @@ -330,9 +331,28 @@ def compatibility_issues( issues.append("model does not declare support for paired perturbation cases") if case is not None: visible = {c for c in case.forcing.columns if not c.startswith("_")} - missing = [v for v in model.needs_forcing if v not in visible] - if missing: - issues.append("forcing does not provide " + ", ".join(missing)) + for what, needed, provided in ( + ("forcing", model.needs_forcing, visible), + ("static", model.needs_static, case.static), + ): + missing = [name for name in needed if name not in provided] + if missing: + issues.append(f"{what} does not provide " + ", ".join(missing)) + # A probe whose verdict rests on a case-supplied input can only judge a + # model that declares it consumes that input. One that estimates its own + # sky or emissivity would be scored against values it never read. The + # adapter smoke test skips this: it asks whether the model can run, not + # whether a verdict can be given. + if check_declared_inputs: + for what, required, declared in ( + ("forcing", probe.requires_forcing, model.needs_forcing + model.uses_forcing), + ("static", probe.requires_static, model.needs_static + model.uses_static), + ): + undeclared = [name for name in required if name not in declared] + if undeclared: + issues.append( + f"model does not declare that it consumes {what} " + ", ".join(undeclared) + ) return issues @@ -382,10 +402,13 @@ def verify_adapter_contract( to the probe's later N/A (INCOMPLETE), not to this smoke test. A probe the model cannot consume is another matter: a step it does not - declare, a forcing it needs and the probe does not generate, a window - that drops a stretch the probe scores. There is nothing to run the - adapter on, so IncompatibleError is raised before it is invoked, on the - same grounds that make `run_probe` call the probe N/A (INCOMPATIBLE). + declare, a forcing or static input it needs and the probe does not + generate, a window that drops a stretch the probe scores. There is + nothing to run the adapter on, so IncompatibleError is raised before it + is invoked, on the same grounds that make `run_probe` call the probe N/A + (INCOMPATIBLE). What the scoring requires the model to have read is not + asked here: that decides whether a verdict can be given, not whether the + adapter can run. The case is cut to the model's evaluation window, as it will be in the real run, so a model that only fits its time budget on the window is @@ -394,7 +417,9 @@ def verify_adapter_contract( # The control variant is the one at the model's own step, so a daily # model is smoke-tested on daily rows rather than on a month of minutes. case = build_case(probe, seed, select_variants(model, probe)[0]) - issues = compatibility_issues(model, probe, case, check_perturbation=False) + issues = compatibility_issues( + model, probe, case, check_perturbation=False, check_declared_inputs=False + ) if issues: raise IncompatibleError(issues) @@ -410,6 +435,7 @@ def verify_adapter_contract( probe, requires_fluxes=model.emits_fluxes, requires_states=model.emits_states, + requires_diagnostics=model.emits_diagnostics, variants=(), criteria=(), ) diff --git a/src/hydroturing/protocol.py b/src/hydroturing/protocol.py index 720aafee..2a625bc3 100644 --- a/src/hydroturing/protocol.py +++ b/src/hydroturing/protocol.py @@ -128,13 +128,16 @@ def stage(io_dir: Path, case: Case, probe: ProbeSpec, model: ModelManifest) -> P # invariant across probes and prevents it identifying the criterion. "fluxes": list(model.emits_fluxes), "states": list(model.emits_states), + "diagnostics": list(model.emits_diagnostics), }, "input": {"forcing": FORCING_FILE, "static": STATIC_FILE}, "output": {"table": RESULT_CSV, "run": RUN_FILE}, "units": {v: UNITS[v] for v in model.emitted if v in UNITS}, "notes": ( "States are absolute storages, not tendencies. The harness " - "differences them itself." + "differences them itself. Row i's ts and rlus are instantaneous " + "values at row i's time, the same instant as row i's rlds, " + "when supplied." ), } request_path = io_dir / REQUEST_FILE diff --git a/src/hydroturing/spec.py b/src/hydroturing/spec.py index efed0d98..8562bcdb 100644 --- a/src/hydroturing/spec.py +++ b/src/hydroturing/spec.py @@ -24,9 +24,14 @@ # positive into the ground at the actual soil surface for `hfg`. Adapters # must correct a deeper-boundary flux for heat storage above that depth; # subsurface storage is not also subtracted from the surface budget. -# `sbl` is a component of `evspsbl`, never an addition to it. -FLUX_VARS = ("pr", "evspsbl", "mrro", "dis", "gwex", "sbl", "hfls", "hfss", "hfg") +# `sbl` is a component of `evspsbl`, never an addition to it. `rlus` is the +# total upward longwave radiation, surface emission plus reflected downward +# longwave, positive away from the surface. +FLUX_VARS = ("pr", "evspsbl", "mrro", "dis", "gwex", "sbl", "hfls", "hfss", "hfg", "rlus") STATE_VARS = ("mrso", "snw", "canopy", "gw", "channel") +# Keep diagnostics out of STATE_VARS: closure sums every reported store, +# and temperature must never be added to water storage. +DIAG_VARS = ("ts",) UNITS = { "pr": "mm day-1", @@ -41,6 +46,9 @@ "hfls": "W m-2", "hfss": "W m-2", "hfg": "W m-2", + "rlus": "W m-2", + # Instantaneous skin temperature; radiation uses kelvin, unlike forcing tas. + "ts": "K", "mrso": "mm", "snw": "mm", "canopy": "mm", @@ -101,6 +109,9 @@ "reference_area_leak", "reference_overshooting", "reference_sublimating", + "reference_radiative", + "reference_air_emitter", + "reference_no_reflection", } @@ -169,10 +180,17 @@ class ProbeSpec: # or a dry-down; such a probe asks for at least a year, and a submitted # model's window is widened to it. min_window_days: int = 0 + # Missing diagnostic outputs cause INCOMPLETE, as for missing fluxes. + requires_diagnostics: tuple[str, ...] = () + # Case-supplied inputs the verdict rests on: forcing columns and + # static.json keys a model must declare it consumes, or it is judged + # against values it never read and is INCOMPATIBLE instead. + requires_forcing: tuple[str, ...] = () + requires_static: tuple[str, ...] = () @property def required_vars(self) -> tuple[str, ...]: - return self.requires_fluxes + self.requires_states + return self.requires_fluxes + self.requires_states + self.requires_diagnostics @property def slug(self) -> str: @@ -265,6 +283,13 @@ class ModelManifest: # FULL_WINDOW for the whole record, or None to take the default for the # kind of model (see DEFAULT_WINDOW_DAYS). window_days: int | str | None = None + # Diagnostics the model reports in addition to its fluxes and states. + emits_diagnostics: tuple[str, ...] = () + # Static inputs the adapter cannot run without. + needs_static: tuple[str, ...] = () + # Optional inputs the adapter consumes whenever the case supplies them. + uses_forcing: tuple[str, ...] = () + uses_static: tuple[str, ...] = () @property def timestep(self) -> str: @@ -276,7 +301,7 @@ def supports_timestep(self, timestep: str) -> bool: @property def emitted(self) -> tuple[str, ...]: - return self.emits_fluxes + self.emits_states + return self.emits_fluxes + self.emits_states + self.emits_diagnostics def missing_for(self, probe: ProbeSpec) -> list[str]: """Variables the probe needs that this model never reports. @@ -366,6 +391,17 @@ def load_probe(path: str | Path) -> ProbeSpec: ) requires = raw.get("requires", {}) + # Refuse a required name no manifest can declare, as load_model refuses an + # unknown emission: a typo would make every model INCOMPLETE, and a flux + # asked for as a state would be met by the flux and never noticed. + unknown = [ + name + for key, known in (("fluxes", FLUX_VARS), ("states", STATE_VARS), ("diagnostics", DIAG_VARS)) + for name in requires.get(key, []) + if name not in known + ] + if unknown: + raise SpecError(f"{spec_file}: unknown variables in requires: {unknown}") return ProbeSpec( id=raw["id"], title=raw["title"], @@ -376,6 +412,9 @@ def load_probe(path: str | Path) -> ProbeSpec: citation=raw.get("citation", ""), requires_fluxes=tuple(requires.get("fluxes", [])), requires_states=tuple(requires.get("states", [])), + requires_diagnostics=tuple(requires.get("diagnostics", [])), + requires_forcing=tuple(requires.get("forcing", [])), + requires_static=tuple(requires.get("static", [])), generator=case["generator"], n_seeds=case["n_seeds"], timestep=case["timestep"], @@ -448,6 +487,7 @@ def load_model(path: str | Path) -> ModelManifest: unknown = [v for v in raw["emits"]["fluxes"] if v not in FLUX_VARS] unknown += [v for v in raw["emits"]["states"] if v not in STATE_VARS] + unknown += [v for v in raw["emits"].get("diagnostics", []) if v not in DIAG_VARS] if unknown: raise SpecError(f"{spec_file}: unknown variables in emits: {unknown}") @@ -463,8 +503,12 @@ def load_model(path: str | Path) -> ModelManifest: timesteps=timesteps, emits_fluxes=tuple(raw["emits"]["fluxes"]), emits_states=tuple(raw["emits"]["states"]), + emits_diagnostics=tuple(raw["emits"].get("diagnostics", [])), runner=runner, needs_forcing=tuple(raw.get("needs_forcing", [])), + needs_static=tuple(raw.get("needs_static", [])), + uses_forcing=tuple(raw.get("uses_forcing", [])), + uses_static=tuple(raw.get("uses_static", [])), supports_perturbation=bool(raw.get("supports", {}).get("perturbation", False)), resources=raw.get("resources", {}), authors=tuple(raw.get("authors", [])), diff --git a/templates/probe.template.yaml b/templates/probe.template.yaml index 2dc4fe6b..ccda50b5 100644 --- a/templates/probe.template.yaml +++ b/templates/probe.template.yaml @@ -52,6 +52,10 @@ description: > #! #! fluxes: pr, evspsbl, mrro, dis (rates) #! states: mrso, snw, canopy (absolute storages, never tendencies) +#! diagnostics: ts (read by a criterion, never integrated) +#! forcing, static: case inputs the verdict rests on; a model must declare +#! it consumes them (needs_* or uses_*) or it is +#! N/A with reason INCOMPATIBLE requires: fluxes: [pr, evspsbl, mrro] states: [mrso, snw, canopy] diff --git a/tests/test_docs_in_sync.py b/tests/test_docs_in_sync.py index 3ed4ebc3..96acf8fe 100644 --- a/tests/test_docs_in_sync.py +++ b/tests/test_docs_in_sync.py @@ -15,7 +15,7 @@ import pytest import yaml -from hydroturing.spec import FLUX_VARS, REPO_ROOT, STATE_VARS, load_model, load_probe +from hydroturing.spec import DIAG_VARS, FLUX_VARS, REPO_ROOT, STATE_VARS, load_model, load_probe PROBES = {p.id: p for p in (load_probe(x) for x in REPO_ROOT.glob("probes/*/*/probe.yaml"))} PROBE_IDS = sorted(PROBES) @@ -197,7 +197,7 @@ def test_site_flowchart_shows_every_probe(): assert html.count(f'"{key}":') == 3, f"{key} lacks a translation in one of the three languages" -@pytest.mark.parametrize("variable", FLUX_VARS + STATE_VARS) +@pytest.mark.parametrize("variable", FLUX_VARS + STATE_VARS + DIAG_VARS) def test_agents_doc_defines_every_variable(variable): assert f"| `{variable}` |" in read("AGENTS.md"), ( f"{variable} is accepted by spec.py but the table in AGENTS.md does not define it" diff --git a/tests/test_harness.py b/tests/test_harness.py index e1a18e7c..7200912c 100644 --- a/tests/test_harness.py +++ b/tests/test_harness.py @@ -169,6 +169,36 @@ def test_missing_forcing_is_reported_as_incompatible(probe): assert "unavailable_driver" in outcome.incompatible[0] +def test_missing_static_is_reported_as_incompatible(probe): + """A static key the adapter would read and the case does not supply is an + incompatibility found before the run, not a KeyError inside the adapter.""" + model = replace(registry.find_model("reference_bucket"), needs_static=("eps",)) + outcome = run_probe(model, probe, [11]) + assert (outcome.verdict, outcome.reason) == (NOT_SCORED, INCOMPATIBLE) + assert outcome.incompatible == ["static does not provide eps"] + + +@pytest.mark.parametrize("kind, name", [("forcing", "rlds"), ("static", "eps")]) +def test_optional_declaration_does_not_cancel_a_required_input(probe, kind, name): + model = replace( + registry.find_model("reference_bucket"), + **{f"needs_{kind}": (name,), f"uses_{kind}": (name,)}, + ) + outcome = run_probe(model, probe, [11]) + assert (outcome.verdict, outcome.reason) == (NOT_SCORED, INCOMPATIBLE) + assert outcome.incompatible == [f"{kind} does not provide {name}"] + + +def test_undeclared_case_inputs_are_reported_as_incompatible(probe): + """A verdict that rests on a case-supplied input can only be given to a + model that says it read that input.""" + rests_on = replace(probe, requires_forcing=("rlds",), requires_static=("eps",)) + outcome = run_probe(registry.find_model("reference_bucket"), rests_on, [11]) + assert (outcome.verdict, outcome.reason) == (NOT_SCORED, INCOMPATIBLE) + assert "forcing rlds" in outcome.incompatible[0] + assert "static eps" in outcome.incompatible[1] + + def test_paired_probe_requires_declared_perturbation_support(probe): paired = replace(probe, variants=("control", "perturbed")) model = replace(registry.find_model("reference_bucket"), supports_perturbation=False) @@ -359,6 +389,93 @@ def test_request_hides_probe_identity_and_generator_seed(probe, tmp_path): assert "spinup_steps" not in request assert request["request"]["fluxes"] == list(model.emits_fluxes) assert request["request"]["states"] == list(model.emits_states) + assert request["request"]["diagnostics"] == list(model.emits_diagnostics) + + +# --- diagnostics ------------------------------------------------------------ +# The third output category is threaded through the probe spec, the manifest, +# the request, the contract check and the verdict. Each join is pinned here so +# that a later change cannot drop one of them silently. + + +def test_a_required_diagnostic_the_model_lacks_is_incomplete(probe): + needs_ts = replace(probe, requires_diagnostics=("ts",)) + model = registry.find_model("reference_bucket") + assert "ts" in needs_ts.required_vars + assert model.missing_for(needs_ts) == ["ts"] + outcome = run_probe(model, needs_ts, [11]) + assert outcome.verdict == NOT_SCORED + assert outcome.reason == INCOMPLETE + assert outcome.missing == ["ts"] + + +def test_manifest_declares_diagnostics_under_their_own_key(tmp_path): + import shutil + + from hydroturing.spec import load_model + + target = tmp_path / "skin_model" + shutil.copytree(registry.MODELS_DIR / "_template", target) + manifest = target / "model.yaml" + manifest.write_text( + manifest.read_text() + .replace("name: _template", "name: skin_model") + .replace(" states: [", " diagnostics: [ts]\n states: [") + ) + model = load_model(target) + assert model.emits_diagnostics == ("ts",) + assert "ts" in model.emitted + assert "ts" not in model.emits_fluxes + model.emits_states + + manifest.write_text(manifest.read_text().replace("diagnostics: [ts]", "diagnostics: [skin]")) + with pytest.raises(SpecError, match="unknown variables"): + load_model(target) + + +def test_required_diagnostics_are_read_from_probe_yaml_and_emitted_ones_reach_the_request(tmp_path): + import shutil + + from hydroturing.spec import load_probe + + source = registry.PROBES_DIR / "mass" / "catchment-closure" + target = tmp_path / "mass" / "catchment-closure" + shutil.copytree(source, target, ignore=shutil.ignore_patterns("__pycache__")) + spec_file = target / "probe.yaml" + spec_file.write_text( + spec_file.read_text().replace( + " states: [mrso, snw, canopy]\n", " states: [mrso, snw, canopy]\n diagnostics: [ts]\n", 1 + ) + ) + needs_ts = load_probe(target) + assert needs_ts.requires_diagnostics == ("ts",) + assert needs_ts.required_vars[-1] == "ts" + + model = replace(registry.find_model("reference_bucket"), emits_diagnostics=("ts",)) + request_path = stage(tmp_path / "io", build_case(needs_ts, 1), needs_ts, model) + request = json.loads(request_path.read_text()) + assert request["request"]["diagnostics"] == ["ts"] + assert request["units"]["ts"] == "K" + + +def test_adapter_verification_rejects_a_declared_but_unreported_diagnostic(probe): + """Declaring `ts` and not writing it is a contract breach, not INCOMPLETE.""" + model = replace(registry.find_model("reference_bucket"), emits_diagnostics=("ts",)) + with pytest.raises(ProtocolError, match=r"missing requested variables: \['ts'\]"): + verify_adapter_contract(model, probe, gate_seeds(probe.id, 1)[0]) + + +def test_a_diagnostic_is_never_counted_as_a_store(probe): + from hydroturing.criteria.base import make_window, reported_states + from hydroturing.protocol import RunResult + + case = build_case(probe, 3) + table = pd.DataFrame({ + "time": case.forcing["time"], "mrso": 100.0, "snw": 0.0, "canopy": 0.0, "ts": 290.0, + }) + run = RunResult(case=case, table=table, meta={}, wall_seconds=0.0) + window = make_window(run, probe) + assert "ts" not in reported_states(window, probe) + assert window.storage(reported_states(window, probe))[0] == pytest.approx(100.0) # --- container isolation ---------------------------------------------------- @@ -869,6 +986,25 @@ def test_unknown_criterion_is_rejected(tmp_path, monkeypatch): load_probe(target) +@pytest.mark.parametrize("key,name", [("diagnostics", "skin"), ("states", "mrro")]) +def test_a_required_name_no_manifest_can_declare_is_rejected(key, name, tmp_path, monkeypatch): + """Reject unknown required outputs and variables in the wrong category.""" + import yaml + + from hydroturing import scaffold + from hydroturing.spec import load_probe + + monkeypatch.setattr(scaffold, "PROBES_DIR", tmp_path / "probes") + target, _ = _scaffold_template(scaffold, tmp_path, "default", "unknown-requirement") + spec_file = target / "probe.yaml" + raw = yaml.safe_load(spec_file.read_text()) + raw["requires"][key] = [name] + spec_file.write_text(yaml.safe_dump(raw, sort_keys=False)) + + with pytest.raises(SpecError, match="unknown variables in requires"): + load_probe(target) + + @pytest.mark.parametrize("kind", _template_kinds()) def test_every_template_scaffolds_into_a_valid_probe(kind, tmp_path, monkeypatch): from hydroturing import scaffold diff --git a/tests/test_radiation_consistency.py b/tests/test_radiation_consistency.py new file mode 100644 index 00000000..dea445a6 --- /dev/null +++ b/tests/test_radiation_consistency.py @@ -0,0 +1,266 @@ +"""Hand-built checks of the radiation identity, tolerance and invalid inputs.""" + +from __future__ import annotations + +import json + +import numpy as np +import pandas as pd +import pytest + +from hydroturing.criteria import get +from hydroturing.criteria.base import FAIL, PASS +from hydroturing.protocol import Case, RunResult + +# CODATA 2018, typed here rather than imported so that the module under test +# cannot vouch for its own constant. +SIGMA = 5.670374419e-8 + + +def upward(ts, rlds, eps): + """The exact identity, the only physics these tests need.""" + return eps * SIGMA * np.asarray(ts, dtype=float) ** 4 + (1.0 - eps) * np.asarray(rlds, dtype=float) + + +def build(ts, rlus=None, rlds=300.0, eps=0.98, spinup=0, static=None): + """A run reporting `ts` and `rlus` under a sky of `rlds`; exact unless told otherwise.""" + ts = np.asarray(ts, dtype=float) + n = len(ts) + rlds = np.broadcast_to(np.asarray(rlds, dtype=float), (n,)) + if rlus is None: + rlus = upward(ts, rlds, eps) + times = pd.date_range("2001-07-01T06:00", periods=n, freq="h") + forcing = pd.DataFrame({"time": times, "rlds": rlds}) + table = pd.DataFrame({"time": times, "ts": ts, "rlus": np.asarray(rlus, dtype=float)}) + case = Case( + probe_id="t", seed=1, forcing=forcing, + static={"eps": eps} if static is None else static, + spinup_steps=spinup, timestep="PT1H", + ) + return RunResult(case=case, table=table, meta={}, wall_seconds=0.0) + + +def score(run, **params): + return get("radiative_identity")(run, None, params) + + +def diurnal(n=48, mean=290.0, amplitude=8.0): + hours = 6.0 + np.arange(n) + return mean + amplitude * np.cos(2.0 * np.pi * (hours - 14.0) / 24.0) + + +def test_exact_identity_passes_at_every_step(): + run = build(diurnal(), rlds=300.0 + 40.0 * np.sin(np.arange(48) / 5.0)) + result = score(run) + assert result.status == PASS + assert result.diagnostics["violating_steps"] == 0 + assert result.diagnostics["scored_steps"] == 48 + assert result.diagnostics["max_abs_residual_w_m2"] < 1e-9 + assert result.value < 1e-9 + + +@pytest.mark.parametrize("ts,error,expected,floor_steps", [ + # At 290 K the exact flux is 399.0 W m-2, so 0.5 percent is about 2.0. + (290.0, 1.9, PASS, 0), + (290.0, -1.9, PASS, 0), + (290.0, 2.1, FAIL, 0), + (290.0, -2.1, FAIL, 0), + # At 180 K the flux is about 64 W m-2 and 0.5 percent of it is below the + # 0.5 W m-2 floor, so the floor is the bound on every step. + (180.0, 0.45, PASS, 24), + (180.0, 0.55, FAIL, 24), +]) +def test_tolerance_is_relative_to_the_reported_flux_with_a_floor(ts, error, expected, floor_steps): + ts = np.full(24, ts) + run = build(ts, rlus=upward(ts, 300.0, 0.98) + error) + result = score(run) + assert result.status == expected + assert result.diagnostics["floor_steps"] == floor_steps + + +def test_the_relative_bound_is_taken_against_the_reported_flux(): + ts = np.full(12, 290.0) + exact = upward(ts, 300.0, 0.98) + # Over-reporting by 0.5 percent of the expected value passes, because the + # bound is 0.5 percent of the larger, reported, flux. + assert score(build(ts, rlus=exact * 1.005)).status == PASS + # Under-reporting by the same amount fails: the bound shrank with the flux. + assert score(build(ts, rlus=exact * 0.995)).status == FAIL + + +@pytest.mark.parametrize("eps", [0.95, 0.98, 0.99]) +def test_reporting_emission_alone_as_upward_longwave_fails(eps): + """Omit reflected sky across the probe's emissivity range. At eps 0.99, + the fixed sky contributes 3 W m-2 of missing reflection.""" + ts = diurnal() + run = build(ts, rlus=eps * SIGMA * ts**4, rlds=300.0, eps=eps) + result = score(run) + assert result.status == FAIL + assert result.diagnostics["violating_steps"] == 48 + assert result.diagnostics["worst_step"]["residual_w_m2"] == pytest.approx(-(1.0 - eps) * 300.0) + + +@pytest.mark.parametrize("offset_k,expected", [ + # Around 290 K under a 300 W m-2 sky, an offset below about 0.4 K sits + # inside 0.5 percent of the flux: the criterion's resolution under these + # conditions, stated rather than hidden. It is coarser at a warmer surface. + (0.3, PASS), + (1.0, FAIL), + (5.0, FAIL), + (-5.0, FAIL), +]) +def test_emitting_at_a_different_temperature_than_reported_fails(offset_k, expected): + """The proposal's first negative control: rlus from the air, ts from the surface.""" + ts = diurnal() + run = build(ts, rlus=upward(ts + offset_k, 300.0, 0.98)) + result = score(run) + assert result.status == expected + if expected == FAIL: + assert result.diagnostics["violating_steps"] > 0 + + +@pytest.mark.parametrize("amplitude_k,expected", [(5.0, PASS), (10.0, FAIL)]) +def test_a_linearised_stefan_boltzmann_law_is_seen_only_on_large_excursions(amplitude_k, expected): + """The proposal's optional third control, reported rather than gated: the + second-order term eps*sigma*6*T0**2*dT**2 is 0.7 W m-2 at 5 K and 2.8 at + 10 K around 288 K, against a bound near 2 W m-2.""" + t0, eps = 288.0, 0.98 + ts = diurnal(mean=t0, amplitude=amplitude_k) + linear = eps * SIGMA * (t0**4 + 4.0 * t0**3 * (ts - t0)) + (1.0 - eps) * 300.0 + result = score(build(ts, rlus=linear, eps=eps)) + assert result.status == expected + second_order = eps * SIGMA * 6.0 * t0**2 * amplitude_k**2 + assert result.diagnostics["max_abs_residual_w_m2"] == pytest.approx(second_order, rel=0.05) + + +def test_spinup_is_not_scored_and_indices_start_at_the_window(): + ts = diurnal(n=72) + rlus = upward(ts, 300.0, 0.98) + rlus[:24] += 50.0 + rlus[24 + 7] += 3.0 + result = score(build(ts, rlus=rlus, spinup=24)) + assert result.status == FAIL + assert result.diagnostics["scored_steps"] == 48 + assert result.diagnostics["violating_steps"] == 1 + assert result.diagnostics["worst_step"]["index"] == 7 + assert result.diagnostics["worst_step"]["time"].startswith("2001-07-02 13:00") + assert result.diagnostics["worst_step"]["residual_w_m2"] == pytest.approx(3.0) + + +@pytest.mark.parametrize("column", ["ts", "rlus"]) +@pytest.mark.parametrize("invalid", [np.nan, np.inf, -np.inf]) +def test_non_finite_model_output_is_a_failure(column, invalid): + run = build(diurnal(n=24)) + run.table.loc[12, column] = invalid + result = score(run) + assert result.status == FAIL + assert result.diagnostics["non_finite_steps"] == 1 + + +def test_finite_temperature_that_overflows_the_identity_fails(): + run = build([1e100], rlus=[500.0]) + result = score(run) + assert result.status == FAIL + assert result.diagnostics["non_finite_calculation_steps"] == 1 + + +def test_large_finite_residuals_keep_the_report_finite(): + run = build([290.0, 290.0], rlus=[1e308, 1e308]) + with np.errstate(over="raise", invalid="raise"): + result = score(run) + assert result.status == FAIL + assert result.diagnostics["mean_abs_residual_w_m2"] == pytest.approx(1e308) + json.dumps(result.diagnostics, allow_nan=False) + + +@pytest.mark.parametrize("temperature", [0.0, -290.0]) +def test_a_temperature_at_or_below_zero_kelvin_is_named_as_a_unit_mistake(temperature): + # Even a matching fourth-power emission must not validate a non-positive K. + ts = np.full(24, temperature) + result = score(build(ts)) + assert result.status == FAIL + assert "kelvin" in result.message + assert result.diagnostics["non_positive_kelvin_steps"] == 24 + + +def test_a_warm_surface_reported_in_celsius_fails_on_the_residual(): + ts = np.full(24, 17.0) + result = score(build(ts, rlus=upward(ts + 273.15, 300.0, 0.98))) + assert result.status == FAIL + assert result.diagnostics["violating_steps"] == 24 + assert "17.00 K" in result.message + + +@pytest.mark.parametrize("setting", [ + {"sigma": 0.0}, {"sigma": -1.0}, {"sigma": np.nan}, {"sigma": np.inf}, + {"rel_tol": -0.01}, {"rel_tol": np.nan}, {"rel_tol": np.inf}, + {"rel_tol": 1e308}, {"sigma": 1e308, "rel_tol": 1e308}, + {"abs_floor": 0.0}, {"abs_floor": -1.0}, {"abs_floor": np.nan}, {"abs_floor": np.inf}, +], ids=lambda s: "=".join(f"{k}{v}" for k, v in s.items())) +def test_an_invalid_setting_is_an_error_rather_than_a_pass(setting): + # About 100 W m-2 of residual, which any sane setting rejects. A NaN or an + # infinite bound would count zero violations and pass it. + ts = np.full(24, 290.0) + run = build(ts, rlus=np.full(24, 500.0)) + assert score(run).status == FAIL + with pytest.raises(ValueError, match="finite"): + score(run, **setting) + + +def test_missing_inputs_and_a_bad_emissivity_are_errors_rather_than_verdicts(): + run = build(diurnal(n=24)) + + run.case.forcing = run.case.forcing.drop(columns="rlds") + with pytest.raises(ValueError, match="forcing column 'rlds'"): + score(run) + + run = build(diurnal(n=24), static={}) + with pytest.raises(ValueError, match="'eps' in the case's static"): + score(run) + + for bad in (0.0, 1.2, -0.5, np.nan, np.inf): + with pytest.raises(ValueError, match="outside"): + score(build(diurnal(n=24), static={"eps": bad})) + + run = build(diurnal(n=24)) + run.case.forcing.loc[3, "rlds"] = np.nan + with pytest.raises(ValueError, match="non-finite"): + score(run) + + run = build(diurnal(n=24)) + run.table = run.table.drop(columns="ts") + with pytest.raises(ValueError, match="needs 'ts'"): + score(run) + + +def test_column_names_emissivity_key_and_constants_are_configurable(): + ts = diurnal(n=24) + sigma = 5.67e-8 + run = build(ts, rlus=0.97 * sigma * ts**4 + 0.03 * 280.0, rlds=280.0, static={"emissivity": 0.97}) + run.case.forcing = run.case.forcing.rename(columns={"rlds": "lw_down"}) + run.table = run.table.rename(columns={"ts": "skin", "rlus": "lw_up"}) + result = score( + run, upward="lw_up", temperature="skin", downward="lw_down", + emissivity="emissivity", sigma=sigma, + ) + assert result.status == PASS + assert result.diagnostics["emissivity"] == 0.97 + assert result.diagnostics["sigma_w_m2_k4"] == sigma + # The two constants differ by 6.6e-5 relative, far inside the tolerance, + # so only an exact residual shows the parameter was used in the arithmetic. + assert result.diagnostics["max_abs_residual_w_m2"] < 1e-9 + + # Loosening either bound changes the verdict where that bound governs, + # and the report says which values were in force. + warm = build(ts, rlus=upward(ts, 300.0, 0.98) + 10.0) + assert score(warm).status == FAIL + loose = score(warm, rel_tol=0.05) + assert loose.status == PASS + assert loose.diagnostics["rel_tol"] == 0.05 + cold = np.full(24, 180.0) + off = build(cold, rlus=upward(cold, 300.0, 0.98) + 3.0) + assert score(off).status == FAIL + lifted = score(off, abs_floor=5.0) + assert lifted.status == PASS + assert lifted.diagnostics["abs_floor_w_m2"] == 5.0 + assert lifted.diagnostics["floor_steps"] == 24 diff --git a/tests/test_radiation_probe.py b/tests/test_radiation_probe.py new file mode 100644 index 00000000..cf835dde --- /dev/null +++ b/tests/test_radiation_probe.py @@ -0,0 +1,189 @@ +"""Radiation forcing, reproducibility, staging and evaluation-window checks.""" + +from __future__ import annotations + +import json +from dataclasses import replace + +import numpy as np +import pandas as pd +import pytest + +from hydroturing import registry +from hydroturing.harness import ( + build_case, + compatibility_issues, + load_generator, + resolve_window_days, + run_probe, + select_window, + verify_adapter_contract, + window_case, +) +from hydroturing.protocol import FORCING_FILE, STATIC_FILE, stage +from hydroturing.scoring import INCOMPATIBLE, NOT_SCORED +from hydroturing.seeds import gate_seeds + +SIGMA = 5.670374419e-8 +COLUMNS = ["time", "pr", "tas", "pet", "rn", "rlds"] + + +@pytest.fixture(scope="module") +def probe(): + return registry.find_probe("energy/radiation-consistency") + + +def test_spec_asks_for_the_identity_and_nothing_else(probe): + # Every extra required variable excludes models from being testable at + # all, so the requirement is pinned to the two outputs the identity needs. + assert probe.requires_fluxes == ("rlus",) + assert probe.requires_states == () + assert probe.requires_diagnostics == ("ts",) + assert [c.name for c in probe.criteria] == ["radiative_identity"] + assert probe.criteria[0].params["emissivity"] == "eps" + # The verdict rests on the case's sky and emissivity, so a model must + # declare that it consumes them. + assert probe.requires_forcing == ("rlds",) + assert probe.requires_static == ("eps",) + + +def test_adapter_verification_asks_only_what_the_adapter_needs(probe): + """The coupled reference reads neither rlds nor eps, so it cannot be + judged here, but its adapter can still be smoke-tested on this case.""" + result = verify_adapter_contract( + registry.find_model("reference_coupled"), probe, gate_seeds(probe.id, 1)[0] + ) + assert len(result.table) == 768 + + +@pytest.mark.parametrize("missing", [("rlds",), ("eps",), ("rlds", "eps")]) +def test_a_model_that_does_not_consume_the_sky_or_the_emissivity_is_not_judged(probe, missing): + """Emitting ts and rlus is not enough: a model with its own downward + longwave or its own emissivity would be scored against values it never + read, so it is INCOMPATIBLE rather than a candidate for VIOLATION.""" + reference = registry.find_model("reference_radiative") + assert compatibility_issues(reference, probe, build_case(probe, 0)) == [] + own_inputs = replace( + reference, + uses_forcing=() if "rlds" in missing else ("rlds",), + uses_static=() if "eps" in missing else ("eps",), + ) + outcome = run_probe(own_inputs, probe, gate_seeds(probe.id, 1)) + assert (outcome.verdict, outcome.reason) == (NOT_SCORED, INCOMPATIBLE) + expected = [] + if "rlds" in missing: + expected.append("model does not declare that it consumes forcing rlds") + if "eps" in missing: + expected.append("model does not declare that it consumes static eps") + assert outcome.incompatible == expected + + +@pytest.mark.parametrize("forcing_kind", ["needs", "uses"]) +@pytest.mark.parametrize("static_kind", ["needs", "uses"]) +def test_required_inputs_accept_mandatory_or_optional_declarations(probe, forcing_kind, static_kind): + model = registry.find_model("reference_radiative") + inputs = dict(needs_forcing=model.needs_forcing, needs_static=(), uses_forcing=(), uses_static=()) + inputs[f"{forcing_kind}_forcing"] += ("rlds",) + inputs[f"{static_kind}_static"] += ("eps",) + model = replace(model, **inputs) + assert compatibility_issues(model, probe, build_case(probe, 0)) == [] + + other = registry.find_probe("energy/surface-energy-closure") + expected = [] + if forcing_kind == "needs": + expected.append("forcing does not provide rlds") + if static_kind == "needs": + expected.append("static does not provide eps") + assert compatibility_issues(model, other, build_case(other, 4242)) == expected + + +def test_generator_reproduces_bytes_and_changes_the_case_with_seed(probe): + generate = load_generator(probe).generate + first, static = generate(0) + repeated, repeated_static = generate(0) + other, other_static = generate(1) + assert first.to_csv(index=False) == repeated.to_csv(index=False) + assert static == repeated_static + for column in ("pr", "tas", "rn", "rlds"): + assert not first[column].equals(other[column]), column + assert static["eps"] != other_static["eps"] + + +@pytest.mark.parametrize("seed", [0, 19]) +def test_case_is_hourly_and_starts_at_dawn(probe, seed): + case = build_case(probe, seed) + assert case.n_steps == 768 + assert case.spinup_steps == 48 + assert case.timestep == "PT1H" + assert list(case.forcing.columns) == COLUMNS + time = pd.to_datetime(case.forcing["time"]) + assert time.diff().iloc[1:].eq(pd.Timedelta(hours=1)).all() + assert time.iloc[0].hour == 6 + assert len(case.after_spinup(case.forcing)) == 720 + + +def test_every_seed_stays_in_the_warm_gray_surface_regime(probe): + """Check input bounds on the gate seeds and twenty additional seeds.""" + generate = load_generator(probe).generate + for seed in [*gate_seeds(probe.id, probe.n_seeds), *range(20)]: + frame, static = generate(seed) + where = f"seed {seed}" + assert 0.95 <= static["eps"] <= 0.99, where + assert static["canopy_capacity_mm"] == 0.0, where + assert np.isfinite(frame[COLUMNS[1:]].to_numpy()).all(), where + + # Warm enough that nothing freezes, which the snow-free boundary needs. + assert frame["tas"].between(15.0, 31.0).all(), where + assert (frame["pr"] >= 0.0).all() and (frame["pet"] > 0.0).all(), where + assert frame["rn"].between(-45.0, 560.0).all(), where + + # The sky radiates as a gray body at the air temperature, between a + # clear and an overcast emissivity, and varies hour to hour. + air_k = frame["tas"].to_numpy() + 273.15 + sky = frame["rlds"].to_numpy() / (SIGMA * air_k**4) + assert np.all((sky > 0.70) & (sky < 0.97)), where + assert frame["rlds"].between(280.0, 470.0).all(), where + assert frame["rlds"].std() > 20.0, where + + # The reflected term the second negative control drops is well above + # the criterion's 0.5 W m-2 floor at every hour, at any emissivity + # the case can draw. Whether it also clears the relative bound is a + # question for the reference models' actual fluxes, not for the input. + assert ((1.0 - static["eps"]) * frame["rlds"] >= 2.5).all(), where + + +def test_the_same_cloud_dims_the_sun_and_brightens_the_sky(probe): + """At noon, lower net radiation must go with a more emissive sky.""" + frame = build_case(probe, 3).forcing + noon = frame[pd.to_datetime(frame["time"]).dt.hour == 12] + air_k = noon["tas"].to_numpy() + 273.15 + sky = noon["rlds"].to_numpy() / (SIGMA * air_k**4) + assert len(noon) == 32 + assert np.corrcoef(noon["rn"].to_numpy(), sky)[0, 1] < -0.5 + + +def test_staging_hands_the_model_the_sky_and_the_emissivity(probe, tmp_path): + case = build_case(probe, 0) + model = registry.find_model("reference_bucket") + request = json.loads(stage(tmp_path, case, probe, model).read_text()) + staged = pd.read_csv(tmp_path / FORCING_FILE) + pd.testing.assert_frame_equal(staged, case.forcing) + static = json.loads((tmp_path / STATIC_FILE).read_text()) + assert static["eps"] == case.static["eps"] + assert request["timestep"] == "PT1H" + assert request["n_steps"] == 768 + assert "ts and rlus are instantaneous values at row i's time" in request["notes"] + assert "the same instant as row i's rlds" in request["notes"] + + +def test_submitted_models_see_every_scored_day(probe): + # A submitted hourly model normally gets seven days; this probe widens + # that to the whole scored record, even when a caller asks for less. + model = replace(registry.find_model("reference_bucket"), name="submitted_skin") + assert resolve_window_days(model, probe) == 30 + assert resolve_window_days(model, probe, override=7) == 30 + case = build_case(probe, 0) + cut = window_case(case, select_window(case, probe, 30)) + assert cut.spinup_steps == 48 + assert cut.window["rows"] == 720 + pd.testing.assert_frame_equal(cut.forcing, case.forcing) diff --git a/tests/test_radiation_references.py b/tests/test_radiation_references.py new file mode 100644 index 00000000..992700d7 --- /dev/null +++ b/tests/test_radiation_references.py @@ -0,0 +1,105 @@ +"""Adapter parity and radiation margins on two gate seeds and emissivity bounds.""" + +from __future__ import annotations + +import tempfile +from dataclasses import replace +from pathlib import Path + +import numpy as np +import pytest + +from hydroturing import registry +from hydroturing.criteria import get +from hydroturing.criteria.base import FAIL, PASS +from hydroturing.harness import build_case, run_probe +from hydroturing.runner import get_runner +from hydroturing.scoring import OK, PASS as PROBE_PASS +from hydroturing.seeds import gate_seeds + +KELVIN = 273.15 +# The design margin: each negative control's worst step must exceed the +# bound by this factor on every seed, so that the verdict does not rest on +# a step that sits at the tolerance. The criterion's own threshold stays 1. +MARGIN = 1.3 + + +@pytest.fixture(scope="module") +def probe(): + return registry.find_probe("energy/radiation-consistency") + + +def run(name, probe, case): + # Read back what the model declares, as verify-adapter does, so that the + # coupled reference can be run on a probe that asks for more than it has. + model = registry.find_model(name) + asks = replace( + probe, requires_fluxes=model.emits_fluxes, requires_states=model.emits_states, + requires_diagnostics=model.emits_diagnostics, + ) + with tempfile.TemporaryDirectory(prefix="ht-radiation-test-") as workdir: + return get_runner(model).run(model, asks, case, Path(workdir)) + + +def score(result, probe): + return get("radiative_identity")(result, probe, dict(probe.criteria[0].params)) + + +def cases(probe, eps=None): + for seed in gate_seeds(probe.id, 2): + case = build_case(probe, seed) + if eps is not None: + case.static["eps"] = eps + yield case + + +def test_positive_control_is_the_coupled_reference_with_a_skin(probe): + for case in cases(probe): + skin = run("reference_radiative", probe, case).table + coupled = run("reference_coupled", probe, case).table + assert set(skin.columns) - set(coupled.columns) == {"rlus", "ts"} + for column in coupled.columns: + assert np.array_equal(skin[column].to_numpy(), coupled[column].to_numpy()), column + ts = skin["ts"].to_numpy() + assert np.isfinite(ts).all() + assert ts.min() > KELVIN, "the snow-free boundary needs a skin that never freezes" + assert ts.max() < KELVIN + 60.0 + + +@pytest.mark.parametrize("name", [ + "reference_radiative", "reference_air_emitter", "reference_no_reflection", +]) +def test_optional_radiation_inputs_do_not_block_surface_energy_closure(name, tmp_path): + probe = registry.find_probe("energy/surface-energy-closure") + model = registry.find_model(name) + outcome = run_probe(model, probe, [4242], workdir=tmp_path) + assert (outcome.verdict, outcome.reason) == (PROBE_PASS, OK), outcome.incompatible + + +@pytest.mark.parametrize("negative", ["reference_air_emitter", "reference_no_reflection"]) +def test_each_negative_control_differs_in_upward_longwave_alone(probe, negative): + for case in cases(probe): + positive = run("reference_radiative", probe, case).table + broken = run(negative, probe, case).table + assert list(broken.columns) == list(positive.columns) + for column in positive.columns: + same = np.array_equal(positive[column].to_numpy(), broken[column].to_numpy()) + assert same == (column != "rlus"), column + + +@pytest.mark.parametrize("eps", [None, 0.95, 0.99], ids=["drawn", "eps=0.95", "eps=0.99"]) +def test_the_probe_separates_the_controls_with_margin(probe, eps): + """At the drawn emissivity and at both ends of the range the case can + draw. The reflected term the second control drops shrinks with (1 - eps), + so 0.99 is where its margin is thinnest.""" + for case in cases(probe, eps): + result = score(run("reference_radiative", probe, case), probe) + assert result.status == PASS + assert result.diagnostics["worst_slack"] < 1e-9 + for negative in ("reference_air_emitter", "reference_no_reflection"): + result = score(run(negative, probe, case), probe) + assert result.status == FAIL, negative + assert result.diagnostics["worst_slack"] >= MARGIN, (negative, case.static["eps"]) + # Dropping the reflected sky is wrong by (1 - eps) * rlds on every + # step, so the second control must violate at every scored step. + assert result.diagnostics["violating_steps"] == result.diagnostics["scored_steps"] diff --git a/tests/test_report_detail.py b/tests/test_report_detail.py index 7440cf4f..34d90655 100644 --- a/tests/test_report_detail.py +++ b/tests/test_report_detail.py @@ -34,6 +34,7 @@ "energy/evaporative-partition": ("partition_shift",), "energy/latent-heat-et-consistency": ("energy_closure", "flux_identity"), "energy/pet-consistency": ("demand_consistency",), + "energy/radiation-consistency": ("radiative_identity",), "energy/surface-energy-closure": ("energy_closure_by_phase",), "mass/antecedent-monotonicity": ("antecedent_monotonicity",), "mass/area-invariance": ("invariance",),