# 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 Streamlit application.
This module defines the interactive control-room experience for the fusion core
simulator. The dashboard exposes five tabs covering equilibrium solve, ignition
performance, nuclear wall loading, power plant balancing, and disruption-shot
replay.
"""
from __future__ import annotations
import numpy as np
import matplotlib.pyplot as plt
import streamlit as st
import importlib
from pathlib import Path
from typing import Any, Final
from scpn_fusion._data_paths import data_root, default_iter_config_path
from scpn_fusion.ui.security_headers import install_tornado_security_headers
install_tornado_security_headers()
try:
FusionKernel: Any = importlib.import_module("scpn_fusion.core._rust_compat").FusionKernel
except (AttributeError, ImportError): # pragma: no cover - platform fallback
FusionKernel = importlib.import_module("scpn_fusion.core.fusion_kernel").FusionKernel
FusionBurnPhysics: Any = importlib.import_module(
"scpn_fusion.core.fusion_ignition_sim"
).FusionBurnPhysics
from scpn_fusion.nuclear.nuclear_wall_interaction import NuclearEngineeringLab
from scpn_fusion.engineering.balance_of_plant import PowerPlantModel
APP_VERSION: Final[str] = "3.4.0"
APP_TITLE: Final[str] = "⚛️ SCPN Fusion Reactor Control Room"
APP_SUBTITLE: Final[str] = "Digital Twin & Engineering Suite"
def _resolve_config_path(config_filename: str = "iter_config.json") -> str:
"""Return the best-effort path to the reactor configuration file.
The default ITER configuration resolves through package data. Custom
filenames first check the active data root, then fall back to the working
directory for local development.
Parameters
----------
config_filename:
Name of the JSON configuration file.
Returns
-------
str
Filesystem path that can be passed to simulation classes.
"""
if config_filename == "iter_config.json":
return str(default_iter_config_path())
data_candidate = data_root() / "validation" / config_filename
if data_candidate.exists():
return str(data_candidate)
return str(Path.cwd() / config_filename)
[docs]
def main() -> None: # pragma: no cover - exercised by Streamlit AppTest.
"""Run the dashboard and render all interactive modules."""
st.set_page_config(page_title="SCPN Fusion Reactor", layout="wide", page_icon="⚛️")
config_path = _resolve_config_path()
st.title(APP_TITLE)
st.markdown(f"### {APP_SUBTITLE} v{APP_VERSION}")
st.sidebar.header("Reactor Parameters")
reactor_size = st.sidebar.slider("Major Radius (m)", 3.0, 9.0, 6.2)
plasma_current = st.sidebar.slider("Target Current (MA)", 1.0, 20.0, 15.0)
aux_heating = st.sidebar.slider("Auxiliary Heating (MW)", 0.0, 100.0, 50.0)
st.sidebar.caption(f"Configured major radius: {reactor_size:.2f} m")
tab1, tab2, tab3, tab4, tab5 = st.tabs(
["Plasma Physics", "Ignition & Q", "Nuclear Engineering", "Power Plant", "Shot Replay"]
)
with tab1:
st.header("Grad-Shafranov Equilibrium")
if st.button("Solve Equilibrium"):
with st.spinner("Solving non-linear MHD equations..."):
kernel = FusionKernel(config_path)
kernel.cfg["physics"]["plasma_current_target"] = float(plasma_current)
kernel.solve_equilibrium()
col1, col2 = st.columns(2)
with col1:
fig, ax = plt.subplots()
ax.contour(kernel.RR, kernel.ZZ, kernel.Psi, levels=20, colors="black")
im = ax.imshow(
kernel.J_phi,
extent=(
float(kernel.R[0]),
float(kernel.R[-1]),
float(kernel.Z[0]),
float(kernel.Z[-1]),
),
origin="lower",
cmap="hot",
alpha=0.6,
)
plt.colorbar(im, label="Current Density")
ax.set_title("Magnetic Flux & Current")
st.pyplot(fig)
with col2:
st.metric("Magnetic Axis Flux", f"{np.max(kernel.Psi):.2f} Wb")
xp, psi_x = kernel.find_x_point(kernel.Psi)
st.metric("X-Point Location", f"R={xp[0]:.2f}, Z={xp[1]:.2f} m")
with tab2:
st.header("Thermonuclear Performance")
if st.button("Run Burn Simulation"):
physics = FusionBurnPhysics(config_path)
physics.solve_equilibrium()
metrics = physics.calculate_thermodynamics(aux_heating)
c1, c2, c3, c4 = st.columns(4)
c1.metric("Fusion Power", f"{metrics['P_fusion_MW']:.1f} MW")
c2.metric("Q-Factor", f"{metrics['Q']:.2f}")
c3.metric("Alpha Heating", f"{metrics['P_alpha_MW']:.1f} MW")
c4.metric("Status", "IGNITION" if metrics["Q"] > 10 else "Driven")
fig, ax = plt.subplots()
ax.bar(
["Aux Heat", "Alpha Heat"],
[metrics["P_aux_MW"], metrics["P_alpha_MW"]],
color=["orange", "red"],
)
ax.bar(["Losses"], [metrics["P_loss_MW"]], color="blue")
ax.set_ylabel("Power (MW)")
ax.set_title("Power Balance")
st.pyplot(fig)
with tab3:
st.header("Nuclear Wall Loading")
if st.button("Calculate Neutron Flux"):
with st.spinner("Ray-tracing neutron paths..."):
lab = NuclearEngineeringLab(config_path)
lab.solve_equilibrium()
Rw, Zw, flux = lab.calculate_neutron_wall_loading()
lifespans, load_mw = lab.analyze_materials(flux)
st.warning(f"Peak Neutron Load: {np.max(load_mw):.2f} MW/m2")
fig, ax = plt.subplots()
sc = ax.scatter(Rw, Zw, c=load_mw, cmap="inferno")
plt.colorbar(sc, label="MW/m2")
ax.set_aspect("equal")
st.pyplot(fig)
st.subheader("Component Lifespan")
st.json(lifespans)
with tab4:
st.header("Balance of Plant (Electricity Generation)")
if st.button("Calculate Grid Output"):
physics = FusionBurnPhysics(config_path)
physics.solve_equilibrium()
plasma_metrics = physics.calculate_thermodynamics(aux_heating)
plant = PowerPlantModel()
plant_metrics = plant.calculate_plant_performance(
plasma_metrics["P_fusion_MW"], aux_heating
)
c1, c2, c3 = st.columns(3)
c1.metric("Gross Electric", f"{plant_metrics['P_gross']:.1f} MWe")
c2.metric("House Load", f"{plant_metrics['P_recirc']:.1f} MWe")
c3.metric("NET TO GRID", f"{plant_metrics['P_net']:.1f} MWe", delta_color="normal")
if plant_metrics["P_net"] > 0:
st.success("✅ SYSTEM IS PRODUCING POWER!")
else:
st.error("❌ SYSTEM IS CONSUMING POWER!")
st.pyplot(plant.plot_sankey_diagram(plant_metrics))
st.subheader("Load Breakdown")
st.json(plant_metrics["breakdown"])
with tab5:
st.header("DIII-D Shot Replay & Disruption Analysis")
disruption_dir = (
data_root() / "validation" / "reference_data" / "diiid" / "disruption_shots"
)
npz_files = sorted(disruption_dir.glob("*.npz")) if disruption_dir.is_dir() else []
if not npz_files:
st.info(f"No disruption shot NPZ files found. Expected location: `{disruption_dir}`")
else:
shot_names = [path.stem for path in npz_files]
selected_shot = st.selectbox("Select Shot", shot_names, index=0)
if st.button("Load & Replay"):
from scpn_fusion.io.tokamak_archive import load_disruption_shot
with st.spinner("Loading shot data..."):
try:
shot_data = load_disruption_shot(
selected_shot, disruption_dir=disruption_dir
)
except Exception as exc: # pragma: no cover - defensive UI path
st.error(f"Failed to load shot: {exc}")
st.stop()
time_s = shot_data["time_s"]
is_disruption = shot_data["is_disruption"]
disruption_idx = shot_data["disruption_time_idx"]
disruption_type = shot_data["disruption_type"]
st.subheader("Shot Metadata")
mc1, mc2, mc3 = st.columns(3)
mc1.metric("Shot Name", selected_shot)
mc2.metric(
"Disruption Type",
disruption_type if is_disruption else "Safe (no disruption)",
)
mc3.metric("Is Disruption", "Yes" if is_disruption else "No")
disruption_time_s = (
float(time_s[disruption_idx])
if is_disruption and 0 <= disruption_idx < len(time_s)
else None
)
st.subheader("Time Series Overview")
fig1, axes = plt.subplots(2, 2, figsize=(12, 8), sharex=True)
signals_2x2 = [
("dBdt_gauss_per_s", "dB/dt (Gauss/s)", "tab:blue"),
("beta_N", r"$\beta_N$", "tab:orange"),
("Ip_MA", "Plasma Current $I_p$ (MA)", "tab:green"),
("q95", "$q_{95}$", "tab:red"),
]
for axis, (key, label, color) in zip(axes.flat, signals_2x2):
axis.plot(time_s, shot_data[key], color=color, linewidth=0.8)
axis.set_ylabel(label)
axis.grid(True, alpha=0.3)
if disruption_time_s is not None:
axis.axvline(
disruption_time_s,
color="red",
linestyle="--",
linewidth=1.2,
alpha=0.8,
label="Disruption",
)
axis.legend(loc="upper right", fontsize=8)
axes[1, 0].set_xlabel("Time (s)")
axes[1, 1].set_xlabel("Time (s)")
fig1.suptitle(f"Shot: {selected_shot}", fontsize=13)
fig1.tight_layout()
st.pyplot(fig1)
plt.close(fig1)
st.subheader("Toroidal Mode Amplitudes")
fig2, ax2 = plt.subplots(figsize=(12, 4))
ax2.plot(time_s, shot_data["n1_amp"], label="n=1 amplitude", linewidth=0.9)
ax2.plot(time_s, shot_data["n2_amp"], label="n=2 amplitude", linewidth=0.9)
ax2.plot(
time_s,
shot_data["locked_mode_amp"],
label="Locked mode amplitude",
linewidth=0.9,
)
if disruption_time_s is not None:
ax2.axvline(
disruption_time_s,
color="red",
linestyle="--",
linewidth=1.2,
alpha=0.8,
label="Disruption",
)
ax2.set_xlabel("Time (s)")
ax2.set_ylabel("Amplitude (a.u.)")
ax2.set_title("Toroidal Mode Structure")
ax2.legend(fontsize=9)
ax2.grid(True, alpha=0.3)
fig2.tight_layout()
st.pyplot(fig2)
plt.close(fig2)
st.subheader("Disruption Risk Score (Sliding Window)")
window_size = 50
n_samples = len(time_s)
risk_scores = np.full(n_samples, np.nan)
from scpn_fusion.control.disruption_predictor import predict_disruption_risk
for i in range(window_size, n_samples):
window_signal = shot_data["dBdt_gauss_per_s"][i - window_size : i]
toroidal_obs = {
"toroidal_n1_amp": float(shot_data["n1_amp"][i]),
"toroidal_n2_amp": float(shot_data["n2_amp"][i]),
"toroidal_n3_amp": float(shot_data["locked_mode_amp"][i]),
}
risk_scores[i] = predict_disruption_risk(window_signal, toroidal_obs)
fig3, ax3 = plt.subplots(figsize=(12, 4))
valid = ~np.isnan(risk_scores)
ax3.plot(time_s[valid], risk_scores[valid], color="darkred", linewidth=1.0)
ax3.axhline(
0.5,
color="orange",
linestyle=":",
linewidth=1.0,
label="Threshold (0.5)",
)
ax3.fill_between(
time_s[valid],
0,
risk_scores[valid],
where=risk_scores[valid] >= 0.5,
color="red",
alpha=0.2,
label="High Risk Region",
)
if disruption_time_s is not None:
ax3.axvline(
disruption_time_s,
color="red",
linestyle="--",
linewidth=1.2,
alpha=0.8,
label="Disruption",
)
ax3.set_xlabel("Time (s)")
ax3.set_ylabel("Risk Score")
ax3.set_ylim(-0.05, 1.05)
ax3.set_title("Disruption Predictor Risk Score vs Time")
ax3.legend(fontsize=9)
ax3.grid(True, alpha=0.3)
fig3.tight_layout()
st.pyplot(fig3)
plt.close(fig3)
st.subheader("Signal Summary Statistics")
stat_keys = [
("dBdt_gauss_per_s", "dB/dt (Gauss/s)"),
("beta_N", "beta_N"),
("Ip_MA", "Ip (MA)"),
("q95", "q95"),
("n1_amp", "n=1 Amplitude"),
("n2_amp", "n=2 Amplitude"),
("locked_mode_amp", "Locked Mode Amp"),
("ne_1e19", "ne (1e19 m^-3)"),
("vertical_position_m", "Vertical Position (m)"),
]
rows = []
for key, label in stat_keys:
arr = shot_data[key]
rows.append(
{
"Signal": label,
"Mean": f"{np.mean(arr):.4f}",
"Max": f"{np.max(arr):.4f}",
"Min": f"{np.min(arr):.4f}",
"Std": f"{np.std(arr):.4f}",
}
)
st.table(rows)
st.sidebar.markdown("---")
st.sidebar.info("SCPN Fusion Core dashboard (2026)")
if __name__ == "__main__":
main()