PlasmaKit

Fusion nuclear engineering, computed with receipts.

PlasmaKit connects plasma conditions to fusion reaction rates, spatially resolved neutron sources, blanket neutronics (TBR, wall load, damage), uncertainty quantification, and the tritium fuel cycle — one Python API spanning the whole plasma → neutron → blanket → tritium chain. Every physics model is cited, every calculation is validated against a published or analytic reference you can re-run yourself, and every result carries a provenance record.

38 validated benchmarks · First-class provenance · Plasma → tritium in one chain

Quick Install

pip install plasmakit
# or
uv add plasmakit

Requires Python ≥ 3.10. Runtime dependencies: NumPy and SciPy. OpenMC is optional (transport only) and installed separately from conda-forge.

30-Second Example

import plasmakit as pk

# A burning plasma
plasma = pk.PlasmaState(
    ion_temperature=15.0,  # keV
    ion_density=1.0e20,  # m^-3
    fuel={"D": 0.5, "T": 0.5},
)

source = pk.NeutronSource(plasma)
source.rate_density()  # 6.87e17 neutrons / m^3 / s
source.mean_energy()  # 14,053 keV (Brysk thermal shift included)
source.power_density().neutron  # W/m^3 carried by neutrons

# Prove the physics is right
pk.validate()  # 38 benchmark cases vs published references

Table of Contents

• Overview
• Validated Accuracy
• Determinism and Provenance
• Plasma States and Reactions
• Spatial Neutron Sources
• Blanket Neutronics
• Tritium Fuel Cycle
• Uncertainty Quantification and Optimization
• End-to-End Worked Example
• API Reference
• Units
• Building from Source
• References
• License
• Changelog

Overview

PlasmaKit is the interface layer of fusion nuclear engineering: it does not replace plasma solvers, Monte Carlo transport codes, or systems codes — it makes them composable. A PlasmaState becomes a spatially resolved neutron source; the source drives an OpenMC blanket calculation; the blanket result feeds a tritium fuel-cycle model; and the uncertainty toolbox wraps any step of that chain.

Key Features

• Fusion reactions — D-T, D-D (both branches), D-³He cross sections and Maxwellian reactivities via the Bosch–Hale (1992) parameterization; Q values and product energies derived from CODATA masses, never hand-copied.
• Neutron spectra — thermally broadened Gaussian (Brysk 1973) spectra with mean shift, FWHM, pdf, and seeded sampling.
• Spatial sources — radial profiles on parameterized Miller flux surfaces (elongation, triangularity, Shafranov shift) or user R-Z fields; exporters to OpenMC ring sources, xarray, and dependency-free VTK.
• Blanket neutronics — layered torus blankets with cited materials; OpenMC coupling returning TBR, neutron wall load, per-layer heating and tritium production, and NRT DPA, each with Monte Carlo uncertainty.
• Tritium fuel cycle — the Abdou compartment model solved exactly by matrix exponential: self-sufficiency, doubling time, required startup inventory.
• UQ and optimization — Sobol QMC propagation, Saltelli sensitivity indices, Metropolis–Hastings estimation, Gaussian-process surrogates, and differential-evolution optimization, all under one fn(**params) convention.
• Validation as a feature — pk.validate() re-runs the entire benchmark registry against published references, and the same registry gates CI.

Validated Accuracy

The exact table pk.validate() prints — every physics model checked against a published value, an analytic closed form, or a hand calculation, in one seeded, deterministic run:

References include Bosch & Hale (Nucl. Fusion 32, 1992), Brysk (Plasma Phys. 15, 1973), Miller et al. (Phys. Plasmas 5, 1998), Norgett–Robinson–Torrens (Nucl. Eng. Des. 33, 1975), Abdou et al. (Fusion Technology 9, 1986), Saltelli et al. (Comput. Phys. Commun. 181, 2010), and analytic closed forms. pk.validate().to_json() produces the same table as a citable machine-readable record.

Determinism and Provenance

• Every stochastic entry point takes a seed — QMC sampling, Sobol indices, MCMC chains, GP training restarts, differential evolution, and spectrum sampling are all bit-reproducible for a fixed seed.
• Every derived result knows how it was made. result.provenance carries the package version, the physics-model identifiers used (e.g. bosch-hale-1992, brysk-1973, nrt-1975), their full citations, the physical inputs, and — for transport — the nuclear-data library, particle count, and RNG seed. provenance.to_json() is a reproducible computational record; blanket results chain the source's provenance so the trail runs unbroken from the Bosch–Hale coefficients to the TBR.

Plasma States and Reactions

The registry derives all energetics from CODATA 2018 nuclide masses at import time — Q values and two-body product energies are computed, not transcribed:

import numpy as np
import plasmakit as pk

# Everything is vectorized: scalar in → float out, array in → array out
temperatures = np.linspace(1.0, 100.0, 500)  # keV
sigma_v = pk.maxwellian_reactivity("DT", temperatures)  # m^3/s, shape (500,)
sigma = pk.cross_section("DT", 64.0)  # m^2 at 64 keV (CM)

spec = pk.neutron_spectrum("DT", ion_temperature=10.0)
spec.mean_energy, spec.fwhm  # keV; FWHM ≈ 177·√T
spec.sample(10_000, rng=np.random.default_rng(0))  # seeded draws

Out-of-range inputs never extrapolate silently — a ValidityRangeWarning names the model and its fitted range.

Spatial Neutron Sources

Real plasmas are not rings. Build a cell-resolved neutron source from radial profiles on shaped flux surfaces, or from 2-D R-Z fields out of an equilibrium code:

profiles = pk.PlasmaProfiles(
    ion_temperature=pk.RadialProfile.parabolic(20.0, 1.0),  # keV, core → edge
    ion_density=pk.RadialProfile.parabolic(1.0e20, 1.0e18),  # m^-3
)
geometry = pk.TokamakGeometry(
    major_radius=6.0, minor_radius=2.0, elongation=1.7, triangularity=0.33
)

source = pk.SpatialNeutronSource.from_profiles(profiles, geometry)
source.total_rate  # neutrons / s
source.total_fusion_power  # W
source.emissivity  # per-cell neutron emissivity field, m^-3 s^-1

source.to_openmc()  # weighted openmc.IndependentSource rings
source.to_xarray()  # labeled Dataset (pip install plasmakit[xarray])
source.to_vtk("source.vtk")  # ParaView/VisIt, no VTK dependency needed

The OpenMC export builds one axisymmetric ring per cell and reaction, each carrying the local Brysk spectrum — the hot core emits harder, wider neutrons than the edge, so wall-load calculations see a realistic source. SpatialNeutronSource.from_rz(r_edges, z_edges, T_field, n_field) accepts equilibrium-code output directly. Cell volumes are validated against analytic torus integrals, and flat profiles reproduce the 0-D results exactly.

Blanket Neutronics

Layered torus blankets with a built-in, literature-cited materials registry:

from plasmakit.materials import beryllium, eurofer97, li4sio4, tungsten

blanket = pk.Blanket(
    layers=(
        pk.Layer("armor", tungsten(), 0.002),
        pk.Layer("first_wall", eurofer97(), 0.02),
        pk.Layer("multiplier", beryllium(), 0.05),
        pk.Layer("breeder", li4sio4(li6_enrichment=0.60), 0.50),
        pk.Layer("shield", eurofer97(), 0.10),
    ),
    major_radius=9.0,
    first_wall_radius=2.9,
)

result = blanket.run_neutronics(plasma, source_rate=7.1e20)  # needs openmc + data

result.tbr  # TallyValue: value ± Monte Carlo std
result.neutron_wall_load  # MW/m^2
result.energy_deposition  # W per layer
result.tritium_production  # atoms/s per layer
result.dpa_per_fpy  # first-wall displacement dose per full-power year

Layers are concentric circular torus shells (Blanket.from_geometry conservatively encloses a shaped Miller LCFS). DPA uses the NRT model (E_d = 40 eV, configurable). OpenMC is not on PyPI — install from conda-forge and point OPENMC_CROSS_SECTIONS at a data library; everything except run_neutronics/to_openmc works without it.

Tritium Fuel Cycle

TBR > 1 is necessary but not sufficient — the plant must keep its inventory alive against decay, recirculation losses, and the huge unburned-fuel loop. TritiumCycle is the Abdou four-compartment linear model (blanket, exhaust, processing, storage) solved exactly by matrix exponential — no ODE tolerances:

cycle = pk.TritiumCycle.from_blanket_result(
    result,
    fractional_burnup=0.02,
    startup_inventory=5.0,  # kg
)

cycle.self_sufficient  # storage accumulating?
cycle.accumulation_rate()  # kg/day at steady state
cycle.doubling_time()  # days to bank a 2nd startup inventory
cycle.required_startup_inventory(days=3650.0)  # kg to survive a 10-year horizon
cycle.simulate(days=365 * 5).inventory("storage")  # kg vs time

TritiumCycle.from_fusion_power(2000.0, tbr=1.15, ...) starts from plant power instead. The model reproduces the classic Abdou result: at low fractional burnup, recirculation losses can defeat even TBR = 1.15.

Uncertainty Quantification and Optimization

One convention everywhere: models are callables invoked as fn(**params), so the same function drives propagation, sensitivity, fitting, surrogates, and optimization.

def fusion_power(ion_temperature, ion_density):
    plasma = pk.PlasmaState(ion_temperature, ion_density, fuel={"D": 0.5, "T": 0.5})
    return pk.fusion_power_density(plasma)


params = {
    "ion_temperature": pk.Distribution.lognormal(mean=15.0, std=2.0),
    "ion_density": pk.Distribution.normal(1.0e20, 5.0e18),
}

result = pk.propagate(fusion_power, params, n_samples=10_000, vectorized=True)
result.mean, result.std, result.percentile(95)

indices = pk.sobol_indices(fusion_power, params)  # which input drives the variance?
indices.first_order["ion_temperature"], indices.total_order["ion_temperature"]

posterior = pk.fit(fusion_power, priors=params, observed=8.0e5, noise_std=5.0e4)
posterior.mean("ion_temperature"), posterior.map_estimate

surrogate = pk.Surrogate.from_function(fusion_power, params, n_train=64)
best = pk.optimize(lambda x, y: (x - 3) ** 2 + y**2, {"x": (0, 5), "y": (-2, 2)})

Through transport, pk.propagate_transport folds the per-run Monte Carlo tally variance into the total: std² = parametric spread + mean tally noise. For expensive objectives, pk.optimize_surrogate runs a one-shot QMC design → GP → optimize → verify loop and reports both the surrogate prediction and the verified true value.

End-to-End Worked Example

The full chain on a DEMO-scale machine (measured with OpenMC 0.15.3 on ENDF/B-VII.1, 10k particles — reproduce with the blanket above):

plasma = pk.PlasmaState(ion_temperature=15.0, ion_density=1.0e20, fuel={"D": 0.5, "T": 0.5})
result = blanket.run_neutronics(plasma, particles=10_000, source_rate=7.1e20)  # ~2 GW fusion
# result.tbr                      -> 1.259 ± 0.006
# result.neutron_wall_load        -> 0.245 MW/m^2 (analytic uncollided: 0.219)

cycle = pk.TritiumCycle.from_blanket_result(
    result,
    fractional_burnup=0.02,
    startup_inventory=5.0,
    extraction_efficiency=0.98,
    processing_loss=1e-4,
)
# cycle.self_sufficient              -> True (accumulating 66.9 g/day)
# cycle.required_startup_inventory() -> 20.3 kg for a 10-year horizon
# cycle.doubling_time()              -> 404 days

Every number above sits inside published DEMO-study ranges, and result.provenance records the nuclear-data library, seeds, and the full model chain that produced it.

API Reference

All names below are importable from the top-level plasmakit namespace.

Core physics

Spatial sources and geometry

Blanket and transport

Tritium fuel cycle

Uncertainty and optimization

Everything else

validate() / BenchmarkReport, Provenance, and TallyValue round out the public API. Package errors derive from plasmakit.errors.PlasmakitError.

Units

One convention, everywhere. No unit-wrapper objects:

Building from Source

git clone git@github.com:alphabench/plasmakit.git
cd plasmakit
uv sync          # runtime + dev dependencies into .venv

Verification Test

uv run ruff check && uv run ruff format --check          # style
uv run mypy                                              # strict typing
uv run pytest                                            # 300+ tests
uv run python -c "import plasmakit as pk; pk.validate()" # 38/38 PASS

For the OpenMC-coupled transport tests (marked transport), see CONTRIBUTING.md — conda-forge OpenMC plus a targeted ~130 MB ENDF/B-VII.1 nuclide library.

References

• H.-S. Bosch and G.M. Hale, Nuclear Fusion 32 (1992) 611 — cross sections and reactivities.
• H. Brysk, Plasma Physics 15 (1973) 611 — neutron spectra.
• R.L. Miller et al., Physics of Plasmas 5 (1998) 973 — flux-surface geometry.
• M.J. Norgett, M.T. Robinson and I.M. Torrens, Nucl. Eng. Des. 33 (1975) 50 — NRT displacement damage.
• P.K. Romano et al., Annals of Nuclear Energy 82 (2015) 90 — OpenMC.
• M.A. Abdou et al., Fusion Technology 9 (1986) 250 and Nuclear Fusion 61 (2021) 013001 — tritium fuel cycle and self-sufficiency.
• A. Saltelli et al., Comput. Phys. Commun. 181 (2010) 259 — Sobol sensitivity estimators.
• C.E. Rasmussen and C.K.I. Williams, Gaussian Processes for Machine Learning, MIT Press (2006).
• R. Storn and K. Price, J. Global Optimization 11 (1997) 341 — differential evolution.

License

MIT — see LICENSE.

Changelog

Canonical history lives in CHANGELOG.md.

v0.1.0 — 2026-08-26

First public release, spanning the full plasma → tritium chain:

• 0-D physics core: PlasmaState, Bosch–Hale cross sections and reactivities, CODATA-derived reaction kinematics, Brysk neutron spectra, NeutronSource with first-class provenance.
• Spatially resolved sources: radial profiles on Miller flux surfaces or R-Z fields, with OpenMC / xarray / dependency-free VTK exporters.
• Blanket neutronics via OpenMC: cited materials registry, layered torus blankets, TBR / wall load / heating / tritium production / NRT DPA with Monte Carlo uncertainties and nuclear-data provenance.
• Uncertainty toolbox: Distribution specs, Sobol-QMC propagate, Saltelli sobol_indices, transport-aware propagation, Metropolis–Hastings fit, Gaussian-process surrogates, and differential-evolution optimize / optimize_surrogate.
• Tritium fuel cycle: exact matrix-exponential Abdou compartment model with self-sufficiency, doubling-time, and startup-inventory analyses.
• 38-case benchmark registry behind plasmakit.validate(), shared with CI.