# SPDX-License-Identifier: AGPL-3.0-or-later
# Commercial license available
# © Concepts 1996–2026 Miroslav Šotek. All rights reserved.
# © Code 2020–2026 Miroslav Šotek. All rights reserved.
# ORCID: 0009-0009-3560-0851
# Contact: www.anulum.li | protoscience@anulum.li
# SCPN Fusion Core — Sawtooth Model (Porcelli trigger + Kadomtsev reconnection)
"""Porcelli-trigger and Kadomtsev-reconnection reduced-order sawtooth contracts."""
from __future__ import annotations
from dataclasses import dataclass
import numpy as np
from numpy.typing import NDArray
from scipy.integrate import trapezoid
FloatArray = NDArray[np.float64]
[docs]
@dataclass
class SawtoothEvent:
"""Record of a single sawtooth crash: timing, radii, and energy release."""
crash_time: float
rho_1: float
rho_mix: float
T_drop: float
seed_energy: float
[docs]
class SawtoothMonitor:
"""Monitors the q-profile for Porcelli-like sawtooth trigger conditions."""
def __init__(self, rho: FloatArray, s_crit: float = 0.1):
self.rho = rho
self.s_crit = s_crit
[docs]
def find_q1_radius(self, q: FloatArray) -> float | None:
"""Find the normalized radius where q=1."""
if q[0] >= 1.0 or np.min(q) >= 1.0:
return None
q_diff = q - 1.0
crossings = np.where(np.diff(np.sign(q_diff)))[0]
if len(crossings) == 0:
return None
idx = crossings[0]
r1, r2 = self.rho[idx], self.rho[idx + 1]
q1, q2 = q[idx], q[idx + 1]
if q1 == q2:
return float(r1)
frac = (1.0 - q1) / (q2 - q1)
rho_1 = r1 + frac * (r2 - r1)
return float(rho_1)
[docs]
def check_trigger(self, q: FloatArray, shear: FloatArray) -> bool:
"""Porcelli-like trigger based on shear at q=1."""
rho_1 = self.find_q1_radius(q)
if rho_1 is None:
return False
idx = np.searchsorted(self.rho, rho_1)
if idx == 0:
s1 = shear[0]
elif idx >= len(self.rho):
s1 = shear[-1]
else:
r1, r2 = self.rho[idx - 1], self.rho[idx]
s_1, s_2 = shear[idx - 1], shear[idx]
if r1 == r2:
s1 = s_1
else:
frac = (rho_1 - r1) / (r2 - r1)
s1 = s_1 + frac * (s_2 - s_1)
return bool(s1 > self.s_crit)
[docs]
def kadomtsev_crash(
rho: FloatArray, T: FloatArray, n: FloatArray, q: FloatArray, R0: float, a: float
) -> tuple[FloatArray, FloatArray, FloatArray, float, float]:
"""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.
"""
monitor = SawtoothMonitor(rho)
rho_1 = monitor.find_q1_radius(q)
if rho_1 is None:
return T.copy(), n.copy(), q.copy(), 0.0, 0.0
# Helical flux proxy: dpsi*/drho = rho (1/q - 1)
integrand = rho * (1.0 / np.maximum(q, 1e-6) - 1.0)
psi_star = np.zeros_like(rho)
for i in range(1, len(rho)):
psi_star[i] = psi_star[i - 1] + 0.5 * (integrand[i - 1] + integrand[i]) * (
rho[i] - rho[i - 1]
)
idx_1 = np.searchsorted(rho, rho_1)
rho_mix = rho[-1]
for i in range(idx_1, len(rho)):
if psi_star[i] <= 0.0:
if i > 0 and psi_star[i - 1] > 0.0:
frac = psi_star[i - 1] / (psi_star[i - 1] - psi_star[i])
rho_mix = rho[i - 1] + frac * (rho[i] - rho[i - 1])
else:
rho_mix = rho[i]
break
idx_mix = np.searchsorted(rho, rho_mix)
if idx_mix == 0:
return T.copy(), n.copy(), q.copy(), rho_1, rho_mix
rho_inner = rho[:idx_mix]
def volume_average(prof: FloatArray) -> float:
prof_inner = prof[:idx_mix]
vol_int = trapezoid(prof_inner * rho_inner, rho_inner)
vol_tot = trapezoid(rho_inner, rho_inner)
if vol_tot == 0:
return float(prof_inner[0])
return float(vol_int / vol_tot)
# Kadomtsev reconnection conserves both particle number and thermal energy
# inside the mixing radius (Kadomtsev 1975; Wesson, Tokamaks, §7.6). Density
# flattens to the particle-conserving volume average, and temperature to the
# energy-conserving pressure average T_mix = <n T> / <n>, so that the mixed
# core leaves 1.5 n T integrated over the mixing volume invariant. Averaging
# T independently as <T> would not conserve thermal energy because
# <n T> != <n> <T> for correlated profiles.
n_mix = volume_average(n)
if n_mix > 0.0:
T_mix = volume_average(n * T) / n_mix
else:
T_mix = volume_average(T)
T_new = T.copy()
n_new = n.copy()
q_new = q.copy()
T_new[:idx_mix] = T_mix
n_new[:idx_mix] = n_mix
q_new[:idx_mix] = 1.01
return T_new, n_new, q_new, rho_1, rho_mix
[docs]
class SawtoothCycler:
"""Tracks and triggers sawtooth crashes."""
def __init__(self, rho: FloatArray, R0: float, a: float, s_crit: float = 0.1):
self.rho = rho
self.R0 = R0
self.a = a
self.monitor = SawtoothMonitor(rho, s_crit)
self.time = 0.0
[docs]
def step(
self, dt: float, q: FloatArray, shear: FloatArray, T: FloatArray, n: FloatArray
) -> SawtoothEvent | None:
"""Advance by dt; trigger Kadomtsev crash if shear at q=1 exceeds s_crit."""
self.time += dt
if self.monitor.check_trigger(q, shear):
T_core_old = T[0]
e_charge = 1.602e-19
T_new, n_new, q_new, rho_1, rho_mix = kadomtsev_crash(
self.rho, T, n, q, self.R0, self.a
)
np.copyto(T, T_new)
np.copyto(n, n_new)
np.copyto(q, q_new)
T_drop = T_core_old - T[0]
# The reconnection conserves total thermal energy across the mixing
# radius, redistributing it from the q<1 core into the annulus. The
# seed energy that feeds downstream reconnection coupling is measured
# as the thermal energy removed from the q=1 core volume.
core_energy_drop = (
1.5
* (n[0] * 1e19)
* (T_drop * 1e3 * e_charge)
* (2 * np.pi**2 * self.R0 * (rho_1 * self.a) ** 2)
)
return SawtoothEvent(
crash_time=self.time,
rho_1=rho_1,
rho_mix=rho_mix,
T_drop=T_drop,
seed_energy=float(core_energy_drop),
)
return None