feat(magicc7): concentration-driven mode for CO2 / CH4 / N2O

changelog/104.docs.md

@@ -0,0 +1 @@
+Added a "Concentration-driven runs" documentation page comparing the MAGICC7, FaIR2 and CICEROSCMPY2 adapters' `RunMode.CONCENTRATION_DRIVEN` support.

changelog/104.feature.md

@@ -0,0 +1 @@
+Added `RunMode.CONCENTRATION_DRIVEN` support to the MAGICC7 adapter, driving the model from `Atmospheric Concentrations|*` inputs (overlaid on the RCMIP3 baseline) instead of emissions. Covers CO2, CH4 and N2O via MAGICC's per-gas concentration flags, plus the F-gas and Montreal-halocarbon groups via the bundled-array `FGAS_FILES_CONC` / `MHALO_FILES_CONC` flags. The concentration-driven decision is made per scenario, so a gas can be concentration-driven in one scenario and emissions-driven in another within the same batch. Requires `rcmip3_bundle_path`.

changelog/104.internal.md

@@ -0,0 +1 @@
+Added an opt-in MAGICC7 validation harness against IPCC AR6 Table 7.SM.4: load the AR6 probabilistic drawnset, run concentration-driven SSP scenarios, rebase GSAT to 1995-2014 and compare per-period percentiles to the published MAGICC7 column. Gated on a local drawnset (not vendored) and the MAGICC binary, so skipped by default.

changelog/109.improvement.md

@@ -0,0 +1 @@
+Sped up the FaIR2 adapter's output extraction by materialising FaIR's forcing array to numpy once and aggregating per ensemble member with integer indexing, instead of repeated per-member xarray orthogonal indexing. For an 842-member native-calibration ensemble this cut a pathway from ~370 s to ~40 s (~9x) and removed the superlinear scaling with member count. Results are numerically identical.

docs/source/concentration-driven.md

@@ -0,0 +1,74 @@
+# Concentration-driven runs
+
+OpenSCM-Runner can drive its three modern adapters — **MAGICC7**, **FaIR2**
+and **CICEROSCMPY2** — from prescribed atmospheric concentrations instead of
+emissions, via {class}`openscm_runner.RunMode`:
+
+```python
+from openscm_runner import RunMode, run
+
+run(
+    climate_models_cfgs={"MAGICC7": ({"core_climatesensitivity": 3,
+                                      "rcmip3_bundle_path": "/path/to/rcmip3"},)},
+    scenarios=scenarios,            # may carry ``Atmospheric Concentrations|*`` rows
+    output_variables=("Surface Air Temperature Change",
+                      "Atmospheric Concentrations|CO2"),
+    mode=RunMode.CONCENTRATION_DRIVEN,
+)
+```
+
+All three adapters share the same contract:
+
+- **Inputs.** Concentrations are supplied as ``Atmospheric Concentrations|<species>``
+  rows in the scenarios DataFrame, overlaid year-by-year (user wins) on a
+  baseline loaded from the canonical RCMIP Phase 3 Zenodo bundle (record
+  `20430630`). Species the user does not supply fall back to the baseline.
+- **`rcmip3_bundle_path`.** Required on the cfg for concentration-driven runs,
+  so the baseline is reproducible across machines.
+- **Outputs.** Each adapter can return ``Atmospheric Concentrations|<species>``
+  for the driven species, so you can check that what comes out matches what was
+  prescribed.
+
+## Capability matrix
+
+| | **MAGICC7** | **FaIR2** | **CICEROSCMPY2** |
+|---|---|---|---|
+| CO2 / CH4 / N2O | ✅ | ✅ | ✅ |
+| F-gases (PFCs / HFCs / SF6 / NF3) | ✅ | ✅ | ✅ |
+| Montreal halocarbons | ✅ | ✅ | ✅ |
+| Requires `rcmip3_bundle_path` | ✅ | ✅ (or conc rows in the scenario) | ✅ |
+| Hybrid baseline + user overlay | ✅ | ✅ | ✅ |
+| Returns concentrations as output | ✅ | ✅ | ✅ |
+| **Mode scoping** | per-scenario, per-gas (WMGHGs) / per-group (F-gas, MHalo) | per-calibration instance (auto-detected) + protocol mixed-mode | adapter-wide global flag |
+| Per-gas mixed mode | WMGHGs independent; F-gas & MHalo all-or-nothing within group | protocol-strict only (CO2 emissions + non-CO2 concentration) | ❌ all-or-nothing |
+| Per-scenario mixing within a batch | ✅ (no batch-consistency constraint) | batch-intersection rule applies in mixed-mode | ❌ global |
+
+**Species coverage is at parity:** all three drive the full WMGHG +
+F-gas + Montreal-halocarbon set.
+
+## Where the adapters differ
+
+The remaining differences are in *how finely modes can be mixed*, and they
+follow each model's native interface rather than any missing feature:
+
+- **MAGICC7** writes a separate config and ``CONC.IN`` files per
+  ``(scenario, model)``, so the decision is genuinely per-scenario and, for
+  CO2/CH4/N2O, per-gas. The F-gas and Montreal-halocarbon groups each share a
+  single ``*_SWITCHFROMCONC2EMIS_YEAR`` flag, so within a group driving is
+  all-or-nothing: the writer fills the whole group from the baseline and warns
+  on any unexpected empty slot. ``HALON1202`` has no RCMIP3 trajectory and
+  falls back to MAGICC's built-in default.
+- **FaIR2** uses a per-species ``input_mode`` that is shared across the whole
+  FaIR *instance* (all scenarios in a calibration). It therefore keeps a
+  batch-consistency rule: its protocol mixed-mode (CO2 emissions-driven,
+  non-CO2 concentration-driven) only engages when every scenario in the batch
+  supplies the concentrations.
+- **CICEROSCMPY2** exposes a single adapter-level ``conc_run`` flag that
+  applies uniformly to every scenario and species in the call.
+
+```{note}
+Because MAGICC7 drives each scenario independently, it has **no**
+batch-consistency requirement: a gas can be concentration-driven in one
+scenario and emissions-driven in another within the same batch. This is
+intentionally more permissive than FaIR2's shared-instance constraint.
+```

docs/source/index.md

@@ -28,6 +28,7 @@
 ```{toctree}
 :caption: Contents
 :maxdepth: 2
+concentration-driven
 idealised-experiments
 notebooks
 development

src/openscm_runner/adapters/fair2_adapter/_output_extractor.py

@@ -158,17 +158,35 @@ def _concentration_unit(species_name: str) -> str:
 # properties_df, available_species) and returns a 1-D numpy array
 # along the time axis (or None if the recipe can't be computed against
 # the species set FaIR actually ran).
+def _forcing_to_numpy(forcing_da):
+    """
+    Materialise FaIR's forcing DataArray as a plain numpy array once.
+
+    Returns ``(forcing_np, species_index)`` where ``forcing_np`` has axis
+    order ``(timebounds, scenario, config, specie)`` and ``species_index``
+    maps species name -> position on the last axis. Extracting/aggregating
+    per (scenario, member) against this numpy array with integer indexing
+    is orders of magnitude faster than repeated xarray ``.sel().isel()``
+    orthogonal indexing (the dominant cost for large ensembles), and is
+    numerically identical (same species order, same summation order).
+    """
+    da = forcing_da.transpose("timebounds", "scenario", "config", "specie")
+    species = list(da["specie"].values)
+    return np.asarray(da.values), {name: i for i, name in enumerate(species)}
+
+
 def _sum_forcing_over(
-    forcing_da, sc_idx: int, member_offset: int, species_to_sum: list[str]
+    forcing_np,
+    species_index: "dict[str, int]",
+    sc_idx: int,
+    member_offset: int,
+    species_to_sum: list[str],
 ) -> np.ndarray:
-    """Helper: sum FaIR's forcing array over a list of species names."""
-    species_present = [s for s in species_to_sum if s in forcing_da["specie"].values]
-    if not species_present:
+    """Sum the forcing array over a list of species names (numpy path)."""
+    idx = [species_index[s] for s in species_to_sum if s in species_index]
+    if not idx:
         return None
-    sliced = forcing_da.sel(specie=species_present).isel(
-        scenario=sc_idx, config=member_offset
-    )
-    return sliced.sum(dim="specie").values
+    return forcing_np[:, sc_idx, member_offset, idx].sum(axis=-1)
 
 
 def _species_by_property(properties_df, predicate) -> list[str]:
@@ -179,20 +197,21 @@ def _species_by_property(properties_df, predicate) -> list[str]:
 
 
 def _build_forcing_aggregations(  # noqa: PLR0912, PLR0915
-    forcing_da, sc_idx: int, member_offset: int, properties_df
+    forcing_np, species_index, species_in_run, sc_idx: int,
+    member_offset: int, properties_df,
 ) -> "dict[str, tuple[np.ndarray, str]]":
     """
     Build the value arrays for the supported forcing aggregations.
 
     Returns a dict keyed by openscm-runner variable name; values are
     ``(values, unit)`` tuples. Aggregations whose species list does
     not intersect FaIR's actual species are omitted, not zero-filled.
+    Operates on the pre-materialised numpy forcing array (see
+    :func:`_forcing_to_numpy`) for speed.
     """
-    species_in_run = list(forcing_da["specie"].values)
-
     def sum_over(species_list: list[str]) -> np.ndarray:
         return _sum_forcing_over(
-            forcing_da, sc_idx, member_offset, species_list
+            forcing_np, species_index, sc_idx, member_offset, species_list
         )
 
     ghgs = _species_by_property(
@@ -230,11 +249,7 @@ def maybe(name: str, values):
         out[name] = (values, forcing_unit)
 
     # Anthropogenic = total forcing minus Solar minus Volcanic
-    total = (
-        forcing_da.isel(scenario=sc_idx, config=member_offset)
-        .sum(dim="specie")
-        .values
-    )
+    total = forcing_np[:, sc_idx, member_offset, :].sum(axis=-1)
     natural_species = [s for s in ("Solar", "Volcanic") if s in species_in_run]
     if natural_species:
         natural = sum_over(natural_species)
@@ -283,13 +298,9 @@ def maybe(name: str, values):
         "Effective Radiative Forcing|Solar": "Solar",
     }
     for openscm_name, fair_specie in single_specie_aliases.items():
-        if fair_specie not in species_in_run:
+        if fair_specie not in species_index:
             continue
-        values = (
-            forcing_da.sel(specie=fair_specie)
-            .isel(scenario=sc_idx, config=member_offset)
-            .values
-        )
+        values = forcing_np[:, sc_idx, member_offset, species_index[fair_specie]]
         maybe(openscm_name, values)
 
     # RCMIP-aligned hierarchical aggregations. Most are aliases for
@@ -367,12 +378,8 @@ def maybe(name: str, values):
     for spec_leaf, spec in (
         ("CO2", "CO2"), ("CH4", "CH4"), ("N2O", "N2O"),
     ):
-        if spec in species_in_run:
-            v = (
-                forcing_da.sel(specie=spec)
-                .isel(scenario=sc_idx, config=member_offset)
-                .values
-            )
+        if spec in species_index:
+            v = forcing_np[:, sc_idx, member_offset, species_index[spec]]
             maybe(f"Effective Radiative Forcing|Anthropogenic|{spec_leaf}", v)
 
     return out
@@ -406,6 +413,9 @@ def extract_outputs(  # noqa: PLR0913, PLR0912, PLR0915
     rows: list = []
 
     species_in_run = list(f.forcing["specie"].values)
+    # Materialise the forcing array once (dominant cost otherwise is
+    # repeated xarray orthogonal indexing per member); aggregate in numpy.
+    forcing_np, forcing_species_index = _forcing_to_numpy(f.forcing)
 
     # Pre-compute aggregations once per (scenario, member) since the
     # recipes share intermediate sums.
@@ -415,7 +425,8 @@ def extract_outputs(  # noqa: PLR0913, PLR0912, PLR0915
 
             aggregations = (
                 _build_forcing_aggregations(
-                    f.forcing, sc_idx, member_offset, properties_df
+                    forcing_np, forcing_species_index, species_in_run,
+                    sc_idx, member_offset, properties_df,
                 )
                 if properties_df is not None
                 else {}