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)

Bases: object

External coil set for free-boundary solve.

Parameters:

• positions (list[tuple[float, float]])

• currents (ndarray[Any, dtype[float64]])

• turns (list[int])

• current_limits (ndarray[Any, dtype[float64]] | None)

• target_flux_points (ndarray[Any, dtype[float64]] | None)

• target_flux_values (ndarray[Any, dtype[float64]] | None)

• x_point_target (ndarray[Any, dtype[float64]] | None)

• x_point_flux_target (float | None)

• divertor_strike_points (ndarray[Any, dtype[float64]] | None)

• divertor_flux_values (ndarray[Any, dtype[float64]] | None)

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)

Bases: FusionKernelFreeBoundaryMixin, 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)

Load reactor configuration from a JSON file.

Parameters:

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

Return type:

None

initialize_grid()

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

Return type:

None

setup_accelerator()

Initialise the optional C++ HPC acceleration bridge.

Return type:

None

calculate_vacuum_field()

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)

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)

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)

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()

Derive the magnetic field components from the solved Psi.

Return type:

None

save_results(filename='equilibrium_nonlinear.npz')

Save the equilibrium state to a compressed NumPy archive.

Parameters:

filename (str) – Output file path.

Return type:

None

phase_sync_step(theta, omega, dt, psi_driver=None, **kwargs)

Evolve oscillator phases using the Kuramoto-Sakaguchi engine.

Return type:

dict[str, Any]

Parameters:

• theta (ndarray[Any, dtype[float64]])

• omega (ndarray[Any, dtype[float64]])

• dt (float)

• psi_driver (float | None)

• kwargs (Any)

phase_sync_step_lyapunov(theta, omega, n_steps, dt, psi_driver=None, **kwargs)

Run a batch of phase sync steps and return Lyapunov stability diagnostics.

Return type:

dict[str, Any]

Parameters:

• theta (ndarray[Any, dtype[float64]])

• omega (ndarray[Any, dtype[float64]])

• n_steps (int)

• dt (float)

• psi_driver (float | None)

• kwargs (Any)

initialize_from_frc(frc_state, kappa=None)

Interpolate a 1D FRC equilibrium onto the 2D Grad-Shafranov grid.

Parameters:

• frc_state (FRCEquilibriumState) – The radial equilibrium state solved by the rigid-rotor engine.

• kappa (float | None) – Separatrix elongation (Z/R). Falls back to config target if omitted.

Return type:

None

FRC Rigid-Rotor Equilibrium

Field-reversed-configuration rigid-rotor equilibrium helpers.

This module implements the Steinhauer no-rotation analytical limit for the FRC rigid-rotor workstream. Rotating BVP support is intentionally fail-closed until the dedicated FUS-C.1 BVP implementation lands.

class scpn_fusion.core.frc_rigid_rotor.RigidRotorFRCInputs(n0, T_i_eV, T_e_eV, theta_dot, R_s, B_ext, delta=None)

Bases: object

Physical inputs for the FRC rigid-rotor analytical equilibrium.

Parameters:

• n0 (float)

• T_i_eV (float)

• T_e_eV (float)

• theta_dot (float)

• R_s (float)

• B_ext (float)

• delta (float | None)

n0: float

T_i_eV: float

T_e_eV: float

theta_dot: float

R_s: float

B_ext: float

delta: float | None = None

class scpn_fusion.core.frc_rigid_rotor.FRCEquilibriumState(rho, psi, psi_normalized, B_z, B_theta, J_theta, p, density_m3, beta, R_null, target_separatrix_radius_m, separatrix_radius_error_m, separatrix_index, field_reversal_passed, s_parameter, energy_J, converged, residual, delta, psi_axis_Wb, psi_separatrix_Wb, psi_normalized_axis_error, psi_normalized_separatrix, psi_normalized_separatrix_error, psi_normalized_residual_linf, psi_normalized_monotonic_passed, psi_normalized_bounds_passed, pressure_balance_ratio, pressure_balance_residual, pressure_balance_residual_linf, pressure_balance_residual_l2, pressure_gradient_analytic_Pa_m, pressure_gradient_residual, pressure_gradient_residual_linf, pressure_gradient_residual_l2, peak_pressure_pa, density_peak_m3, input_density_m3, central_density_residual_m3, central_density_relative_error, beta_peak, beta_separatrix_average, particle_line_density_m1, separatrix_pressure_energy_J_m, separatrix_magnetic_deficit_energy_J_m, separatrix_energy_closure_relative_error, input_thermal_pressure_pa, thermal_pressure_ratio, flux_derivative_residual, flux_derivative_residual_linf, flux_derivative_residual_l2, ampere_residual, ampere_residual_linf, ampere_residual_l2, peak_j_theta_A_m2, separatrix_bz_gradient_T_m, separatrix_expected_bz_gradient_T_m, separatrix_gradient_relative_error, separatrix_current_density_A_m2, separatrix_expected_current_density_A_m2, separatrix_current_density_relative_error, sheet_current_integral_A_m, expected_sheet_current_integral_A_m, sheet_current_integral_relative_error, force_balance_residual, force_balance_residual_linf, force_balance_residual_l2, model)

Bases: object

Radial FRC equilibrium state returned by solve_frc_equilibrium().

Parameters:

• rho (ndarray[Any, dtype[float64]])

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

• psi_normalized (ndarray[Any, dtype[float64]])

• B_z (ndarray[Any, dtype[float64]])

• B_theta (ndarray[Any, dtype[float64]])

• J_theta (ndarray[Any, dtype[float64]])

• p (ndarray[Any, dtype[float64]])

• density_m3 (ndarray[Any, dtype[float64]])

• beta (ndarray[Any, dtype[float64]])

• R_null (float)

• target_separatrix_radius_m (float)

• separatrix_radius_error_m (float)

• separatrix_index (int)

• field_reversal_passed (bool)

• s_parameter (float)

• energy_J (float)

• converged (bool)

• residual (float)

• delta (float)

• psi_axis_Wb (float)

• psi_separatrix_Wb (float)

• psi_normalized_axis_error (float)

• psi_normalized_separatrix (float)

• psi_normalized_separatrix_error (float)

• psi_normalized_residual_linf (float)

• psi_normalized_monotonic_passed (bool)

• psi_normalized_bounds_passed (bool)

• pressure_balance_ratio (float)

• pressure_balance_residual (ndarray[Any, dtype[float64]])

• pressure_balance_residual_linf (float)

• pressure_balance_residual_l2 (float)

• pressure_gradient_analytic_Pa_m (ndarray[Any, dtype[float64]])

• pressure_gradient_residual (ndarray[Any, dtype[float64]])

• pressure_gradient_residual_linf (float)

• pressure_gradient_residual_l2 (float)

• peak_pressure_pa (float)

• density_peak_m3 (float)

• input_density_m3 (float)

• central_density_residual_m3 (float)

• central_density_relative_error (float)

• beta_peak (float)

• beta_separatrix_average (float)

• particle_line_density_m1 (float)

• separatrix_pressure_energy_J_m (float)

• separatrix_magnetic_deficit_energy_J_m (float)

• separatrix_energy_closure_relative_error (float)

• input_thermal_pressure_pa (float)

• thermal_pressure_ratio (float)

• flux_derivative_residual (ndarray[Any, dtype[float64]])

• flux_derivative_residual_linf (float)

• flux_derivative_residual_l2 (float)

• ampere_residual (ndarray[Any, dtype[float64]])

• ampere_residual_linf (float)

• ampere_residual_l2 (float)

• peak_j_theta_A_m2 (float)

• separatrix_bz_gradient_T_m (float)

• separatrix_expected_bz_gradient_T_m (float)

• separatrix_gradient_relative_error (float)

• separatrix_current_density_A_m2 (float)

• separatrix_expected_current_density_A_m2 (float)

• separatrix_current_density_relative_error (float)

• sheet_current_integral_A_m (float)

• expected_sheet_current_integral_A_m (float)

• sheet_current_integral_relative_error (float)

• force_balance_residual (ndarray[Any, dtype[float64]])

• force_balance_residual_linf (float)

• force_balance_residual_l2 (float)

• model (str)

rho: ndarray[Any, dtype[float64]]

psi: ndarray[Any, dtype[float64]]

psi_normalized: ndarray[Any, dtype[float64]]

B_z: ndarray[Any, dtype[float64]]

B_theta: ndarray[Any, dtype[float64]]

J_theta: ndarray[Any, dtype[float64]]

p: ndarray[Any, dtype[float64]]

density_m3: ndarray[Any, dtype[float64]]

beta: ndarray[Any, dtype[float64]]

R_null: float

target_separatrix_radius_m: float

separatrix_radius_error_m: float

separatrix_index: int

field_reversal_passed: bool

s_parameter: float

energy_J: float

converged: bool

residual: float

delta: float

psi_axis_Wb: float

psi_separatrix_Wb: float

psi_normalized_axis_error: float

psi_normalized_separatrix: float

psi_normalized_separatrix_error: float

psi_normalized_residual_linf: float

psi_normalized_monotonic_passed: bool

psi_normalized_bounds_passed: bool

pressure_balance_ratio: float

pressure_balance_residual: ndarray[Any, dtype[float64]]

pressure_balance_residual_linf: float

pressure_balance_residual_l2: float

pressure_gradient_analytic_Pa_m: ndarray[Any, dtype[float64]]

pressure_gradient_residual: ndarray[Any, dtype[float64]]

pressure_gradient_residual_linf: float

pressure_gradient_residual_l2: float

peak_pressure_pa: float

density_peak_m3: float

input_density_m3: float

central_density_residual_m3: float

central_density_relative_error: float

beta_peak: float

beta_separatrix_average: float

particle_line_density_m1: float

separatrix_pressure_energy_J_m: float

separatrix_magnetic_deficit_energy_J_m: float

separatrix_energy_closure_relative_error: float

input_thermal_pressure_pa: float

thermal_pressure_ratio: float

flux_derivative_residual: ndarray[Any, dtype[float64]]

flux_derivative_residual_linf: float

flux_derivative_residual_l2: float

ampere_residual: ndarray[Any, dtype[float64]]

ampere_residual_linf: float

ampere_residual_l2: float

peak_j_theta_A_m2: float

separatrix_bz_gradient_T_m: float

separatrix_expected_bz_gradient_T_m: float

separatrix_gradient_relative_error: float

separatrix_current_density_A_m2: float

separatrix_expected_current_density_A_m2: float

separatrix_current_density_relative_error: float

sheet_current_integral_A_m: float

expected_sheet_current_integral_A_m: float

sheet_current_integral_relative_error: float

force_balance_residual: ndarray[Any, dtype[float64]]

force_balance_residual_linf: float

force_balance_residual_l2: float

model: str

class scpn_fusion.core.frc_rigid_rotor.FRCValidationReport(finite, monotonic_grid, null_error_m, target_separatrix_radius_m, field_reversal_passed, psi_axis_Wb, psi_separatrix_Wb, psi_normalized_axis_error, psi_normalized_separatrix, psi_normalized_separatrix_error, psi_normalized_residual_linf, psi_normalized_monotonic_passed, psi_normalized_bounds_passed, psi_normalized_passed, pressure_peak_error_m, edge_field_error_T, pressure_balance_ratio, pressure_balance_residual_linf, pressure_balance_residual_l2, pressure_balance_passed, pressure_gradient_residual_linf, pressure_gradient_residual_l2, pressure_gradient_passed, thermal_pressure_ratio, density_peak_m3, input_density_m3, central_density_residual_m3, central_density_relative_error, density_consistency_passed, beta_peak, beta_separatrix_average, particle_line_density_m1, beta_limit_passed, separatrix_pressure_energy_J_m, separatrix_magnetic_deficit_energy_J_m, separatrix_energy_closure_relative_error, energy_inventory_passed, flux_derivative_residual_linf, flux_derivative_residual_l2, flux_closure_passed, ampere_residual_linf, ampere_residual_l2, ampere_closure_passed, separatrix_bz_gradient_T_m, separatrix_expected_bz_gradient_T_m, separatrix_gradient_relative_error, separatrix_current_density_A_m2, separatrix_expected_current_density_A_m2, separatrix_current_density_relative_error, sheet_current_integral_A_m, expected_sheet_current_integral_A_m, sheet_current_integral_relative_error, sheet_current_passed, current_sheet_passed, force_balance_residual_linf, force_balance_residual_l2, force_balance_passed, passed)

Bases: object

Validation diagnostics for the analytical FRC state.

Parameters:

• finite (bool)

• monotonic_grid (bool)

• null_error_m (float)

• target_separatrix_radius_m (float)

• field_reversal_passed (bool)

• psi_axis_Wb (float)

• psi_separatrix_Wb (float)

• psi_normalized_axis_error (float)

• psi_normalized_separatrix (float)

• psi_normalized_separatrix_error (float)

• psi_normalized_residual_linf (float)

• psi_normalized_monotonic_passed (bool)

• psi_normalized_bounds_passed (bool)

• psi_normalized_passed (bool)

• pressure_peak_error_m (float)

• edge_field_error_T (float)

• pressure_balance_ratio (float)

• pressure_balance_residual_linf (float)

• pressure_balance_residual_l2 (float)

• pressure_balance_passed (bool)

• pressure_gradient_residual_linf (float)

• pressure_gradient_residual_l2 (float)

• pressure_gradient_passed (bool)

• thermal_pressure_ratio (float)

• density_peak_m3 (float)

• input_density_m3 (float)

• central_density_residual_m3 (float)

• central_density_relative_error (float)

• density_consistency_passed (bool)

• beta_peak (float)

• beta_separatrix_average (float)

• particle_line_density_m1 (float)

• beta_limit_passed (bool)

• separatrix_pressure_energy_J_m (float)

• separatrix_magnetic_deficit_energy_J_m (float)

• separatrix_energy_closure_relative_error (float)

• energy_inventory_passed (bool)

• flux_derivative_residual_linf (float)

• flux_derivative_residual_l2 (float)

• flux_closure_passed (bool)

• ampere_residual_linf (float)

• ampere_residual_l2 (float)

• ampere_closure_passed (bool)

• separatrix_bz_gradient_T_m (float)

• separatrix_expected_bz_gradient_T_m (float)

• separatrix_gradient_relative_error (float)

• separatrix_current_density_A_m2 (float)

• separatrix_expected_current_density_A_m2 (float)

• separatrix_current_density_relative_error (float)

• sheet_current_integral_A_m (float)

• expected_sheet_current_integral_A_m (float)

• sheet_current_integral_relative_error (float)

• sheet_current_passed (bool)

• current_sheet_passed (bool)

• force_balance_residual_linf (float)

• force_balance_residual_l2 (float)

• force_balance_passed (bool)

• passed (bool)

finite: bool

monotonic_grid: bool

null_error_m: float

target_separatrix_radius_m: float

field_reversal_passed: bool

psi_axis_Wb: float

psi_separatrix_Wb: float

psi_normalized_axis_error: float

psi_normalized_separatrix: float

psi_normalized_separatrix_error: float

psi_normalized_residual_linf: float

psi_normalized_monotonic_passed: bool

psi_normalized_bounds_passed: bool

psi_normalized_passed: bool

pressure_peak_error_m: float

edge_field_error_T: float

pressure_balance_ratio: float

pressure_balance_residual_linf: float

pressure_balance_residual_l2: float

pressure_balance_passed: bool

pressure_gradient_residual_linf: float

pressure_gradient_residual_l2: float

pressure_gradient_passed: bool

thermal_pressure_ratio: float

density_peak_m3: float

input_density_m3: float

central_density_residual_m3: float

central_density_relative_error: float

density_consistency_passed: bool

beta_peak: float

beta_separatrix_average: float

particle_line_density_m1: float

beta_limit_passed: bool

separatrix_pressure_energy_J_m: float

separatrix_magnetic_deficit_energy_J_m: float

separatrix_energy_closure_relative_error: float

energy_inventory_passed: bool

flux_derivative_residual_linf: float

flux_derivative_residual_l2: float

flux_closure_passed: bool

ampere_residual_linf: float

ampere_residual_l2: float

ampere_closure_passed: bool

separatrix_bz_gradient_T_m: float

separatrix_expected_bz_gradient_T_m: float

separatrix_gradient_relative_error: float

separatrix_current_density_A_m2: float

separatrix_expected_current_density_A_m2: float

separatrix_current_density_relative_error: float

sheet_current_integral_A_m: float

expected_sheet_current_integral_A_m: float

sheet_current_integral_relative_error: float

sheet_current_passed: bool

current_sheet_passed: bool

force_balance_residual_linf: float

force_balance_residual_l2: float

force_balance_passed: bool

passed: bool

scpn_fusion.core.frc_rigid_rotor.rotating_frc_bvp_acceptance_status()

Return the fail-closed acceptance status for the unresolved rotating BVP.

This is a machine-readable claim boundary for B.8/FUS-C.1. It deliberately reports the current accepted no-rotation contract and the exact external reference still required before nonzero theta_dot support can be certified.

Return type:

dict[str, object]

scpn_fusion.core.frc_rigid_rotor.ion_gyroradius_m(T_i_eV, B_T, *, mass_amu=2.014)

Return thermal ion gyroradius in metres using sqrt(2 m_i T_i) / (e B).

Return type:

float

Parameters:

• T_i_eV (float)

• B_T (float)

• mass_amu (float)

scpn_fusion.core.frc_rigid_rotor.solve_frc_equilibrium(inputs, rho_grid, *, solver='numpy', tolerance=1e-10, max_iter=200)

Solve the fail-closed Steinhauer no-rotation rigid-rotor FRC equilibrium.

The accepted production path covers the analytical theta_dot == 0 limit on a finite, strictly increasing radial grid that starts at the magnetic axis and extends beyond the requested separatrix radius. Rotating theta_dot != 0 equilibria intentionally raise NotImplementedError until the full BVP derivation and reference parity gate are verified.

Return type:

FRCEquilibriumState

Parameters:

• inputs (RigidRotorFRCInputs)

• rho_grid (ndarray[Any, dtype[float64]])

• solver (Literal['numpy', 'rust'])

• tolerance (float)

• max_iter (int)

scpn_fusion.core.frc_rigid_rotor.frc_no_rotation_jax_observables(rho_normalized_grid, *, n0, T_i_eV, T_e_eV, R_s, B_ext, delta, mass_amu=2.014)

Return differentiable observables for the accepted no-rotation FRC contract.

The independent grid is x = r / R_s. Keeping the grid normalised makes gradients with respect to R_s well-defined because the separatrix interval remains the fixed domain 0 <= x <= 1. The implemented equations are the same Steinhauer no-rotation field, cylindrical flux primitive, magnetic-pressure-balance profile, and Eq. 27 s integral used by the NumPy and Rust solver paths.

This helper intentionally does not implement the rotating rigid-rotor BVP.

Return type:

dict[str, Any]

Parameters:

• rho_normalized_grid (ndarray[Any, dtype[float64]])

• n0 (Any)

• T_i_eV (float)

• T_e_eV (float)

• R_s (Any)

• B_ext (Any)

• delta (Any | None)

• mass_amu (float)

scpn_fusion.core.frc_rigid_rotor.null_radius(state)

Return the interpolated radius where the axial field reverses.

Return type:

float

Parameters:

state (FRCEquilibriumState)

scpn_fusion.core.frc_rigid_rotor.s_parameter(state, mass_amu=2.014)

Return the Steinhauer Eq. 27 s parameter carried by this FRC state.

Return type:

float

Parameters:

• state (FRCEquilibriumState)

• mass_amu (float)

scpn_fusion.core.frc_rigid_rotor.force_balance_residual(state)

Return radial dp/dr - (J x B)_r force-balance residual in Pa/m.

Return type:

ndarray[Any, dtype[float64]]

Parameters:

state (FRCEquilibriumState)

scpn_fusion.core.frc_rigid_rotor.ampere_residual(state)

Return radial Ampere closure residual mu_0 J_theta + dB_z/dr in T/m.

Return type:

ndarray[Any, dtype[float64]]

Parameters:

state (FRCEquilibriumState)

scpn_fusion.core.frc_rigid_rotor.flux_derivative_residual(state)

Return cylindrical flux closure residual dpsi/dr - r B_z in T m.

Return type:

ndarray[Any, dtype[float64]]

Parameters:

state (FRCEquilibriumState)

scpn_fusion.core.frc_rigid_rotor.psi_normalized_profile(state)

Return separatrix-normalised cylindrical flux psi_N with psi_N(R_s)=1.

Return type:

ndarray[Any, dtype[float64]]

Parameters:

state (FRCEquilibriumState)

scpn_fusion.core.frc_rigid_rotor.pressure_balance_residual(state)

Return local pressure-balance residual p + B_z^2/(2 mu_0) - B_ext^2/(2 mu_0).

Return type:

ndarray[Any, dtype[float64]]

Parameters:

state (FRCEquilibriumState)

scpn_fusion.core.frc_rigid_rotor.pressure_gradient_residual(state)

Return finite-grid dp/dr minus analytical magnetic-pressure gradient.

Return type:

ndarray[Any, dtype[float64]]

Parameters:

state (FRCEquilibriumState)

scpn_fusion.core.frc_rigid_rotor.density_profile(state)

Return solved density profile n(r) = p(r) / ((T_i + T_e) e) in m^-3.

Return type:

ndarray[Any, dtype[float64]]

Parameters:

state (FRCEquilibriumState)

scpn_fusion.core.frc_rigid_rotor.beta_profile(state)

Return local plasma beta profile relative to the external axial-field pressure.

Return type:

ndarray[Any, dtype[float64]]

Parameters:

state (FRCEquilibriumState)

scpn_fusion.core.frc_rigid_rotor.validate_equilibrium(state, *, tolerance=1e-06, flux_tolerance=0.02, pressure_balance_tolerance=0.02, pressure_gradient_tolerance=0.02, density_tolerance=0.02, beta_limit_tolerance=0.02, energy_inventory_tolerance=1e-10, ampere_tolerance=0.02, current_sheet_tolerance=0.02, sheet_current_tolerance=0.02, psi_normalized_tolerance=1e-10, force_balance_tolerance=None)

Validate finite values, null placement, pressure peaking, and optional force balance.

Return type:

FRCValidationReport

Parameters:

• state (FRCEquilibriumState)

• tolerance (float)

• flux_tolerance (float)

• pressure_balance_tolerance (float)

• pressure_gradient_tolerance (float)

• density_tolerance (float)

• beta_limit_tolerance (float)

• energy_inventory_tolerance (float)

• ampere_tolerance (float)

• current_sheet_tolerance (float)

• sheet_current_tolerance (float)

• psi_normalized_tolerance (float)

• force_balance_tolerance (float | None)

Free-Boundary Kernel Adapters

Free-boundary and coil-optimisation helpers for FusionKernel.

scpn_fusion.core.fusion_kernel_free_boundary.build_mutual_inductance_matrix(kernel, coils, obs_points)

Build the coil-to-point flux-response matrix M[coil, point].

Return type:

ndarray[Any, dtype[float64]]

Parameters:

• kernel (Any)

• coils (CoilSet)

• obs_points (FloatArray)

scpn_fusion.core.fusion_kernel_free_boundary.compute_external_flux(kernel, coils)

Sum external coil flux on the kernel’s full (Z, R) grid.

Return type:

ndarray[Any, dtype[float64]]

Parameters:

• kernel (Any)

• coils (CoilSet)

scpn_fusion.core.fusion_kernel_free_boundary.green_function(R_src, Z_src, R_obs, Z_obs)

Evaluate the axisymmetric circular-filament poloidal-flux Green’s function.

Parameters are cylindrical source and observation coordinates in metres. The return value is the flux linkage contribution per ampere-turn. The singular self-observation limit is regularised to zero because this helper is used for external coil-to-grid coupling rather than coil self-inductance.

Return type:

float

Parameters:

• R_src (float)

• Z_src (float)

• R_obs (float)

• Z_obs (float)

scpn_fusion.core.fusion_kernel_free_boundary.interp_psi(kernel, R_pt, Z_pt)

Interpolate the kernel flux grid at an arbitrary cylindrical point.

Return type:

float

Parameters:

• kernel (Any)

• R_pt (float)

• Z_pt (float)

scpn_fusion.core.fusion_kernel_free_boundary.optimize_coil_currents(kernel, coils, target_flux, tikhonov_alpha=0.0001)

Find bounded coil currents that reproduce target flux at control points.

Return type:

ndarray[Any, dtype[float64]]

Parameters:

• kernel (Any)

• coils (CoilSet)

• target_flux (ndarray[Any, dtype[float64]])

• tikhonov_alpha (float)

scpn_fusion.core.fusion_kernel_free_boundary.reconstruct_boundary_flux_from_coils(kernel, coils, *, boundary_points, limiter_points=None, axis_point=None, x_points=None, target_flux=None)

Reconstruct free-boundary contour flux directly from coil Green functions.

This is the forward free-boundary vacuum contract: it evaluates external coil flux on a named boundary/limiter contour using the same circular filament Green’s function as the grid vacuum solve. If target_flux is provided, diagnostics report pointwise residuals without fitting a scale or replaying Dirichlet boundary values. Optional limiter, magnetic-axis, and X-point metadata are sampled with the same response matrix convention so diagnostics cover topology surfaces instead of only the LCFS contour.

Return type:

dict[str, Any]

Parameters:

• kernel (Any)

• coils (CoilSet)

• boundary_points (ndarray[Any, dtype[float64]])

• limiter_points (ndarray[Any, dtype[float64]] | None)

• axis_point (ndarray[Any, dtype[float64]] | None)

• x_points (ndarray[Any, dtype[float64]] | None)

• target_flux (ndarray[Any, dtype[float64]] | None)

scpn_fusion.core.fusion_kernel_free_boundary.resolve_shape_target_flux(kernel, coils)

Resolve shape-control target flux values at configured control points.

Return type:

ndarray[Any, dtype[float64]]

Parameters:

• kernel (Any)

• coils (CoilSet)

scpn_fusion.core.fusion_kernel_free_boundary.solve_free_boundary(kernel, coils, max_outer_iter=20, tol=0.0001, optimize_shape=False, tikhonov_alpha=0.0001, limiter_points=None, axis_point=None, x_points=None)

Solve the Grad-Shafranov system with externally driven coil boundary flux.

Each outer iteration applies coil-generated boundary flux, runs the kernel equilibrium solve, and optionally re-optimises coil currents against the configured shape-control points. The result reports the outer iteration count, the final maximum flux-grid change, and the final coil currents.

Return type:

dict[str, Any]

Parameters:

• kernel (Any)

• coils (CoilSet)

• max_outer_iter (int)

• tol (float)

• optimize_shape (bool)

• tikhonov_alpha (float)

• limiter_points (ndarray[Any, dtype[float64]] | None)

• axis_point (ndarray[Any, dtype[float64]] | None)

• x_points (ndarray[Any, dtype[float64]] | None)

Free-boundary adapter methods mixed into FusionKernel.

class scpn_fusion.core.fusion_kernel_free_boundary_mixin.FusionKernelFreeBoundaryMixin

Bases: object

Attach free-boundary operations to FusionKernel.

The implementation lives in scpn_fusion.core.fusion_kernel_free_boundary so the core kernel remains focused on grid state and equilibrium solving. This mixin preserves the existing FusionKernel method API while routing each operation to the validated adapter function used by tests and tooling.

build_coilset_from_config()

Build a validated free-boundary coil set from self.cfg.

Return type:

CoilSet

reconstruct_coil_currents_from_magnetic_probes(coils, *, flux_points=None, flux_measurements=None, b_probe_points=None, b_probe_directions=None, b_probe_measurements=None, measurement_sigma=None, tikhonov_alpha=1e-06)

Fit free-boundary coil currents from magnetic diagnostic measurements.

Return type:

dict[str, Any]

Parameters:

• coils (CoilSet)

• flux_points (ndarray[Any, dtype[float64]] | None)

• flux_measurements (ndarray[Any, dtype[float64]] | None)

• b_probe_points (ndarray[Any, dtype[float64]] | None)

• b_probe_directions (list[str] | tuple[str, ...] | None)

• b_probe_measurements (ndarray[Any, dtype[float64]] | None)

• measurement_sigma (ndarray[Any, dtype[float64]] | None)

• tikhonov_alpha (float)

reconstruct_boundary_flux_from_coils(coils, *, boundary_points, limiter_points=None, axis_point=None, x_points=None, target_flux=None)

Reconstruct boundary-contour vacuum flux from coil Green functions.

Return type:

dict[str, Any]

Parameters:

• coils (CoilSet)

• boundary_points (ndarray[Any, dtype[float64]])

• limiter_points (ndarray[Any, dtype[float64]] | None)

• axis_point (ndarray[Any, dtype[float64]] | None)

• x_points (ndarray[Any, dtype[float64]] | None)

• target_flux (ndarray[Any, dtype[float64]] | None)

optimize_coil_currents(coils, target_flux, tikhonov_alpha=0.0001)

Fit bounded coil currents for the requested target flux vector.

Return type:

ndarray[Any, dtype[float64]]

Parameters:

• coils (CoilSet)

• target_flux (ndarray[Any, dtype[float64]])

• tikhonov_alpha (float)

solve_free_boundary(coils, max_outer_iter=20, tol=0.0001, optimize_shape=False, tikhonov_alpha=0.0001, limiter_points=None, axis_point=None, x_points=None)

Run the free-boundary outer loop with optional coil optimisation.

Return type:

dict[str, Any]

Parameters:

• coils (CoilSet)

• max_outer_iter (int)

• tol (float)

• optimize_shape (bool)

• tikhonov_alpha (float)

• limiter_points (ndarray[Any, dtype[float64]] | None)

• axis_point (ndarray[Any, dtype[float64]] | None)

• x_points (ndarray[Any, dtype[float64]] | None)

Configuration-to-CoilSet construction for free-boundary kernels.

scpn_fusion.core.fusion_kernel_coilset_config.build_coilset_from_config(kernel)

Build a validated CoilSet.

The configuration contract is intentionally strict because these values feed the free-boundary least-squares and Green’s-function paths directly. Coil positions must be finite cylindrical (R, Z) coordinates with R > 0; turns must be positive integers; optional current limits must be finite and positive; and every optional target or diagnostic vector must have a length consistent with its associated point array.

Return type:

CoilSet

Parameters:

kernel (Any)

Fusion Kernel Solver Runtime

Shared numeric helpers for FusionKernel solver paths.

scpn_fusion.core.fusion_kernel_numerics.sanitize_numeric_array(arr, *, cap=1e+250)

Return finite array with values clipped to [-cap, cap].

Return type:

ndarray[Any, dtype[float64]]

Parameters:

• arr (ndarray[Any, dtype[float64]])

• cap (float)

scpn_fusion.core.fusion_kernel_numerics.stable_rms(arr)

Compute RMS without overflow from squaring very large magnitudes.

Return type:

float

Parameters:

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

Iterative linear/nonlinear GS sub-solvers for FusionKernel.

class scpn_fusion.core.fusion_kernel_iterative_solver.FusionKernelIterativeSolverMixin

Bases: object

Mixin providing Jacobi, SOR, and multigrid sub-solvers for GS solves.

Newton and top-level equilibrium dispatch routines for FusionKernel.

class scpn_fusion.core.fusion_kernel_newton_solver.FusionKernelNewtonSolverMixin

Bases: object

Mixin providing Newton-Kantorovich equilibrium dispatch for FusionKernel.

solve_equilibrium(preserve_initial_state=False, boundary_flux=None)

Run the full Picard-iteration equilibrium solver.

Iterates: topology analysis -> source update -> elliptic solve -> under-relaxation until the residual drops below the configured convergence threshold or the maximum iteration count is reached.

When solver.solver_method is "anderson", Anderson acceleration is applied every few Picard steps to speed up convergence.

Returns:

{"psi": FloatArray, "converged": bool, "iterations": int, "residual": float, "residual_history": list[float], "gs_residual": float, "gs_residual_best": float, "gs_residual_history": list[float], "wall_time_s": float, "solver_method": str}

Return type:

dict[str, Any]

Parameters:

• preserve_initial_state (bool, optional) – Keep current interior self.Psi values and only enforce boundary conditions before the iterative solve. Default is False.

• boundary_flux (FloatArray | None, optional) – Explicit boundary map to enforce during solve. Must have shape (NZ, NR) when provided.

Raises:

RuntimeError – If the solver produces NaN or Inf and solver.fail_on_diverge is enabled in the active configuration.

Shared Newton-solver option parsing and telemetry helpers.

class scpn_fusion.core.fusion_kernel_newton_runtime.GmresTelemetry(nonconverged_count=0, breakdown_count=0, last_info=0)

Bases: object

Track GMRES non-convergence and breakdown events during Newton solves.

Parameters:

• nonconverged_count (int)

• breakdown_count (int)

• last_info (int)

nonconverged_count: int = 0

breakdown_count: int = 0

last_info: int = 0

record(*, info, iter_idx, fail_on_breakdown, nonconverged_budget)

Record one GMRES status code and enforce configured failure budgets.

Return type:

None

Parameters:

• info (int)

• iter_idx (int)

• fail_on_breakdown (bool)

• nonconverged_budget (int | None)

class scpn_fusion.core.fusion_kernel_newton_runtime.NewtonDispatchConfig(max_iter, tol, picard_alpha, fail_on_diverge, require_gs_residual, gs_tol, mu0, warmup_steps, newton_alpha, use_newton_line_search, line_search_c, max_backtracks, gmres_preconditioner_mode, ilu_drop_tol, ilu_fill_factor, gmres_fail_on_breakdown, gmres_nonconverged_budget)

Bases: object

Validated Newton-dispatch options extracted from kernel configuration.

Parameters:

• max_iter (int)

• tol (float)

• picard_alpha (float)

• fail_on_diverge (bool)

• require_gs_residual (bool)

• gs_tol (float)

• mu0 (float)

• warmup_steps (int)

• newton_alpha (float)

• use_newton_line_search (bool)

• line_search_c (float)

• max_backtracks (int)

• gmres_preconditioner_mode (str)

• ilu_drop_tol (float)

• ilu_fill_factor (float)

• gmres_fail_on_breakdown (bool)

• gmres_nonconverged_budget (int | None)

max_iter: int

tol: float

picard_alpha: float

fail_on_diverge: bool

require_gs_residual: bool

gs_tol: float

mu0: float

warmup_steps: int

newton_alpha: float

use_newton_line_search: bool

line_search_c: float

max_backtracks: int

gmres_preconditioner_mode: str

ilu_drop_tol: float

ilu_fill_factor: float

gmres_fail_on_breakdown: bool

gmres_nonconverged_budget: int | None

scpn_fusion.core.fusion_kernel_newton_runtime.parse_newton_dispatch_config(cfg)

Parse and validate Newton-dispatch options from kernel config.

Return type:

NewtonDispatchConfig

Parameters:

cfg (dict[str, Any])

Shared runtime helpers for FusionKernel solver mixins.

scpn_fusion.core.fusion_kernel_solver_runtime.apply_gs_operator(kernel, v)

Apply the discrete GS* operator to array v.

Return type:

ndarray[Any, dtype[float64]]

Parameters:

• kernel (Any)

• v (ndarray[Any, dtype[float64]])

scpn_fusion.core.fusion_kernel_solver_runtime.compute_gs_residual(kernel, Source)

Compute the GS residual r = L*[psi] - Source on interior points.

Return type:

ndarray[Any, dtype[float64]]

Parameters:

• kernel (Any)

• Source (ndarray[Any, dtype[float64]])

scpn_fusion.core.fusion_kernel_solver_runtime.compute_gs_residual_rms(kernel, Source)

Return RMS GS residual over interior points.

Return type:

float

Parameters:

• kernel (Any)

• Source (ndarray[Any, dtype[float64]])

scpn_fusion.core.fusion_kernel_solver_runtime.compute_profile_jacobian(kernel, Psi_axis, Psi_boundary, mu0)

Compute dJ_phi/dpsi as a 2D diagonal scaling field.

Return type:

ndarray[Any, dtype[float64]]

Parameters:

• kernel (Any)

• Psi_axis (float)

• Psi_boundary (float)

• mu0 (float)

scpn_fusion.core.fusion_kernel_solver_runtime.solve_newton_linear_system(*, kernel, J_op, rhs, diag_term, n_interior, nz, nr, gmres_preconditioner_mode, ilu_drop_tol, ilu_fill_factor, iter_idx)

Solve J * delta = rhs with optional GMRES preconditioning.

Return type:

tuple[ndarray[Any, dtype[float64]], int]

Parameters:

• kernel (Any)

• J_op (Any)

• rhs (ndarray[Any, dtype[float64]])

• diag_term (ndarray[Any, dtype[float64]])

• n_interior (int)

• nz (int)

• nr (int)

• gmres_preconditioner_mode (str)

• ilu_drop_tol (float)

• ilu_fill_factor (float)

• iter_idx (int)

scpn_fusion.core.fusion_kernel_solver_runtime.solve_via_rust_multigrid(kernel, preserve_initial_state=False, boundary_flux=None)

Delegate the full equilibrium solve to the Rust multigrid backend.

Return type:

dict[str, Any]

Parameters:

• kernel (Any)

• preserve_initial_state (bool)

• boundary_flux (ndarray[Any, dtype[float64]] | 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 include Lao et al., Nucl. Fusion 25 (1985) 1611, FreeQDSK documentation, and the ITER IMAS equilibrium IDS.

The parser expects the standard Fortran fixed-width layout: five values per line using (5e16.9), a 48-character header description followed by idum, nw and nh, the 20-scalar equilibrium block, five 1-D profile arrays, the row-major psirz flux map, and boundary/limiter (R, Z) point lists.

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_convention='raw_canonical', source_convention_adapter='not_applied', source_convention_adapter_pass=False, source_convention_metadata=<factory>)

Bases: object

Container for all data in a G-EQDSK file.

Parameters:

• description (str)

• nw (int)

• nh (int)

• rdim (float)

• zdim (float)

• rcentr (float)

• rleft (float)

• zmid (float)

• rmaxis (float)

• zmaxis (float)

• simag (float)

• sibry (float)

• bcentr (float)

• current (float)

• 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]])

• source_convention (str)

• source_convention_adapter (str)

• source_convention_adapter_pass (bool)

• source_convention_metadata (dict[str, object])

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]]

source_convention: str = 'raw_canonical'

source_convention_adapter: str = 'not_applied'

source_convention_adapter_pass: bool = False

source_convention_metadata: dict[str, object]

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)

Map raw ψ to normalised ψ_N.

Return type:

ndarray[Any, dtype[float64]]

Parameters:

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

to_config(name='eqdsk')

Convert to a FusionKernel JSON config with GEQDSK shape metadata.

GEQDSK files do not contain an external-coil set, so the returned coils list remains empty. They do contain plasma-boundary and limiter contours; those are exported into free_boundary so native free-boundary workflows can use the EFIT boundary as an isoflux shape contract instead of silently discarding it.

Return type:

dict[str, object]

Parameters:

name (str)

scpn_fusion.core.eqdsk.read_geqdsk(path, *, source_convention_mode='raw_canonical')

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.

• source_convention_mode ({"raw_canonical", "public_sparc_named_adapter"}) –

  raw_canonical keeps parsed ffprime and pprime in raw file units (default, strict mode).

  public_sparc_named_adapter applies a documented, source-aware adapter only for explicitly recognized public SPARC files.

Returns:

Parsed equilibrium data.

Return type:

GEqdsk

scpn_fusion.core.eqdsk.write_geqdsk(eq, path)

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

Newton-Raphson PF-coil force-balance adjustment for equilibrium configs.

class scpn_fusion.core.force_balance.ForceBalanceSolver(config_path)

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)

Iteratively adjust paired PF currents until radial force is balanced.

Return type:

dict[str, float | int | bool]

Parameters:

• max_iterations (int)

• jacobian_floor (float)

save_config(output_path=None)

Write the current balanced kernel configuration to disk.

Parameters:

output_path (str | Path | None)

Compact Reactor Optimiser

Reduced compact-reactor design-space search and reporting utilities.

class scpn_fusion.core.compact_reactor_optimizer.CompactReactorArchitect

Bases: object

Optimise compact fusion reactor geometry for minimum size.

plasma_physics_model(R, a, B0)

Estimate fusion power, plasma current, and volume for a compact tokamak design.

Return type:

tuple[float, float, float]

Parameters:

• R (float)

• a (float)

• B0 (float)

radial_build_constraints(R, a, B0)

Check simple inboard radial-build and high-field coil feasibility constraints.

Return type:

tuple[bool, float]

Parameters:

• R (float)

• a (float)

• B0 (float)

visualize_space(designs, label='')

Write a scatter plot of feasible compact-reactor designs for review.

Return type:

None

Parameters:

• designs (list[dict[str, float]])

• label (str)

calculate_economics(d)

Calculate the cost of electricity (CoE) in $/MWh.

Uses the Sheffield model CoE ~ (C_cap * F_cap + C_om) / (8760 * P_net * f_avail).

Return type:

tuple[float, float]

Parameters:

d (dict[str, float])

report_design(d)

Print the engineering and economic summary for a selected reactor design.

Return type:

None

Parameters:

d (dict[str, float])

find_minimum_reactor(target_power_MW=1.0, use_temhd=True)

Search grid candidates and report the smallest design meeting the power target.

Return type:

None

Parameters:

• target_power_MW (float)

• use_temhd (bool)

Global Design Scanner

Monte Carlo reactor-design scans with heat-exhaust and magnet constraints.

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)

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)

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)

Sample the global reactor design space and return valid design metrics.

run_compact_scan(n_samples=2000, seed=42)

Run the compact-reactor search envelope with deterministic sampling.

analyze_pareto(df)

Print and plot the viable Pareto region for a design-scan dataframe.

Integrated Transport Solver

Top-level integrated transport solver orchestration and factory helpers.

This module wires together transport-model and runtime mixins into the public IntegratedTransportSolver interface. It also exposes conservative fallback paths (Python implementations when Rust bindings are unavailable) and preserves legacy compatibility symbols required by downstream callers.

Docstring scope in this module is intentionally high-level because most detailed numerical contracts are documented in the extracted model/runtime mixins.

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)

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)

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). solver must provide the TransportSolver-compatible Ti, Te and evolve_profiles API.

Returns the estimated error norm.

Return type:

float

Parameters:

• solver (Any)

• P_aux (float)

• enforce_numerical_recovery (bool)

• max_numerical_recoveries (int | None)

adapt_dt(error)

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

Bases: RuntimeError, FusionCoreError

Raised when a physics constraint is violated.

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

Bases: TransportSolverInitializationMixin, 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)

__init__(config_path, *, multi_ion=False)

Build a transport solver with deterministic defaults and requested mode.

Parameters:

• config_path (str | Path)

• multi_ion (bool)

Return type:

None

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)

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)

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

RF heating workflows with simplified ICRH/ECRH ray-tracing examples.

class scpn_fusion.core.rf_heating.RFHeatingSystem(config_path)

Bases: object

Simulate Ion Cyclotron Resonance Heating (ICRH).

Uses ray-tracing to track EM waves launched from the antenna and absorbed at the resonance layer.

Parameters:

config_path (str)

get_plasma_params(R, Z)

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

Return type:

tuple[float, float, float, float]

Parameters:

• R (float)

• Z (float)

dispersion_relation(R, Z, k_R, k_Z)

Calculate the local dispersion relation D(omega, k) = 0.

Includes warm-plasma thermal (finite-Larmor-radius) corrections for ICRH.

Return type:

float

Parameters:

• R (float)

• Z (float)

• k_R (float)

• k_Z (float)

ray_equations(state, t)

Evaluate the Hamiltonian ray-tracing equations.

dr/dt = dD/dk and dk/dt = -dD/dr for the local dispersion relation D.

Return type:

list[float]

Parameters:

• state (ndarray[Any, dtype[float64]])

• t (float)

trace_rays(n_rays=10)

Trace ICRH rays from an outboard antenna and return trajectories.

Parameters:

n_rays (int) – Number of launch rays.

Return type:

tuple[list[ndarray[Any, dtype[float64]]], float]

Returns:

• trajectories (list of ndarray) – Ray coordinates and wave-vectors as ODE solution arrays.

• B_res – Resonant magnetic field value in Tesla for the configured frequency.

plot_heating(trajectories, B_res)

Render and persist a diagnostic plot of ray paths and resonance layer.

Return type:

None

Parameters:

• trajectories (list[ndarray[Any, dtype[float64]]])

• B_res (float)

compute_power_deposition(trajectories, P_rf_mw=20.0, n_radial_bins=50)

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.

Return type:

tuple[ndarray[Any, dtype[float64]], ndarray[Any, dtype[float64]], float]

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)

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()

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)

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[Any, dtype[float64]], ndarray[Any, dtype[float64]], 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:

• P_ecrh_mw (float)

• n_radial_bins (int)

• T_e_keV (float)

• n_e (float)

• launch_angle_deg (float)

MHD Sawtooth

Sawtooth instability toy solver and diagnostic visualizations.

class scpn_fusion.core.mhd_sawtooth.ReducedMHD(nr=100)

Bases: object

Simulate 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)

__init__(nr=100)

Create a reduced-MHD 1D m=1, n=1 model on a radial grid.

Parameters:

nr (int)

Return type:

None

laplacian(f, m=1)

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

Return type:

ndarray[Any, dtype[complex128]]

Parameters:

• f (ndarray[Any, dtype[complex128]])

• m (int)

step(dt=0.01)

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)

Solve Del^2 phi = U for phi using the Thomas Algorithm (O(N)).

Optimized for tridiagonal Laplacian in cylindrical coordinates.

Return type:

ndarray[Any, dtype[complex128]]

Parameters:

U (ndarray[Any, dtype[complex128]])

scpn_fusion.core.mhd_sawtooth.run_sawtooth_sim()

Run a simple sawtooth simulation and emit diagnostic figures as PNG files.

Return type:

None

Hall-MHD Discovery

Reduced Hall-MHD discovery workflow for zonal-flow and tearing diagnostics.

The implementation is intentionally lightweight and deterministic so that CI and local benchmark tasks can exercise a full discovery-style simulation path:

• semi-spectral Hall-MHD dynamics in periodic geometry,

• short parameter sweeps for eta/nu response grids,

• bisection-style tearing threshold search, and

• automated visual output for inspection.

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

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)

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)

Compute the 2D Poisson bracket [A, B] in spectral space.

The bracket uses FFT-domain derivatives and returns the transformed commutator dxA * dyB - dyA * dxB.

Parameters:

• A_k – Real- or complex-valued spectral field A(kx, ky).

• B_k – Real- or complex-valued spectral field B(kx, ky).

Returns:

Fourier representation of the Poisson bracket [A, B].

dynamics(phi, psi)

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()

Advance one reduced Hall-MHD pseudo-time step (RK2).

Returns:

(total_energy, zonal_energy) for the post-step state in spectral coordinates. total_energy is the total potential energy proxy from spectral coefficients; zonal_energy accumulates non-zero ky=0 modes as a zonal-flow metric.

parameter_sweep(eta_range, nu_range, n_steps=5, sim_steps=200)

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)

Bisection search for marginal tearing stability threshold.

scpn_fusion.core.hall_mhd_discovery.run_discovery_sim()

Run the standalone Hall-MHD discovery demo and emit figure artifacts.

Side effects:

• prints periodic progress snapshots,

• writes Hall_MHD_Discovery.png,

• writes Hall_MHD_Structure.png.

Stability Analyser

Reduced-order equilibrium force and MHD stability diagnostics.

class scpn_fusion.core.stability_analyzer.StabilityAnalyzer(config_path)

Bases: object

Linear eigenvalue stability analysis of plasma rigid-body motion.

Analyses the n=0 mode (radial and vertical displacement) via a force-balance stiffness matrix.

Parameters:

config_path (str)

get_vacuum_field_at(R, Z)

Interpolate vacuum-field components and decay index at position (R, Z).

Returns (Bz, Br, n_index) where Bz = 1/R dPsi/dR, Br = -1/R dPsi/dZ and n_index is the field decay index.

Return type:

tuple[float, float, float]

Parameters:

• R (float)

• Z (float)

calculate_forces(R, Z, Ip)

Compute radial and vertical forces acting on the plasma ring.

Returns (F_R, F_Z, n_index) where F_R = F_hoop + F_Lorentz_R and F_Z = F_Lorentz_Z.

Return type:

tuple[float, float, float]

Parameters:

• R (float)

• Z (float)

• Ip (float)

analyze_stability(R_target=6.2, Z_target=0.0)

Run force-balance linearization and print eigenvalue stability summary.

Return type:

None

Parameters:

• R_target (float)

• Z_target (float)

analyze_mhd_stability(R0=6.2, a=2.0, B0=5.3, Ip_MA=15.0, transport_solver=None)

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)

Generate and persist a radial-force contour landscape around (R0, Z0).

Return type:

None

Parameters:

• R0 (float)

• Z0 (float)

MHD Stability Criteria

MHD stability analysis suite — seven criteria.

1. Mercier — interchange stability (D_M >= 0)

2. Ballooning — first-stability boundary (Connor-Hastie-Taylor)

3. Kruskal-Shafranov — external kink (q_edge > 1)

4. Troyon — normalised beta limit (beta_N < g)

5. NTM — neoclassical tearing mode seeding threshold

6. RWM — resistive wall mode (beta_N between no-wall and ideal-wall limits)

7. Peeling-ballooning — coupled pedestal stability (ELM boundary)

Criteria 4–7 live in stability_mhd_extended.py and are re-exported here for backward compatibility.

References

• Freidberg, Ideal MHD, Cambridge (2014), Ch. 12

• Connor, Hastie & Taylor, Phys. Rev. Lett. 40:396 (1978)

• Kruskal & Schwarzschild, Proc. R. Soc. Lond. A 223:348 (1954)

• Troyon et al., Plasma Phys. Control. Fusion 26:209 (1984)

• La Haye, Phys. Plasmas 13:055501 (2006)

• Snyder et al., Phys. Plasmas 9:2037 (2002)

• Snyder et al., Nucl. Fusion 51:103016 (2011)

class scpn_fusion.core.stability_mhd.QProfile(rho, q, shear, alpha_mhd, q_min, q_min_rho, q_edge)

Bases: object

Safety-factor profile and derived quantities.

Parameters:

• rho (ndarray[Any, dtype[float64]])

• q (ndarray[Any, dtype[float64]])

• shear (ndarray[Any, dtype[float64]])

• alpha_mhd (ndarray[Any, dtype[float64]])

• q_min (float)

• q_min_rho (float)

• q_edge (float)

rho: ndarray[Any, dtype[float64]]

q: ndarray[Any, dtype[float64]]

shear: ndarray[Any, dtype[float64]]

alpha_mhd: ndarray[Any, dtype[float64]]

q_min: float

q_min_rho: float

q_edge: float

class scpn_fusion.core.stability_mhd.MercierResult(rho, D_M, stable, first_unstable_rho)

Bases: object

Mercier interchange stability result.

Parameters:

• rho (ndarray[Any, dtype[float64]])

• D_M (ndarray[Any, dtype[float64]])

• stable (ndarray[Any, dtype[bool_]])

• first_unstable_rho (float | None)

rho: ndarray[Any, dtype[float64]]

D_M: ndarray[Any, dtype[float64]]

stable: ndarray[Any, dtype[bool_]]

first_unstable_rho: float | None

class scpn_fusion.core.stability_mhd.BallooningResult(rho, s, alpha, alpha_crit, stable, margin)

Bases: object

First-stability ballooning boundary result.

Parameters:

• rho (ndarray[Any, dtype[float64]])

• s (ndarray[Any, dtype[float64]])

• alpha (ndarray[Any, dtype[float64]])

• alpha_crit (ndarray[Any, dtype[float64]])

• stable (ndarray[Any, dtype[bool_]])

• margin (ndarray[Any, dtype[float64]])

rho: ndarray[Any, dtype[float64]]

s: ndarray[Any, dtype[float64]]

alpha: ndarray[Any, dtype[float64]]

alpha_crit: ndarray[Any, dtype[float64]]

stable: ndarray[Any, dtype[bool_]]

margin: ndarray[Any, dtype[float64]]

class scpn_fusion.core.stability_mhd.KruskalShafranovResult(q_edge, stable, margin)

Bases: object

External kink stability result (Kruskal-Shafranov criterion).

Parameters:

• q_edge (float)

• stable (bool)

• margin (float)

q_edge: float

stable: bool

margin: float

class scpn_fusion.core.stability_mhd.TroyonResult(beta_N, beta_N_crit_nowall, beta_N_crit_wall, stable_nowall, stable_wall, margin_nowall)

Bases: object

Troyon normalised-beta-limit result.

Parameters:

• beta_N (float)

• beta_N_crit_nowall (float)

• beta_N_crit_wall (float)

• stable_nowall (bool)

• stable_wall (bool)

• margin_nowall (float)

beta_N: float

beta_N_crit_nowall: float

beta_N_crit_wall: float

stable_nowall: bool

stable_wall: bool

margin_nowall: float

class scpn_fusion.core.stability_mhd.NTMResult(rho, delta_prime, j_bs_drive, w_marginal, ntm_unstable, most_unstable_rho)

Bases: object

Neoclassical tearing mode seeding analysis result.

Parameters:

• rho (ndarray[Any, dtype[float64]])

• delta_prime (ndarray[Any, dtype[float64]])

• j_bs_drive (ndarray[Any, dtype[float64]])

• w_marginal (ndarray[Any, dtype[float64]])

• ntm_unstable (ndarray[Any, dtype[bool_]])

• most_unstable_rho (float | None)

rho: ndarray[Any, dtype[float64]]

delta_prime: ndarray[Any, dtype[float64]]

j_bs_drive: ndarray[Any, dtype[float64]]

w_marginal: ndarray[Any, dtype[float64]]

ntm_unstable: ndarray[Any, dtype[bool_]]

most_unstable_rho: float | None

class scpn_fusion.core.stability_mhd.RWMResult(beta_N, beta_N_crit_nowall, beta_N_crit_wall, stable, mode_growth_rate)

Bases: object

Resistive Wall Mode stability result.

Parameters:

• beta_N (float)

• beta_N_crit_nowall (float)

• beta_N_crit_wall (float)

• stable (bool)

• mode_growth_rate (float)

beta_N: float

beta_N_crit_nowall: float

beta_N_crit_wall: float

stable: bool

mode_growth_rate: float

class scpn_fusion.core.stability_mhd.PeelingBallooningResult(j_edge_norm, alpha_edge_norm, stability_distance, stable, elm_type)

Bases: object

Coupled peeling-ballooning stability at the H-mode pedestal.

Snyder et al., Phys. Plasmas 9:2037 (2002); Nucl. Fusion 51:103016 (2011).

Parameters:

• j_edge_norm (float)

• alpha_edge_norm (float)

• stability_distance (float)

• stable (bool)

• elm_type (str | None)

j_edge_norm: float

alpha_edge_norm: float

stability_distance: float

stable: bool

elm_type: str | None

class scpn_fusion.core.stability_mhd.StabilitySummary(mercier, ballooning, kruskal_shafranov, troyon, ntm, rwm, peeling_ballooning, n_criteria_checked, n_criteria_stable, overall_stable)

Bases: object

Combined result from all MHD stability criteria.

Parameters:

• mercier (MercierResult)

• ballooning (BallooningResult)

• kruskal_shafranov (KruskalShafranovResult)

• troyon (TroyonResult | None)

• ntm (NTMResult | None)

• rwm (RWMResult | None)

• peeling_ballooning (PeelingBallooningResult | None)

• n_criteria_checked (int)

• n_criteria_stable (int)

• overall_stable (bool)

mercier: MercierResult

ballooning: BallooningResult

kruskal_shafranov: KruskalShafranovResult

troyon: TroyonResult | None

ntm: NTMResult | None

rwm: RWMResult | None

peeling_ballooning: PeelingBallooningResult | None

n_criteria_checked: int

n_criteria_stable: int

overall_stable: bool

scpn_fusion.core.stability_mhd.compute_q_profile(rho, ne, Ti, Te, R0, a, B0, Ip_MA, kappa=1.0, delta=0.0)

Compute the safety-factor profile from a shape-aware approximation.

Uses a parabolic current profile and a Uckan-style geometric correction for elongation (kappa) and triangularity (delta).

Parameters:

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

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

• Ti (array — ion/electron temperature [keV])

• Te (array — ion/electron temperature [keV])

• R0 (float — major radius [m])

• a (float — minor radius [m])

• B0 (float — toroidal field on axis [T])

• Ip_MA (float — total plasma current [MA])

• kappa (float — elongation)

• delta (float — triangularity)

Return type:

QProfile

scpn_fusion.core.stability_mhd.mercier_stability(qp)

Evaluate interchange stability using the Suydam approximation.

This is a low-order cylindrical proxy for the full Mercier criterion, not a complete toroidal Mercier evaluation.

D_M = s^2 / 4 - alpha_mhd Stable where D_M >= 0.

Return type:

MercierResult

Parameters:

qp (QProfile)

scpn_fusion.core.stability_mhd.ballooning_stability(qp)

Evaluate the first ballooning stability boundary.

The critical normalised pressure gradient (Connor-Hastie-Taylor 1978):

alpha_crit(s) = s*(1 - s/2) for s < 1

= 0.6*s for s >= 1

Stable where alpha <= alpha_crit.

Parameters:

qp (QProfile)

Return type:

BallooningResult

scpn_fusion.core.stability_mhd.kruskal_shafranov_stability(qp)

Evaluate the Kruskal-Shafranov external kink stability criterion.

The plasma is stable against the m=1/n=1 external kink mode when the edge safety factor satisfies q(edge) > 1.

Parameters:

qp (QProfile)

Return type:

KruskalShafranovResult

References

Kruskal & Schwarzschild, Proc. R. Soc. Lond. A 223:348 (1954) Shafranov, Sov. Phys. Tech. Phys. 15:175 (1970)

scpn_fusion.core.stability_mhd.run_full_stability_check(qp, beta_t=None, Ip_MA=None, a=None, B0=None, R0=None, j_bs=None, j_total=None, j_edge=None, p_ped_Pa=None, kappa=1.7, delta=0.3)

Run all available MHD stability criteria and return a summary.

Mercier, ballooning, and Kruskal-Shafranov are always evaluated. Troyon requires beta_t, Ip_MA, a, and B0. NTM requires j_bs, j_total, and a. Peeling-ballooning requires j_edge, p_ped_Pa, R0, a, B0.

Parameters:

• qp (QProfile — pre-computed safety-factor profile)

• beta_t (float, optional — total toroidal beta (dimensionless))

• Ip_MA (float, optional — plasma current [MA])

• a (float, optional — minor radius [m])

• B0 (float, optional — toroidal field on axis [T])

• R0 (float, optional — major radius [m])

• j_bs (array, optional — bootstrap current density [A/m^2])

• j_total (array, optional — total current density [A/m^2])

• j_edge (float, optional — edge parallel current density [A/m^2])

• p_ped_Pa (float, optional — pedestal-top pressure [Pa])

• kappa (float — elongation (default 1.7))

• delta (float — triangularity (default 0.3))

Return type:

StabilitySummary

Extended MHD stability criteria: Troyon, NTM, RWM, peeling-ballooning.

Separated from stability_mhd.py for module-size compliance. All result dataclasses remain in stability_mhd.py.

References

• Troyon et al., Plasma Phys. Control. Fusion 26:209 (1984)

• La Haye, Phys. Plasmas 13:055501 (2006)

• Snyder et al., Phys. Plasmas 9:2037 (2002)

• Snyder et al., Nucl. Fusion 51:103016 (2011)

scpn_fusion.core.stability_mhd_extended.troyon_beta_limit(beta_t, Ip_MA, a, B0, g_nowall=2.8, g_wall=3.5)

Evaluate the Troyon normalised-beta limit.

The normalised beta is defined as:

\[\beta_N = \frac{\beta_t}{I_N} \quad\text{where}\quad I_N = \frac{I_p\,[\text{MA}]}{a\,[\text{m}]\,B_0\,[\text{T}]}\]

The factor of 100 converts fractional beta_t into the conventional percent-based normalisation (units: % m T / MA).

Stability requires beta_N < g, where g ~ 2.8 (no-wall) or g ~ 3.5 (ideal-wall).

Parameters:

• beta_t (float) – Total toroidal beta (dimensionless, e.g. 0.025).

• Ip_MA (float) – Plasma current [MA].

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

• B0 (float) – Toroidal magnetic field on axis [T].

• g_nowall (float) – Troyon coefficient without wall (default 2.8).

• g_wall (float) – Troyon coefficient with ideal wall (default 3.5).

Return type:

TroyonResult

References

Troyon et al., Plasma Phys. Control. Fusion 26:209 (1984)

scpn_fusion.core.stability_mhd_extended.ntm_stability(qp, j_bs, j_total, a, r_s_delta_prime=-2.0)

Reduced-order neoclassical tearing mode (NTM) seeding analysis.

The modified Rutherford equation for island width w is (reduced-order):

\[\tau_R\,\frac{dw}{dt} = r_s\,\Delta'(w) + \frac{j_\text{bs}}{j_\phi}\,\frac{a}{w}\]

An NTM is potentially unstable at a given radius when the bootstrap drive exceeds the classical stabilisation. The marginal island width (below which the island shrinks) is:

\[w_\text{marg} = -\frac{j_\text{bs} / j_\phi}{r_s \Delta'}\,a\]

Parameters:

• qp (QProfile) – Safety-factor profile providing the rational-surface locations.

• j_bs (array) – Bootstrap current density [A/m^2].

• j_total (array) – Total current density [A/m^2].

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

• r_s_delta_prime (float) – Classical tearing stability index (negative = classically stable, default -2.0).

Return type:

NTMResult

References

La Haye, Phys. Plasmas 13:055501 (2006) Sauter et al., Phys. Plasmas 4:1654 (1997)

scpn_fusion.core.stability_mhd_extended.rwm_stability(beta_N, g_nowall=2.8, g_wall=3.5)

Evaluate Resistive Wall Mode (RWM) stability.

The RWM occurs when beta_N exceeds the no-wall limit but remains below the ideal-wall limit. The mode grows on the resistive time-scale of the wall (tau_wall).

Parameters:

• beta_N (float) – Normalised beta [% m T / MA].

• g_nowall (float) – Troyon no-wall limit (default 2.8).

• g_wall (float) – Troyon ideal-wall limit (default 3.5).

Return type:

RWMResult

scpn_fusion.core.stability_mhd_extended.peeling_ballooning_stability(qp, j_edge, p_ped_Pa, R0, a, B0, kappa=1.7, delta=0.3)

Coupled peeling-ballooning stability for H-mode pedestal.

The instability boundary in normalised (j, alpha) space is approximately elliptical with shaping corrections:

\[\left(\frac{j_\parallel}{j_{\rm crit}}\right)^2 + \left(\frac{\alpha}{\alpha_{\rm crit}}\right)^2 = 1\]

Parameters:

• qp (QProfile) – Pre-computed safety-factor profile.

• j_edge (float) – Edge parallel current density [A/m^2].

• p_ped_Pa (float) – Pedestal-top pressure [Pa].

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

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

• B0 (float) – Toroidal field on axis [T].

• kappa (float) – Elongation (default 1.7).

• delta (float) – Triangularity (default 0.3).

Return type:

PeelingBallooningResult

References

Snyder et al., Phys. Plasmas 9:2037 (2002) Snyder et al., Nucl. Fusion 51:103016 (2011)

Runaway Electron Dynamics

Runaway electron utility contracts for DREAM-style reduced-order workflows.

class scpn_fusion.core.runaway_electrons.RunawayParams(ne_20, Te_keV, E_par, Z_eff, B0, R0, a=2.0)

Bases: object

Scalar plasma parameters used by reduced-order runaway balances.

Parameters:

• ne_20 (float)

• Te_keV (float)

• E_par (float)

• Z_eff (float)

• B0 (float)

• R0 (float)

• a (float)

ne_20: float

Te_keV: float

E_par: float

Z_eff: float

B0: float

R0: float

a: float = 2.0

class scpn_fusion.core.runaway_electrons.DreamFluidBalance(dreicer_source, avalanche_source, loss_source, total_source, runaway_fraction, growth_time_s)

Bases: object

DREAM-style fluid runaway density balance for a single plasma state.

Parameters:

• dreicer_source (float)

• avalanche_source (float)

• loss_source (float)

• total_source (float)

• runaway_fraction (float)

• growth_time_s (float)

dreicer_source: float

avalanche_source: float

loss_source: float

total_source: float

runaway_fraction: float

growth_time_s: float

scpn_fusion.core.runaway_electrons.dreicer_field(ne_20, Te_keV, coulomb_log=15.0)

Dreicer field E_D = n_e e^3 lnΛ / (4π ε₀² T_e).

Connor & Hastie, Nucl. Fusion 15, 415 (1975).

Return type:

float

Parameters:

• ne_20 (float)

• Te_keV (float)

• coulomb_log (float)

scpn_fusion.core.runaway_electrons.critical_field(ne_20, coulomb_log=15.0)

Critical field E_c = n_e e^3 lnΛ / (4π ε₀² m_e c²).

Rosenbluth & Putvinski, Nucl. Fusion 37, 1355 (1997). E_D / E_c = m_e c² / T_e ≈ 51 at T_e = 10 keV.

Return type:

float

Parameters:

• ne_20 (float)

• coulomb_log (float)

scpn_fusion.core.runaway_electrons.dreicer_generation_rate(params, coulomb_log=15.0)

Primary (seed) runaway electrons generation rate [m^-3 s^-1].

Connor & Hastie (1975).

Return type:

float

Parameters:

• params (RunawayParams)

• coulomb_log (float)

scpn_fusion.core.runaway_electrons.avalanche_growth_rate(params, n_RE, coulomb_log=15.0)

Avalanche multiplication rate [m^-3 s^-1].

Rosenbluth & Putvinski, Nucl. Fusion 37, 1355 (1997), Eq. 66. dn_RE/dt = n_RE (E/E_c - 1) / τ_c where τ_c already includes lnΛ.

Return type:

float

Parameters:

• params (RunawayParams)

• n_RE (float)

• coulomb_log (float)

scpn_fusion.core.runaway_electrons.hot_tail_seed(Te_pre_keV, Te_post_keV, ne_20, quench_time_ms, *, vc_vte_ref=4.0, quench_exponent=0.2)

Seed RE density from thermal quench hot-tail mechanism [m^-3].

Smith, H.M. et al., Phys. Plasmas 15, 072502 (2008). Faster TQ → lower v_c/v_te → more seed electrons.

Return type:

float

Parameters:

• Te_pre_keV (float)

• Te_post_keV (float)

• ne_20 (float)

• quench_time_ms (float)

• vc_vte_ref (float)

• quench_exponent (float)

scpn_fusion.core.runaway_electrons.dream_fluid_density_balance(params, n_RE, *, loss_time_s=inf, max_runaway_fraction=1.0, coulomb_log=15.0)

Evaluate the scalar density balance used by DREAM-style fluid runs.

The contract is dn_RE/dt = S_Dreicer + gamma_avalanche n_RE - n_RE/tau_loss. It is a fluid benchmark contract, not a kinetic DREAM distribution solver.

Return type:

DreamFluidBalance

Parameters:

• params (RunawayParams)

• n_RE (float)

• loss_time_s (float)

• max_runaway_fraction (float)

• coulomb_log (float)

class scpn_fusion.core.runaway_electrons.RunawayEvolution(params)

Bases: object

Time-domain reduced-order solver for runaway electron density evolution.

Parameters:

params (RunawayParams)

balance(n_RE, E_par, *, loss_time_s=inf, max_runaway_fraction=1.0)

Return the instantaneous DREAM-style fluid density balance.

Return type:

DreamFluidBalance

Parameters:

• n_RE (float)

• E_par (float)

• loss_time_s (float)

• max_runaway_fraction (float)

step(dt, n_RE, E_par, *, loss_time_s=inf, max_runaway_fraction=1.0)

Advance one explicit time step of the reduced-order runaway balance.

Return type:

float

Parameters:

• dt (float)

• n_RE (float)

• E_par (float)

• loss_time_s (float)

• max_runaway_fraction (float)

evolve(n_RE_0, E_par_profile, t_span, dt)

Integrate runaway density over a profile-driven E∥(t) history.

Return type:

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

Parameters:

• n_RE_0 (float)

• E_par_profile (Callable[[float], float])

• t_span (tuple[float, float])

• dt (float)

current_fraction(n_RE, I_p_MA)

Fraction of total plasma current carried by REs (v ≈ c).

Return type:

float

Parameters:

• n_RE (float)

• I_p_MA (float)

class scpn_fusion.core.runaway_electrons.RunawayMitigationAssessment

Bases: object

Collection of reduced-order mitigation formulas and guard checks.

static required_density_for_suppression(E_par, Z_eff, coulomb_log=15.0)

Density [10^20 m^-3] needed to make E_c > E_par.

Return type:

float

Parameters:

• E_par (float)

• Z_eff (float)

• coulomb_log (float)

static maximum_re_energy(B0, R0, *, ne_20=1.0, Z_eff=1.5, pitch_angle_rad=0.35)

Maximum RE energy [MeV] from orbit and synchrotron limits.

Return type:

float

Parameters:

• B0 (float)

• R0 (float)

• ne_20 (float)

• Z_eff (float)

• pitch_angle_rad (float)

static wall_heat_load(n_RE, E_max_MeV, A_wet, volume=800.0, *, mean_energy_MeV=None)

Deposited runaway-electron beam energy density [MJ/m^2].

Return type:

float

Parameters:

• n_RE (float)

• E_max_MeV (float)

• A_wet (float)

• volume (float)

• mean_energy_MeV (float | None)

Neural Equilibrium

PCA + MLP surrogate for Grad-Shafranov equilibrium reconstruction.

Maps the canonical 12-feature equilibrium descriptor to PCA coefficients and then to ψ(R,Z), providing ~1000× speedup over the full Picard iteration.

Training modes:

1. From FusionKernel — Generate training data by perturbing reactor parameters 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, shape, and scalar equilibrium parameters 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)

Bases: object

Configuration for the neural equilibrium model.

Parameters:

• n_components (int)

• hidden_sizes (tuple[int, ...])

• n_input_features (int)

• grid_shape (tuple[int, int])

• lambda_gs (float)

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)

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)

• test_mse (float)

• test_max_error (float)

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

scpn_fusion.core.neural_equilibrium.iter_surrogate_artifact_status(*, report_path=PosixPath('/home/runner/work/scpn-fusion-core/scpn-fusion-core/validation/reports/iter_surrogate_weight_validation.json'), weights_path=PosixPath('/home/runner/work/scpn-fusion-core/scpn-fusion-core/weights/neural_equilibrium_iter_v1.npz'), high_fidelity_report_path=PosixPath('/home/runner/work/scpn-fusion-core/scpn-fusion-core/validation/reports/iter_surrogate_training_report.json'))

Return the bounded acceptance status for the tracked ITER surrogate artifact.

Return type:

dict[str, object]

Parameters:

• report_path (Path)

• weights_path (Path)

• high_fidelity_report_path (Path)

class scpn_fusion.core.neural_equilibrium.SimpleMLP(layer_sizes, seed=42)

Bases: object

Feedforward MLP with ReLU hidden layers and linear output.

Parameters:

• layer_sizes (list[int])

• seed (int)

forward(x)

Compute a forward pass through all MLP layers.

Parameters:

x (NDArray) – Input feature batch with shape (batch, n_features).

Returns:

Predicted latent coefficients.

Return type:

NDArray

predict(x)

Alias for forward().

Parameters:

x (NDArray) – Input feature batch.

Returns:

Network output with shape (batch, n_outputs).

Return type:

NDArray

class scpn_fusion.core.neural_equilibrium.MinimalPCA(n_components=20)

Bases: object

Minimal PCA via SVD, no sklearn required.

Parameters:

n_components (int)

fit(X)

Fit principal-component basis on centred data.

Parameters:

X (NDArray) – Data matrix where rows are samples and columns are flattened coordinates.

Return type:

MinimalPCA

transform(X)

Project samples into PCA coefficient space.

Parameters:

X (NDArray) – Input matrix with the same feature dimension as fitted data.

Return type:

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

inverse_transform(Z)

Reconstruct flattened fields from PCA coefficients.

Parameters:

Z (NDArray) – Coefficient matrix.

Return type:

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

fit_transform(X)

Fit PCA model and immediately transform the training matrix.

Parameters:

X (NDArray) – Data matrix to fit and project.

Return type:

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

class scpn_fusion.core.neural_equilibrium.NeuralEquilibriumAccelerator(config=None)

Bases: object

PCA + MLP surrogate for Grad-Shafranov equilibrium.

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

Parameters:

config (NeuralEqConfig | None)

evaluate_surrogate(X_test, Y_test_raw)

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

Return type:

dict[str, float]

Parameters:

• X_test (ndarray[Any, dtype[_ScalarType_co]])

• Y_test_raw (ndarray[Any, dtype[_ScalarType_co]])

train_from_geqdsk(geqdsk_paths, n_perturbations=25, seed=42)

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:

• geqdsk_paths (list[Path])

• n_perturbations (int)

• seed (int)

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)

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'))

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'))

Load model from .npz.

Return type:

None

Parameters:

path (str | Path)

benchmark(features, n_runs=100)

Time inference over n_runs and return stats.

Return type:

dict[str, float]

Parameters:

• features (ndarray[Any, dtype[_ScalarType_co]])

• n_runs (int)

Neural Transport

Compatibility façade for the neural transport surrogate.

The implementation is split across focused helper modules while this module retains the historical import surface used by callers and tests.

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)

Bases: object

Stored weights for a variable-depth feedforward MLP.

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)

• gb_scale (bool)

• gated (bool)

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]]

First hidden-layer weight matrix, or an empty matrix when absent.

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

First hidden-layer bias vector, or an empty vector when absent.

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

Second hidden-layer weight matrix, or an empty matrix when absent.

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

Second hidden-layer bias vector, or an empty vector when absent.

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

Third hidden-layer weight matrix, or an empty matrix when absent.

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

Third hidden-layer bias vector, or an empty vector when absent.

property n_layers: int

Number of stored weight layers.

class scpn_fusion.core.neural_transport.NeuralTransportModel(weights_path=None)

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 (Optional[str | Path])

predict(inp)

Predict turbulent transport fluxes for given local parameters.

Return type:

TransportFluxes

Parameters:

inp (TransportInputs)

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

Predict transport coefficients on the full radial profile.

Return type:

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

Parameters:

• rho (ndarray[Any, dtype[float64]])

• te (ndarray[Any, dtype[float64]])

• ti (ndarray[Any, dtype[float64]])

• ne (ndarray[Any, dtype[float64]])

• q_profile (ndarray[Any, dtype[float64]])

• s_hat_profile (ndarray[Any, dtype[float64]])

• r_major (float)

• a_minor (float)

• b_toroidal (float)

scpn_fusion.core.neural_transport.NeuralTransportSurrogate

alias of NeuralTransportModel

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)

Bases: object

Predicted turbulent transport fluxes.

Parameters:

• chi_e (float)

• chi_i (float)

• d_e (float)

• channel (str)

• chi_e_itg (float)

• chi_e_tem (float)

• chi_e_etg (float)

• chi_i_itg (float)

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

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)

Bases: object

Local plasma parameters at a single radial location.

All quantities are in SI / conventional tokamak units.

Parameters:

• rho (float)

• te_kev (float)

• ti_kev (float)

• ne_19 (float)

• grad_te (float)

• grad_ti (float)

• grad_ne (float)

• q (float)

• s_hat (float)

• beta_e (float)

• r_major_m (float)

• a_minor_m (float)

• b_tesla (float)

• z_eff (float)

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

scpn_fusion.core.neural_transport.critical_gradient_model(inp, *, stiffness=2.0)

Reduced multichannel gyrokinetic closure used as analytic fallback.

Return type:

TransportFluxes

Parameters:

• inp (TransportInputs)

• stiffness (float)

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)

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:

• rho (ndarray[Any, dtype[float64]])

• te (ndarray[Any, dtype[float64]])

• ti (ndarray[Any, dtype[float64]])

• ne (ndarray[Any, dtype[float64]])

• q_profile (ndarray[Any, dtype[float64]])

• s_hat_profile (ndarray[Any, dtype[float64]])

• r_major (float)

• a_minor (float)

• b_toroidal (float)

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)

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)

Advance one deterministic spectral turbulence step and return the field.

Return type:

ndarray

Parameters:

• dt (float)

• damping (float)

class scpn_fusion.core.fno_turbulence_suppressor.FNO_Controller(weights_path=None, *, allow_legacy=False)

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)

Load explicitly enabled legacy JAX-FNO weights from an .npz file.

Return type:

None

Parameters:

path (str)

predict_and_suppress(field)

Predict a bounded suppression command and postprocessed field response.

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)

Run the seeded turbulence suppression loop and return deterministic metrics.

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).

class scpn_fusion.core.fno_training.MultiLayerFNO(modes=12, width=32, n_layers=4, seed=42)

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:

• modes (int)

• width (int)

• n_layers (int)

• seed (int)

forward_with_hidden(x_field)

Return the projected field and final hidden representation.

Return type:

Tuple[ndarray[Any, dtype[float64]], ndarray[Any, dtype[float64]]]

Parameters:

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

forward(x_field)

Evaluate the FNO field-to-field surrogate for one input field.

Return type:

ndarray[Any, dtype[float64]]

Parameters:

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

save_weights(path)

Serialise FNO architecture metadata and NumPy weights to path.

Return type:

None

Parameters:

path (str | Path)

load_weights(path)

Load FNO architecture metadata and NumPy weights from path.

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)

Train MultiLayerFNO with pure NumPy.

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

Return type:

Dict[str, object]

Parameters:

• n_samples (int)

• epochs (int)

• lr (float)

• modes (int)

• width (int)

• save_path (str | Path)

• batch_size (int)

• seed (int)

• patience (int)

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)

Compatibility wrapper over extracted multi-regime training runtime.

Return type:

Dict[str, object]

Parameters:

• n_samples (int)

• epochs (int)

• lr (float)

• modes (int)

• width (int)

• save_path (str | Path)

• batch_size (int)

• seed (int)

• patience (int)

• regime_weights (Dict[str, float] | None)

Turbulence Oracle

Turbulence oracle for reduced-order edge-chaos forecasting.

class scpn_fusion.core.turbulence_oracle.DriftWavePhysics(N=64, seed=None)

Bases: object

2D Hasegawa-Wakatani solver for plasma edge turbulence.

Variables: phi (electrostatic potential), n (density fluctuation). Pass seed for reproducible stochastic initial conditions.

bracket(A_k, B_k)

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

static spectral_dissipation_multiplier(k2)

Return the Hasegawa-Wakatani fourth-order hyperviscous multiplier.

step()

Runge-Kutta 4th Order Step with Stability Clamp

class scpn_fusion.core.turbulence_oracle.OracleESN(input_dim, reservoir_size=500, spectral_radius=0.95, seed=None)

Bases: object

Echo State Network (Reservoir Computing). Specialized for predicting chaotic time series. Pass seed for reproducible reservoir construction.

train(inputs, targets)

Ridge Regression training

predict(u_current, steps=50)

Generate closed-loop ESN predictions from an initial state.

scpn_fusion.core.turbulence_oracle.run_turbulence_oracle()

Run an end-to-end Drift-Wave + ESN prediction and save diagnostic plot.

Fusion Ignition Simulator

Zero-dimensional ignition and dynamic burn calculations for equilibrium states.

exception scpn_fusion.core.fusion_ignition_sim.BurnPhysicsError

Bases: RuntimeError, FusionCoreError

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

class scpn_fusion.core.fusion_ignition_sim.FusionBurnPhysics(config_path)

Bases: FusionKernel

Extend the Grad-Shafranov solver with thermonuclear burn physics.

Computes fusion power, alpha heating, and the Q factor from the magnetic equilibrium.

Parameters:

config_path (str)

bosch_hale_dt(T_keV)

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

Accepts a scalar temperature or a temperature profile array, matching the flux-grid usage in calculate_thermodynamics().

Return type:

float | ndarray[Any, dtype[float64]]

Parameters:

T_keV (float | ndarray[Any, dtype[float64]])

calculate_thermodynamics(P_aux_MW=50.0)

Map the magnetic equilibrium to thermodynamics and fusion power.

Parameters:

P_aux_MW (float) – External heating power (NBI/ECRH) in megawatts.

Return type:

dict[str, float]

scpn_fusion.core.fusion_ignition_sim.run_ignition_experiment()

Run the standalone auxiliary-power ignition scan and write its plot.

Return type:

None

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)

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)

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)

Compute the 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()

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)

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)

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

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

Return type:

dict[str, object]

Parameters:

• R0 (float)

• a (float)

• B_t (float)

• I_p (float)

• kappa (float)

Divertor Thermal Simulation

Reduced divertor heat-exhaust model with tungsten and liquid-lithium lanes.

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

Bases: object

Simulate Heat Exhaust in a Compact Fusion Reactor.

Compares Solid Tungsten vs. Liquid Lithium Vapor Shielding.

Parameters:

• P_sol_MW (float)

• R_major (float)

• B_pol (float)

solve_2point_transport(expansion_factor=10.0, f_rad=0.5)

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:

• expansion_factor (float)

• f_rad (float)

calculate_heat_load(expansion_factor=10.0)

Calculate Peak Heat Flux using 2-Point Model Physics.

Return type:

float

Parameters:

expansion_factor (float)

simulate_tungsten()

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)

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)

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:

• flow_velocity_m_s (float)

• channel_length_m (float)

• channel_half_gap_m (float)

• density_kg_m3 (float)

• viscosity_pa_s (float)

• conductivity_s_m (float)

estimate_evaporation_rate(surface_temp_c, flow_velocity_m_s)

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)

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()

Run the standalone divertor comparison and write the diagnostic plot.

Return type:

None

Sandpile Fusion Reactor

Toy SOC/turbulence sandpile model and HJB-style control demo contracts.

class scpn_fusion.core.sandpile_fusion_reactor.TokamakSandpile(size=100)

Bases: object

Model self-organized criticality in plasma turbulence.

Implements a modified Bak-Tang-Wiesenfeld (BTW) sandpile model to simulate self-organized criticality (SOC). Reference: ‘Running Sandpile’ models for the L-H transition.

Parameters:

size (int)

drive()

Add energy (sand) to the core.

Return type:

None

relax(suppression_strength=0.0)

Run one avalanche relaxation step.

If a local slope exceeds the critical gradient, redistribute energy. suppression_strength is 0.0 to 1.0 (effect of HJB control).

Return type:

int

Parameters:

suppression_strength (float)

calculate_profile()

Reconstruct the height profile (temperature) from gradients.

Return type:

ndarray[Any, dtype[int64]]

class scpn_fusion.core.sandpile_fusion_reactor.HJB_Avalanche_Controller

Bases: object

Heuristic closed-loop controller for avalanche suppression.

Observed state: current avalanche size. Action: increase/decrease magnetic shear (Z_crit). Reward: maximise core height (confinement) minus action cost.

act(last_avalanche_size, current_core_temp)

Return a bounded suppression action from the latest avalanche size.

Parameters:

• last_avalanche_size (int) – Observed avalanche size from previous step.

• current_core_temp (float) – Current core pseudo-temperature proxy used to bias control action.

Returns:

Shear/control signal in the range [0, 1].

Return type:

float

scpn_fusion.core.sandpile_fusion_reactor.run_sandpile_simulation()

Run the sandpile control demonstration and write a PNG diagnostic report.

Return type:

None

Warm Dense Matter Engine

Whole-device-model driver for coupled transport-equilibrium-wall simulation loops.

class scpn_fusion.core.wdm_engine.WholeDeviceModel(config_path)

Bases: object

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

thomas_fermi_pressure(n_e_m3, T_eV)

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)

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)

Run the whole-device discharge timeline and collect time-series state.

Return type:

list[dict[str, float | str]]

plot_results(history)

Plot discharge evolution and export WDM summary figure.

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)

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)

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)

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[Any, dtype[float64]])

• P_fusion_percentiles (ndarray[Any, dtype[float64]])

• Q_percentiles (ndarray[Any, dtype[float64]])

• n_samples (int)

tau_E: float

P_fusion: float

Q: float

tau_E_sigma: float

P_fusion_sigma: float

Q_sigma: float

tau_E_percentiles: ndarray[Any, dtype[float64]]

P_fusion_percentiles: ndarray[Any, dtype[float64]]

Q_percentiles: ndarray[Any, dtype[float64]]

n_samples: int = 0

scpn_fusion.core.uncertainty.ipb98_tau_e(scenario, params=None)

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.

Return type:

float | ndarray[Any, dtype[float64]]

Parameters:

Ti_keV (float | ndarray[Any, dtype[float64]])

scpn_fusion.core.uncertainty.fusion_power_from_tau(scenario, tau_E)

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 the 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:

• scenario (PlasmaScenario)

• tau_E (float)

scpn_fusion.core.uncertainty.quantify_uncertainty(scenario, n_samples=10000, seed=None)

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)

Bases: object

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

Parameters:

• config_path (Optional[str | Path])

• kernel (Optional[FusionKernel])

• equilibrium_3d (Optional[VMECStyleEquilibrium3D])

• solve_equilibrium (bool)

generate_plasma_surface(resolution_toroidal=60, resolution_poloidal=60, radial_steps=512)

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[numpy.ndarray, numpy.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)

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

Return type:

VMECStyleEquilibrium3D

Parameters:

• lcfs_resolution (int)

• radial_steps (int)

• nfp (int)

• toroidal_modes (list[FourierMode3D] | None)

build_stellarator_w7x_like_equilibrium(*, lcfs_resolution=96, radial_steps=512, nfp=5, edge_ripple=0.08, vertical_ripple=0.05)

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(*, n_toroidal_samples=72)

Generate external coil centreline geometry as sampled 3D loops.

Return type:

list[dict[str, object]]

Parameters:

n_toroidal_samples (int)

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)

Create a reduced 3D field-line tracer.

Return type:

FieldLineTracer3D

Parameters:

• rotational_transform (float)

• helical_coupling_scale (float)

• radial_coupling_scale (float)

• nfp (int)

• toroidal_modes (list[FourierMode3D] | None)

• lcfs_resolution (int)

• radial_steps (int)

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)

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

Return type:

tuple[FieldLineTrace3D, dict[float, ndarray[Any, dtype[float64]]]]

Parameters:

• rho0 (float)

• theta0 (float)

• phi0 (float)

• toroidal_turns (int)

• steps_per_turn (int)

• phi_planes (list[float] | None)

• rotational_transform (float)

• helical_coupling_scale (float)

• radial_coupling_scale (float)

• nfp (int)

• toroidal_modes (list[FourierMode3D] | None)

export_obj(vertices, faces, filename='plasma.obj')

Write a triangular plasma mesh to Wavefront OBJ format.

Return type:

Path

Parameters:

• vertices (ndarray[Any, dtype[float64]])

• faces (ndarray[Any, dtype[int64]])

• filename (str | Path)

export_preview_png(vertices, faces, filename='plasma_preview.png', *, dpi=140)

Render mesh preview to a PNG using matplotlib.

Return type:

Path

Parameters:

• vertices (ndarray[Any, dtype[float64]])

• faces (ndarray[Any, dtype[int64]])

• filename (str | Path)

• dpi (int)

scpn_fusion.core.geometry_3d.main(argv=None)

Run the command-line OBJ mesh export workflow.

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)

Bases: object

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

Parameters:

• m (int)

• n (int)

• r_cos (float)

• r_sin (float)

• z_cos (float)

• z_sin (float)

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)

Bases: object

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

Parameters:

• r_axis (float)

• z_axis (float)

• a_minor (float)

• kappa (float)

• triangularity (float)

• nfp (int)

• modes (Iterable[FourierMode3D] | None)

classmethod from_axisymmetric_lcfs(lcfs_points, *, r_axis, z_axis, nfp=1, modes=None)

Infer baseline shaping from traced axisymmetric LCFS points.

Return type:

VMECStyleEquilibrium3D

Parameters:

• lcfs_points (ndarray[Any, dtype[float64]])

• r_axis (float)

• z_axis (float)

• nfp (int)

• modes (Iterable[FourierMode3D] | None)

to_vmec_like_dict()

Export a deterministic VMEC-like boundary payload.

Return type:

dict[str, object]

classmethod from_vmec_like_dict(payload)

Rebuild reduced equilibrium from VMEC-like boundary payload.

Return type:

VMECStyleEquilibrium3D

Parameters:

payload (dict[str, object])

flux_to_cylindrical(rho, theta, phi)

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

Return type:

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

Parameters:

• rho (float | ndarray[Any, dtype[float64]])

• theta (float | ndarray[Any, dtype[float64]])

• phi (float | ndarray[Any, dtype[float64]])

flux_to_cartesian(rho, theta, phi)

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

Return type:

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

Parameters:

• rho (float | ndarray[Any, dtype[float64]])

• theta (float | ndarray[Any, dtype[float64]])

• phi (float | ndarray[Any, dtype[float64]])

sample_surface(*, rho=1.0, resolution_toroidal=60, resolution_poloidal=60)

Sample a triangulated flux surface from native 3D coordinates.

Return type:

tuple[ndarray[Any, dtype[float64]], ndarray[Any, dtype[int64]]]

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)

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)

• non_decreasing_steps (int)

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)

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:

• eq (VMECStyleEquilibrium3D)

• b0_tesla (float)

• r0_major (float)

• p0_pa (float)

• j0_ma_m2 (float)

• pressure_exp (float)

• current_exp (float)

compute_force_residual(n_rho=12, n_theta=24, n_phi=16)

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:

• n_rho (int)

• n_theta (int)

• n_phi (int)

solve(*, max_iterations=200, tolerance=0.0001, learning_rate=0.01, armijo_c=0.0001, n_rho=12, n_theta=24, n_phi=16)

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:

• max_iterations (int)

• tolerance (float)

• learning_rate (float)

• armijo_c (float)

• n_rho (int)

• n_theta (int)

• n_phi (int)

3D Field-Line Tracing

Reduced 3D field-line tracing and Poincare diagnostics.

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

Bases: object

Sampled field-line trajectory in native and Cartesian coordinates.

Parameters:

• rho (ndarray[Any, dtype[float64]])

• theta (ndarray[Any, dtype[float64]])

• phi (ndarray[Any, dtype[float64]])

• xyz (ndarray[Any, dtype[float64]])

rho: ndarray[Any, dtype[float64]]

theta: ndarray[Any, dtype[float64]]

phi: ndarray[Any, dtype[float64]]

xyz: ndarray[Any, dtype[float64]]

class scpn_fusion.core.fieldline_3d.PoincareSection3D(phi_plane, xyz, rz)

Bases: object

Poincare section points for one toroidal cut plane.

Parameters:

• phi_plane (float)

• xyz (ndarray[Any, dtype[float64]])

• rz (ndarray[Any, dtype[float64]])

phi_plane: float

xyz: ndarray[Any, dtype[float64]]

rz: ndarray[Any, dtype[float64]]

class scpn_fusion.core.fieldline_3d.ToroidalAsymmetryObservables3D(n1_amp, n2_amp, n3_amp, z_n1_amp, asymmetry_index, radial_spread)

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)

n1_amp: float

n2_amp: float

n3_amp: float

z_n1_amp: float

asymmetry_index: float

radial_spread: float

as_dict()

Return controller-facing toroidal asymmetry observables by stable key.

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)

Bases: object

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

Parameters:

• equilibrium (VMECStyleEquilibrium3D)

• rotational_transform (float | Callable[[float], float])

• helical_coupling_scale (float)

• radial_coupling_scale (float)

trace_line(*, rho0=0.95, theta0=0.0, phi0=0.0, toroidal_turns=16, steps_per_turn=256)

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

Return type:

FieldLineTrace3D

Parameters:

• rho0 (float)

• theta0 (float)

• phi0 (float)

• toroidal_turns (int)

• steps_per_turn (int)

poincare_section(trace, *, phi_plane=0.0)

Compute Poincare intersection points for one toroidal plane.

Return type:

PoincareSection3D

Parameters:

• trace (FieldLineTrace3D)

• phi_plane (float)

poincare_map(trace, *, phi_planes)

Compute Poincare sections for multiple toroidal planes.

Return type:

dict[float, PoincareSection3D]

Parameters:

• trace (FieldLineTrace3D)

• phi_planes (Iterable[float])

toroidal_asymmetry_observables(trace)

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)

Bases: object

Latency summary for multigrid and SNN benchmark lanes.

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)

Bases: object

Equilibrium-kernel latency summary with nominal and faulted runs.

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)

Bases: object

Reduced runtime bridge with deterministic benchmark estimates.

Parameters:

seed (int)

static available_equilibrium_backends()

Return equilibrium latency backends available in the current runtime.

Return type:

tuple[str, ...]

benchmark(*, backend, trials=64, grid_size=64)

Benchmark deterministic multigrid and SNN kernels for one backend.

Return type:

RuntimeBenchmark

Parameters:

• backend (str)

• trials (int)

• grid_size (int)

benchmark_pair(*, trials=64, grid_size=64)

Benchmark CPU and GPU-sim lanes and return estimated speedups.

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)

Measure equilibrium-kernel latency under nominal and injected-fault inputs.

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)

Bases: object

Synthetic core-turbulence dataset bundle.

Parameters:

• features (ndarray[Any, dtype[float64]])

• chi_i (ndarray[Any, dtype[float64]])

features: ndarray[Any, dtype[float64]]

chi_i: ndarray[Any, dtype[float64]]

class scpn_fusion.core.gyro_swin_surrogate.SpeedBenchmark(gene_proxy_s_per_sample, surrogate_s_per_sample, speedup)

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)

Compute the reference synthetic turbulence law used as the CI benchmark target.

Return type:

ndarray[Any, dtype[float64]]

Parameters:

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

scpn_fusion.core.gyro_swin_surrogate.generate_synthetic_gyrokinetic_dataset(seed, samples)

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

Return type:

TurbulenceDataset

Parameters:

• seed (int)

• samples (int)

class scpn_fusion.core.gyro_swin_surrogate.GyroSwinLikeSurrogate(hidden_dim=64, ridge=0.0005, seed=42)

Bases: object

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

Parameters:

• hidden_dim (int)

• ridge (float)

• seed (int)

fit(features, targets)

Fit random-feature ridge weights to turbulence targets.

Return type:

None

Parameters:

• features (ndarray[Any, dtype[float64]])

• targets (ndarray[Any, dtype[float64]])

predict(features)

Predict non-negative ion heat diffusivity for feature rows.

Return type:

ndarray[Any, dtype[float64]]

Parameters:

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

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

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[Any, dtype[float64]]

Parameters:

• features (ndarray[Any, dtype[float64]])

• iterations (int)

scpn_fusion.core.gyro_swin_surrogate.rmse_percent(y_true, y_pred)

Return RMSE as a percentage of mean absolute reference value.

Return type:

float

Parameters:

• y_true (ndarray[Any, dtype[float64]])

• y_pred (ndarray[Any, dtype[float64]])

scpn_fusion.core.gyro_swin_surrogate.benchmark_speedup(features, surrogate, min_baseline_s=0.15, min_surrogate_s=0.02)

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

Return type:

SpeedBenchmark

Parameters:

• features (ndarray[Any, dtype[float64]])

• surrogate (GyroSwinLikeSurrogate)

• min_baseline_s (float)

• min_surrogate_s (float)

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)

Bases: object

Synthetic magnetic-shadow training features and target fractions.

Parameters:

• features (ndarray[Any, dtype[float64]])

• shadow_fraction (ndarray[Any, dtype[float64]])

features: ndarray[Any, dtype[float64]]

shadow_fraction: ndarray[Any, dtype[float64]]

scpn_fusion.core.heat_ml_shadow_surrogate.synthetic_shadow_reference(features)

Synthetic reference law for divertor magnetic-shadow fraction.

Return type:

ndarray[Any, dtype[float64]]

Parameters:

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

scpn_fusion.core.heat_ml_shadow_surrogate.generate_shadow_dataset(seed, samples)

Generate a deterministic synthetic HEAT-ML shadow training dataset.

Return type:

ShadowDataset

Parameters:

• seed (int)

• samples (int)

class scpn_fusion.core.heat_ml_shadow_surrogate.HeatMLShadowSurrogate(ridge=0.0001)

Bases: object

Compact polynomial HEAT-ML surrogate with deterministic ridge fit.

Parameters:

ridge (float)

fit(features, target)

Fit ridge-regularised polynomial weights to shadow fractions.

Return type:

None

Parameters:

• features (ndarray[Any, dtype[float64]])

• target (ndarray[Any, dtype[float64]])

fit_synthetic(seed=42, samples=2048)

Fit the surrogate on a deterministic synthetic reference dataset.

Return type:

None

Parameters:

• seed (int)

• samples (int)

predict_shadow_fraction(features)

Predict clipped magnetic-shadow fractions for divertor feature rows.

Return type:

ndarray[Any, dtype[float64]]

Parameters:

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

predict_divertor_flux(q_div_baseline_w_m2, features)

Apply predicted magnetic-shadow attenuation to baseline divertor flux.

Return type:

ndarray[Any, dtype[float64]]

Parameters:

• q_div_baseline_w_m2 (ndarray[Any, dtype[float64]] | float)

• features (ndarray[Any, dtype[float64]])

scpn_fusion.core.heat_ml_shadow_surrogate.rmse_percent(y_true, y_pred)

Return RMSE as a percentage of mean absolute reference magnitude.

Return type:

float

Parameters:

• y_true (ndarray[Any, dtype[float64]])

• y_pred (ndarray[Any, dtype[float64]])

scpn_fusion.core.heat_ml_shadow_surrogate.benchmark_inference_seconds(model, samples=200000)

Measure vectorised shadow-fraction inference time for synthetic samples.

Return type:

float

Parameters:

• model (HeatMLShadowSurrogate)

• samples (int)

Rust Compatibility Layer

Backward-compatibility layer for the optional Rust acceleration backend.

Imports from scpn_fusion_rs when available, falling back to pure-Python implementations otherwise.

Usage:

from scpn_fusion.core._rust_compat import FusionKernel, RUST_BACKEND

class scpn_fusion.core._rust_compat.RustAcceleratedKernel(config_path)

Bases: object

Drop-in wrapper around the Rust PyFusionKernel.

Mirrors the Python FusionKernel attribute interface (.Psi, .R, .Z, .RR, .ZZ, .cfg, etc.) and delegates the equilibrium solve to Rust for ~20x speedup while keeping all attribute accesses compatible with downstream code.

Parameters:

config_path (str | Path)

solve_equilibrium()

Solve Grad-Shafranov equilibrium via Rust backend.

Return type:

Any

compute_b_field()

Compute magnetic field components from Psi gradient.

Return type:

None

find_x_point(Psi)

Locate the null point (B=0) using local minimisation.

Matches the Python FusionKernel.find_x_point() interface.

Return type:

tuple[tuple[float, float], float]

Parameters:

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

calculate_thermodynamics(p_aux_mw)

Calculate thermodynamics via Rust backend.

Return type:

Any

Parameters:

p_aux_mw (float)

calculate_vacuum_field()

Compute vacuum field with Python reference implementation.

Return type:

ndarray[Any, dtype[float64]]

set_solver_method(method)

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

Return type:

None

Parameters:

method (str)

solver_method()

Get current solver method name.

Return type:

str

save_results(filename='equilibrium_nonlinear.npz')

Save current state to .npz file.

Return type:

None

Parameters:

filename (str)

scpn_fusion.core._rust_compat.rust_shafranov_bv(*args, **kwargs)

Return type:

Any

Parameters:

• args (Any)

• kwargs (Any)

scpn_fusion.core._rust_compat.rust_solve_coil_currents(*args, **kwargs)

Return type:

Any

Parameters:

• args (Any)

• kwargs (Any)

scpn_fusion.core._rust_compat.rust_measure_magnetics(*args, **kwargs)

Return type:

Any

Parameters:

• args (Any)

• kwargs (Any)

scpn_fusion.core._rust_compat.rust_simulate_tearing_mode(steps, seed=None)

Return type:

Any

Parameters:

• steps (int)

• seed (int | None)

class scpn_fusion.core._rust_compat.RustSnnPool(n_neurons=50, gain=10.0, window_size=20, *, allow_numpy_fallback=True, seed=42)

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)

Process error through SNN pool and return scalar control output.

Return type:

float

Parameters:

error (float)

property n_neurons: int

Number of neurons in the active pool backend.

property gain: float

Controller gain used by the active pool backend.

property backend: str

Name of the active SNN pool backend.

class scpn_fusion.core._rust_compat.RustSnnController(target_r=6.2, target_z=0.0, *, allow_numpy_fallback=True, seed=42)

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)

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

Return type:

tuple[float, float]

Parameters:

• measured_r (float)

• measured_z (float)

property target_r: float

Target major-radius position passed to the active backend.

property target_z: float

Target vertical position passed to the active backend.

property backend: str

Name of the active SNN controller backend.

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)

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:

• source (ndarray[Any, dtype[float64]])

• psi_bc (ndarray[Any, dtype[float64]])

• r_min (float)

• r_max (float)

• z_min (float)

• z_max (float)

• nr (int)

• nz (int)

• tol (float)

• max_cycles (int)

Research Bridges

Bridge harness between plasma-equilibrium outputs and Opentrons protocol synthesis.

class scpn_fusion.core.lazarus_bridge.LazarusBridge(config_path)

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()

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

generate_protocol(resonance_score)

Creates an Opentrons Python Protocol based on the physics state.

run_bridge_simulation()

Run end-to-end bridge scenario from plasma equilibrium to protocol output.

visualize_bridge(score)

Save a symbolic bridge visualization for the current resonance score.

Bridge that maps tokamak control dynamics into VIBRANA sonification settings.

class scpn_fusion.core.vibrana_bridge.VibranaFusionBridge(config_path)

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)

The Transduction Function: Physics -> Acoustics

run_sonification_session(duration_steps=100)

Execute an autonomous sonification run and collect audio-control telemetry.

visualize_soundscape(log)

Generate and save a 2-panel summary of sonification dynamics.

Runtime integration wrapper for the optional local Quantum Lab scripts.

scpn_fusion.core.quantum_bridge.run_quantum_suite(*, base_path=None, script_timeout_seconds=1800.0)

Execute the staged Quantum Lab workflow scripts in sequence.

Parameters:

• base_path (str | Path | None) – Optional explicit path to the Quantum Lab directory. If omitted, the default repository-local path is used.

• script_timeout_seconds (float) – Per-script timeout in seconds. Must be finite and strictly positive.

Returns:

Execution report with success flag, resolved base path, and script names.

Return type:

dict[str, object]

Raises:

• FileNotFoundError – If the Quantum Lab directory or required scripts are missing.

• RuntimeError – If any script fails or times out.

• ValueError – If the timeout value is invalid.

Current Diffusion

One-dimensional current-diffusion and safety-factor profile evolution tools.

class scpn_fusion.core.current_diffusion.FluxEvolutionTrajectory(time_s, rho, psi, hall_drive, resistive_loss, source, damping_rate, source_increment, damping_decrement, update_residual, dt_s)

Bases: object

Time history for the non-adiabatic MIF/FRC flux carrier equation.

Arrays use SI units. psi has shape (n_steps + 1, n_rho) and stores the evolved poloidal-flux profile. hall_drive stores R_null E_theta, resistive_loss stores eta J_theta, and damping_rate stores 1 / tau_psi at the same time/grid points. source_increment and damping_decrement have shape (n_steps, n_rho) and close the exact discrete balance psi[n+1] = psi[n] - damping_decrement[n] + source_increment[n].

Parameters:

• time_s (ndarray[Any, dtype[float64]])

• rho (ndarray[Any, dtype[float64]])

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

• hall_drive (ndarray[Any, dtype[float64]])

• resistive_loss (ndarray[Any, dtype[float64]])

• source (ndarray[Any, dtype[float64]])

• damping_rate (ndarray[Any, dtype[float64]])

• source_increment (ndarray[Any, dtype[float64]])

• damping_decrement (ndarray[Any, dtype[float64]])

• update_residual (ndarray[Any, dtype[float64]])

• dt_s (float)

time_s: ndarray[Any, dtype[float64]]

rho: ndarray[Any, dtype[float64]]

psi: ndarray[Any, dtype[float64]]

hall_drive: ndarray[Any, dtype[float64]]

resistive_loss: ndarray[Any, dtype[float64]]

source: ndarray[Any, dtype[float64]]

damping_rate: ndarray[Any, dtype[float64]]

source_increment: ndarray[Any, dtype[float64]]

damping_decrement: ndarray[Any, dtype[float64]]

update_residual: ndarray[Any, dtype[float64]]

dt_s: float

scpn_fusion.core.current_diffusion.solve_flux_evolution_nonadiabatic(rho, psi0, *, tau_psi_fn, R_null_t, E_theta_t, eta_spitzer_fn, J_theta_t, dt, n_steps)

Evolve the pulsed non-adiabatic flux constraint.

The implemented carrier is the Ono-style flattened FRC flux equation

dpsi/dt = -psi / tau_psi + R_null E_theta - eta J_theta.

Each time step samples the source at both endpoints and advances the local linear damping term analytically using the midpoint damping rate. This is exact for constant damping/source, second-order for smooth time-dependent drive terms, and reduces exactly to constant psi when Hall drive, resistive loss, and damping are all zero.

Return type:

FluxEvolutionTrajectory

Parameters:

• rho (ndarray[Any, dtype[float64]])

• psi0 (ndarray[Any, dtype[float64]])

• tau_psi_fn (Callable[[...], ndarray[Any, dtype[float64]] | float])

• R_null_t (Callable[[float], float])

• E_theta_t (Callable[[float, ndarray[Any, dtype[float64]]], ndarray[Any, dtype[float64]] | float])

• eta_spitzer_fn (Callable[[ndarray[Any, dtype[float64]]], ndarray[Any, dtype[float64]] | float])

• J_theta_t (Callable[[float, ndarray[Any, dtype[float64]]], ndarray[Any, dtype[float64]] | float])

• dt (float)

• n_steps (int)

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

Sauter neoclassical parallel resistivity [Ohm-m].

Return type:

float

Parameters:

• Te_keV (float)

• ne_19 (float)

• Z_eff (float)

• epsilon (float)

• q (float)

• R0 (float)

scpn_fusion.core.current_diffusion.q_from_psi(rho, psi, R0, a, B0)

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

Return type:

ndarray[Any, dtype[float64]]

Parameters:

• rho (ndarray[Any, dtype[float64]])

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

• R0 (float)

• a (float)

• B0 (float)

scpn_fusion.core.current_diffusion.resistive_diffusion_time(a, eta)

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

Return type:

float

Parameters:

• a (float)

• eta (float)

class scpn_fusion.core.current_diffusion.CurrentDiffusionSolver(rho, R0, a, B0)

Bases: object

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

Parameters:

• rho (FloatArray)

• R0 (float)

• a (float)

• B0 (float)

step(dt, Te, ne, Z_eff, j_bs, j_cd, j_ext=None)

Advance poloidal flux psi by one timestep dt.

Return type:

ndarray[Any, dtype[float64]]

Parameters:

• dt (float)

• Te (ndarray[Any, dtype[float64]])

• ne (ndarray[Any, dtype[float64]])

• Z_eff (float)

• j_bs (ndarray[Any, dtype[float64]])

• j_cd (ndarray[Any, dtype[float64]])

• j_ext (ndarray[Any, dtype[float64]] | None)

Current Drive Sources

Reduced ECCD, NBI, and LHCD current-drive source models.

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

Bases: object

Electron Cyclotron Current Drive (ECCD).

Parameters:

• P_ec_MW (float)

• rho_dep (float)

• sigma_rho (float)

• eta_cd (float)

P_absorbed(rho)

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

Return type:

ndarray[Any, dtype[float64]]

Parameters:

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

j_cd(rho, ne_19, Te_keV)

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

Return type:

ndarray[Any, dtype[float64]]

Parameters:

• rho (ndarray[Any, dtype[float64]])

• ne_19 (ndarray[Any, dtype[float64]])

• Te_keV (ndarray[Any, dtype[float64]])

class scpn_fusion.core.current_drive.NBISource(P_nbi_MW, E_beam_keV, rho_tangency, sigma_rho=0.15)

Bases: object

Neutral Beam Injection (NBI).

Parameters:

• P_nbi_MW (float)

• E_beam_keV (float)

• rho_tangency (float)

• sigma_rho (float)

P_heating(rho)

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

Return type:

ndarray[Any, dtype[float64]]

Parameters:

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

j_cd(rho, ne_19, Te_keV, Ti_keV)

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

Return type:

ndarray[Any, dtype[float64]]

Parameters:

• rho (ndarray[Any, dtype[float64]])

• ne_19 (ndarray[Any, dtype[float64]])

• Te_keV (ndarray[Any, dtype[float64]])

• Ti_keV (ndarray[Any, dtype[float64]])

class scpn_fusion.core.current_drive.LHCDSource(P_lh_MW, rho_dep, sigma_rho, eta_cd=0.15)

Bases: object

Lower Hybrid Current Drive (LHCD).

Parameters:

• P_lh_MW (float)

• rho_dep (float)

• sigma_rho (float)

• eta_cd (float)

P_absorbed(rho)

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

Return type:

ndarray[Any, dtype[float64]]

Parameters:

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

j_cd(rho, ne_19, Te_keV)

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

Return type:

ndarray[Any, dtype[float64]]

Parameters:

• rho (ndarray[Any, dtype[float64]])

• ne_19 (ndarray[Any, dtype[float64]])

• Te_keV (ndarray[Any, dtype[float64]])

class scpn_fusion.core.current_drive.CurrentDriveMix(a=1.0)

Bases: object

Combines multiple current drive sources.

Parameters:

a (float)

add_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)

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

Return type:

ndarray[Any, dtype[float64]]

Parameters:

• rho (ndarray[Any, dtype[float64]])

• ne (ndarray[Any, dtype[float64]])

• Te (ndarray[Any, dtype[float64]])

• Ti (ndarray[Any, dtype[float64]])

total_heating_power(rho)

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

Return type:

ndarray[Any, dtype[float64]]

Parameters:

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

total_driven_current(rho, ne, Te, Ti, *, elongation=1.0)

Integrate total driven current [A] over circular or elongated area.

Return type:

float

Parameters:

• rho (ndarray[Any, dtype[float64]])

• ne (ndarray[Any, dtype[float64]])

• Te (ndarray[Any, dtype[float64]])

• Ti (ndarray[Any, dtype[float64]])

• elongation (float | ndarray[Any, dtype[float64]])

Sawtooth Dynamics

Porcelli-trigger and Kadomtsev-reconnection reduced-order sawtooth contracts.

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

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)

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)

Bases: object

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

Parameters:

• rho (FloatArray)

• s_crit (float)

find_q1_radius(q)

Find the normalized radius where q=1.

Return type:

float | None

Parameters:

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

check_trigger(q, shear)

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

Return type:

bool

Parameters:

• q (ndarray[Any, dtype[float64]])

• shear (ndarray[Any, dtype[float64]])

scpn_fusion.core.sawtooth.kadomtsev_crash(rho, T, n, q, R0, a)

Apply Kadomtsev reconnection inside the q=1 mixing radius.

Returns (T_new, n_new, q_new, rho_1, rho_mix). The mixing flattens the density and temperature inside rho_mix so that both the particle content and the thermal energy 1.5 n T integrated over the mixing volume are conserved (Kadomtsev 1975; Wesson, Tokamaks, §7.6). The reconnected core safety factor is reset just above unity.

Return type:

tuple[ndarray[Any, dtype[float64]], ndarray[Any, dtype[float64]], ndarray[Any, dtype[float64]], float, float]

Parameters:

• rho (ndarray[Any, dtype[float64]])

• T (ndarray[Any, dtype[float64]])

• n (ndarray[Any, dtype[float64]])

• q (ndarray[Any, dtype[float64]])

• R0 (float)

• a (float)

class scpn_fusion.core.sawtooth.SawtoothCycler(rho, R0, a, s_crit=0.1)

Bases: object

Tracks and triggers sawtooth crashes.

Parameters:

• rho (FloatArray)

• R0 (float)

• a (float)

• s_crit (float)

step(dt, q, shear, T, n)

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

Return type:

SawtoothEvent | None

Parameters:

• dt (float)

• q (ndarray[Any, dtype[float64]])

• shear (ndarray[Any, dtype[float64]])

• T (ndarray[Any, dtype[float64]])

• n (ndarray[Any, dtype[float64]])

NTM Island Dynamics

Neoclassical tearing mode search/suppression helpers and island dynamics.

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

Bases: object

Represents a q=m/n rational surface.

Parameters:

• rho (float)

• r_s (float)

• m (int)

• n (int)

• q (float)

• shear (float)

rho: float

r_s: float

m: int

n: int

q: float

shear: float

scpn_fusion.core.ntm_dynamics.eccd_stabilization_factor(d_cd, w)

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

Return type:

float

Parameters:

• d_cd (float)

• w (float)

scpn_fusion.core.ntm_dynamics.find_rational_surfaces(q, rho, a, m_max=5, n_max=3)

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

Return type:

list[RationalSurface]

Parameters:

• q (ndarray[Any, dtype[float64]])

• rho (ndarray[Any, dtype[float64]])

• a (float)

• m_max (int)

• n_max (int)

class scpn_fusion.core.ntm_dynamics.NTMIslandDynamics(r_s, m, n, a, R0, B0, Delta_prime_0=None)

Bases: object

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

Parameters:

• r_s (float)

• m (int)

• n (int)

• a (float)

• R0 (float)

• B0 (float)

• Delta_prime_0 (float | None)

delta_prime_model(w)

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)

Evaluate the RHS of the Modified Rutherford Equation.

Return type:

float

Parameters:

• w (float)

• j_bs (float)

• j_phi (float)

• j_cd (float)

• eta (float)

• w_d (float)

• w_pol (float)

• d_cd (float)

evolve(w0, t_span, dt, j_bs, j_phi, j_cd, eta, w_d=0.001, w_pol=0.0005, d_cd=0.05)

Integrate w(t) using RK4.

Return type:

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

Parameters:

• w0 (float)

• t_span (tuple[float, float])

• dt (float)

• j_bs (float)

• j_phi (float)

• j_cd (float)

• eta (float)

• w_d (float)

• w_pol (float)

• d_cd (float)

class scpn_fusion.core.ntm_dynamics.NTMController(w_onset=0.02, w_target=0.005)

Bases: object

Monitors and controls NTM islands via ECCD.

Parameters:

• w_onset (float)

• w_target (float)

step(w, rho_rs, max_power=20.0)

Return requested ECCD power in MW.

Return type:

float

Parameters:

• w (float)

• rho_rs (float)

• max_power (float)

SOL Two-Point Model

SOL two-point closure models using Eich scaling and simple Spitzer-Harm balance.

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

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)

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)

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

Return type:

float

Parameters:

• P_SOL_MW (float)

• R0 (float)

• B_pol (float)

• epsilon (float)

scpn_fusion.core.sol_model.peak_target_heat_flux(P_SOL_MW, R0, lambda_q_m, f_expansion=5.0, alpha_deg=3.0)

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

Return type:

float

Parameters:

• P_SOL_MW (float)

• R0 (float)

• lambda_q_m (float)

• f_expansion (float)

• alpha_deg (float)

class scpn_fusion.core.sol_model.TwoPointSOL(R0, a, q95, B_pol, kappa=1.0)

Bases: object

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

Parameters:

• R0 (float)

• a (float)

• q95 (float)

• B_pol (float)

• kappa (float)

solve(P_SOL_MW, n_u_19, f_rad=0.0)

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

Return type:

SOLSolution

Parameters:

• P_SOL_MW (float)

• n_u_19 (float)

• f_rad (float)