Energy Confinement Time From the Current API¶

This notebook keeps the analytical comparison, but uses direct fusdb objects and focused Reactor.run(...) calls throughout.

In [ ]:

from pathlib import Path
import math
import sys
import warnings

import matplotlib.pyplot as plt
import numpy as np

root = Path.cwd().resolve()
while root != root.parent and not (root / "src" / "fusdb").is_dir():
    root = root.parent

src_path = str(root / "src")
if src_path not in sys.path:
    sys.path.insert(0, src_path)

warnings.filterwarnings("ignore", message="Multiple registered relations define outputs.*")

from fusdb import Reactor, Variable
from fusdb.registry import KEV_TO_J
from fusdb.relations.geometry.plasma_geometry import (
    kappa_ipb_from_volume,
    sauter_plasma_cross_sectional_surface,
    sauter_plasma_volume,
)

reactor_path = root / "reactors" / "DEMO_2022" / "reactor.yaml"
print(f"Repository root: {root}")
print(f"Reactor path: {reactor_path}")

from pathlib import Path import math import sys import warnings import matplotlib.pyplot as plt import numpy as np root = Path.cwd().resolve() while root != root.parent and not (root / "src" / "fusdb").is_dir(): root = root.parent src_path = str(root / "src") if src_path not in sys.path: sys.path.insert(0, src_path) warnings.filterwarnings("ignore", message="Multiple registered relations define outputs.*") from fusdb import Reactor, Variable from fusdb.registry import KEV_TO_J from fusdb.relations.geometry.plasma_geometry import ( kappa_ipb_from_volume, sauter_plasma_cross_sectional_surface, sauter_plasma_volume, ) reactor_path = root / "reactors" / "DEMO_2022" / "reactor.yaml" print(f"Repository root: {root}") print(f"Reactor path: {reactor_path}")

In [ ]:

reactor_inputs = Reactor.from_yaml(reactor_path)
reactor_name = reactor_inputs.name or reactor_path.stem

R = reactor_inputs.get_variable("R").value

a = reactor_inputs.get_variable("a").value
A = reactor_inputs.get_variable("A").value or (R / a)
I_p = reactor_inputs.get_variable("I_p").value
B0 = reactor_inputs.get_variable("B0").value
H98_y2 = reactor_inputs.get_variable("H98_y2").value
afuel = reactor_inputs.get_variable("afuel").value

n_avg = reactor_inputs.get_variable("n_avg")
f_GW = reactor_inputs.get_variable("f_GW").value
if n_avg is None:
    n_GW_input = 1e20 * (I_p / 1e6) / (math.pi * a**2)
    n_avg_value = f_GW * n_GW_input
else:
    n_avg_value = n_avg.value

n_la = reactor_inputs.get_variable("n_la")
n_la_value = n_la.value if n_la is not None and n_la.value is not None else n_avg_value

T_avg_value = reactor_inputs.get_variable("T_avg").value
T_e_value = reactor_inputs.get_variable("T_e")
T_i_value = reactor_inputs.get_variable("T_i")
T_e_scalar = float(np.mean(np.asarray(T_e_value.value))) if T_e_value is not None and T_e_value.value is not None else T_avg_value
T_i_scalar = float(np.mean(np.asarray(T_i_value.value))) if T_i_value is not None and T_i_value.value is not None else T_avg_value

n_e_value = reactor_inputs.get_variable("n_e")
n_i_value = reactor_inputs.get_variable("n_i")
n_e_scalar = float(np.mean(np.asarray(n_e_value.value))) if n_e_value is not None and n_e_value.value is not None else n_avg_value
n_i_scalar = float(np.mean(np.asarray(n_i_value.value))) if n_i_value is not None and n_i_value.value is not None else n_avg_value

kappa = reactor_inputs.get_variable("kappa")
kappa_95 = reactor_inputs.get_variable("kappa_95")
delta = reactor_inputs.get_variable("delta")
delta_95 = reactor_inputs.get_variable("delta_95")
squareness = reactor_inputs.get_variable("squareness")

kappa_value = kappa.value if kappa is not None and kappa.value is not None else 1.12 * kappa_95.value
delta_value = delta.value if delta is not None and delta.value is not None else 1.5 * delta_95.value
squareness_value = squareness.value if squareness is not None and squareness.value is not None else 0.0

eps_value = a / R
S_phi_value = float(sauter_plasma_cross_sectional_surface.func(a, kappa_value, delta_value, squareness_value))
V_p_value = float(sauter_plasma_volume.func(R, delta_value, eps_value, S_phi_value))
kappa_ipb_value = float(kappa_ipb_from_volume.func(V_p_value, R, a))

p_th_value = (n_e_scalar * T_e_scalar + n_i_scalar * T_i_scalar) * KEV_TO_J
W_th_manual = 1.5 * p_th_value * V_p_value

print(f"{reactor_name}: manual inputs ready")
print(f"  R = {R:.3f} m, a = {a:.3f} m, A = {A:.3f}")
print(f"  n_avg = {n_avg_value:.3e} m^-3, n_la = {n_la_value:.3e} m^-3")
print(f"  T_e = {T_e_scalar:.3f} keV, T_i = {T_i_scalar:.3f} keV")
print(f"  V_p = {V_p_value:.3f} m^3, kappa_ipb = {kappa_ipb_value:.3f}")
print(f"  W_th = {W_th_manual/1e9:.4f} GJ")
        

reactor_inputs = Reactor.from_yaml(reactor_path) reactor_name = reactor_inputs.name or reactor_path.stem R = reactor_inputs.get_variable("R").value a = reactor_inputs.get_variable("a").value A = reactor_inputs.get_variable("A").value or (R / a) I_p = reactor_inputs.get_variable("I_p").value B0 = reactor_inputs.get_variable("B0").value H98_y2 = reactor_inputs.get_variable("H98_y2").value afuel = reactor_inputs.get_variable("afuel").value n_avg = reactor_inputs.get_variable("n_avg") f_GW = reactor_inputs.get_variable("f_GW").value if n_avg is None: n_GW_input = 1e20 * (I_p / 1e6) / (math.pi * a**2) n_avg_value = f_GW * n_GW_input else: n_avg_value = n_avg.value n_la = reactor_inputs.get_variable("n_la") n_la_value = n_la.value if n_la is not None and n_la.value is not None else n_avg_value T_avg_value = reactor_inputs.get_variable("T_avg").value T_e_value = reactor_inputs.get_variable("T_e") T_i_value = reactor_inputs.get_variable("T_i") T_e_scalar = float(np.mean(np.asarray(T_e_value.value))) if T_e_value is not None and T_e_value.value is not None else T_avg_value T_i_scalar = float(np.mean(np.asarray(T_i_value.value))) if T_i_value is not None and T_i_value.value is not None else T_avg_value n_e_value = reactor_inputs.get_variable("n_e") n_i_value = reactor_inputs.get_variable("n_i") n_e_scalar = float(np.mean(np.asarray(n_e_value.value))) if n_e_value is not None and n_e_value.value is not None else n_avg_value n_i_scalar = float(np.mean(np.asarray(n_i_value.value))) if n_i_value is not None and n_i_value.value is not None else n_avg_value kappa = reactor_inputs.get_variable("kappa") kappa_95 = reactor_inputs.get_variable("kappa_95") delta = reactor_inputs.get_variable("delta") delta_95 = reactor_inputs.get_variable("delta_95") squareness = reactor_inputs.get_variable("squareness") kappa_value = kappa.value if kappa is not None and kappa.value is not None else 1.12 * kappa_95.value delta_value = delta.value if delta is not None and delta.value is not None else 1.5 * delta_95.value squareness_value = squareness.value if squareness is not None and squareness.value is not None else 0.0 eps_value = a / R S_phi_value = float(sauter_plasma_cross_sectional_surface.func(a, kappa_value, delta_value, squareness_value)) V_p_value = float(sauter_plasma_volume.func(R, delta_value, eps_value, S_phi_value)) kappa_ipb_value = float(kappa_ipb_from_volume.func(V_p_value, R, a)) p_th_value = (n_e_scalar * T_e_scalar + n_i_scalar * T_i_scalar) * KEV_TO_J W_th_manual = 1.5 * p_th_value * V_p_value print(f"{reactor_name}: manual inputs ready") print(f" R = {R:.3f} m, a = {a:.3f} m, A = {A:.3f}") print(f" n_avg = {n_avg_value:.3e} m^-3, n_la = {n_la_value:.3e} m^-3") print(f" T_e = {T_e_scalar:.3f} keV, T_i = {T_i_scalar:.3f} keV") print(f" V_p = {V_p_value:.3f} m^3, kappa_ipb = {kappa_ipb_value:.3f}") print(f" W_th = {W_th_manual/1e9:.4f} GJ")

In [ ]:

alpha_P = -0.69

gamma = (
    H98_y2
    * 0.0562
    * (I_p / 1e6) ** 0.93
    * B0 ** 0.15
    * (n_la_value / 1e19) ** 0.41
    * R ** 1.97
    * kappa_ipb_value ** 0.78
    * A ** (-0.58)
    * afuel ** 0.19
    * (1e6) ** 0.69
)

P_loss_manual = (W_th_manual / gamma) ** (1.0 / (1.0 + alpha_P)) if gamma > 0.0 else np.inf
tau_E_manual = W_th_manual / P_loss_manual

tau_E_scaling = (
    H98_y2
    * 0.0562
    * (I_p / 1e6) ** 0.93
    * B0 ** 0.15
    * (n_la_value / 1e19) ** 0.41
    * (P_loss_manual / 1e6) ** (-0.69)
    * R ** 1.97
    * kappa_ipb_value ** 0.78
    * A ** (-0.58)
    * afuel ** 0.19
)

print("Analytical solution:")
print(f"  gamma = {gamma:.6e}")
print(f"  P_loss = {P_loss_manual/1e6:.2f} MW")
print(f"  tau_E (W_th/P_loss) = {tau_E_manual:.4f} s")
print(f"  tau_E (scaling)     = {tau_E_scaling:.4f} s")
print(f"  |difference|        = {abs(tau_E_manual - tau_E_scaling):.2e} s")
        

alpha_P = -0.69 gamma = ( H98_y2 * 0.0562 * (I_p / 1e6) ** 0.93 * B0 ** 0.15 * (n_la_value / 1e19) ** 0.41 * R ** 1.97 * kappa_ipb_value ** 0.78 * A ** (-0.58) * afuel ** 0.19 * (1e6) ** 0.69 ) P_loss_manual = (W_th_manual / gamma) ** (1.0 / (1.0 + alpha_P)) if gamma > 0.0 else np.inf tau_E_manual = W_th_manual / P_loss_manual tau_E_scaling = ( H98_y2 * 0.0562 * (I_p / 1e6) ** 0.93 * B0 ** 0.15 * (n_la_value / 1e19) ** 0.41 * (P_loss_manual / 1e6) ** (-0.69) * R ** 1.97 * kappa_ipb_value ** 0.78 * A ** (-0.58) * afuel ** 0.19 ) print("Analytical solution:") print(f" gamma = {gamma:.6e}") print(f" P_loss = {P_loss_manual/1e6:.2f} MW") print(f" tau_E (W_th/P_loss) = {tau_E_manual:.4f} s") print(f" tau_E (scaling) = {tau_E_scaling:.4f} s") print(f" |difference| = {abs(tau_E_manual - tau_E_scaling):.2e} s")

In [ ]:

reactor_solver = Reactor.from_yaml(reactor_path)
reactor_solver.relation_include = ("Energy confinement balance",)

# Fix the analytical W_th and P_loss as solver inputs, then let the
# energy-confinement balance relation reconcile tau_E = W_th / P_loss.
for name, fixed_value in (("W_th", W_th_manual), ("P_loss", P_loss_manual)):
    variable = reactor_solver.get_variable(name)
    if variable is None:
        reactor_solver.add_variable(Variable(name=name, value=fixed_value, fixed=True))
    else:
        variable.set_input(fixed_value)
        variable.fixed = True

tau_E_variable = reactor_solver.get_variable("tau_E")
if tau_E_variable is not None:
    tau_E_variable.set_input(None)
    tau_E_variable.fixed = False

solver_result = reactor_solver.run("reconcile")
P_loss_run = solver_result["variables"]["P_loss"].value
tau_E_run = solver_result["variables"]["tau_E"].value
balance_status = solver_result["relation_status"].get("Energy confinement balance", {})

print(f"Energy confinement balance satisfied: {balance_status.get('verified', False)}")
print(f"  P_loss = {P_loss_run/1e6:.2f} MW")
print(f"  W_th   = {W_th_manual/1e9:.4f} GJ (manual seed)")
print(f"  tau_E  = {tau_E_run:.4f} s")

reactor_solver = Reactor.from_yaml(reactor_path) reactor_solver.relation_include = ("Energy confinement balance",) # Fix the analytical W_th and P_loss as solver inputs, then let the # energy-confinement balance relation reconcile tau_E = W_th / P_loss. for name, fixed_value in (("W_th", W_th_manual), ("P_loss", P_loss_manual)): variable = reactor_solver.get_variable(name) if variable is None: reactor_solver.add_variable(Variable(name=name, value=fixed_value, fixed=True)) else: variable.set_input(fixed_value) variable.fixed = True tau_E_variable = reactor_solver.get_variable("tau_E") if tau_E_variable is not None: tau_E_variable.set_input(None) tau_E_variable.fixed = False solver_result = reactor_solver.run("reconcile") P_loss_run = solver_result["variables"]["P_loss"].value tau_E_run = solver_result["variables"]["tau_E"].value balance_status = solver_result["relation_status"].get("Energy confinement balance", {}) print(f"Energy confinement balance satisfied: {balance_status.get('verified', False)}") print(f" P_loss = {P_loss_run/1e6:.2f} MW") print(f" W_th = {W_th_manual/1e9:.4f} GJ (manual seed)") print(f" tau_E = {tau_E_run:.4f} s")

In [ ]:

print("Comparison: analytical vs focused Reactor.run()")
print(f"  tau_E difference = {tau_E_manual - tau_E_run:.4e} s")

print("Comparison: analytical vs focused Reactor.run()") print(f" tau_E difference = {tau_E_manual - tau_E_run:.4e} s")

In [ ]:

P_loss_range = np.linspace(0.5 * P_loss_manual, 1.5 * P_loss_manual, 200)
tau_E_energy = W_th_manual / P_loss_range
tau_E_scaling_range = (
    H98_y2
    * 0.0562
    * (I_p / 1e6) ** 0.93
    * B0 ** 0.15
    * (n_la_value / 1e19) ** 0.41
    * (P_loss_range / 1e6) ** (-0.69)
    * R ** 1.97
    * kappa_ipb_value ** 0.78
    * A ** (-0.58)
    * afuel ** 0.19
)

fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(P_loss_range / 1e6, tau_E_energy, linewidth=2, label="Energy balance")
ax.plot(P_loss_range / 1e6, tau_E_scaling_range, linewidth=2, label="IPB98(y,2) scaling")
ax.plot(P_loss_manual / 1e6, tau_E_manual, "o", markersize=8, label="Analytical intersection")
ax.plot(P_loss_run / 1e6, tau_E_run, "x", markersize=10, markeredgewidth=2, label="Focused Reactor.run()")
ax.set_xlabel("P_loss [MW]")
ax.set_ylabel("tau_E [s]")
ax.set_title(f"{reactor_name}: confinement-time consistency")
ax.grid(True, alpha=0.3)
ax.legend()
plt.show()
        

P_loss_range = np.linspace(0.5 * P_loss_manual, 1.5 * P_loss_manual, 200) tau_E_energy = W_th_manual / P_loss_range tau_E_scaling_range = ( H98_y2 * 0.0562 * (I_p / 1e6) ** 0.93 * B0 ** 0.15 * (n_la_value / 1e19) ** 0.41 * (P_loss_range / 1e6) ** (-0.69) * R ** 1.97 * kappa_ipb_value ** 0.78 * A ** (-0.58) * afuel ** 0.19 ) fig, ax = plt.subplots(figsize=(10, 6)) ax.plot(P_loss_range / 1e6, tau_E_energy, linewidth=2, label="Energy balance") ax.plot(P_loss_range / 1e6, tau_E_scaling_range, linewidth=2, label="IPB98(y,2) scaling") ax.plot(P_loss_manual / 1e6, tau_E_manual, "o", markersize=8, label="Analytical intersection") ax.plot(P_loss_run / 1e6, tau_E_run, "x", markersize=10, markeredgewidth=2, label="Focused Reactor.run()") ax.set_xlabel("P_loss [MW]") ax.set_ylabel("tau_E [s]") ax.set_title(f"{reactor_name}: confinement-time consistency") ax.grid(True, alpha=0.3) ax.legend() plt.show()