Your First Simulation

This hands-on tutorial walks you through your first SCPN-Fusion-Core simulation, from installation to plotting results. By the end you will have solved a Grad-Shafranov equilibrium, computed transport profiles, and run a basic stability check.

Prerequisite: Fusion Engineering 101 (for physics context).

Installation

pip install scpn-fusion

For the Rust-accelerated backend (10–50x faster):

pip install scpn-fusion-rs

Verify the installation:

import scpn_fusion
print(scpn_fusion.__version__)  # should print 3.9.3 or later

Step 1: Solve an Equilibrium

Create an ITER-like equilibrium with 65x65 grid resolution:

from scpn_fusion.core import FusionKernel

kernel = FusionKernel(
    R0=6.2,       # major radius [m]
    a=2.0,        # minor radius [m]
    B0=5.3,       # toroidal field [T]
    Ip=15e6,      # plasma current [A]
    kappa=1.7,    # elongation
    delta=0.33,   # triangularity
    grid_size=65,
)

result = kernel.solve()
print(f"Converged in {result.iterations} iterations")
print(f"Beta_p = {result.beta_p:.3f}")
print(f"li = {result.li:.3f}")

The solver computes the poloidal flux \(\psi(R, Z)\) and derives all equilibrium quantities: flux surfaces, q-profile, pressure, and current density.

Step 2: Visualise Flux Surfaces

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(5, 8))
ax.contour(kernel.R_grid, kernel.Z_grid, result.psi, levels=20)
ax.set_xlabel("R [m]")
ax.set_ylabel("Z [m]")
ax.set_aspect("equal")
ax.set_title("Poloidal Flux Surfaces")
plt.tight_layout()
plt.show()

You should see nested elliptical contours centred on the magnetic axis, with a separatrix at the outermost closed surface.

Step 3: Inspect the q-Profile

fig, ax = plt.subplots()
ax.plot(result.rho, result.q_profile)
ax.set_xlabel(r"$\rho$ (normalised radius)")
ax.set_ylabel("q")
ax.set_title("Safety Factor Profile")
ax.axhline(1.0, color="red", linestyle="--", label="q = 1 (sawtooth)")
ax.axhline(1.5, color="orange", linestyle="--", label="q = 3/2 (NTM)")
ax.axhline(2.0, color="green", linestyle="--", label="q = 2 (NTM)")
ax.legend()
plt.show()

The q-profile typically starts near 1 on axis and rises to 3–5 at the edge. Rational surfaces where \(q = m/n\) are marked: these are the locations of MHD instabilities.

Step 4: Run an MHD Stability Check

from scpn_fusion.core import run_full_stability_check

stability = run_full_stability_check(
    q_profile=result.q_profile,
    rho=result.rho,
    pressure=result.pressure_profile,
    R0=6.2, a=2.0, B0=5.3, Ip=15e6, kappa=1.7,
)

print(f"Mercier stable: {stability.mercier.stable}")
print(f"Kruskal-Shafranov: q_edge = {stability.kruskal_shafranov.q_edge:.2f}")
print(f"Troyon beta_N limit: {stability.troyon.beta_N_limit:.2f}")
print(f"Overall stable: {stability.stable}")

Step 5: Compute Transport

Run the 1.5D radial transport solver for a single timestep:

from scpn_fusion.core.integrated_transport_solver import IntegratedTransportSolver
import numpy as np

nr = 50
rho = np.linspace(0, 1, nr)

solver = IntegratedTransportSolver(
    rho=rho,
    R0=6.2,
    a=2.0,
    B0=5.3,
)

# Parabolic initial profiles
Te = 10.0 * (1 - rho**2)   # keV
ne = 10.0 * (1 - rho**2)   # 10^19 m^-3

# Heating source: 50 MW NBI, Gaussian deposition
P_heat = 50e6 * np.exp(-((rho - 0.3) / 0.15)**2)

Te_new = solver.step_energy(Te, ne, P_heat, dt=0.01)

fig, ax = plt.subplots()
ax.plot(rho, Te, label="Before")
ax.plot(rho, Te_new, label="After (10 ms)")
ax.set_xlabel(r"$\rho$")
ax.set_ylabel("Te [keV]")
ax.legend()
plt.show()

Step 6: Compute Confinement Scaling

Check how your scenario compares to the ITER H-mode database:

from scpn_fusion.core import ipb98y2_tau_e, compute_h_factor

tau_scaling = ipb98y2_tau_e(
    Ip_MA=15.0, B_T=5.3, n_19=10.0, P_MW=50.0,
    R_m=6.2, epsilon=2.0/6.2, kappa=1.7, M=2.5,
)
print(f"IPB98(y,2) tau_E = {tau_scaling:.2f} s")

tau_actual = 3.7  # measured or simulated
H = compute_h_factor(tau_actual, tau_scaling)
print(f"H-factor = {H:.2f}")

An H-factor of 1.0 means the plasma confines exactly as the scaling predicts; ITER targets \(H \geq 1.0\).

Step 7: Read Real Experimental Data

Load a SPARC GEQDSK file and compare:

from scpn_fusion.core import read_geqdsk

eq = read_geqdsk("validation/reference_data/sparc/lmode_vv.geqdsk")
print(f"R0 = {eq.R0:.3f} m, a = {eq.a:.3f} m")
print(f"B_T = {eq.B_T:.2f} T, Ip = {eq.Ip/1e6:.1f} MA")

What’s Next?

You’ve completed the basics. Depending on your interests: