D0FUS (Design 0-dimensional for Fusion Systems) is a Python tokamak systems code for fast 0D/1D design-space exploration, covering plasma physics, superconducting magnet engineering, and techno-economic assessment. It is developed at CEA-IRFM.

The codebase weighs about 18 000 lines of pure Python distributed across five functional modules and a visualisation library, totalling 235 documented functions.

• Pure Python, NumPy/SciPy only. No compilation, no Makefile, runs in Spyder, Jupyter or as a batch job.
• Two fidelity levels (Academic and Refined) selectable through a single GlobalConfig field, sharing the same interface and solver.
• Three execution modes (RUN, SCAN, OPTIMIZATION) auto-detected from the input file syntax.
• Library mode: every physical and engineering function is callable in isolation, outside the solver loop.
• Distinctive engineering features: radially graded TF coils, REBCO Jc scaling laws, three TF mechanical configurations (bucking, wedging, plug).

Three installation options are available. Option A is recommended for most users (no prior Python knowledge required).

This is the simplest path. Spyder comes with its own bundled Python, with no separate Python installation needed.

Go to spyder-ide.org and download the installer for your operating system (Windows, macOS, or Linux). Run the installer and follow the on-screen instructions. Default settings are fine.

If you have git installed, open a terminal (macOS/Linux) or Command Prompt (Windows) and run:

git clone https://github.com/IRFM/D0FUS.git

Otherwise, go to github.com/IRFM/D0FUS, click Code → Download ZIP, and extract the archive somewhere on your computer (e.g. C:\Users\you\D0FUS\ on Windows or ~/D0FUS/ on macOS/Linux).

Launch Spyder from your applications menu (or desktop shortcut). You should see a layout with a code editor on the left, a file explorer on the top right, and an IPython console on the bottom right.

Go to File → Open… and open D0FUS.py from the D0FUS/ folder. Spyder will automatically set the working directory to D0FUS/. You can confirm this by checking the folder path displayed in the toolbar at the top.

In the IPython console (bottom-right panel), type the following and press Enter:

%pip install -r requirements.txt

The %pip magic command installs packages directly into Spyder's internal Python. Since the Spyder standalone bundles its own pip, no prior installation is needed. Since the working directory is already set to D0FUS/, no path specification is needed. Wait for all packages to finish installing.

If %pip fails, try: !pip install -r requirements.txt. If pip itself is missing (unlikely with the standalone), reinstalling Spyder usually fixes it.

Press F5 (or click the green ▶ Run file button). A file-picker dialog will appear. Navigate to D0FUS_INPUTS/ and select 1_run_ITER.txt. D0FUS will run and print results to the IPython console. Figures will open automatically.

Miniforge is a free, community-maintained conda distribution with no commercial restrictions. Recommended if you manage multiple Python projects.

Install Miniforge, then open a Miniforge Prompt (Windows) or terminal (macOS/Linux):

conda create -n d0fus python=3.11
conda activate d0fus
conda install pip
conda install spyder
git clone https://github.com/IRFM/D0FUS.git
cd D0FUS
pip install -r requirements.txt
spyder

Once Spyder is open, follow Steps 5–6 from Option A. Alternatively, run directly from the terminal:

python D0FUS.py
# When prompted, enter the path to the input file:
# D0FUS_INPUTS/1_run_ITER.txt

For use without a GUI, or to import D0FUS as a library in your own scripts:

pip install d0fus

To run the ITER example case from a script:

from D0FUS_EXE import D0FUS_run

results = D0FUS_run.run("D0FUS_INPUTS/1_run_ITER.txt")
print(f"Q = {results['Q']:.1f},  COE = {results['COE']:.0f} EUR/MWh")

Note: input files (D0FUS_INPUTS/) are not included in the PyPI package. Clone the repository separately to access them.

D0FUS/
├── D0FUS_BIB/                          # Core library modules
│   ├── D0FUS_import.py                     # Centralised imports
│   ├── D0FUS_parameterization.py           # Physical constants, GlobalConfig dataclass
│   ├── D0FUS_physical_functions.py         # Plasma physics (profiles, bootstrap, q, RE…)
│   ├── D0FUS_radial_build_functions.py     # Engineering (TF/CS/CICC/quench)
│   ├── D0FUS_cost_functions.py             # Techno-economic models (Sheffield, Whyte)
│   ├── D0FUS_cost_data.py                  # Reference cost data and currency conversions
│   └── D0FUS_figures.py                    # Full figure catalogue (32 plot functions)
│
├── D0FUS_INPUTS/                       # Input parameter files
│   ├── 1_run_ITER.txt                      # ITER Q=10 reference case (RUN mode)
│   ├── 2_scan_ITER.txt                     # 2D scan around ITER point (SCAN mode)
│   └── 3_genetic_ITER.txt                  # Genetic optimisation around ITER (OPTIMIZATION mode)
│
├── D0FUS_OUTPUTS/                      # Generated outputs (auto-created)
│   ├── Run_D0FUS_YYYYMMDD_HHMMSS/         # Single run results + figures
│   ├── Scan_D0FUS_YYYYMMDD_HHMMSS/        # 2D scan maps
│   └── Genetic_D0FUS_YYYYMMDD_HHMMSS/     # Optimisation results
│
├── D0FUS_EXE/                          # Execution modules
│   ├── D0FUS_run.py                        # Single design point
│   ├── D0FUS_scan.py                       # 2D parameter scan
│   └── D0FUS_genetic.py                    # Genetic algorithm optimisation
│
├── D0FUS.py                            # Main entry point
├── requirements.txt
└── README.md

D0FUS detects the execution mode from the input file syntax:

RUN mode evaluates a single tokamak configuration and returns 99 scalar outputs grouped into 8 families (plasma parameters, magnetic fields, power balance, radial build, current profile, techno-economics, RE indicators, diagnostics). Because several physical quantities are mutually coupled, a convergence loop is unavoidable: in pulsed mode the solver reduces to a scalar root-finding on the helium ash fraction $f_\alpha$, solved in 8 to 12 iterations by Brent's method; in steady-state mode the auxiliary power becomes $P_\mathrm{aux} = P_\mathrm{fus}/Q$ and the system becomes two-dimensional, solved by Powell's hybrid method on $(f_\alpha, Q)$. Post-convergence calculations cover TF mechanical sizing, magnetic flux budget and CS sizing, divertor heat loads, L-H power threshold, the techno-economic assessment, and runaway electron diagnostics.

SCAN mode generates 2D maps over two parameters, visualising feasibility regions bounded by stability limits (Greenwald, Troyon, kink) and engineering constraints. Iso-contours of any output quantity registered in OUTPUT_REGISTRY can be overlaid. Grid points are independent and evaluated in parallel across all available CPU cores using joblib with the loky backend.

OPTIMIZATION mode uses a DEAP genetic algorithm to navigate the design space adaptively, concentrating evaluations in the most promising regions. Four fitness objectives are available:

Constraints (Greenwald, normalised beta, kink safety factor, optional capital cost ceiling) are enforced via soft penalties: a design that just barely fails one of the limits stays in the running, which helps the search converge toward the feasible boundary rather than fleeing it.

Most physics and engineering models are available at two fidelity levels, selectable through a single GlobalConfig field. Both share the same interface, the same GlobalConfig object, the same solver infrastructure, and have comparable execution times.

Every function in D0FUS_BIB can be called in isolation, outside the solver loop. This makes the code equally useful as a library for standalone analyses and as an integrated systems code:

from D0FUS_BIB.D0FUS_radial_build_functions import J_non_Cu_REBCO

Jc = J_non_Cu_REBCO(B=20, T=4.2)  # A/m^2, at 20 T and 4.2 K

Comparing two bootstrap models, evaluating a critical current density at a specific field and temperature, or registering a new scaling law all reduce to writing or calling one function with the standard signature.

On a modern 14-core laptop:

These execution times make interactive design-space exploration entirely feasible without access to a computing cluster.

All user-adjustable parameters are gathered into a single typed dataclass GlobalConfig (115 fields organised into 15 categories), each with a physically motivated default value inspired by ITER and EU-DEMO. When an input file is provided, only the specified parameters are overwritten, so a complete tokamak calculation can be set up in a few lines:

R0 = 7
Bmax_TF = 14
Supra_choice = REBCO

All unspecified parameters silently take their default values.

Options for Option_Kappa: Wenninger, Stambaugh, Freidberg, Manual.

*Only relevant when Operation_mode = Pulsed.

Note on CD models: only CD_source = Academic is fully validated. Technology-specific models (LHCD, ECCD, NBCD, Multi) are under active development and should not be used for production runs.

These parameters control the post-convergence RE indicator computation. They do not enter the main physics solver.

Set cost_model = None to skip cost computation entirely.

RUN mode: fixed values only:

R0 = 9
Bmax_TF = 13
Supra_choice = Nb3Sn

See D0FUS_INPUTS/1_run_ITER.txt for a complete example.

SCAN mode: exactly 2 scan parameters with [min, max, n_points]:

R0 = [3, 9, 25]
a  = [1, 3, 25]

See D0FUS_INPUTS/2_scan_ITER.txt for a complete example.

OPTIMIZATION mode: 2+ parameters with [min, max]:

R0     = [3, 9]
a      = [1, 3]
Bmax_TF = [10, 16]
fitness_objective = COE
C_invest_max = 20000

See D0FUS_INPUTS/3_genetic_ITER.txt for a complete example.

D0FUS_OUTPUTS/Run_D0FUS_YYYYMMDD_HHMMSS/
├── input_parameters.txt        # Copy of input configuration
├── output_results.txt          # Complete calculation results (99 scalars, 8 families)
└── figures/                    # 10 run-specific PNG figures (150 dpi)
    ├── 01_cross_section.png
    ├── 02_miller_surfaces.png
    ├── 03_shaping_profiles.png
    ├── 04_kinetic_profiles.png
    ├── 05_q_profile.png
    ├── 06_radiation_profile.png
    ├── 07_TF_side_view.png
    ├── 08_CICC_TF.png
    ├── 09_CS_cross_section.png
    └── 10_CICC_CS.png

Output quantities include:

• Plasma parameters: Ip, n̄_e, T̄, β_N, Q, τ_E, q₉₅, f_bs
• Magnetic fields: B0, B_CS, B_pol
• Power balance: P_fus, P_CD, P_rad, P_sep, P_elec
• Radial build: TF coil thickness, CS inner/outer radius, first-wall area
• Current profile: j_Ohm(ρ), j_CD(ρ), j_bs(ρ), q(ρ) (self-consistent iteration)
• Techno-economics: COE [EUR/MWh], C_invest [M EUR]
• RE indicators: I_RE_seed, I_RE_aval, f_RE/Ip, E_RE_kin

D0FUS_OUTPUTS/Scan_D0FUS_YYYYMMDD_HHMMSS/
├── scan_parameters.txt
└── scan_map_[iso]_[bg].png     # 2D map (300 dpi)

Available scan output quantities (OUTPUT_REGISTRY):

D0FUS_OUTPUTS/Genetic_D0FUS_YYYYMMDD_HHMMSS/
├── optimization_config.txt     # Bounds, algorithm settings
├── optimization_results.txt    # Best solution, COE, C_invest, Hall of Fame
└── convergence_plot.png        # Fitness evolution over generations

plot_run() generates 10 run-specific figures. plot_all() generates the full 32-figure catalogue, including:

D0FUS supports two plasma geometry models, selected via Plasma_geometry:

The Miller model computes V′(ρ) numerically from the Jacobian of the (R, Z) flux-surface coordinates, with κ(0) = 1 and δ(0) = 0 enforced on-axis (Ball & Parra 2015). Volume integrals (P_fus, ⟨nT⟩, W_th, P_rad, bootstrap current) are weighted by V′(ρ)/V when in Refined mode.

Two q(ρ) models are available:

• Analytical (f_q_profile): parameterised as j(ρ) ∝ (1 − ρ²)^α_J, with α_J derived automatically from the resistive diffusion timescale τ_R.
• Self-consistent (f_q_profile_selfconsistent): Picard iteration on j_total(ρ) = j_Ohm(ρ) + j_CD(ρ) + j_bs(ρ). Produces a decomposed current profile visualisation.

The q₉₅ boundary value is computed via the Option_q95 selector:

Four models are supported via Bootstrap_choice:

Scaling_Law options: IPB98(y,2), ITPA20, ITPA20-IL, DS03, L-mode, L-mode OK, ITER89-P.

Line radiation from seeded or intrinsic impurities is modelled via the Mavrin (2018) cooling function tabulation (validated against ADAS). Supported species: W, Ar, Ne, C, N, Kr. Core radiation (ρ < rho_rad_core) is subtracted from the power balance; edge radiation reduces the divertor load P_sep.

Post-convergence RE indicators are computed by compute_RE_indicators() using:

• Hot-tail seed: Smith & Verwichte (2008), exponentially sensitive to τ_TQ
• Avalanche amplification: Breizman et al. (2019)
• Coulomb logarithm: single relativistic lnΛ ≈ 15.3 via _coulomb_log_relativistic()
• Literature benchmarks: Martín-Solís et al. (2017)

Outputs include I_RE_seed, I_RE_aval, f_RE/Ip, and kinetic energy E_RE_kin. These are indicators for comparative design ranking, not quantitative predictions. Configurable via tau_TQ, Te_final_eV, and pellet_dilution.

The TF coil inboard leg is shaped as a Princeton-D contour (Gralnick & Tenney 1976). Conductor sizing follows a three-level helium fraction hierarchy distinguishing the cooling pipe (f_He_pipe), interstitial void (f_void, LTS only), and active SC fraction. Three structural models are available via Radial_build_model: a simplified analytical model (Academic), the D0FUS stress model (Refined), and the multi-layer CIRCE0D solver (CIRCE). Three TF mechanical configurations are supported via Choice_Buck_Wedg (bucking, wedging, plug), and the TF winding pack can be radially graded with $\alpha(R)$ profiles obtained by inward integration with Picard iteration.

Critical current density scaling laws are implemented for NbTi, Nb₃Sn (ITER strand parameterisation), and REBCO (Senatore 2024 / Fujikura 2019 datasets). Quench protection is sized via the Maddock hot-spot criterion, with dump time computed from stored magnetic energy and conductor current.

Capital investment and COE are computed using the Sheffield & Milora (2016) volume-based cost scaling (2010 USD, converted to 2025 EUR). Component costs cover the SC coil set, blanket, shield, auxiliary heating, heat transfer system, balance of plant, buildings, and annual O&M. A simplified surface-proportional model (Whyte 2024) is also available for cross-checks. The cost calculation is performed after convergence and does not feed back into the physics.

Contributions are welcome. Please contact:

• Email: [email protected]

This project is licensed under the CeCILL-C License, a French free software license compatible with the GNU LGPL.

See the LICENSE file for details.

© 2025 CEA/IRFM

(back to top)