Fix krypton calculation

.github/workflows/workflow_actions.yml

@@ -10,14 +10,17 @@ on:
 
 env:
   FORCE_JAVASCRIPT_ACTIONS_TO_NODE24: true
+  # Use non-interactive matplotlib backend, to prevent the workflow from
+  # trying to make interactive windows.
+  MPLBACKEND: Agg
 
 concurrency:
   group: ${{ github.workflow }}-${{ github.ref }}
   cancel-in-progress: true
 
 jobs:
   # --- JOB 1: RUN TESTS ---
-  test:
+  pytest:
     name: Test (Python ${{ matrix.python-version }})
     runs-on: ubuntu-latest
     strategy:
@@ -41,26 +44,48 @@ jobs:
         run: poetry install --with dev
 
       - name: Run Pytest
-        env:
-          MPLBACKEND: Agg
         run: poetry run pytest
 
       - name: Smoke test package installation
-        env:
-          MPLBACKEND: Agg
         run: |
           poetry build
           python -m venv test_env
           source test_env/bin/activate
           pip install dist/*.whl
           radas_config -o ./new_config.yaml
           radas -s hydrogen -c ./new_config.yaml
+
+  # --- JOB 2: RUN RADAS COMMAND, CHECK IT WORKS WITH DEFAULT CONFIG ---
+  test_radas:
+    name: Test radas command
+    runs-on: ubuntu-latest
+    if: (github.event_name == 'release' || github.event_name == 'workflow_dispatch')
+    permissions:
+      contents: write
+
+    steps:
+      - uses: actions/checkout@v5
+
+      - name: Install Poetry
+        run: pipx install "poetry>=2,<3"
+
+      - name: Set up Python
+        uses: actions/setup-python@v6
+        with:
+          python-version: '3.12'
+          cache: 'poetry'
+
+      - name: Install dependencies
+        run: poetry install
 
-  # --- JOB 2: BUILD RELEASE ARTIFACTS ---
+      - name: Run radas command
+        run: poetry run radas -vvv
+
+  # --- JOB 3: BUILD RELEASE ARTIFACTS ---
   build_release:
     name: Build Release
-    needs: test
-    if: startsWith(github.ref, 'refs/tags')
+    needs: [pytest, test_radas]
+    if: github.event_name == 'release'
     runs-on: ubuntu-latest
     steps:
       - uses: actions/checkout@v5
@@ -78,11 +103,11 @@ jobs:
           path: dist/
           retention-days: 1
 
-  # --- JOB 3: PUBLISH TO PYPI ---
+  # --- JOB 4: PUBLISH TO PYPI ---
   publish:
     name: Publish to PyPI
     needs: build_release
-    if: startsWith(github.ref, 'refs/tags')
+    if: github.event_name == 'release'
     runs-on: ubuntu-latest
     environment:
       name: pypi-publish

pyproject.toml

@@ -5,7 +5,7 @@ name = "radas"
 # used a date-based format (e.g., 2025.6.2), package managers will treat
 # Semantic Versioning (i.e. 1.0.0) as an "older" release. The epoch ensures that
 # the new versioning takes precedence over all legacy YYYY.MM.version releases.
-version = "1!1.0.0"
+version = "1!1.0.1"
 description = "Plasma radiated power calculated using OpenADAS"
 authors = ["Commonwealth Fusion Systems"]
 readme = "README.md"

radas/adas_interface/download_adas_datasets.py

@@ -32,7 +32,7 @@ def download_species_data(
         year_key = f"{year}"[-2:]
         dataset_prefix = dataset_config["prefix"].lower()
 
-        output_filename = data_file_dir / f"{species_name}_{dataset_type}.dat"
+        output_filename = data_file_dir / f"{species_name}_{dataset_type}_{year_key}.dat"
         query_path = f"{url_base}/download/{reader_class}/{dataset_prefix}{year_key}/{dataset_prefix}{year_key}_{species_key}.dat"
 
         if not output_filename.exists():

radas/adas_interface/read_adf11_file.py

@@ -14,15 +14,15 @@ def load_library(library_name: str, filepath: Path):
 
 
 def read_adf11_file(
-    data_file_dir, species_name, dataset_type
+    data_file_dir, species_name, year, dataset_type
 ) -> dict:
     """Open and read an ADF11 OpenADAS file.
 
     Uses the format specification from https://www.adas.ac.uk/man/appxa-11.pdf
     """
 
-
-    filename = data_file_dir / f"{species_name}_{dataset_type}.dat"
+    year_key = f"{year}"[-2:]
+    filename = data_file_dir / f"{species_name}_{dataset_type}_{year_key}.dat"
     if not filename.exists():
         raise FileNotFoundError(f"{filename} does not exist.")
 

radas/cli.py

@@ -4,6 +4,7 @@
 from pathlib import Path
 from functools import partial
 from typing import Optional
+import contextlib
 
 from .shared import open_yaml_file, default_config_file
 from .adas_interface.download_adas_datasets import download_species_data
@@ -14,7 +15,6 @@
 from .time_evolution import calculate_time_evolution
 from .unit_handling import convert_units, ureg
 from .mavrin_reference import compare_radas_to_mavrin
-from .interpolate_rates import interpolate_dataset
 
 
 @click.command()
@@ -69,14 +69,22 @@ def run_radas_cli(
         verbose=verbose,
         debug=debug,
     )
-    try:
-        from ipdb import launch_ipdb_on_exception
-
-        with launch_ipdb_on_exception():
+
+    if debug:
+        with _post_mortem_debugger():
             run_radas(**kwargs)
-    except ModuleNotFoundError:
+    else:
         run_radas(**kwargs)
 
+def _post_mortem_debugger():
+    """Context manager that drops into ipdb on unhandled exceptions, or a no-op if ipdb is not installed."""
+    try:
+        from ipdb import launch_ipdb_on_exception
+    except ModuleNotFoundError:
+        print("Warning: --debug set but ipdb is not installed; "
+                "install with `pip install ipdb` for post-mortem debugging.")
+        return contextlib.nullcontext()
+    return launch_ipdb_on_exception()
 
 def run_radas(
     directory: Path,
@@ -126,22 +134,9 @@ def run_radas(
                 (species_name in species) or (species == ("all",))
             ):
                 datasets[species_name] = read_rate_coeff(
-                    data_file_dir, species_name, configuration
+                    data_file_dir, species_name, configuration, verbose=verbose,
                 )
 
-        if ("electron_density_resolution" in configuration["globals"]) or ("electron_temp_resolution") in configuration["globals"]:
-            if verbose:
-                print("Interpolating rate coefficients")
-
-            new_datasets = dict()
-            for species_name, dataset in datasets.items():
-                electron_density_resolution = configuration["globals"].get("electron_density_resolution", dataset.sizes["dim_electron_density"])
-                electron_temp_resolution = configuration["globals"].get("electron_temp_resolution", dataset.sizes["dim_electron_temp"])
-                new_datasets[species_name] = interpolate_dataset(dataset,
-                                                                 electron_density_resolution = electron_density_resolution,
-                                                                 electron_temp_resolution = electron_temp_resolution)
-            datasets = new_datasets
-
         output_dir.mkdir(exist_ok=True, parents=True)
         if not debug:
             with mp.Pool() as pool:
@@ -150,11 +145,15 @@ def run_radas(
                         species_name: datasets[species_name] for species_name in species
                     }
 
+                # Sort by atomic number, to process heavier elements first since they take longer
+                # N.b. there will be a race condition, so it might not appear these start first
+                sorted_datasets = dict(sorted(datasets.items(), key=lambda item: item[1].atomic_number, reverse=True))
+
                 pool.map(
                     partial(
                         run_radas_computation, output_dir=output_dir, verbose=verbose
                     ),
-                    [(ds) for ds in datasets.values()],
+                    [(ds) for ds in sorted_datasets.values()],
                 )
         else:
             for ds in datasets.values():
@@ -199,6 +198,9 @@ def run_radas_computation(dataset: xr.Dataset, output_dir: Path, verbose: int):
     output_dir.mkdir(exist_ok=True)
     dataset.pint.dequantify().to_netcdf(output_dir / f"{dataset.species_name}.nc")
 
+    if verbose:
+        print(f"Finished computation for {dataset.species_name}")
+
 
 @click.command()
 @click.option(

radas/config.yaml

@@ -10,14 +10,13 @@ globals:
 
   # electron density (ne) * residence time (tau) (in m^-3 s)
   ne_tau:
-    value: [0.5e+17]
+    value: [0.5e+16, 0.5e+17, 0.5e+18]
     units: "m^-3 s"
 
   # Number of log-spaced points to use for electron density and electron temp.
   # Interpolation is used to map from raw data to points.
-  # To disable interpolation, do not provide these inputs.
-  electron_density_resolution: 50
-  electron_temp_resolution: 100
+  electron_density_resolution: 20
+  electron_temp_resolution: 80
 
 data_file_config:
   adf11:
@@ -459,4 +458,4 @@ species:
     atomic_number: 81
   lead:
     atomic_symbol: "Pb"
-    atomic_number: 82
+    atomic_number: 82

radas/interpolate_rates.py

@@ -1,78 +1,86 @@
-"""Routines to interpolate a dataset of rate coefficients to higher resolution."""
+"""Routines for log-log interpolation of rate coefficients with boundary clipping."""
 import xarray as xr
 import numpy as np
 from scipy.interpolate import RectBivariateSpline
 from numpy.typing import NDArray
+import warnings
 
-def interpolate_array(array: xr.DataArray, new_electron_density: NDArray[np.floating], new_electron_temp: NDArray[np.floating]) -> xr.DataArray:
-    """Interpolate array onto new values for the electron density and electron temp.
+def is_significantly_below(requested, limit):
+    return requested < limit and not np.isclose(requested, limit)
+
+def is_significantly_above(requested, limit):
+    return requested > limit and not np.isclose(requested, limit)
+
+def interpolate_array(
+    array: xr.DataArray, 
+    new_electron_density: NDArray[np.floating], 
+    new_electron_temp: NDArray[np.floating]
+) -> xr.DataArray:
+    """
+    Interpolate rate coefficients onto a new density/temperature grid in log-log space.
 
-    The interpolation is performed for logarithmic values.
+    Uses nearest-neighbor extrapolation by clipping out-of-bounds coordinates to 
+    the original grid edges.
     """
     units = array.pint.units
     array = array.pint.dequantify().squeeze()
 
+    # Handle zero-value edge cases (log of zero is undefined)
     if np.allclose(array, 0.0, atol=0.0, rtol=1e-6):
-        # If all values of the array are zero, return a zero array.
-        return xr.DataArray(np.zeros((np.size(new_electron_temp), np.size(new_electron_density))),
+        return xr.DataArray(
+            np.zeros((np.size(new_electron_temp), np.size(new_electron_density))),
             coords=dict(dim_electron_temp=new_electron_temp, dim_electron_density=new_electron_density)
         ) * units
-    elif np.any(array <= 0.0):
-        # If only some of the values of the array are zero, raise an error.
-        raise NotImplementedError("Cannot handle zero-valued entries in non-zero rate coefficients.")
 
+    if np.any(array <= 0.0):
+        raise NotImplementedError("Cannot log-interpolate rate coefficients containing zeros.")
+
+    # Check if extrapolation is needed and raise a warning if this is the case.
+    out_of_bounds_msg = []
+
+    req_dens_min, req_dens_max = new_electron_density.min(), new_electron_density.max()
+    grid_dens_min, grid_dens_max = array.dim_electron_density.min(), array.dim_electron_density.max()
+
+    if is_significantly_below(req_dens_min, grid_dens_min) or is_significantly_above(req_dens_max, grid_dens_max):
+        out_of_bounds_msg.append(
+            f"Density requested [{req_dens_min:.2e}, {req_dens_max:.2e}] "
+            f"exceeds grid [{grid_dens_min:.2e}, {grid_dens_max:.2e}]."
+        )
+
+    # Check Temperature Bounds
+    req_temp_min, req_temp_max = new_electron_temp.min(), new_electron_temp.max()
+    grid_temp_min, grid_temp_max = array.dim_electron_temp.min(), array.dim_electron_temp.max()
+
+    if is_significantly_below(req_temp_min, grid_temp_min) or is_significantly_above(req_temp_max, grid_temp_max):
+        out_of_bounds_msg.append(
+            f"Temperature requested [{req_temp_min:.2e}, {req_temp_max:.2e}] "
+            f"exceeds grid [{grid_temp_min:.2e}, {grid_temp_max:.2e}]."
+        )
+
+    if out_of_bounds_msg:
+        full_msg = "Nearest-neighbour extrapolation used for off-grid values: " + " ".join(out_of_bounds_msg)
+        warnings.warn(full_msg, RuntimeWarning)
+
+    # ------------------------------------
+
+    # Prepare original grid and data in log10 space
     x = np.log10(array.dim_electron_density)
     y = np.log10(array.dim_electron_temp)
     z = np.log10(array.transpose("dim_electron_density", "dim_electron_temp").pint.magnitude)
 
+    # Transform target coordinates to log10
     x_interp = np.log10(new_electron_density)
     y_interp = np.log10(new_electron_temp)
-    z_interp = np.power(10, RectBivariateSpline(x, y, z)(x_interp, y_interp, grid=True).T)
-
-    return xr.DataArray(z_interp,
-        coords=dict(dim_electron_temp=new_electron_temp, dim_electron_density=new_electron_density)
-    ) * units
 
-def interpolate_dataset(dataset: xr.Dataset, electron_density_resolution: int, electron_temp_resolution: int) -> xr.Dataset:
-    """Interpolate all rate coefficients in a dataset."""
-    new_electron_density = np.logspace(
-        np.log10(dataset["dim_electron_density"].min().item()),
-        np.log10(dataset["dim_electron_density"].max().item()),
-        num = electron_density_resolution
-    )
+    # Force nearest-neighbor extrapolation by clipping points to the grid domain
+    x_clipped = np.clip(x_interp, x.min().item(), x.max().item())
+    y_clipped = np.clip(y_interp, y.min().item(), y.max().item())
 
-    new_electron_temp = np.logspace(
-        np.log10(dataset["dim_electron_temp"].min().item()),
-        np.log10(dataset["dim_electron_temp"].max().item()),
-        num = electron_temp_resolution
-    )
+    # Perform spline interpolation and revert from log space
+    z_interp_log = RectBivariateSpline(x, y, z)(x_clipped, y_clipped, grid=True)
+    z_interp = np.power(10, z_interp_log.T)
 
-    new_dataset = xr.Dataset().assign_attrs(dataset.attrs)
-
-    for key, array in dataset.items():
-        if key in [
-            "electron_density",
-            "electron_temp",
-        ]:
-            # Don't copy in the coordinate arrays which we'll interpolate
-            continue
-        elif key in [
-            "ne_tau"
-        ]:
-            # Directly copy in the coordinate arrays which we'll leave unchanged
-            new_dataset[key] = array
-        elif array.ndim == 0:
-            # Directly copy in scalar arrays
-            new_dataset[key] = array
-        elif (("dim_electron_density" in array.coords)
-              and ("dim_electron_temp" in array.coords)
-              and ("dim_charge_state" in array.coords)):
-            # For each charge state, interpolate the rate coefficient
-            new_dataset[key] = array.groupby("dim_charge_state").map(interpolate_array, args=(new_electron_density, new_electron_temp))
-        else:
-            raise NotImplementedError(f"Could not process array '{key}' with coords {array.coords}")
-
-    new_dataset["electron_density"] = xr.DataArray(new_electron_density, dims="dim_electron_density") * dataset["electron_density"].pint.units
-    new_dataset["electron_temp"] = xr.DataArray(new_electron_temp, dims="dim_electron_temp") * dataset["electron_temp"].pint.units
-
-    return new_dataset
+    return xr.DataArray(
+        z_interp,
+        coords=dict(dim_electron_temp=new_electron_temp, dim_electron_density=new_electron_density)
+    ) * units