SCPN Fusion Core — Compact Reactor Design Search

This notebook demonstrates the Minimum Viable Reactor (MVR-0.96) optimizer, which performs multi-objective design-space exploration to find the smallest tokamak configuration achieving ignition.

What you'll learn:

1. How the plasma physics model computes fusion power from geometry and magnetic field
2. How radial build constraints enforce HTS magnet limits
3. How the optimizer scans design space and finds the minimum viable reactor
4. How to visualize the design space and interpret results

---

import numpy as np
import matplotlib.pyplot as plt
import time

# Add the repo root to path if running from examples/
import sys, os
sys.path.insert(0, os.path.join(os.path.dirname(os.getcwd()), 'src'))

from scpn_fusion.core.compact_reactor_optimizer import CompactReactorArchitect

1. Initialize the Optimizer

The CompactReactorArchitect encodes:

• HTS magnet limits: J_crit = 1500 MA/m², B_max = 30 T
• Shielding: λ_shield = 10 cm
• Neutron fluence: 5×10²² n/m² lifetime limit

architect = CompactReactorArchitect()
print(f"Critical current density: {architect.J_crit_base} MA/m²")
print(f"Max coil field: {architect.B_max_coil} T")
print(f"Shield thickness: {architect.lambda_shield} m")

2. Explore the Physics Model

The core physics model computes fusion power from three parameters:

• R — Major radius (m)
• a — Minor radius (m)
• B0 — Toroidal field at plasma center (T)

Let's sweep R and B0 to see how fusion power scales.

# Sweep major radius and field
radii = np.linspace(0.5, 4.0, 50)
fields = np.linspace(5.0, 20.0, 50)
RR, BB = np.meshgrid(radii, fields)
PP = np.zeros_like(RR)

A = 2.5  # aspect ratio
for i in range(len(fields)):
    for j in range(len(radii)):
        a = radii[j] / A
        P_fus, _, _ = architect.plasma_physics_model(radii[j], a, fields[i])
        PP[i, j] = P_fus

fig, ax = plt.subplots(figsize=(10, 7))
c = ax.contourf(RR, BB, np.log10(PP + 1e-3), levels=20, cmap='plasma')
plt.colorbar(c, label='log₁₀(P_fusion / MW)')
ax.set_xlabel('Major Radius R (m)', fontsize=12)
ax.set_ylabel('Toroidal Field B₀ (T)', fontsize=12)
ax.set_title('Fusion Power Scaling (Aspect Ratio = 2.5)', fontsize=14)
plt.tight_layout()
plt.show()

3. Run the MVR Optimizer

The optimizer scans (R, B0, A) space with engineering constraints:

• Radial build must accommodate HTS coils within J_crit limits
• Divertor heat flux < 100 MW/m² (TEMHD liquid divertor)
• First-wall neutron flux < 5 MW/m²

Target: find the smallest R that delivers ≥ 5 MW fusion power.

t0 = time.perf_counter()
architect.find_minimum_reactor(target_power_MW=5.0, use_temhd=True)
dt = time.perf_counter() - t0
print(f"\nOptimization completed in {dt:.3f} s")

4. Compare TEMHD vs Solid Divertor

The TEMHD (thermoelectric MHD) liquid metal divertor allows 10× higher heat flux limits, enabling significantly more compact designs.

print("=" * 60)
print("COMPARISON: TEMHD vs Solid Divertor")
print("=" * 60)

print("\n--- WITH TEMHD liquid divertor (q_div < 100 MW/m²) ---")
architect.find_minimum_reactor(target_power_MW=5.0, use_temhd=True)

print("\n--- WITH solid divertor (q_div < 10 MW/m²) ---")
architect.find_minimum_reactor(target_power_MW=5.0, use_temhd=False)

5. Parameter Sensitivity

How does the minimum viable radius change with target fusion power?

# Simplified sweep: find minimum R for each target power
targets = [1.0, 2.0, 5.0, 10.0, 50.0, 100.0, 500.0]
min_radii_temhd = []
min_radii_solid = []

for target in targets:
    best_R_t, best_R_s = 999.0, 999.0
    for R in np.linspace(0.3, 5.0, 100):
        for B0 in np.linspace(5.0, 20.0, 30):
            for A in [2.0, 2.5, 3.0]:
                a = R / A
                P_fus, _, _ = architect.plasma_physics_model(R, a, B0)
                if P_fus < target:
                    continue
                ok, B_coil = architect.radial_build_constraints(R, a, B0)
                if not ok:
                    continue
                f_rad = 0.90
                P_sep = (0.2 * P_fus + 5.0) * (1.0 - f_rad)
                lambda_q = 0.63 * (B0**(-1.19)) * 1e-3
                Area_div = 2 * np.pi * R * lambda_q * 20.0
                q_div = P_sep / Area_div
                Area_wall = 4 * np.pi**2 * R * a
                q_wall = (0.8 * P_fus) / Area_wall
                if q_wall < 5.0:
                    if q_div < 100.0 and R < best_R_t:
                        best_R_t = R
                    if q_div < 10.0 and R < best_R_s:
                        best_R_s = R
    min_radii_temhd.append(best_R_t if best_R_t < 999 else None)
    min_radii_solid.append(best_R_s if best_R_s < 999 else None)

fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(targets, min_radii_temhd, 'bo-', label='TEMHD divertor', linewidth=2, markersize=8)
ax.plot(targets, min_radii_solid, 'rs--', label='Solid divertor', linewidth=2, markersize=8)
ax.set_xscale('log')
ax.set_xlabel('Target Fusion Power (MW)', fontsize=12)
ax.set_ylabel('Minimum Major Radius R (m)', fontsize=12)
ax.set_title('Compact Reactor Scaling: TEMHD vs Solid Divertor', fontsize=14)
ax.legend(fontsize=12)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

Performance Benchmarks

Timing the key computations in this notebook:

1. Single physics model evaluation (plasma_physics_model)
2. Full MVR optimizer sweep (find_minimum_reactor)
3. Radial build constraint check (radial_build_constraints)

import timeit

# 1. Single plasma physics model evaluation
def bench_physics_model():
    architect.plasma_physics_model(1.85, 0.74, 12.0)

t_physics = timeit.repeat(bench_physics_model, number=1000, repeat=5)
print(f"plasma_physics_model (1000 calls):")
print(f"  Mean: {np.mean(t_physics)*1000:.1f} ms +/- {np.std(t_physics)*1000:.1f} ms")
print(f"  Per call: {np.mean(t_physics):.3f} us")

# 2. Single radial build constraint check
def bench_radial_build():
    architect.radial_build_constraints(1.85, 0.74, 12.0)

t_radial = timeit.repeat(bench_radial_build, number=1000, repeat=5)
print(f"\nradial_build_constraints (1000 calls):")
print(f"  Mean: {np.mean(t_radial)*1000:.1f} ms +/- {np.std(t_radial)*1000:.1f} ms")
print(f"  Per call: {np.mean(t_radial):.3f} us")

# 3. Full MVR optimizer sweep (find_minimum_reactor)
t_mvr = timeit.repeat(
    lambda: architect.find_minimum_reactor(target_power_MW=5.0, use_temhd=True),
    number=1, repeat=3
)
print(f"\nfind_minimum_reactor (full sweep):")
print(f"  Mean: {np.mean(t_mvr)*1000:.1f} ms +/- {np.std(t_mvr)*1000:.1f} ms")

# 4. 50x50 parameter space scan (R, B0 grid from Section 2)
def bench_param_scan():
    A = 2.5
    for i in range(50):
        for j in range(50):
            a = radii[j] / A
            architect.plasma_physics_model(radii[j], a, fields[i])

t_scan = timeit.repeat(bench_param_scan, number=1, repeat=3)
print(f"\n50x50 parameter scan (2500 evaluations):")
print(f"  Mean: {np.mean(t_scan)*1000:.1f} ms +/- {np.std(t_scan)*1000:.1f} ms")

Summary

Key takeaways from the compact reactor search:

1. TEMHD divertors enable ~2× smaller reactors by relaxing heat flux limits
2. HTS magnets (REBCO, 30 T) are the key enabling technology for compact designs
3. Fusion power scales roughly as R² × B⁴, making high-field compact designs attractive
4. The MVR-0.96 achieves ignition-class performance at R ~ 1 m scale

For deeper analysis, see:

• docs/COMPACT_REACTOR_FINDINGS.md
• SCPN_FUSION_CORE_COMPREHENSIVE_STUDY.md (Section: Compact Reactor Optimization)