scpn_fusion.core – Physics

The core subpackage contains the plasma physics engines: equilibrium solvers, transport models, stability analysers, turbulence surrogates, heating models, and geometry tools.

Equilibrium Solver

Non-linear free-boundary Grad-Shafranov equilibrium solver.

Solves the Grad-Shafranov equation for toroidal plasma equilibrium using Picard iteration with under-relaxation. Supports both L-mode (linear) and H-mode (mtanh pedestal) pressure/current profiles.

The solver can optionally offload the inner elliptic solve to a compiled C++ library via HPCBridge, or to the Rust multigrid backend when available.

class scpn_fusion.core.fusion_kernel.CoilSet(positions=<factory>, currents=<factory>, turns=<factory>, current_limits=None, target_flux_points=None, target_flux_values=None, x_point_target=None, x_point_flux_target=None, divertor_strike_points=None, divertor_flux_values=None)[source]

Bases: object

External coil set for free-boundary solve.

Parameters:
positions

Coil centre coordinates [m].

Type:

list of (R, Z) tuples

currents

Current per coil [A].

Type:

NDArray

turns

Number of turns per coil.

Type:

list of int

current_limits

Per-coil maximum absolute current [A]. Shape (n_coils,). When set, optimize_coil_currents enforces these bounds.

Type:

NDArray or None

target_flux_points

Points (R, Z) on the desired separatrix for shape optimisation. Shape (n_pts, 2).

Type:

NDArray or None

target_flux_values

Optional target flux values [Wb] at target_flux_points. Shape (n_pts,). When omitted, solve_free_boundary uses an isoflux target inferred from current interpolation.

Type:

NDArray or None

x_point_target

Target X-point position (R, Z) [m]. Shape (2,).

Type:

NDArray or None

x_point_flux_target

Target flux value [Wb] at the X-point.

Type:

float or None

divertor_strike_points

Divertor strike-point positions (R, Z) [m]. Shape (n, 2).

Type:

NDArray or None

divertor_flux_values

Target flux values [Wb] at divertor strike points. Shape (n,).

Type:

NDArray or None

positions: list[tuple[float, float]]
currents: ndarray[Any, dtype[float64]]
turns: list[int]
current_limits: ndarray[Any, dtype[float64]] | None = None
target_flux_points: ndarray[Any, dtype[float64]] | None = None
target_flux_values: ndarray[Any, dtype[float64]] | None = None
x_point_target: ndarray[Any, dtype[float64]] | None = None
x_point_flux_target: float | None = None
divertor_strike_points: ndarray[Any, dtype[float64]] | None = None
divertor_flux_values: ndarray[Any, dtype[float64]] | None = None
class scpn_fusion.core.fusion_kernel.FusionKernel(config_path)[source]

Bases: FusionKernelNewtonSolverMixin, FusionKernelIterativeSolverMixin

Non-linear free-boundary Grad-Shafranov equilibrium solver.

Parameters:

config_path (str | Path) – Path to a JSON configuration file describing the reactor geometry, coil set, physics parameters and solver settings.

Psi

Poloidal flux on the (Z, R) grid.

Type:

FloatArray

J_phi

Toroidal current density on the (Z, R) grid.

Type:

FloatArray

B_R, B_Z

Radial and vertical magnetic field components (set after solve).

Type:

FloatArray

load_config(path)[source]

Load reactor configuration from a JSON file.

Parameters:

path (str | Path) – Filesystem path to the configuration JSON.

Return type:

None

initialize_grid()[source]

Build the computational (R, Z) mesh from the loaded config.

Return type:

None

setup_accelerator()[source]

Initialise the optional C++ HPC acceleration bridge.

Return type:

None

calculate_vacuum_field()[source]

Compute the vacuum poloidal flux from the external coil set.

Uses elliptic integrals (toroidal Green’s function) for each coil.

Returns:

Vacuum flux Psi_vac on the (NZ, NR) grid.

Return type:

FloatArray

find_x_point(Psi)[source]

Locate the X-point (magnetic null) in the lower divertor region.

Parameters:

Psi (FloatArray) – Poloidal flux array on the (NZ, NR) grid.

Returns:

((R_x, Z_x), Psi_x) — position and flux value at the X-point.

Return type:

tuple[tuple[float, float], float]

static mtanh_profile(psi_norm, params)[source]

Evaluate a modified-tanh pedestal profile (vectorised).

Parameters:

  • psi_norm (FloatArray) – Normalised poloidal flux (0 at axis, 1 at separatrix).

  • params (dict[str, float]) – Profile shape parameters with keys ped_top, ped_width, ped_height, core_alpha.

Returns:

Profile value; zero outside the plasma region.

Return type:

FloatArray

update_plasma_source_nonlinear(Psi_axis, Psi_boundary)[source]

Compute the toroidal current density J_phi from the GS source.

Uses J_phi = R p'(psi) + FF'(psi) / (mu0 R) with either L-mode (linear) or H-mode (mtanh) profiles, then renormalises to match the target plasma current.

Parameters:

  • Psi_axis (float) – Poloidal flux at the magnetic axis (O-point).

  • Psi_boundary (float) – Poloidal flux at the separatrix (X-point or limiter).

Returns:

Updated J_phi on the (NZ, NR) grid.

Return type:

FloatArray

compute_b_field()[source]

Derive the magnetic field components from the solved Psi.

Return type:

None

optimize_coil_currents(coils, target_flux, tikhonov_alpha=0.0001)[source]

Find coil currents that best reproduce target_flux at the control points.

Solves the bounded linear least-squares problem:

min_I || M^T I - target_flux ||^2 + alpha * ||I||^2 s.t. -I_max <= I <= I_max (per coil)

where M is the mutual-inductance matrix.

Parameters:

  • coils (CoilSet) – Must have target_flux_points set (shape (n_pts, 2)).

  • target_flux (FloatArray, shape (n_pts,)) – Desired poloidal flux at each control point.

  • tikhonov_alpha (float) – Regularisation strength to penalise large currents.

Returns:

Optimised coil currents [A].

Return type:

FloatArray, shape (n_coils,)

solve_free_boundary(coils, max_outer_iter=20, tol=0.0001, optimize_shape=False, tikhonov_alpha=0.0001)[source]

Free-boundary GS solve with external coil currents.

Iterates between updating boundary flux from coils and solving the internal GS equation. When optimize_shape=True and the coil set has target_flux_points, an additional outer loop optimises the coil currents to match the desired plasma boundary shape.

Parameters:

  • coils (CoilSet) – External coil set.

  • max_outer_iter (int) – Maximum free-boundary iterations.

  • tol (float) – Convergence tolerance on max abs(delta psi).

  • optimize_shape (bool) – When True, run coil-current optimisation at each outer step.

  • tikhonov_alpha (float) – Tikhonov regularisation for coil optimisation.

Returns:

{"outer_iterations": int, "final_diff": float, "coil_currents": NDArray}

Return type:

dict

save_results(filename='equilibrium_nonlinear.npz')[source]

Save the equilibrium state to a compressed NumPy archive.

Parameters:

filename (str) – Output file path.

Return type:

None

GEQDSK I/O

Reader and writer for the G-EQDSK (EFIT) equilibrium file format.

The G-EQDSK format is the de facto standard for tokamak equilibrium data exchange. It stores the poloidal flux map ψ(R,Z) on a uniform (R,Z) grid together with 1-D profile arrays and scalar quantities.

References

  • Fortran fixed-width: 5 values per line, (5e16.9)

  • Header line: 48-char description, 3 ints (idum, nw, nh)

  • Scalars block (20 values): rdim, zdim, rcentr, rleft, zmid, rmaxis, zmaxis, simag, sibry, bcentr, current, …

  • 1-D arrays (each nw values): fpol, pres, ffprime, pprime, qpsi

  • 2-D array: psirz (nh × nw values, row-major)

  • Boundary & limiter point counts + (R,Z) pairs

class scpn_fusion.core.eqdsk.GEqdsk(description='', nw=0, nh=0, rdim=0.0, zdim=0.0, rcentr=0.0, rleft=0.0, zmid=0.0, rmaxis=0.0, zmaxis=0.0, simag=0.0, sibry=0.0, bcentr=0.0, current=0.0, fpol=<factory>, pres=<factory>, ffprime=<factory>, pprime=<factory>, qpsi=<factory>, psirz=<factory>, rbdry=<factory>, zbdry=<factory>, rlim=<factory>, zlim=<factory>)[source]

Bases: object

Container for all data in a G-EQDSK file.

Parameters:
description: str = ''
nw: int = 0
nh: int = 0
rdim: float = 0.0
zdim: float = 0.0
rcentr: float = 0.0
rleft: float = 0.0
zmid: float = 0.0
rmaxis: float = 0.0
zmaxis: float = 0.0
simag: float = 0.0
sibry: float = 0.0
bcentr: float = 0.0
current: float = 0.0
fpol: ndarray[Any, dtype[float64]]
pres: ndarray[Any, dtype[float64]]
ffprime: ndarray[Any, dtype[float64]]
pprime: ndarray[Any, dtype[float64]]
qpsi: ndarray[Any, dtype[float64]]
psirz: ndarray[Any, dtype[float64]]
rbdry: ndarray[Any, dtype[float64]]
zbdry: ndarray[Any, dtype[float64]]
rlim: ndarray[Any, dtype[float64]]
zlim: ndarray[Any, dtype[float64]]
property r: ndarray[Any, dtype[float64]]

1-D array of R grid values.

property z: ndarray[Any, dtype[float64]]

1-D array of Z grid values.

property psi_norm: ndarray[Any, dtype[float64]]

Normalised poloidal flux ψ_N ∈ [0, 1] (axis=0, boundary=1).

psi_to_norm(psi)[source]

Map raw ψ to normalised ψ_N.

Return type:

ndarray[Any, dtype[float64]]

Parameters:

psi (ndarray[Any, dtype[float64]])

to_config(name='eqdsk')[source]

Convert to FusionKernel JSON config dict (approximate).

Return type:

dict[str, object]

Parameters:

name (str)

scpn_fusion.core.eqdsk.read_geqdsk(path)[source]

Read a G-EQDSK file and return a GEqdsk container.

Handles both whitespace-separated and fixed-width Fortran formats, including the common case where values run together without spaces (e.g. 2.385E+00-1.216E+01).

Parameters:

path (str or Path) – Path to the G-EQDSK file.

Returns:

Parsed equilibrium data.

Return type:

GEqdsk

scpn_fusion.core.eqdsk.write_geqdsk(eq, path)[source]

Write a GEqdsk to a G-EQDSK file.

Parameters:

  • eq (GEqdsk) – Equilibrium data.

  • path (str or Path) – Output file path.

Return type:

None

Force Balance

class scpn_fusion.core.force_balance.ForceBalanceSolver(config_path)[source]

Bases: object

Automatic Newton-Raphson Solver to find Coil Currents that result in ZERO Radial Force on the plasma at the target Radius. This guarantees static equilibrium at R_target.

solve_for_equilibrium(target_R=6.2, target_Z=0.0, *, max_iterations=10, jacobian_floor=1e-12)[source]
Parameters:

  • max_iterations (int)

  • jacobian_floor (float)

save_config(output_path=None)[source]
Parameters:

output_path (str | Path | None)

Compact Reactor Optimiser

class scpn_fusion.core.compact_reactor_optimizer.CompactReactorArchitect[source]

Bases: object

Optimizes Fusion Reactor Geometry for MINIMUM SIZE.

plasma_physics_model(R, a, B0)[source]
radial_build_constraints(R, a, B0)[source]
visualize_space(designs, label='')[source]
calculate_economics(d)[source]

Calculates Cost of Electricity (CoE) [$/MWh]. Sheffield model: CoE ~ (C_cap * F_cap + C_om) / (8760 * P_net * f_avail)

report_design(d)[source]
find_minimum_reactor(target_power_MW=1.0, use_temhd=True)[source]

Global Design Scanner

class scpn_fusion.core.global_design_scanner.GlobalDesignExplorer(base_config_path, *, divertor_flux_cap_mw_m2=45.0, zeff_cap=0.4, hts_peak_cap_t=21.0)[source]

Bases: object

Monte Carlo Design Space Explorer. Searches for the Pareto Frontier of Fusion Reactors. Objectives: Maximize Q, Minimize Radius (Cost), Minimize Wall Load.

evaluate_design(R_maj, B_field, I_plasma)[source]

Runs a full-stack simulation for a specific design point.

run_scan(n_samples=2000, *, r_bounds=(2.0, 9.0), b_bounds=(4.0, 12.0), i_bounds=(5.0, 25.0), seed=None, q95_min=3.0)[source]
run_compact_scan(n_samples=2000, seed=42)[source]
analyze_pareto(df)[source]

Integrated Transport Solver

class scpn_fusion.core.integrated_transport_solver.AdaptiveTimeController(dt_init=0.01, dt_min=1e-05, dt_max=1.0, tol=0.001, safety=0.9)[source]

Bases: object

Richardson-extrapolation adaptive time controller for CN transport.

Compares one full CN step vs. two half-steps to estimate the local truncation error, then uses a PI controller to adjust dt.

Parameters:

  • dt_init (float — initial time step [s])

  • dt_min (float — minimum allowed dt)

  • dt_max (float — maximum allowed dt)

  • tol (float — target local error tolerance)

  • safety (float — safety factor (< 1) for step adjustment)

estimate_error(solver, P_aux, *, enforce_numerical_recovery=False, max_numerical_recoveries=None)[source]

Estimate local error via Richardson extrapolation.

Takes one full CN step of size dt and two half-steps of size dt/2, then compares. The solver state is advanced by the half-step result (more accurate).

Returns the estimated error norm.

Return type:

float

Parameters:
adapt_dt(error)[source]

Adjust dt using a PI controller.

dt *= min(2, safety * (tol/err)^(0.7/p) * (err_prev/err)^(0.4/p))

Return type:

None

Parameters:

error (float)

exception scpn_fusion.core.integrated_transport_solver.PhysicsError[source]

Bases: RuntimeError, FusionCoreError

Raised when a physics constraint is violated.

class scpn_fusion.core.integrated_transport_solver.TransportSolver(config_path, *, multi_ion=False)[source]

Bases: TransportSolverModelMixin, TransportSolverRuntimeMixin, FusionKernel

1.5D Integrated Transport Code. Solves Heat and Particle diffusion equations on flux surfaces, coupled self-consistently with the 2D Grad-Shafranov equilibrium.

When multi_ion=True, the solver evolves separate D/T fuel densities, He-ash transport with pumping (configurable tau_He), independent electron temperature Te, coronal-equilibrium tungsten radiation (Pütterich et al. 2010), and per-cell Bremsstrahlung.

Parameters:

  • config_path (str | Path)

  • multi_ion (bool)

scpn_fusion.core.integrated_transport_solver.IntegratedTransportSolver

alias of TransportSolver

scpn_fusion.core.integrated_transport_solver.chang_hinton_chi_profile(rho, T_i, n_e_19, q, R0, a, B0, A_ion=2.0, Z_eff=1.5)[source]

Chang-Hinton (1982) neoclassical ion thermal diffusivity profile [m²/s].

Parameters:

  • rho (array — normalised radius [0,1])

  • T_i (array — ion temperature [keV])

  • n_e_19 (array — electron density [10^19 m^-3])

  • q (array — safety factor profile)

  • R0 (float — major radius [m])

  • a (float — minor radius [m])

  • B0 (float — toroidal field [T])

  • A_ion (float — ion mass number (default 2 = deuterium))

  • Z_eff (float — effective charge)

Returns:

chi_nc

Return type:

array — neoclassical chi_i [m²/s]

scpn_fusion.core.integrated_transport_solver.calculate_sauter_bootstrap_current_full(rho, Te, Ti, ne, q, R0, a, B0, Z_eff=1.5)[source]

Full Sauter bootstrap current model (Sauter et al., Phys. Plasmas 6, 1999).

Parameters:

  • rho (array — normalised radius [0,1])

  • Te (array — electron temperature [keV])

  • Ti (array — ion temperature [keV])

  • ne (array — electron density [10^19 m^-3])

  • q (array — safety factor profile)

  • R0 (float — major radius [m])

  • a (float — minor radius [m])

  • B0 (float — toroidal field [T])

  • Z_eff (float — effective charge)

Returns:

j_bs

Return type:

array — bootstrap current density [A/m^2]

RF Heating

class scpn_fusion.core.rf_heating.RFHeatingSystem(config_path)[source]

Bases: object

Simulates Ion Cyclotron Resonance Heating (ICRH). Uses Ray-Tracing to track EM waves launching from the antenna and absorbing at the resonance layer.

get_plasma_params(R, Z)[source]

Returns B_mod, density, and derivatives at (R,Z).

dispersion_relation(R, Z, k_R, k_Z)[source]

Calculates the local dispersion D(w, k) = 0. Harden with Warm Plasma thermal corrections for ICRH.

ray_equations(state, t)[source]

Hamiltonian Ray Tracing equations. dr/dt = dD/dk dk/dt = -dD/dr

trace_rays(n_rays=10)[source]
plot_heating(trajectories, B_res)[source]
compute_power_deposition(trajectories, P_rf_mw=20.0, n_radial_bins=50)[source]

Compute radial power deposition profile from ray absorption.

Uses cyclotron damping: the imaginary part of the wave vector causes exponential power decay along each ray path. The absorption coefficient is proportional to exp(-(omega - n*Omega_ci)^2 / (k_par * v_thi)^2).

Parameters:

  • trajectories (list of ndarray) – Ray paths from trace_rays().

  • P_rf_mw (float) – Total injected RF power (MW).

  • n_radial_bins (int) – Number of radial bins for the deposition profile.

Returns:

  • rho_bins (ndarray) – Normalised radius bin centres (0–1).

  • P_dep_mw_m3 (ndarray) – Power deposition density (MW/m^3) per radial bin.

  • absorption_efficiency (float) – Fraction of injected power absorbed (0–1).

class scpn_fusion.core.rf_heating.ECRHHeatingSystem(b0_tesla=5.3, r0_major=6.2, freq_ghz=170.0, harmonic=1)[source]

Bases: object

Electron Cyclotron Resonance Heating (ECRH) with ray-tracing absorption.

Operates at 170 GHz (ITER ECRH frequency). Resonance occurs where omega = n * Omega_ce (electron cyclotron frequency), with n=1 (fundamental) or n=2 (second harmonic).

Parameters:

  • b0_tesla (float) – On-axis toroidal magnetic field (T).

  • r0_major (float) – Major radius (m).

  • freq_ghz (float) – ECRH frequency (GHz). Default: 170.

  • harmonic (int) – Cyclotron harmonic number (1 or 2).

resonance_radius()[source]

Major radius where n*Omega_ce = omega.

Return type:

float

compute_deposition(P_ecrh_mw=20.0, n_radial_bins=50, T_e_keV=20.0, n_e=1e+20, launch_angle_deg=0.0)[source]

Compute ECRH power deposition profile.

Uses Gaussian deposition centred at the resonance layer with width determined by the Doppler broadening of the electron cyclotron resonance.

Return type:

tuple[ndarray, ndarray, float]

Returns:

  • rho_bins (ndarray) – Normalised radius (0–1).

  • P_dep_mw_m3 (ndarray) – Power deposition density.

  • absorption_efficiency (float) – Fraction absorbed (0–1).

Parameters:

MHD Sawtooth

class scpn_fusion.core.mhd_sawtooth.ReducedMHD(nr=100)[source]

Bases: object

Simulates the internal m=1, n=1 Kink Mode (Sawtooth Instability). Solves Reduced MHD equations in cylindrical geometry. Demonstrates Magnetic Reconnection (Kadomtsev Model).

Parameters:

nr (int)

laplacian(f, m=1)[source]

Radial Laplacian: 1/r d/dr (r df/dr) - m^2/r^2 f

Return type:

ndarray

Parameters:
step(dt=0.01)[source]

Time integration (Semi-Implicit) Equations: 1. d(W_11)/dt = [J_eq, psi_11] + [J_11, psi_eq] + Dissipation 2. d(psi_11)/dt = B_parallel * phi_11 + eta * J_11

Return type:

tuple[float, bool]

Parameters:

dt (float)

solve_poisson(U)[source]

Solves Del^2 phi = U for phi using the Thomas Algorithm (O(N)). Optimized for tridiagonal Laplacian in cylindrical coordinates.

Return type:

ndarray

Parameters:

U (ndarray)

scpn_fusion.core.mhd_sawtooth.run_sawtooth_sim()[source]
Return type:

None

Hall-MHD Discovery

scpn_fusion.core.hall_mhd_discovery.spitzer_resistivity(T_e_eV, Z_eff=1.0, ln_lambda=17.0)[source]

Spitzer resistivity [Ohm*m]. eta = 1.65e-9 * Z_eff * ln_lambda / T_e^1.5

class scpn_fusion.core.hall_mhd_discovery.HallMHD(N=64, eta=None, nu=None)[source]

Bases: object

3D-like Reduced Hall-MHD with magnetic flutter and shear flows.

Fields: phi (stream function), psi (magnetic flux), U (vorticity), J (current density).

poisson_bracket(A_k, B_k)[source]
dynamics(phi, psi)[source]

Reduced MHD Equations: d(U)/dt = [phi, U] + [J, psi] + nu*del^4 U d(psi)/dt = [phi, psi] + rho_s^2 [J, psi] - eta*J where U = del^2 phi, J = del^2 psi

step()[source]
parameter_sweep(eta_range, nu_range, n_steps=5, sim_steps=200)[source]

Run grid of simulations varying eta and nu. Returns dict with growth rates.

find_tearing_threshold(eta_range=(1e-06, 0.01), n_bisect=10, sim_steps=500)[source]

Bisection search for marginal tearing stability threshold.

scpn_fusion.core.hall_mhd_discovery.run_discovery_sim()[source]

Stability Analyser

class scpn_fusion.core.stability_analyzer.StabilityAnalyzer(config_path)[source]

Bases: object

Performs Linear Stability Analysis (Eigenvalue analysis) of the plasma rigid-body motion (n=0 mode).

get_vacuum_field_at(R, Z)[source]

Interpolates Bz and Br from vacuum coils at position (R,Z). Bz = 1/R * dPsi/dR Br = -1/R * dPsi/dZ

calculate_forces(R, Z, Ip)[source]

Calculates radial and vertical forces acting on the plasma ring. F_R = F_Hoop + F_Lorentz_R F_Z = F_Lorentz_Z

analyze_stability(R_target=6.2, Z_target=0.0)[source]
analyze_mhd_stability(R0=6.2, a=2.0, B0=5.3, Ip_MA=15.0, transport_solver=None)[source]

Run Mercier and ballooning stability analysis.

Uses either profiles from transport_solver (if provided) or default parabolic profiles.

Parameters:

  • R0 (float — major radius [m])

  • a (float — minor radius [m])

  • B0 (float — toroidal field [T])

  • Ip_MA (float — plasma current [MA])

  • transport_solver (TransportSolver, optional) – If given, uses its rho/ne/Ti/Te profiles.

Return type:

dict with keys q_profile, mercier, ballooning

plot_stability_landscape(R0, Z0)[source]

Neural Equilibrium

PCA + MLP surrogate for Grad-Shafranov equilibrium reconstruction.

Maps coil currents (or profile parameters) → PCA coefficients → ψ(R,Z) with ~1000× speedup over the full Picard iteration.

Training modes:

  1. From FusionKernel — Generate training data by perturbing coil currents and running the GS solver. Requires a valid config JSON.

  2. From SPARC GEQDSKs — Train on real equilibrium data from CFS. Uses the GEQDSK’s profile parameters (p’, FF’, I_p) as input features and ψ(R,Z) as targets. No coil model needed.

Status: Reduced-order surrogate. Not on the critical control path. Use for rapid design-space exploration and batch equilibrium sweeps.

class scpn_fusion.core.neural_equilibrium.NeuralEqConfig(n_components=20, hidden_sizes=(128, 64, 32), n_input_features=12, grid_shape=(129, 129), lambda_gs=0.1)[source]

Bases: object

Configuration for the neural equilibrium model.

Parameters:
n_components: int = 20
hidden_sizes: tuple[int, ...] = (128, 64, 32)
n_input_features: int = 12
grid_shape: tuple[int, int] = (129, 129)
lambda_gs: float = 0.1
class scpn_fusion.core.neural_equilibrium.TrainingResult(n_samples, n_components, explained_variance, final_loss, train_time_s, weights_path, val_loss=nan, test_mse=nan, test_max_error=nan)[source]

Bases: object

Training summary.

Parameters:
n_samples: int
n_components: int
explained_variance: float
final_loss: float
train_time_s: float
weights_path: str
val_loss: float = nan
test_mse: float = nan
test_max_error: float = nan
class scpn_fusion.core.neural_equilibrium.SimpleMLP(layer_sizes, seed=42)[source]

Bases: object

Feedforward MLP with ReLU hidden layers and linear output.

Parameters:
forward(x)[source]
Return type:

ndarray[Any, dtype[TypeVar(_ScalarType_co, bound= generic, covariant=True)]]

Parameters:

x (ndarray[Any, dtype[_ScalarType_co]])

predict(x)[source]
Return type:

ndarray[Any, dtype[TypeVar(_ScalarType_co, bound= generic, covariant=True)]]

Parameters:

x (ndarray[Any, dtype[_ScalarType_co]])

class scpn_fusion.core.neural_equilibrium.MinimalPCA(n_components=20)[source]

Bases: object

Minimal PCA via SVD, no sklearn required.

Parameters:

n_components (int)

fit(X)[source]
Return type:

MinimalPCA

Parameters:

X (ndarray[Any, dtype[_ScalarType_co]])

transform(X)[source]
Return type:

ndarray[Any, dtype[TypeVar(_ScalarType_co, bound= generic, covariant=True)]]

Parameters:

X (ndarray[Any, dtype[_ScalarType_co]])

inverse_transform(Z)[source]
Return type:

ndarray[Any, dtype[TypeVar(_ScalarType_co, bound= generic, covariant=True)]]

Parameters:

Z (ndarray[Any, dtype[_ScalarType_co]])

fit_transform(X)[source]
Return type:

ndarray[Any, dtype[TypeVar(_ScalarType_co, bound= generic, covariant=True)]]

Parameters:

X (ndarray[Any, dtype[_ScalarType_co]])

class scpn_fusion.core.neural_equilibrium.NeuralEquilibriumAccelerator(config=None)[source]

Bases: object

PCA + MLP surrogate for Grad-Shafranov equilibrium.

Can be trained from SPARC GEQDSK files (preferred) or from a FusionKernel config with coil perturbations.

Parameters:

config (NeuralEqConfig | None)

evaluate_surrogate(X_test, Y_test_raw)[source]

Evaluate on test set. Returns dict with mse, max_error, gs_residual.

Return type:

dict[str, float]

Parameters:
train_from_geqdsk(geqdsk_paths, n_perturbations=25, seed=42)[source]

Train on real SPARC GEQDSK equilibria with perturbations.

For each GEQDSK file, generates n_perturbations by scaling p’/FF’ profiles, yielding n_files × n_perturbations training pairs.

Return type:

TrainingResult

Parameters:
Input features (12-dim):
[I_p, B_t, R_axis, Z_axis, pprime_scale, ffprime_scale,

simag, sibry, kappa, delta_upper, delta_lower, q95]

Output: flattened ψ(R,Z) → PCA coefficients

Uses 70/15/15 train/val/test split with val-loss early stopping (patience=20) and combined MSE + lambda_gs * GS residual loss.

predict(features)[source]

Predict ψ(R,Z) from input features.

Parameters:

features (NDArray) – Shape (n_features,) or (batch, n_features).

Returns:

Shape (nh, nw) or (batch, nh, nw).

Return type:

NDArray

save_weights(path=PosixPath('/home/runner/work/scpn-fusion-core/scpn-fusion-core/weights/neural_equilibrium_sparc.npz'))[source]

Save model to .npz (no pickle dependency).

Return type:

None

Parameters:

path (str | Path)

load_weights(path=PosixPath('/home/runner/work/scpn-fusion-core/scpn-fusion-core/weights/neural_equilibrium_sparc.npz'))[source]

Load model from .npz.

Return type:

None

Parameters:

path (str | Path)

benchmark(features, n_runs=100)[source]

Time inference over n_runs and return stats.

Return type:

dict[str, float]

Parameters:

Neural Transport

Neural-network surrogate for turbulent transport coefficients.

Replaces the simple critical-gradient transport model with a trained MLP that reproduces gyrokinetic-level predictions at millisecond inference speeds. When no trained weights are available the module falls back to an analytic critical-gradient model, so existing code keeps working without any neural network dependency.

The architecture follows the QLKNN paradigm (van de Plassche et al., Phys. Plasmas 27, 022310, 2020): a small MLP maps local plasma parameters to turbulent fluxes (heat, particle) across ITG/TEM/ETG channels.

Training data

The module is designed for the QLKNN-10D public dataset:

Download the dataset and run the training recipe in docs/NEURAL_TRANSPORT_TRAINING.md to produce an .npz weight file that this module loads at construction time.

References

class scpn_fusion.core.neural_transport.TransportInputs(rho=0.5, te_kev=5.0, ti_kev=5.0, ne_19=5.0, grad_te=6.0, grad_ti=6.0, grad_ne=2.0, q=1.5, s_hat=0.8, beta_e=0.01, r_major_m=6.2, a_minor_m=2.0, b_tesla=5.3, z_eff=1.0)[source]

Bases: object

Local plasma parameters at a single radial location.

All quantities are in SI / conventional tokamak units.

Parameters:

  • rho (float) – Normalised toroidal flux coordinate (0 = axis, 1 = edge).

  • te_kev (float) – Electron temperature [keV].

  • ti_kev (float) – Ion temperature [keV].

  • ne_19 (float) – Electron density [10^19 m^-3].

  • grad_te (float) – Normalised electron temperature gradient R/L_Te.

  • grad_ti (float) – Normalised ion temperature gradient R/L_Ti.

  • grad_ne (float) – Normalised electron density gradient R/L_ne.

  • q (float) – Safety factor.

  • s_hat (float) – Magnetic shear s = (r/q)(dq/dr).

  • beta_e (float) – Electron beta (kinetic pressure / magnetic pressure).

  • r_major_m (float) – Major radius used in reduced gyro-Bohm normalization [m].

  • a_minor_m (float) – Minor radius used for inverse-aspect-ratio effects [m].

  • b_tesla (float) – Toroidal magnetic field used in reduced gyro-Bohm normalization [T].

  • z_eff (float) – Effective charge used for collisionality estimate.

rho: float = 0.5
te_kev: float = 5.0
ti_kev: float = 5.0
ne_19: float = 5.0
grad_te: float = 6.0
grad_ti: float = 6.0
grad_ne: float = 2.0
q: float = 1.5
s_hat: float = 0.8
beta_e: float = 0.01
r_major_m: float = 6.2
a_minor_m: float = 2.0
b_tesla: float = 5.3
z_eff: float = 1.0
class scpn_fusion.core.neural_transport.TransportFluxes(chi_e=0.0, chi_i=0.0, d_e=0.0, channel='stable', chi_e_itg=0.0, chi_e_tem=0.0, chi_e_etg=0.0, chi_i_itg=0.0)[source]

Bases: object

Predicted turbulent transport fluxes.

Parameters:

  • chi_e (float) – Electron thermal diffusivity [m^2/s].

  • chi_i (float) – Ion thermal diffusivity [m^2/s].

  • d_e (float) – Particle diffusivity [m^2/s].

  • channel (str) – Dominant instability channel (“ITG”, “TEM”, “ETG”, or “stable”).

  • chi_e_itg (float) – Electron transport contributions from the reduced ITG/TEM/ETG channels.

  • chi_e_tem (float) – Electron transport contributions from the reduced ITG/TEM/ETG channels.

  • chi_e_etg (float) – Electron transport contributions from the reduced ITG/TEM/ETG channels.

  • chi_i_itg (float) – Ion transport contribution from the reduced ITG channel.

chi_e: float = 0.0
chi_i: float = 0.0
d_e: float = 0.0
channel: str = 'stable'
chi_e_itg: float = 0.0
chi_e_tem: float = 0.0
chi_e_etg: float = 0.0
chi_i_itg: float = 0.0
scpn_fusion.core.neural_transport.critical_gradient_model(inp, *, stiffness=2.0)[source]

Reduced multichannel gyrokinetic closure used as analytic fallback.

This is still a reduced-order model, but it separates three dominant turbulence branches: :rtype: TransportFluxes

  • ITG: ion-driven, mostly ion heat transport with a smaller electron tail

  • TEM: trapped-electron-mode, electron-dominant with trapped-particle drive

  • ETG: electron-temperature-gradient branch, high-k electron heat channel

Parameters:
Return type:

TransportFluxes

scpn_fusion.core.neural_transport.reduced_gyrokinetic_profile_model(rho, te, ti, ne, q_profile, s_hat_profile, *, r_major=6.2, a_minor=2.0, b_toroidal=5.3)[source]

Evaluate the reduced ITG/TEM/ETG closure across a radial profile.

Return type:

tuple[ndarray[Any, dtype[float64]], ndarray[Any, dtype[float64]], ndarray[Any, dtype[float64]], dict[str, Any]]

Parameters:
class scpn_fusion.core.neural_transport.MLPWeights(layers_w=<factory>, layers_b=<factory>, input_mean=<factory>, input_std=<factory>, output_scale=<factory>, log_transform=False, gb_scale=False, gated=False)[source]

Bases: object

Stored weights for a variable-depth feedforward MLP.

Architecture: input(10) → hidden1 → [hidden2 → …] → output(3) Activation: GELU on hidden layers, softplus on output (ensures chi > 0).

Supports 2+ layer architectures. The number of layers is auto-detected from the .npz keys at load time (w1/b1, w2/b2, …). Legacy 3-layer weights (w1→w2→w3) load without modification.

Parameters:
layers_w: list[ndarray[Any, dtype[float64]]]
layers_b: list[ndarray[Any, dtype[float64]]]
input_mean: ndarray[Any, dtype[float64]]
input_std: ndarray[Any, dtype[float64]]
output_scale: ndarray[Any, dtype[float64]]
log_transform: bool = False
gb_scale: bool = False
gated: bool = False
property w1: ndarray[Any, dtype[float64]]
property b1: ndarray[Any, dtype[float64]]
property w2: ndarray[Any, dtype[float64]]
property b2: ndarray[Any, dtype[float64]]
property w3: ndarray[Any, dtype[float64]]
property b3: ndarray[Any, dtype[float64]]
property n_layers: int
class scpn_fusion.core.neural_transport.NeuralTransportModel(weights_path=None)[source]

Bases: object

Neural transport surrogate with analytic fallback.

On construction, attempts to load MLP weights from weights_path. If loading fails (file missing, wrong format), the model transparently falls back to critical_gradient_model().

Parameters:

weights_path (str or Path, optional) – Path to a .npz file containing MLP weights. The file must contain arrays w1, b1, ..., wN, bN, input_mean, input_std, output_scale. The number of layers N is auto-detected.

predict(inp)[source]

Predict turbulent transport fluxes for given local parameters.

Uses the neural MLP if weights are loaded, otherwise falls back to the analytic critical-gradient model.

Parameters:

inp (TransportInputs) – Local plasma parameters at a single radial point.

Returns:

Predicted heat and particle diffusivities.

Return type:

TransportFluxes

predict_profile(rho, te, ti, ne, q_profile, s_hat_profile, r_major=6.2, a_minor=2.0, b_toroidal=5.3)[source]

Predict transport coefficients on the full radial profile.

Computes normalised gradients from the profile arrays via central finite differences, then evaluates the surrogate at each radial point. When the MLP is loaded the entire profile is evaluated in a single batched forward pass (no Python loop).

Parameters:

  • rho (FloatArray) – Normalised radius grid (0 to 1), shape (N,).

  • te (FloatArray) – Electron/ion temperature profiles [keV], shape (N,).

  • ti (FloatArray) – Electron/ion temperature profiles [keV], shape (N,).

  • ne (FloatArray) – Electron density profile [10^19 m^-3], shape (N,).

  • q_profile (FloatArray) – Safety factor profile, shape (N,).

  • s_hat_profile (FloatArray) – Magnetic shear profile, shape (N,).

  • r_major (float) – Major radius [m] for gradient normalisation.

  • a_minor (float) – Minor radius [m].

  • b_toroidal (float)

Returns:

chi_e, chi_i, d_e – Transport coefficient profiles, each shape (N,).

Return type:

FloatArray

scpn_fusion.core.neural_transport.NeuralTransportSurrogate

alias of NeuralTransportModel

FNO Turbulence Suppressor

Legacy turbulence suppression compatibility module.

Default execution uses a deterministic reduced-order compatibility backend. The historical JAX-FNO backend is available only via explicit opt-in: allow_legacy=True or SCPN_ENABLE_LEGACY_FNO=1.

class scpn_fusion.core.fno_turbulence_suppressor.SpectralTurbulenceGenerator(size=64, *, seed=None, rng=None)[source]

Bases: object

Generates synthetic ITG turbulence with Fourier-space drift-wave dynamics.

Parameters:

  • size (int)

  • seed (Optional[int])

  • rng (Optional[np.random.Generator])

step(dt=0.01, damping=0.0)[source]
Return type:

ndarray

Parameters:
class scpn_fusion.core.fno_turbulence_suppressor.FNO_Controller(weights_path=None, *, allow_legacy=False)[source]

Bases: object

Predict turbulence suppression with explicit legacy opt-in semantics.

Default backend is deterministic compatibility mode. Legacy JAX-FNO can be enabled with allow_legacy=True or SCPN_ENABLE_LEGACY_FNO=1.

Parameters:

  • weights_path (Optional[str])

  • allow_legacy (bool)

load_weights(path)[source]
Return type:

None

Parameters:

path (str)

predict_and_suppress(field)[source]
Return type:

Tuple[float, ndarray]

Parameters:

field (ndarray)

scpn_fusion.core.fno_turbulence_suppressor.run_fno_simulation(time_steps=200, weights_path=None, *, seed=42, allow_legacy=False, save_plot=True, output_path='FNO_Turbulence_Result.png', verbose=True)[source]
Return type:

dict[str, Any]

Parameters:

  • time_steps (int)

  • weights_path (str | None)

  • seed (int)

  • allow_legacy (bool)

  • save_plot (bool)

  • output_path (str)

  • verbose (bool)

FNO Training

Pure-NumPy training for a multi-layer Fourier Neural Operator turbulence model (LEGACY).

Note

As of v3.6.0, this module is superseded by the JAX-accelerated version in fno_jax_training.py, which provides 100x faster training and higher accuracy (~0.001 loss).

scpn_fusion.core.fno_training.gelu(x)[source]

Fast GeLU approximation.

Return type:

ndarray

Parameters:

x (ndarray)

class scpn_fusion.core.fno_training.AdamOptimizer(beta1=0.9, beta2=0.999, eps=1e-08)[source]

Bases: object

Minimal Adam optimizer for NumPy arrays.

Parameters:
step(params, grads, lr)[source]
Return type:

None

Parameters:
class scpn_fusion.core.fno_training.MultiLayerFNO(modes=12, width=32, n_layers=4, seed=42)[source]

Bases: object

Multi-layer FNO model: Input [N,N] -> Lift (1->width) -> 4x FNO layers -> Project (width->1) -> [N,N].

Training routine updates the project head with Adam while keeping the spectral backbone fixed. This keeps the implementation NumPy-only and fast enough for iterative dataset generation.

Parameters:
forward_with_hidden(x_field)[source]
Return type:

Tuple[ndarray, ndarray]

Parameters:

x_field (ndarray)

forward(x_field)[source]
Return type:

ndarray

Parameters:

x_field (ndarray)

save_weights(path)[source]
Return type:

None

Parameters:

path (str | Path)

load_weights(path)[source]
Return type:

None

Parameters:

path (str | Path)

scpn_fusion.core.fno_training.train_fno(n_samples=10000, epochs=500, lr=0.001, modes=12, width=32, save_path=PosixPath('/home/runner/work/scpn-fusion-core/scpn-fusion-core/weights/fno_turbulence.npz'), batch_size=8, seed=42, patience=50)[source]

Train MultiLayerFNO with pure NumPy.

Returns a history dictionary with loss curves and saved model metadata.

Return type:

Dict[str, object]

Parameters:
scpn_fusion.core.fno_training.train_fno_multi_regime(n_samples=10000, epochs=500, lr=0.001, modes=12, width=32, save_path=PosixPath('/home/runner/work/scpn-fusion-core/scpn-fusion-core/weights/fno_turbulence_sparc.npz'), batch_size=8, seed=42, patience=50, regime_weights=None)[source]

Compatibility wrapper over extracted multi-regime training runtime.

Return type:

Dict[str, object]

Parameters:

Turbulence Oracle

class scpn_fusion.core.turbulence_oracle.DriftWavePhysics(N=64)[source]

Bases: object

2D Hasegawa-Wakatani solver for plasma edge turbulence.

Variables: phi (electrostatic potential), n (density fluctuation).

bracket(A_k, B_k)[source]

Calculates Poisson Bracket [A, B] with de-aliasing

step()[source]

Runge-Kutta 4th Order Step with Stability Clamp

class scpn_fusion.core.turbulence_oracle.OracleESN(input_dim, reservoir_size=500, spectral_radius=0.95)[source]

Bases: object

Echo State Network (Reservoir Computing). Specialized for predicting chaotic time series.

train(inputs, targets)[source]

Ridge Regression training

predict(u_current, steps=50)[source]
scpn_fusion.core.turbulence_oracle.run_turbulence_oracle()[source]

Fusion Ignition Simulator

exception scpn_fusion.core.fusion_ignition_sim.BurnPhysicsError[source]

Bases: RuntimeError, FusionCoreError

Raised when strict 0-D burn physics contracts are violated.

class scpn_fusion.core.fusion_ignition_sim.FusionBurnPhysics(config_path)[source]

Bases: FusionKernel

Extends the Grad-Shafranov Solver with Thermonuclear Physics. Calculates Fusion Power, Alpha Heating, and Q-Factor.

bosch_hale_dt(T_keV)[source]

D-T <sigma v> [m^3/s]. Bosch & Hale, NF 32 (1992) 611.

calculate_thermodynamics(P_aux_MW=50.0)[source]

Maps Magnetic Equilibrium -> Thermodynamics -> Fusion Power. P_aux_MW: External Heating Power (NBI/ECRH) in MegaWatts.

scpn_fusion.core.fusion_ignition_sim.run_ignition_experiment()[source]
class scpn_fusion.core.fusion_ignition_sim.DynamicBurnModel(R0=6.2, a=2.0, B_t=5.3, I_p=15.0, kappa=1.7, n_e20=1.0, M_eff=2.5, Z_eff=1.65)[source]

Bases: object

Self-consistent dynamic burn model with ITER98y2 confinement scaling.

Solves the coupled ODEs for plasma energy and temperature evolution:

dW/dt = P_alpha + P_aux - W/tau_E - P_rad - P_brems

with:

  • Bosch-Hale D-T reactivity

  • ITER IPB98(y,2) confinement scaling

  • Bremsstrahlung radiation losses

  • Impurity radiation (Z_eff-dependent)

  • He ash accumulation and dilution

  • H-mode power threshold (Martin scaling)

This enables finding validated Q>=10 operating points.

Parameters:

  • R0 (float) – Major radius (m).

  • a (float) – Minor radius (m).

  • B_t (float) – Toroidal field (T).

  • I_p (float) – Plasma current (MA).

  • kappa (float) – Elongation.

  • n_e20 (float) – Line-averaged electron density (10^20 m^-3).

  • M_eff (float) – Effective ion mass (amu). D-T = 2.5.

  • Z_eff (float) – Effective charge.

static bosch_hale_dt(T_keV)[source]

D-T <sigma v> [m^3/s]. Bosch & Hale, NF 32 (1992) 611.

Return type:

float

Parameters:

T_keV (float)

iter98y2_tau_e(P_loss_mw)[source]

ITER IPB98(y,2) energy confinement scaling (s).

Return type:

float

Parameters:

P_loss_mw (float)

tau_E = 0.0562 * I_p^0.93 * B_t^0.15 * n_e19^0.41 * P^-0.69 *

R^1.97 * (a/R)^0.58 * kappa^0.78 * M^0.19

h_mode_threshold_mw()[source]

H-mode power threshold (Martin 2008 scaling).

P_thr = 0.0488 * n_e20^0.717 * B_t^0.803 * S^0.941 where S is the plasma surface area.

Return type:

float

simulate(P_aux_mw=50.0, T_initial_keV=5.0, duration_s=100.0, dt_s=0.01, f_he_initial=0.02, tau_he_factor=5.0, pumping_efficiency=0.8, *, enforce_temperature_limit=False, max_temperature_clamp_events=None, warn_on_temperature_cap=True, emit_repeated_temperature_warnings=False)[source]

Run dynamic burn simulation.

Parameters:

  • P_aux_mw (float) – Auxiliary heating power (MW).

  • T_initial_keV (float) – Initial ion temperature (keV).

  • duration_s (float) – Simulation duration (s).

  • dt_s (float) – Time step (s).

  • f_he_initial (float) – Initial helium fraction.

  • tau_he_factor (float) – He confinement / energy confinement ratio.

  • pumping_efficiency (float) – He pumping efficiency (0–1).

  • enforce_temperature_limit (bool) – When True, raise BurnPhysicsError instead of clamping temperatures above 25 keV.

  • max_temperature_clamp_events (int or None) – Optional cap on number of clamp events. If exceeded, raises BurnPhysicsError.

  • warn_on_temperature_cap (bool) – Emit warning when a cap event occurs.

  • emit_repeated_temperature_warnings (bool) – Emit warning for every cap event when True; otherwise warn once.

Return type:

dict with time-histories and final metrics including Q factor.

static find_q10_operating_point(R0=6.2, a=2.0, B_t=5.3, I_p=15.0, kappa=1.7)[source]

Scan auxiliary power to find Q>=10 operating point.

Returns the scan results including the optimal P_aux for Q=10.

Return type:

dict

Parameters:

Divertor Thermal Simulation

class scpn_fusion.core.divertor_thermal_sim.DivertorLab(P_sol_MW=50.0, R_major=2.1, B_pol=2.0)[source]

Bases: object

Simulates Heat Exhaust in a Compact Fusion Reactor. Compares Solid Tungsten vs. Liquid Lithium Vapor Shielding.

Parameters:
solve_2point_transport(expansion_factor=10.0, f_rad=0.5)[source]

Two-Point Model (2PM) for SOL Transport. Balances upstream pressure with target flux constraints. T_u = (7/2 * L_c * q_par / kappa0)^(2/7) n_u determines if we are in sheath-limited or conduction-limited regime.

Return type:

tuple[float, float]

Parameters:
calculate_heat_load(expansion_factor=10.0)[source]

Calculates Peak Heat Flux using 2-Point Model Physics.

Return type:

float

Parameters:

expansion_factor (float)

simulate_tungsten()[source]

1D Thermal limit of Tungsten Monoblock. Simple conduction model: T_surf = q * d / k + T_coolant

Return type:

tuple[float, str]

simulate_lithium_vapor(*, relaxation=0.7, max_iter=50, tol=0.1)[source]

Self-Consistent Vapor Shielding Physics.

Parameters:

  • relaxation (float) – Under-relaxation factor in (0, 1) for iterative convergence.

  • max_iter (int) – Maximum Picard iterations.

  • tol (float) – Absolute temperature convergence tolerance [°C].

Return type:

tuple[float, float, float]

calculate_mhd_pressure_loss(flow_velocity_m_s, channel_length_m=1.2, channel_half_gap_m=0.012, density_kg_m3=510.0, viscosity_pa_s=0.0025, conductivity_s_m=800000.0)[source]

Reduced TEMHD pressure-loss model using a Hartmann-flow correction.

Returns pressure-loss summary for the provided channel flow speed.

Return type:

dict[str, float]

Parameters:
estimate_evaporation_rate(surface_temp_c, flow_velocity_m_s)[source]

Velocity-dependent lithium evaporation estimate [kg m^-2 s^-1].

Return type:

float

Parameters:

  • surface_temp_c (float)

  • flow_velocity_m_s (float)

simulate_temhd_liquid_metal(flow_velocity_m_s, expansion_factor=15.0)[source]

Reduced TEMHD divertor state including MHD pressure loss and evaporation.

Return type:

dict[str, object]

Parameters:

  • flow_velocity_m_s (float)

  • expansion_factor (float)

scpn_fusion.core.divertor_thermal_sim.run_divertor_sim()[source]
Return type:

None

Sandpile Fusion Reactor

class scpn_fusion.core.sandpile_fusion_reactor.TokamakSandpile(size=100)[source]

Bases: object

Implements a modified Bak-Tang-Wiesenfeld (BTW) sandpile model to simulate self-organized criticality (SOC) in plasma turbulence.

Reference: ‘Running Sandpile’ models for L-H transition.

drive()[source]

Adds energy (sand) to the core.

relax(suppression_strength=0.0)[source]

The Avalanche Step. If slope > critical, redistribute energy. suppression_strength: 0.0 to 1.0 (Effect of HJB control)

calculate_profile()[source]

Reconstructs the Height profile (Temperature) from gradients.

class scpn_fusion.core.sandpile_fusion_reactor.HJB_Avalanche_Controller[source]

Bases: object

A heuristic closed-loop controller for avalanche suppression. Observed State: Current Avalanche Size. Action: Increase/Decrease Magnetic Shear (Z_crit). Reward: Maximize Core Height (Confinement) - Action Cost.

act(last_avalanche_size, current_core_temp)[source]
scpn_fusion.core.sandpile_fusion_reactor.run_sandpile_simulation()[source]

Warm Dense Matter Engine

class scpn_fusion.core.wdm_engine.WholeDeviceModel(config_path)[source]

Bases: object

SCPN WDM: Coupled Multi-Physics Simulation. Loops: Equilibrium <-> Transport <-> Wall <-> Radiation

thomas_fermi_pressure(n_e_m3, T_eV)[source]

Hardened Thomas-Fermi Equation of State (EOS) heuristic. Accounts for electron degeneracy pressure in the WDM regime. P_total = P_ideal + P_deg

calculate_redeposition_fraction(T_edge_eV, B_field_T)[source]

Estimates the fraction of sputtered atoms that are promptly redeposited. Redeposition fraction f_redep ~ 1 - (lambda_ion / rho_L) For heavy impurities (W) in high B-field, redeposition can exceed 90%.

run_discharge(duration_sec=10.0)[source]
plot_results(history)[source]

Uncertainty Quantification

Bayesian uncertainty quantification for fusion performance predictions.

Provides error bars on confinement time, fusion power, and Q-factor by Monte Carlo sampling over scaling-law parameter uncertainties. Uses the IPB98(y,2) H-mode confinement scaling as the baseline model.

Full-chain propagation (equilibrium -> transport -> fusion) lives in uncertainty_full_chain.

References

  • ITER Physics Basis, Nucl. Fusion 39 (1999) 2175

  • Verdoolaege et al., Nucl. Fusion 61 (2021) 076006

class scpn_fusion.core.uncertainty.PlasmaScenario(I_p, B_t, P_heat, n_e, R, A, kappa, M=2.5)[source]

Bases: object

Input plasma parameters for a confinement prediction.

Parameters:
I_p: float
B_t: float
P_heat: float
n_e: float
R: float
A: float
kappa: float
M: float = 2.5
class scpn_fusion.core.uncertainty.UQResult(tau_E, P_fusion, Q, tau_E_sigma, P_fusion_sigma, Q_sigma, tau_E_percentiles=<factory>, P_fusion_percentiles=<factory>, Q_percentiles=<factory>, n_samples=0)[source]

Bases: object

Uncertainty-quantified prediction result.

Parameters:
tau_E: float
P_fusion: float
Q: float
tau_E_sigma: float
P_fusion_sigma: float
Q_sigma: float
tau_E_percentiles: ndarray
P_fusion_percentiles: ndarray
Q_percentiles: ndarray
n_samples: int = 0
scpn_fusion.core.uncertainty.ipb98_tau_e(scenario, params=None)[source]

Compute IPB98(y,2) confinement time for given plasma parameters.

Parameters:

  • scenario (PlasmaScenario) – Plasma parameters.

  • params (dict, optional) – Scaling law coefficients. Defaults to IPB98 central values.

Return type:

float — confinement time in seconds.

scpn_fusion.core.uncertainty.bosch_hale_reactivity(Ti_keV)

D-T fusion reactivity <sigma v> [m^3/s].

Bosch & Hale, Nuclear Fusion 32 (1992) 611, Table IV. Valid for 0.2 <= T <= 100 keV. Accepts scalar or array.

scpn_fusion.core.uncertainty.fusion_power_from_tau(scenario, tau_E)[source]

Estimate fusion power from cross-section integrated thermal reactivity. P_fus = n_D * n_T * <sigma v> * V * E_fusion

Includes alpha self-heating iteration: in a burning plasma, 20% of fusion energy (3.5 MeV alphas out of 17.6 MeV total) heats the plasma, raising Ti and thus reactivity. Three fixed-point iterations converge for ITER-class Q~10 scenarios.

Return type:

float

Parameters:
scpn_fusion.core.uncertainty.quantify_uncertainty(scenario, n_samples=10000, seed=None)[source]

Monte Carlo uncertainty quantification for fusion performance.

Samples scaling-law coefficients from their Gaussian posteriors and propagates through the confinement and fusion power models.

Parameters:

  • scenario (PlasmaScenario) – Plasma parameters (held fixed).

  • n_samples (int) – Number of Monte Carlo samples (default 10,000).

  • seed (int, optional) – Random seed for reproducibility.

Return type:

UQResult — central estimates + error bars + percentiles.

Geometry 3D

3D flux-surface mesh generation utilities.

This module extracts an approximate last-closed flux surface (LCFS) from the 2D Grad-Shafranov equilibrium and revolves it toroidally into a triangular mesh.

class scpn_fusion.core.geometry_3d.Reactor3DBuilder(config_path=None, *, kernel=None, equilibrium_3d=None, solve_equilibrium=True)[source]

Bases: object

Build a 3D toroidal mesh from a solved 2D equilibrium.

Parameters:
generate_plasma_surface(resolution_toroidal=60, resolution_poloidal=60, radial_steps=512)[source]

Generate a triangulated toroidal surface mesh.

Returns:

(vertices, faces), where vertices are shaped (N, 3) and faces are shaped (M, 3) with zero-based indices.

Return type:

tuple[np.ndarray, np.ndarray]

Parameters:

  • resolution_toroidal (int)

  • resolution_poloidal (int)

  • radial_steps (int)

build_vmec_like_equilibrium(*, lcfs_resolution=96, radial_steps=512, nfp=1, toroidal_modes=None)[source]

Infer a reduced VMEC-like 3D equilibrium from solved 2D LCFS.

Return type:

VMECStyleEquilibrium3D

Parameters:
build_stellarator_w7x_like_equilibrium(*, lcfs_resolution=96, radial_steps=512, nfp=5, edge_ripple=0.08, vertical_ripple=0.05)[source]

Build a reduced W7-X-like non-axisymmetric equilibrium extension.

Return type:

VMECStyleEquilibrium3D

Parameters:

  • lcfs_resolution (int)

  • radial_steps (int)

  • nfp (int)

  • edge_ripple (float)

  • vertical_ripple (float)

generate_coil_meshes()[source]

Return compact metadata for external coil geometry placeholders.

Return type:

list[dict[str, float | str]]

create_fieldline_tracer(*, rotational_transform=0.45, helical_coupling_scale=0.08, radial_coupling_scale=0.0, nfp=1, toroidal_modes=None, lcfs_resolution=96, radial_steps=512)[source]

Create a reduced 3D field-line tracer.

Return type:

FieldLineTracer3D

Parameters:
generate_poincare_map(*, rho0=0.95, theta0=0.0, phi0=0.0, toroidal_turns=20, steps_per_turn=256, phi_planes=None, rotational_transform=0.45, helical_coupling_scale=0.08, radial_coupling_scale=0.0, nfp=1, toroidal_modes=None)[source]

Trace one field line and return Poincare map in (R, Z).

Return type:

tuple[FieldLineTrace3D, dict[float, ndarray]]

Parameters:
export_obj(vertices, faces, filename='plasma.obj')[source]
Return type:

Path

Parameters:
export_preview_png(vertices, faces, filename='plasma_preview.png', *, dpi=140)[source]

Render mesh preview to a PNG using matplotlib.

Return type:

Path

Parameters:
scpn_fusion.core.geometry_3d.main(argv=None)[source]
Return type:

int

Parameters:

argv (list[str] | None)

3D Equilibrium

VMEC-like reduced 3D equilibrium interface and flux coordinates.

This module provides a compact Fourier-parameterized equilibrium surrogate that maps native flux coordinates (rho, theta, phi) into cylindrical and Cartesian geometry. It is intentionally reduced-order (not a full VMEC solver), but exposes a compatible interface for non-axisymmetric n != 0 shaping.

class scpn_fusion.core.equilibrium_3d.FourierMode3D(m, n, r_cos=0.0, r_sin=0.0, z_cos=0.0, z_sin=0.0)[source]

Bases: object

Single VMEC-like Fourier harmonic for 3D boundary shaping.

Parameters:
m: int
n: int
r_cos: float = 0.0
r_sin: float = 0.0
z_cos: float = 0.0
z_sin: float = 0.0
class scpn_fusion.core.equilibrium_3d.VMECStyleEquilibrium3D(*, r_axis, z_axis, a_minor, kappa=1.0, triangularity=0.0, nfp=1, modes=None)[source]

Bases: object

Reduced VMEC-like 3D equilibrium in native flux coordinates.

Parameters:
classmethod from_axisymmetric_lcfs(lcfs_points, *, r_axis, z_axis, nfp=1, modes=None)[source]

Infer baseline shaping from traced axisymmetric LCFS points.

Return type:

VMECStyleEquilibrium3D

Parameters:
to_vmec_like_dict()[source]

Export a deterministic VMEC-like boundary payload.

Return type:

dict[str, object]

classmethod from_vmec_like_dict(payload)[source]

Rebuild reduced equilibrium from VMEC-like boundary payload.

Return type:

VMECStyleEquilibrium3D

Parameters:

payload (dict[str, object])

flux_to_cylindrical(rho, theta, phi)[source]

Map native flux coordinates (rho, theta, phi) to (R, Z, phi).

Return type:

tuple[ndarray, ndarray, ndarray]

Parameters:
flux_to_cartesian(rho, theta, phi)[source]

Map native flux coordinates (rho, theta, phi) to (x, y, z).

Return type:

tuple[ndarray, ndarray, ndarray]

Parameters:
sample_surface(*, rho=1.0, resolution_toroidal=60, resolution_poloidal=60)[source]

Sample a triangulated flux surface from native 3D coordinates.

Return type:

tuple[ndarray, ndarray]

Parameters:

  • rho (float)

  • resolution_toroidal (int)

  • resolution_poloidal (int)

class scpn_fusion.core.equilibrium_3d.ForceBalanceResult(converged, iterations, residual_norm, initial_residual, modes, force_residual_history, armijo_reject_count=0, non_decreasing_steps=0)[source]

Bases: object

Result of reduced-order 3D force-balance solve.

Parameters:
converged: bool
iterations: int
residual_norm: float
initial_residual: float
modes: list[FourierMode3D]
force_residual_history: list[float]
armijo_reject_count: int = 0
non_decreasing_steps: int = 0
class scpn_fusion.core.equilibrium_3d.ForceBalance3D(eq, *, b0_tesla=5.3, r0_major=6.2, p0_pa=500000.0, j0_ma_m2=1.0, pressure_exp=2.0, current_exp=1.5)[source]

Bases: object

Reduced-order 3D force-balance solver via spectral variational method.

Takes a 2D Grad-Shafranov base equilibrium (Psi, pressure, current profiles) and iterates 3D Fourier mode coefficients to minimize the MHD force residual:

||J x B - nabla p||^2

evaluated on a (rho, theta, phi) sampling grid. This is NOT a full 3D VMEC solver but provides a genuine force-balance closure by coupling the Fourier boundary parameterization to the MHD equilibrium condition.

The method is a spectral Galerkin projection: compute the force residual on a 3D sampling grid, project onto each Fourier mode via inner product, and update the mode amplitudes by gradient descent with Armijo line search.

Parameters:
compute_force_residual(n_rho=12, n_theta=24, n_phi=16)[source]

Compute ||J x B - nabla p||^2 on a (rho, theta, phi) grid.

Returns the volume-averaged L2 norm of the force residual.

Return type:

float

Parameters:
solve(*, max_iterations=200, tolerance=0.0001, learning_rate=0.01, armijo_c=0.0001, n_rho=12, n_theta=24, n_phi=16)[source]

Iterate Fourier coefficients to minimise force residual.

Uses gradient descent with Armijo back-tracking line search. If the equilibrium has no modes, a default set (m=0..2, n=0..1) is seeded.

Return type:

ForceBalanceResult

Parameters:

3D Field-Line Tracing

Reduced 3D field-line tracing and Poincare diagnostics.

class scpn_fusion.core.fieldline_3d.FieldLineTrace3D(rho, theta, phi, xyz)[source]

Bases: object

Sampled field-line trajectory in native and Cartesian coordinates.

Parameters:
rho: ndarray
theta: ndarray
phi: ndarray
xyz: ndarray
class scpn_fusion.core.fieldline_3d.PoincareSection3D(phi_plane, xyz, rz)[source]

Bases: object

Poincare section points for one toroidal cut plane.

Parameters:
phi_plane: float
xyz: ndarray
rz: ndarray
class scpn_fusion.core.fieldline_3d.ToroidalAsymmetryObservables3D(n1_amp, n2_amp, n3_amp, z_n1_amp, asymmetry_index, radial_spread)[source]

Bases: object

Reduced toroidal asymmetry indicators for instability-aware control.

Parameters:
n1_amp: float
n2_amp: float
n3_amp: float
z_n1_amp: float
asymmetry_index: float
radial_spread: float
as_dict()[source]
Return type:

dict[str, float]

class scpn_fusion.core.fieldline_3d.FieldLineTracer3D(equilibrium, *, rotational_transform=0.45, helical_coupling_scale=0.08, radial_coupling_scale=0.0)[source]

Bases: object

Trace reduced field lines over a VMEC-like 3D equilibrium surface.

Parameters:
trace_line(*, rho0=0.95, theta0=0.0, phi0=0.0, toroidal_turns=16, steps_per_turn=256)[source]

Trace one reduced field line in (rho, theta, phi) coordinates.

Return type:

FieldLineTrace3D

Parameters:
poincare_section(trace, *, phi_plane=0.0)[source]

Compute Poincare intersection points for one toroidal plane.

Return type:

PoincareSection3D

Parameters:
poincare_map(trace, *, phi_planes)[source]

Compute Poincare sections for multiple toroidal planes.

Return type:

dict[float, PoincareSection3D]

Parameters:
toroidal_asymmetry_observables(trace)[source]

Extract reduced toroidal asymmetry observables from a field-line trace.

Return type:

ToroidalAsymmetryObservables3D

Parameters:

trace (FieldLineTrace3D)

GPU Runtime Bridge

Deterministic CPU vs GPU-sim bridge for multigrid and SNN inference lanes.

class scpn_fusion.core.gpu_runtime.RuntimeBenchmark(backend, multigrid_p95_ms_est, snn_p95_ms_est, multigrid_mean_ms_wall, snn_mean_ms_wall)[source]

Bases: object

Parameters:

  • backend (str)

  • multigrid_p95_ms_est (float)

  • snn_p95_ms_est (float)

  • multigrid_mean_ms_wall (float)

  • snn_mean_ms_wall (float)

backend: str
multigrid_p95_ms_est: float
snn_p95_ms_est: float
multigrid_mean_ms_wall: float
snn_mean_ms_wall: float
class scpn_fusion.core.gpu_runtime.EquilibriumLatencyBenchmark(backend, trials, grid_size, fault_runs, p95_ms_est, mean_ms_wall, p95_ms_wall, fault_p95_ms_est, fault_mean_ms_wall, fault_p95_ms_wall, sub_ms_target_pass, latency_spike_over_10ms)[source]

Bases: object

Parameters:

  • backend (str)

  • trials (int)

  • grid_size (int)

  • fault_runs (int)

  • p95_ms_est (float)

  • mean_ms_wall (float)

  • p95_ms_wall (float)

  • fault_p95_ms_est (float)

  • fault_mean_ms_wall (float)

  • fault_p95_ms_wall (float)

  • sub_ms_target_pass (bool)

  • latency_spike_over_10ms (bool)

backend: str
trials: int
grid_size: int
fault_runs: int
p95_ms_est: float
mean_ms_wall: float
p95_ms_wall: float
fault_p95_ms_est: float
fault_mean_ms_wall: float
fault_p95_ms_wall: float
sub_ms_target_pass: bool
latency_spike_over_10ms: bool
class scpn_fusion.core.gpu_runtime.GPURuntimeBridge(seed=42)[source]

Bases: object

Reduced runtime bridge with deterministic benchmark estimates.

Parameters:

seed (int)

static available_equilibrium_backends()[source]
Return type:

tuple[str, ...]

benchmark(*, backend, trials=64, grid_size=64)[source]
Return type:

RuntimeBenchmark

Parameters:
benchmark_pair(*, trials=64, grid_size=64)[source]
Return type:

dict[str, float | dict[str, float]]

Parameters:

  • trials (int)

  • grid_size (int)

benchmark_equilibrium_latency(*, backend='auto', trials=128, grid_size=64, iterations=4, fault_runs=10, sensor_noise_std=0.015, bit_flips_per_run=3, seed=42)[source]
Return type:

EquilibriumLatencyBenchmark

Parameters:

  • backend (str)

  • trials (int)

  • grid_size (int)

  • iterations (int)

  • fault_runs (int)

  • sensor_noise_std (float)

  • bit_flips_per_run (int)

  • seed (int)

Gyro-Swin Surrogate

Deterministic GyroSwin-like surrogate for core-turbulence benchmarking.

This module intentionally uses synthetic data and an offline, zero-dependency training path (NumPy only). It is scoped for CI and reproducibility while providing a measurable speed/accuracy comparison against a slower “GENE-like proxy” baseline.

class scpn_fusion.core.gyro_swin_surrogate.TurbulenceDataset(features, chi_i)[source]

Bases: object

Synthetic core-turbulence dataset bundle.

Parameters:
features: ndarray
chi_i: ndarray
class scpn_fusion.core.gyro_swin_surrogate.SpeedBenchmark(gene_proxy_s_per_sample, surrogate_s_per_sample, speedup)[source]

Bases: object

Per-sample timing comparison against a slower proxy solver.

Parameters:

  • gene_proxy_s_per_sample (float)

  • surrogate_s_per_sample (float)

  • speedup (float)

gene_proxy_s_per_sample: float
surrogate_s_per_sample: float
speedup: float
scpn_fusion.core.gyro_swin_surrogate.synthetic_core_turbulence_target(features)[source]

Reference synthetic turbulence law used as the CI benchmark target.

Return type:

ndarray

Parameters:

features (ndarray)

scpn_fusion.core.gyro_swin_surrogate.generate_synthetic_gyrokinetic_dataset(seed, samples)[source]

Generate deterministic synthetic “JET/ITPA-like” core turbulence samples.

Return type:

TurbulenceDataset

Parameters:
class scpn_fusion.core.gyro_swin_surrogate.GyroSwinLikeSurrogate(hidden_dim=64, ridge=0.0005, seed=42)[source]

Bases: object

Lightweight random-feature surrogate emulating a transformer-style map.

Parameters:
fit(features, targets)[source]
Return type:

None

Parameters:
predict(features)[source]
Return type:

ndarray

Parameters:

features (ndarray)

scpn_fusion.core.gyro_swin_surrogate.gene_proxy_predict(features, iterations=800)[source]

Slow “GENE-like” reference proxy for speed benchmarking.

The proxy intentionally performs per-sample iterative updates to emulate heavier solver cost, while remaining deterministic and bounded for CI.

Return type:

ndarray

Parameters:
scpn_fusion.core.gyro_swin_surrogate.rmse_percent(y_true, y_pred)[source]
Return type:

float

Parameters:
scpn_fusion.core.gyro_swin_surrogate.benchmark_speedup(features, surrogate, min_baseline_s=0.15, min_surrogate_s=0.02)[source]

Measure per-sample speedup of surrogate vs slow baseline proxy.

Return type:

SpeedBenchmark

Parameters:

Heat ML Shadow Surrogate

Deterministic HEAT-ML surrogate for magnetic-shadow divertor optimization.

class scpn_fusion.core.heat_ml_shadow_surrogate.ShadowDataset(features, shadow_fraction)[source]

Bases: object

Parameters:
features: ndarray
shadow_fraction: ndarray
scpn_fusion.core.heat_ml_shadow_surrogate.synthetic_shadow_reference(features)[source]

Synthetic reference law for divertor magnetic-shadow fraction.

Return type:

ndarray

Parameters:

features (ndarray)

scpn_fusion.core.heat_ml_shadow_surrogate.generate_shadow_dataset(seed, samples)[source]
Return type:

ShadowDataset

Parameters:
class scpn_fusion.core.heat_ml_shadow_surrogate.HeatMLShadowSurrogate(ridge=0.0001)[source]

Bases: object

Compact polynomial HEAT-ML surrogate with deterministic ridge fit.

Parameters:

ridge (float)

fit(features, target)[source]
Return type:

None

Parameters:
fit_synthetic(seed=42, samples=2048)[source]
Return type:

None

Parameters:
predict_shadow_fraction(features)[source]
Return type:

ndarray

Parameters:

features (ndarray)

predict_divertor_flux(q_div_baseline_w_m2, features)[source]
Return type:

ndarray

Parameters:
scpn_fusion.core.heat_ml_shadow_surrogate.rmse_percent(y_true, y_pred)[source]
Return type:

float

Parameters:
scpn_fusion.core.heat_ml_shadow_surrogate.benchmark_inference_seconds(model, samples=200000)[source]
Return type:

float

Parameters:

Rust Compatibility Layer

class scpn_fusion.core._rust_compat.RustAcceleratedKernel(config_path)[source]

Bases: object

Drop-in wrapper around Rust PyFusionKernel that mirrors the Python FusionKernel attribute interface (.Psi, .R, .Z, .RR, .ZZ, .cfg, etc.).

Delegates equilibrium solve to Rust for ~20x speedup while keeping all attribute accesses compatible with downstream code.

solve_equilibrium()[source]

Solve Grad-Shafranov equilibrium via Rust backend.

compute_b_field()[source]

Compute magnetic field components from Psi gradient.

find_x_point(Psi)[source]

Locate the null point (B=0) using local minimization. Matches Python FusionKernel.find_x_point() interface.

calculate_thermodynamics(p_aux_mw)[source]

Calculate thermodynamics via Rust backend.

calculate_vacuum_field()[source]

Compute vacuum field with Python reference implementation.

set_solver_method(method)[source]

Set inner linear solver: ‘sor’ or ‘multigrid’.

Return type:

None

Parameters:

method (str)

solver_method()[source]

Get current solver method name.

Return type:

str

save_results(filename='equilibrium_nonlinear.npz')[source]

Save current state to .npz file.

scpn_fusion.core._rust_compat.rust_shafranov_bv(*args, **kwargs)[source]
scpn_fusion.core._rust_compat.rust_solve_coil_currents(*args, **kwargs)[source]
scpn_fusion.core._rust_compat.rust_measure_magnetics(*args, **kwargs)[source]
scpn_fusion.core._rust_compat.rust_simulate_tearing_mode(*args, **kwargs)[source]
class scpn_fusion.core._rust_compat.RustSnnPool(n_neurons=50, gain=10.0, window_size=20, *, allow_numpy_fallback=True, seed=42)[source]

Bases: object

Compatibility wrapper for Rust SpikingControllerPool.

Uses the Rust implementation when available and falls back to a deterministic NumPy LIF population otherwise.

Parameters:

  • n_neurons (int) – Number of LIF neurons per sub-population (positive/negative).

  • gain (float) – Output scaling factor.

  • window_size (int) – Sliding window length for rate-code averaging.

  • allow_numpy_fallback (bool) – When False, raise ImportError if Rust extension is unavailable.

  • seed (int) – Seed used by deterministic NumPy compatibility backend.

step(error)[source]

Process error through SNN pool and return scalar control output.

Return type:

float

Parameters:

error (float)

property n_neurons: int
property gain: float
property backend: str
class scpn_fusion.core._rust_compat.RustSnnController(target_r=6.2, target_z=0.0, *, allow_numpy_fallback=True, seed=42)[source]

Bases: object

Compatibility wrapper for Rust NeuroCyberneticController.

Uses the Rust implementation when available and falls back to paired deterministic NumPy LIF pools otherwise.

Parameters:

  • target_r (float) – Target major-radius position [m].

  • target_z (float) – Target vertical position [m].

  • allow_numpy_fallback (bool) – When False, raise ImportError if Rust extension is unavailable.

  • seed (int) – Seed used by deterministic NumPy compatibility backend.

step(measured_r, measured_z)[source]

Process measured (R, Z) position and return (ctrl_R, ctrl_Z).

Return type:

tuple[float, float]

Parameters:
property target_r: float
property target_z: float
property backend: str
scpn_fusion.core._rust_compat.rust_multigrid_vcycle(source, psi_bc, r_min, r_max, z_min, z_max, nr, nz, tol=1e-06, max_cycles=500)[source]

Call Rust multigrid V-cycle if available, else return None.

Return type:

tuple of (psi, residual, n_cycles, converged), or None when Rust is unavailable.

Parameters:

Experimental Bridges

class scpn_fusion.core.lazarus_bridge.LazarusBridge(config_path)[source]

Bases: object

Connects Layer 3 (Fusion Energy) to Layer 6 (Biological Regeneration).

Hypothesis: The stable fusion plasma generates a specific electromagnetic frequency spectrum (Alpha/Theta resonance) which can imprint structural information into biological substrates (Water Memory / DNA).

Mechanism: 1. Monitor Fusion Stability (Energy Metric). 2. If Metric > Threshold (Golden Ratio Resonance), trigger synthesis. 3. Generate Opentrons Protocol to mix Reagents (TERT/SIRT6) in specific ratios derived from plasma shape.

calculate_bio_resonance()[source]

Maps Plasma Geometry to Biological Efficacy. Metric = (Elongation / 1.618) * (Triangularity / 0.3) * confinement_time

generate_protocol(resonance_score)[source]

Creates an Opentrons Python Protocol based on the physics state.

run_bridge_simulation()[source]
visualize_bridge(score)[source]
class scpn_fusion.core.vibrana_bridge.VibranaFusionBridge(config_path)[source]

Bases: object

Layer 3 (Fusion) -> Layer 5 (VIBRANA Audio) Bridge.

“The Song of the Sun”

Mechanism: Real-time Sonification of Tokamak Physics. - Plasma Current (Ip) -> Carrier Frequency (Pitch) - Stability Error (dR) -> Chaos Intensity (Dissonance) - Confinement (Tau) -> Binaural Beat (Brain State) - Magnetic Flux (Psi) -> Harmonic Series (Timbre)

This provides an intuitive “Auditory Dashboard” for the operator. Stable Plasma = Harmonic, Angelic Sound (Solfeggio/Golden Ratio). Unstable Plasma = Chaotic, Dissonant Screeching (Lorenz Attractor).

map_physics_to_audio(t, Ip, err_R, err_Z, psi_matrix)[source]

The Transduction Function: Physics -> Acoustics

run_sonification_session(duration_steps=100)[source]
visualize_soundscape(log)[source]
scpn_fusion.core.quantum_bridge.run_quantum_suite(*, base_path=None, script_timeout_seconds=1800.0)[source]
Return type:

dict[str, object]

Parameters:

  • base_path (str | Path | None)

  • script_timeout_seconds (float)

Current Diffusion

scpn_fusion.core.current_diffusion.neoclassical_resistivity(Te_keV, ne_19, Z_eff, epsilon, q=1.0, R0=2.0)[source]

Sauter neoclassical parallel resistivity [Ohm-m].

Return type:

float

Parameters:
scpn_fusion.core.current_diffusion.q_from_psi(rho, psi, R0, a, B0)[source]

q(rho) = -rho a^2 B0 / (R0 dpsi/drho), L’Hopital at axis.

Return type:

ndarray

Parameters:
scpn_fusion.core.current_diffusion.resistive_diffusion_time(a, eta)[source]

tau_R = mu0 a^2 / eta [seconds].

Return type:

float

Parameters:
class scpn_fusion.core.current_diffusion.CurrentDiffusionSolver(rho, R0, a, B0)[source]

Bases: object

Implicit Crank-Nicolson solver for 1D poloidal flux diffusion.

Parameters:
step(dt, Te, ne, Z_eff, j_bs, j_cd, j_ext=None)[source]

Advance poloidal flux psi by one timestep dt.

Return type:

ndarray

Parameters:

Current Drive Sources

class scpn_fusion.core.current_drive.ECCDSource(P_ec_MW, rho_dep, sigma_rho, eta_cd=0.03)[source]

Bases: object

Electron Cyclotron Current Drive (ECCD).

Parameters:
P_absorbed(rho)[source]

Absorbed power density profile [W/m^3].

Return type:

ndarray

Parameters:

rho (ndarray)

j_cd(rho, ne_19, Te_keV)[source]

Driven current density profile [A/m^2].

Return type:

ndarray

Parameters:
class scpn_fusion.core.current_drive.NBISource(P_nbi_MW, E_beam_keV, rho_tangency, sigma_rho=0.15)[source]

Bases: object

Neutral Beam Injection (NBI).

Parameters:
P_heating(rho)[source]

Heating power deposition profile [W/m^3].

Return type:

ndarray

Parameters:

rho (ndarray)

j_cd(rho, ne_19, Te_keV, Ti_keV)[source]

Driven current density profile [A/m^2].

Return type:

ndarray

Parameters:
class scpn_fusion.core.current_drive.LHCDSource(P_lh_MW, rho_dep, sigma_rho, eta_cd=0.15)[source]

Bases: object

Lower Hybrid Current Drive (LHCD).

Parameters:
P_absorbed(rho)[source]

Absorbed power density profile [W/m^3], Gaussian in rho.

Return type:

ndarray

Parameters:

rho (ndarray)

j_cd(rho, ne_19, Te_keV)[source]

Driven current density [A/m^2]. j_cd = eta_cd * P_abs / (ne * Te).

Return type:

ndarray

Parameters:
class scpn_fusion.core.current_drive.CurrentDriveMix(a=1.0)[source]

Bases: object

Combines multiple current drive sources.

Parameters:

a (float)

add_source(source)[source]

Register a current drive source (ECCD, NBI, or LHCD).

Return type:

None

Parameters:

source (ECCDSource | NBISource | LHCDSource)

total_j_cd(rho, ne, Te, Ti)[source]

Sum driven current density [A/m^2] from all registered sources.

Return type:

ndarray

Parameters:
total_heating_power(rho)[source]

Sum heating power density [W/m^3] from all registered sources.

Return type:

ndarray

Parameters:

rho (ndarray)

total_driven_current(rho, ne, Te, Ti)[source]

Integrate total driven current [A] assuming circular cross-section.

Return type:

float

Parameters:

Sawtooth Dynamics

class scpn_fusion.core.sawtooth.SawtoothEvent(crash_time, rho_1, rho_mix, T_drop, seed_energy)[source]

Bases: object

Record of a single sawtooth crash: timing, radii, and energy release.

Parameters:
crash_time: float
rho_1: float
rho_mix: float
T_drop: float
seed_energy: float
class scpn_fusion.core.sawtooth.SawtoothMonitor(rho, s_crit=0.1)[source]

Bases: object

Monitors the q-profile for Porcelli-like sawtooth trigger conditions.

Parameters:

  • rho (np.ndarray)

  • s_crit (float)

find_q1_radius(q)[source]

Find the normalized radius where q=1.

Return type:

float | None

Parameters:

q (ndarray)

check_trigger(q, shear)[source]

Porcelli-like trigger based on shear at q=1.

Return type:

bool

Parameters:
scpn_fusion.core.sawtooth.kadomtsev_crash(rho, T, n, q, R0, a)[source]

Apply Kadomtsev reconnection. Returns (T_new, n_new, q_new, rho_1, rho_mix).

Return type:

tuple[ndarray, ndarray, ndarray, float, float]

Parameters:
class scpn_fusion.core.sawtooth.SawtoothCycler(rho, R0, a, s_crit=0.1)[source]

Bases: object

Tracks and triggers sawtooth crashes.

Parameters:
step(dt, q, shear, T, n)[source]

Advance by dt; trigger Kadomtsev crash if shear at q=1 exceeds s_crit.

Return type:

SawtoothEvent | None

Parameters:

NTM Island Dynamics

class scpn_fusion.core.ntm_dynamics.RationalSurface(rho, r_s, m, n, q, shear)[source]

Bases: object

Represents a q=m/n rational surface.

Parameters:
rho: float
r_s: float
m: int
n: int
q: float
shear: float
scpn_fusion.core.ntm_dynamics.eccd_stabilization_factor(d_cd, w)[source]

ECCD stabilization efficiency: f = (w/d_cd) exp(-w^2/(4 d_cd^2)).

Return type:

float

Parameters:
scpn_fusion.core.ntm_dynamics.find_rational_surfaces(q, rho, a, m_max=5, n_max=3)[source]

Locate radii where q(rho) = m/n.

Return type:

list[RationalSurface]

Parameters:
class scpn_fusion.core.ntm_dynamics.NTMIslandDynamics(r_s, m, n, a, R0, B0, Delta_prime_0=None)[source]

Bases: object

Solves the Modified Rutherford Equation (MRE) for an NTM island.

Parameters:
delta_prime_model(w)[source]

Classical stability index with finite-width correction.

Return type:

float

Parameters:

w (float)

dw_dt(w, j_bs, j_phi, j_cd, eta, w_d=0.001, w_pol=0.0005, d_cd=0.05)[source]

Evaluate the RHS of the Modified Rutherford Equation.

Return type:

float

Parameters:
evolve(w0, t_span, dt, j_bs, j_phi, j_cd, eta, w_d=0.001, w_pol=0.0005, d_cd=0.05)[source]

Integrate w(t) using RK4.

Return type:

tuple[ndarray, ndarray]

Parameters:
class scpn_fusion.core.ntm_dynamics.NTMController(w_onset=0.02, w_target=0.005)[source]

Bases: object

Monitors and controls NTM islands via ECCD.

Parameters:
step(w, rho_rs, max_power=20.0)[source]

Returns requested ECCD power in MW.

Return type:

float

Parameters:

SOL Two-Point Model

class scpn_fusion.core.sol_model.SOLSolution(T_upstream_eV, T_target_eV, n_target_19, q_parallel_MW_m2, lambda_q_mm)[source]

Bases: object

Two-point SOL solution: upstream/target conditions and heat-flux width.

Parameters:
T_upstream_eV: float
T_target_eV: float
n_target_19: float
q_parallel_MW_m2: float
lambda_q_mm: float
scpn_fusion.core.sol_model.eich_heat_flux_width(P_SOL_MW, R0, B_pol, epsilon)[source]

Eich scaling: lambda_q [mm] = 1.35 P^-0.02 R0^0.04 Bpol^-0.92 eps^0.42.

Return type:

float

Parameters:
scpn_fusion.core.sol_model.peak_target_heat_flux(P_SOL_MW, R0, lambda_q_m, f_expansion=5.0, alpha_deg=3.0)[source]

Peak target heat flux [MW/m^2].

Return type:

float

Parameters:
class scpn_fusion.core.sol_model.TwoPointSOL(R0, a, q95, B_pol, kappa=1.0)[source]

Bases: object

Spitzer-Harm two-point model with Eich heat-flux width scaling.

Parameters:
solve(P_SOL_MW, n_u_19, f_rad=0.0)[source]

Solve for target temperature, density, and parallel heat flux.

Return type:

SOLSolution

Parameters: