Frontier Physics — Phase 3 Modules (v0.16.0)¶

Demonstrates the ten new physics modules added in scpn-control v0.16.0:

  1. Gyrokinetic transport — ITG/TEM/ETG dispersion, quasilinear fluxes
  2. Ballooning stability — s-alpha diagram, marginal stability
  3. Current diffusion — neoclassical resistive evolution
  4. Current drive — ECCD + NBI source profiles
  5. NTM dynamics — modified Rutherford equation
  6. Sawtooth — Kadomtsev reconnection crash model
  7. SOL model — two-point divertor heat flux
  8. Integrated scenario — ITER baseline full-pulse simulation

Dependencies: numpy, scipy, matplotlib (no GPU, no optional deps)

In [1]:
import sys
from pathlib import Path

sys.path.insert(0, str(Path(".").resolve().parent / "src"))

import numpy as np
import matplotlib.pyplot as plt

1. Gyrokinetic Transport — ITG Spectrum & Quasilinear Fluxes¶

The GyrokineticTransportModel computes ion and electron thermal diffusivities from a local gyrokinetic dispersion relation (ITG/TEM/ETG branches).

In [2]:
from scpn_control.core.gyrokinetic_transport import (
    GyrokineticsParams,
    GyrokineticTransportModel,
    compute_spectrum,
    quasilinear_fluxes,
)

params = GyrokineticsParams(
    R_L_Ti=6.0,  # R/L_Ti — ion temperature gradient
    R_L_Te=6.0,  # R/L_Te — electron temperature gradient
    R_L_ne=2.0,  # R/L_ne — density gradient
    q=1.5,
    s_hat=0.8,  # magnetic shear
    alpha_MHD=0.1,
    Te_Ti=1.0,
    Z_eff=1.5,
    nu_star=0.05,
    beta_e=0.01,
)

spec = compute_spectrum(params, n_modes=16)
fluxes = quasilinear_fluxes(params, spec)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))
ax1.plot(spec.k_y, spec.gamma_linear, "bo-", ms=4)
ax1.set_xlabel(r"$k_\theta \rho_s$")
ax1.set_ylabel(r"$\gamma$ [a.u.]")
ax1.set_title("Linear growth rate spectrum")
ax1.grid(True, alpha=0.3)

ax2.bar(["chi_i", "chi_e", "D_e"], [fluxes.chi_i, fluxes.chi_e, fluxes.D_e])
ax2.set_ylabel(r"Transport coefficient [m$^2$/s]")
ax2.set_title("Quasilinear fluxes")
plt.tight_layout()
plt.show()

print(f"chi_i = {fluxes.chi_i:.3f} m^2/s")
print(f"chi_e = {fluxes.chi_e:.3f} m^2/s")
print(f"D_e   = {fluxes.D_e:.3f} m^2/s")
No description has been provided for this image 
chi_i = 43.957 m^2/s
chi_e = 0.000 m^2/s
D_e   = 0.000 m^2/s

Radial transport profile¶

GyrokineticTransportModel.evaluate_profile computes diffusivities across the full radius.

In [3]:
model = GyrokineticTransportModel(n_modes=16)
rho = np.linspace(0.1, 0.9, 20)

profiles = {
    "Te_keV": 5.0 * (1 - rho**2),
    "Ti_keV": 5.0 * (1 - rho**2),
    "ne_19": 8.0 * (1 - 0.5 * rho**2),
    "q": 1.0 + 2.0 * rho**2,
    "R0": 6.2,
    "a": 2.0,
}

chi_i, chi_e, D_e = model.evaluate_profile(rho, profiles)

plt.figure(figsize=(8, 4))
plt.plot(rho, chi_i, "r-", lw=2, label=r"$\chi_i$")
plt.plot(rho, chi_e, "b--", lw=2, label=r"$\chi_e$")
plt.plot(rho, D_e, "g:", lw=2, label=r"$D_e$")
plt.xlabel(r"$\rho$")
plt.ylabel(r"m$^2$/s")
plt.title("Gyrokinetic transport profile")
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()
No description has been provided for this image

2. Ballooning Stability — s-alpha Diagram¶

The BallooningEquation solves the 1-D ballooning eigenvalue problem on a field line. find_marginal_stability locates the critical pressure gradient for each value of magnetic shear.

In [4]:
from scpn_control.core.ballooning_solver import (
    BallooningEquation,
    find_marginal_stability,
)

# Single-point solve
eq_stable = BallooningEquation(s=0.5, alpha=0.2)
res_stable = eq_stable.solve()
print(f"s=0.5, alpha=0.2 -> stable={res_stable.is_stable}")

eq_unstable = BallooningEquation(s=1.0, alpha=1.5)
res_unstable = eq_unstable.solve()
print(f"s=1.0, alpha=1.5 -> stable={res_unstable.is_stable}")

# Marginal stability curve
s_range = np.linspace(0.1, 2.0, 30)
alpha_crit = np.array([find_marginal_stability(s) for s in s_range])

plt.figure(figsize=(8, 5))
plt.plot(s_range, alpha_crit, "k-", lw=2)
plt.fill_between(s_range, 0, alpha_crit, alpha=0.15, color="green", label="Stable")
plt.fill_between(s_range, alpha_crit, 2.0, alpha=0.15, color="red", label="Unstable")
plt.xlabel("Magnetic shear s")
plt.ylabel(r"Normalized pressure gradient $\alpha$")
plt.title("s-alpha ballooning stability diagram")
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

s=0.5, alpha=0.2 -> stable=True
s=1.0, alpha=1.5 -> stable=False
No description has been provided for this image

3. Current Diffusion + Current Drive¶

The CurrentDiffusionSolver evolves the poloidal flux profile via Crank-Nicolson diffusion with neoclassical resistivity. ECCDSource and NBISource provide non-inductive current drive.

In [5]:
from scpn_control.core.current_diffusion import CurrentDiffusionSolver
from scpn_control.core.current_drive import ECCDSource, NBISource, CurrentDriveMix

nr = 50
rho = np.linspace(0, 1, nr)
R0, a, B0 = 6.2, 2.0, 5.3

cd_solver = CurrentDiffusionSolver(rho, R0, a, B0)

eccd = ECCDSource(P_ec_MW=10.0, rho_dep=0.4, sigma_rho=0.05)
nbi = NBISource(P_nbi_MW=33.0, E_beam_keV=1000.0, rho_tangency=0.2)
mix = CurrentDriveMix(a=a)
mix.add_source(eccd)
mix.add_source(nbi)

Te = 10.0 * (1 - rho**2) + 0.5
Ti = 9.0 * (1 - rho**2) + 0.5
ne = 8.0 * (1 - 0.3 * rho**2)

j_cd = mix.total_j_cd(rho, ne, Te, Ti)
j_bs = 0.3 * ne * np.sqrt(rho) * (1 - rho)

# Evolve 20 diffusion steps
psi_history = [cd_solver.psi.copy()]
for _ in range(20):
    cd_solver.step(dt=0.5, Te=Te, ne=ne, Z_eff=1.5, j_bs=j_bs, j_cd=j_cd)
    psi_history.append(cd_solver.psi.copy())

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))
ax1.plot(rho, j_cd, "r-", lw=2, label="j_CD (ECCD+NBI)")
ax1.plot(rho, j_bs, "b--", lw=2, label="j_BS (bootstrap)")
ax1.set_xlabel(r"$\rho$")
ax1.set_ylabel("Current density [a.u.]")
ax1.set_title("Non-inductive current sources")
ax1.legend()
ax1.grid(True, alpha=0.3)

for i in [0, 5, 10, 20]:
    ax2.plot(rho, psi_history[i], label=f"t={i * 0.5:.1f}s")
ax2.set_xlabel(r"$\rho$")
ax2.set_ylabel(r"$\psi$ [Wb]")
ax2.set_title("Flux profile evolution")
ax2.legend()
ax2.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
No description has been provided for this image

4. NTM Dynamics — Modified Rutherford Equation¶

The NTMIslandDynamics class integrates the modified Rutherford equation for neoclassical tearing mode (NTM) island evolution.

In [6]:
from scpn_control.core.ntm_dynamics import NTMIslandDynamics

ntm = NTMIslandDynamics(r_s=0.6, m=2, n=1, a=2.0, R0=6.2, B0=5.3)

# Island evolution without ECCD stabilization
dt = 0.01
w = 0.01  # initial seed island width [m]
w_no_eccd = [w]
for _ in range(500):
    dw = ntm.dw_dt(w, j_bs=1e5, j_phi=5e5, j_cd=0.0, eta=1e-7)
    w = max(w + dw * dt, 1e-4)
    w_no_eccd.append(w)

# Island evolution with ECCD stabilization
w = 0.01
w_eccd = [w]
for _ in range(500):
    dw = ntm.dw_dt(w, j_bs=1e5, j_phi=5e5, j_cd=3e5, eta=1e-7, d_cd=0.03)
    w = max(w + dw * dt, 1e-4)
    w_eccd.append(w)

t = np.arange(len(w_no_eccd)) * dt
plt.figure(figsize=(8, 4))
plt.plot(t, np.array(w_no_eccd) * 100, "r-", lw=2, label="No ECCD")
plt.plot(t, np.array(w_eccd) * 100, "b-", lw=2, label="With ECCD")
plt.xlabel("Time [s]")
plt.ylabel("Island width [cm]")
plt.title("2/1 NTM island evolution")
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()
No description has been provided for this image

5. Sawtooth — Kadomtsev Crash¶

The kadomtsev_crash function applies the Kadomtsev reconnection model to flatten temperature and density profiles inside the mixing radius.

In [7]:
from scpn_control.core.sawtooth import kadomtsev_crash

rho = np.linspace(0, 1, 200)
T_pre = 10.0 * (1 - rho**2) ** 1.5 + 0.5
n_pre = 8.0 * (1 - 0.3 * rho**2)
q_pre = 0.85 + 1.5 * rho**2  # q(0) < 1 triggers sawtooth

T_post, n_post, q_post, rho_1, rho_mix = kadomtsev_crash(
    rho,
    T_pre,
    n_pre,
    q_pre,
    R0=6.2,
    a=2.0,
)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))
ax1.plot(rho, T_pre, "r--", lw=1.5, label="Pre-crash")
ax1.plot(rho, T_post, "b-", lw=2, label="Post-crash")
ax1.axvline(rho_1, color="gray", ls=":", label=f"q=1 at rho={rho_1:.2f}")
ax1.axvline(rho_mix, color="orange", ls=":", label=f"rho_mix={rho_mix:.2f}")
ax1.set_xlabel(r"$\rho$")
ax1.set_ylabel("T [keV]")
ax1.set_title("Temperature profile")
ax1.legend(fontsize=9)
ax1.grid(True, alpha=0.3)

ax2.plot(rho, q_pre, "r--", lw=1.5, label="Pre-crash")
ax2.plot(rho, q_post, "b-", lw=2, label="Post-crash")
ax2.axhline(1.0, color="gray", ls=":", alpha=0.5)
ax2.set_xlabel(r"$\rho$")
ax2.set_ylabel("q")
ax2.set_title("Safety factor profile")
ax2.legend(fontsize=9)
ax2.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
No description has been provided for this image

6. SOL Model — Two-Point Divertor Heat Flux¶

The TwoPointSOL model computes upstream and target plasma conditions from power crossing the separatrix. The eich_heat_flux_width function implements the Eich scaling (Eich et al., NF 2013).

In [8]:
from scpn_control.core.sol_model import (
    TwoPointSOL,
    eich_heat_flux_width,
)

sol = TwoPointSOL(R0=6.2, a=2.0, q95=3.0, B_pol=1.0, kappa=1.7)

P_SOL_range = np.linspace(10, 100, 20)
T_target = []
q_par = []

for P in P_SOL_range:
    soln = sol.solve(P_SOL_MW=P, n_u_19=5.0)
    T_target.append(soln.T_target_eV)
    q_par.append(soln.q_parallel_MW_m2)

lambda_q = eich_heat_flux_width(P_SOL_MW=50.0, R0=6.2, B_pol=1.0, epsilon=2.0 / 6.2)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))
ax1.plot(P_SOL_range, T_target, "r-", lw=2)
ax1.set_xlabel("P_SOL [MW]")
ax1.set_ylabel("T_target [eV]")
ax1.set_title("Target temperature vs power")
ax1.grid(True, alpha=0.3)

ax2.plot(P_SOL_range, q_par, "b-", lw=2)
ax2.set_xlabel("P_SOL [MW]")
ax2.set_ylabel(r"$q_\parallel$ [MW/m$^2$]")
ax2.set_title("Parallel heat flux to divertor")
ax2.grid(True, alpha=0.3)
plt.suptitle(f"Eich heat flux width: lambda_q = {lambda_q:.2f} mm", fontsize=11)
plt.tight_layout()
plt.show()
No description has been provided for this image

7. Integrated Scenario — ITER 15 MA Baseline¶

The IntegratedScenarioSimulator couples transport, current diffusion, current drive, sawteeth, NTMs, and SOL into a single time-stepping loop.

In [9]:
import dataclasses
from scpn_control.core.integrated_scenario import (
    IntegratedScenarioSimulator,
    iter_baseline_scenario,
)

# Run a short 10s window of ITER baseline
config = dataclasses.replace(iter_baseline_scenario(), t_end=10.0, dt=0.5)

sim = IntegratedScenarioSimulator(config)
states = sim.run()

times = [s.time for s in states]
Te_axis = [s.Te[0] for s in states]
q_axis = [s.q[0] for s in states]

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 6), sharex=True)
ax1.plot(times, Te_axis, "r-", lw=2)
ax1.set_ylabel("T_e(0) [keV]")
ax1.set_title("ITER 15 MA -- Integrated scenario")
ax1.grid(True, alpha=0.3)

ax2.plot(times, q_axis, "b-", lw=2)
ax2.axhline(1.0, color="gray", ls=":", alpha=0.5)
ax2.set_xlabel("Time [s]")
ax2.set_ylabel("q(0)")
ax2.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

print(f"Simulated {len(states)} steps, t = {times[0]:.1f} to {times[-1]:.1f} s")
print(f"NTM islands: {states[-1].ntm_island_widths}")
No description has been provided for this image 
Simulated 20 steps, t = 0.5 to 10.0 s
NTM islands: {}

Summary: scpn-control v0.16.0 adds first-principles models for gyrokinetic transport, ideal MHD stability, current profile evolution, NTM/sawtooth instabilities, and scrape-off layer heat exhaust. These modules compose into the IntegratedScenarioSimulator for full-pulse tokamak simulation.

See docs/api.md for the complete API reference.