1#!/usr/bin/env python3
2"""
3Tokamak Stability Analysis
4
5Analysis of MHD stability in a tokamak plasma confinement device.
6Computes safety factor profiles, checks Kruskal-Shafranov criterion,
7and estimates tearing mode stability.
8
9Key results:
10- Safety factor q(r) profile for different current distributions
11- Kruskal-Shafranov limit: q(a) > 2-3 for stability
12- Tearing mode stability parameter Δ'
13- Rational surfaces where resonant modes can develop
14
15Physics:
16- Safety factor: q = (r B_φ) / (R₀ B_θ) ≈ (r/R₀)(B_φ/B_θ)
17- Kruskal-Shafranov: Prevents external kink modes
18- Tearing modes: m,n modes at rational surfaces q = m/n
19
20Author: Claude
21Date: 2026-02-14
22"""
23
24import numpy as np
25import matplotlib.pyplot as plt
26from scipy.integrate import odeint, cumtrapz
27from scipy.interpolate import interp1d
28from typing import Tuple, Callable
29
30
31class TokamakEquilibrium:
32 """
33 Tokamak equilibrium and stability analysis.
34
35 Simple cylindrical model with:
36 - Major radius R₀
37 - Minor radius a
38 - Toroidal field B_φ(r)
39 - Poloidal field B_θ(r) from current profile j(r)
40 """
41
42 def __init__(self, R0: float = 1.0, a: float = 0.3,
43 B0: float = 2.0, I_p: float = 1.0e6):
44 """
45 Initialize tokamak parameters.
46
47 Args:
48 R0: Major radius (m)
49 a: Minor radius (m)
50 B0: Toroidal field at R0 (T)
51 I_p: Total plasma current (A)
52 """
53 self.R0 = R0
54 self.a = a
55 self.B0 = B0
56 self.I_p = I_p
57
58 # Constants
59 self.mu0 = 4 * np.pi * 1e-7 # Permeability of free space
60
61 def safety_factor(self, r: np.ndarray, j_profile: Callable) -> np.ndarray:
62 """
63 Compute safety factor q(r).
64
65 In cylindrical approximation:
66 q(r) ≈ (r * B_φ) / (R₀ * B_θ(r))
67
68 where B_θ(r) is computed from current profile via Ampere's law:
69 B_θ(r) = (μ₀/2πr) ∫₀ʳ j(r') 2πr' dr' = (μ₀/r) ∫₀ʳ j(r') r' dr'
70
71 Args:
72 r: Radial coordinate array
73 j_profile: Current density function j(r)
74
75 Returns:
76 Safety factor q(r)
77 """
78 # Compute current density
79 j = j_profile(r)
80
81 # Integrate to get enclosed current I(r)
82 # I(r) = ∫₀ʳ j(r') 2πr' dr'
83 integrand = j * r
84 I_enclosed = 2 * np.pi * cumtrapz(integrand, r, initial=0)
85
86 # Poloidal magnetic field from Ampere's law
87 # B_θ(r) = μ₀ I(r) / (2π r)
88 B_theta = np.zeros_like(r)
89 B_theta[1:] = self.mu0 * I_enclosed[1:] / (2 * np.pi * r[1:])
90 B_theta[0] = 0 # On axis
91
92 # Toroidal field (approximately constant in simple model)
93 B_phi = self.B0 * np.ones_like(r)
94
95 # Safety factor
96 q = np.zeros_like(r)
97 q[1:] = (r[1:] * B_phi[1:]) / (self.R0 * B_theta[1:])
98
99 # On-axis limit: q(0) ≈ 2 B₀ / (μ₀ R₀ j₀)
100 if j[0] > 0:
101 q[0] = 2 * B_phi[0] / (self.mu0 * self.R0 * j[0])
102 else:
103 q[0] = q[1]
104
105 return q
106
107 def kruskal_shafranov_criterion(self, q_a: float) -> Tuple[bool, str]:
108 """
109 Check Kruskal-Shafranov criterion for external kink stability.
110
111 Criterion: q(a) > q_crit ≈ 2-3
112
113 More precisely:
114 - q(a) > 2/(m-1) for m-th mode stability
115 - For m=2 (most dangerous): q(a) > 2
116 - Empirically, q(a) > 3 provides good stability margin
117
118 Args:
119 q_a: Safety factor at plasma edge
120
121 Returns:
122 (stable, message)
123 """
124 if q_a > 3.0:
125 return True, f"Stable (q(a)={q_a:.2f} > 3)"
126 elif q_a > 2.0:
127 return True, f"Marginally stable (q(a)={q_a:.2f}, 2 < q < 3)"
128 else:
129 return False, f"Unstable to kink (q(a)={q_a:.2f} < 2)"
130
131 def find_rational_surfaces(self, r: np.ndarray, q: np.ndarray,
132 m_max: int = 5, n: int = 1) -> list:
133 """
134 Find rational surfaces where q = m/n.
135
136 These are locations where resonant MHD modes can develop.
137
138 Args:
139 r: Radial coordinate
140 q: Safety factor profile
141 m_max: Maximum poloidal mode number
142 n: Toroidal mode number (typically n=1)
143
144 Returns:
145 List of (m, n, r_res) tuples
146 """
147 rational_surfaces = []
148
149 # Interpolate q(r) for finding roots
150 q_interp = interp1d(r, q, kind='cubic', bounds_error=False, fill_value='extrapolate')
151
152 for m in range(2, m_max + 1):
153 q_res = m / n
154
155 # Check if resonance exists in domain
156 if q.min() < q_res < q.max():
157 # Find radius where q(r) = m/n
158 # Binary search
159 r_search = r[(q > q_res * 0.9) & (q < q_res * 1.1)]
160 if len(r_search) > 0:
161 # Refine using interpolation
162 r_fine = np.linspace(r_search.min(), r_search.max(), 1000)
163 q_fine = q_interp(r_fine)
164 idx = np.argmin(np.abs(q_fine - q_res))
165 r_res = r_fine[idx]
166
167 rational_surfaces.append((m, n, r_res))
168
169 return rational_surfaces
170
171 def tearing_mode_delta_prime(self, r: np.ndarray, q: np.ndarray,
172 m: int, n: int) -> float:
173 """
174 Estimate tearing mode stability parameter Δ'.
175
176 Δ' measures the jump in dψ/dr across the rational surface.
177 - Δ' > 0: Unstable (tearing mode grows)
178 - Δ' < 0: Stable (current profile resists tearing)
179
180 Simplified estimate using current gradient:
181 Δ' ≈ (2μ₀ R₀ / B_φ) * (dj/dr) / q' at q = m/n
182
183 This is a crude approximation; real calculation requires
184 solving the Grad-Shafranov equation and matching inner/outer solutions.
185
186 Args:
187 r: Radial coordinate
188 q: Safety factor
189 m, n: Mode numbers
190
191 Returns:
192 Δ' estimate (normalized)
193 """
194 # Find rational surface
195 q_res = m / n
196 rational_surfaces = self.find_rational_surfaces(r, q, m_max=m, n=n)
197
198 if not any(surf[0] == m for surf in rational_surfaces):
199 return 0.0 # No resonance
200
201 r_res = next(surf[2] for surf in rational_surfaces if surf[0] == m)
202
203 # Compute q' = dq/dr at resonance
204 dq_dr = np.gradient(q, r)
205 q_prime_interp = interp1d(r, dq_dr, kind='linear')
206 q_prime_res = q_prime_interp(r_res)
207
208 # Simplified Δ' estimate
209 # Positive if current profile is peaked (dj/dr < 0 at r_res)
210 # This is a very rough proxy
211 delta_prime = -1.0 / (r_res * q_prime_res + 1e-10) # Normalized
212
213 return delta_prime
214
215
216# =============================================================================
217# Current Profile Models
218# =============================================================================
219
220def parabolic_current(alpha: float = 2.0):
221 """
222 Parabolic current profile: j(r) = j₀ (1 - (r/a)^α)
223
224 Args:
225 alpha: Peaking parameter (α=1: linear, α=2: parabolic)
226
227 Returns:
228 Function j(r) normalized to integrate to 1
229 """
230 def j_func(r: np.ndarray, a: float = 0.3) -> np.ndarray:
231 rho = r / a
232 j = (1 - rho**alpha)
233 j[r > a] = 0
234 # Normalize
235 norm = np.trapz(j * r, r) * 2 * np.pi
236 return j / (norm + 1e-10)
237
238 return j_func
239
240
241def hollow_current(r_peak: float = 0.5, width: float = 0.1):
242 """
243 Hollow current profile (off-axis peak).
244
245 j(r) = j₀ exp(-(r - r_peak)² / (2*width²))
246
247 Hollow current profiles can be unstable to tearing modes.
248
249 Returns:
250 Function j(r)
251 """
252 def j_func(r: np.ndarray, a: float = 0.3) -> np.ndarray:
253 j = np.exp(-(r - r_peak * a)**2 / (2 * (width * a)**2))
254 j[r > a] = 0
255 # Normalize
256 norm = np.trapz(j * r, r) * 2 * np.pi
257 return j / (norm + 1e-10)
258
259 return j_func
260
261
262def bootstrap_current(alpha: float = 1.5, beta: float = 2.0):
263 """
264 Bootstrap current profile (typical in advanced tokamaks).
265
266 Peaked off-axis due to pressure gradient.
267
268 j(r) = j₀ (r/a)^α (1 - (r/a)^β)
269
270 Returns:
271 Function j(r)
272 """
273 def j_func(r: np.ndarray, a: float = 0.3) -> np.ndarray:
274 rho = r / a
275 j = rho**alpha * (1 - rho**beta)
276 j[r > a] = 0
277 j[r == 0] = 0
278 # Normalize
279 norm = np.trapz(j * r, r) * 2 * np.pi
280 return j / (norm + 1e-10)
281
282 return j_func
283
284
285# =============================================================================
286# Analysis and Visualization
287# =============================================================================
288
289def compare_current_profiles():
290 """Compare different current profiles and their q-profiles."""
291 print("=" * 70)
292 print("Tokamak Current Profiles and Safety Factor Analysis")
293 print("=" * 70)
294
295 # Setup
296 tokamak = TokamakEquilibrium(R0=1.0, a=0.3, B0=2.0, I_p=1.0e6)
297 r = np.linspace(0, tokamak.a, 200)
298
299 # Current profiles to compare
300 profiles = {
301 'Parabolic (α=2)': parabolic_current(alpha=2.0),
302 'Parabolic (α=1)': parabolic_current(alpha=1.0),
303 'Hollow': hollow_current(r_peak=0.6, width=0.15),
304 'Bootstrap': bootstrap_current(alpha=1.5, beta=2.0)
305 }
306
307 fig, axes = plt.subplots(2, 2, figsize=(14, 10))
308
309 # Results storage
310 results = {}
311
312 for idx, (name, j_func) in enumerate(profiles.items()):
313 print(f"\n{name}:")
314 print("-" * 40)
315
316 # Compute profiles
317 j = j_func(r, tokamak.a)
318 q = tokamak.safety_factor(r, lambda x: j_func(x, tokamak.a))
319
320 # Store results
321 results[name] = {'r': r, 'j': j, 'q': q}
322
323 # Safety factor at edge
324 q_a = q[-1]
325 print(f" q(0) = {q[0]:.2f}")
326 print(f" q(a) = {q_a:.2f}")
327
328 # Kruskal-Shafranov
329 stable, msg = tokamak.kruskal_shafranov_criterion(q_a)
330 print(f" Kruskal-Shafranov: {msg}")
331
332 # Rational surfaces
333 rational = tokamak.find_rational_surfaces(r, q, m_max=5, n=1)
334 print(f" Rational surfaces q=m/n:")
335 for m, n, r_res in rational:
336 delta_prime = tokamak.tearing_mode_delta_prime(r, q, m, n)
337 stability = "unstable" if delta_prime > 0 else "stable"
338 print(f" m/n={m}/{n} at r={r_res:.3f} m (r/a={r_res/tokamak.a:.2f}), "
339 f"Δ'={delta_prime:.2f} ({stability})")
340
341 # Current density profile
342 ax = axes[0, 0]
343 ax.plot(r / tokamak.a, j / j.max(), label=name, linewidth=2)
344
345 # Safety factor profile
346 ax = axes[0, 1]
347 ax.plot(r / tokamak.a, q, label=name, linewidth=2)
348
349 # Finalize current density plot
350 ax = axes[0, 0]
351 ax.set_xlabel('r/a (normalized radius)')
352 ax.set_ylabel('j(r) / j_max (normalized)')
353 ax.set_title('Current Density Profiles')
354 ax.legend()
355 ax.grid(True, alpha=0.3)
356
357 # Finalize safety factor plot
358 ax = axes[0, 1]
359 ax.set_xlabel('r/a')
360 ax.set_ylabel('q(r)')
361 ax.set_title('Safety Factor Profiles')
362 ax.axhline(y=2, color='r', linestyle='--', label='q=2 (KS limit)')
363 ax.axhline(y=3, color='orange', linestyle='--', label='q=3')
364 # Mark rational surfaces
365 for m in [2, 3, 4]:
366 ax.axhline(y=m, color='gray', linestyle=':', alpha=0.5, linewidth=0.8)
367 ax.text(0.02, m, f'q={m}', fontsize=8, color='gray')
368 ax.set_ylim([0, 8])
369 ax.legend()
370 ax.grid(True, alpha=0.3)
371
372 # Detailed view: Parabolic vs Hollow
373 ax = axes[1, 0]
374 for name in ['Parabolic (α=2)', 'Hollow']:
375 r_plot = results[name]['r']
376 q_plot = results[name]['q']
377 ax.plot(r_plot / tokamak.a, q_plot, label=name, linewidth=2.5)
378
379 # Mark rational surfaces
380 rational = tokamak.find_rational_surfaces(r_plot, q_plot, m_max=5, n=1)
381 for m, n, r_res in rational:
382 ax.plot(r_res / tokamak.a, m/n, 'o', markersize=8)
383 ax.text(r_res / tokamak.a + 0.02, m/n, f'{m}/{n}', fontsize=9)
384
385 ax.set_xlabel('r/a')
386 ax.set_ylabel('q(r)')
387 ax.set_title('Rational Surfaces (m/n resonances)')
388 ax.axhline(y=2, color='r', linestyle='--', alpha=0.5)
389 ax.set_ylim([1, 6])
390 ax.legend()
391 ax.grid(True, alpha=0.3)
392
393 # Shear profile: s = (r/q) dq/dr
394 ax = axes[1, 1]
395 for name in ['Parabolic (α=2)', 'Bootstrap']:
396 r_plot = results[name]['r']
397 q_plot = results[name]['q']
398
399 # Magnetic shear
400 dq_dr = np.gradient(q_plot, r_plot)
401 shear = (r_plot / (q_plot + 1e-10)) * dq_dr
402
403 ax.plot(r_plot / tokamak.a, shear, label=name, linewidth=2)
404
405 ax.set_xlabel('r/a')
406 ax.set_ylabel('Magnetic shear s')
407 ax.set_title('Magnetic Shear Profile')
408 ax.axhline(y=0, color='k', linestyle='-', alpha=0.3)
409 ax.legend()
410 ax.grid(True, alpha=0.3)
411
412 plt.tight_layout()
413 plt.savefig('/opt/projects/01_Personal/03_Study/examples/Numerical_Simulation/10_projects/tokamak_stability.png',
414 dpi=150, bbox_inches='tight')
415 plt.close()
416 print("\nFigure saved: tokamak_stability.png")
417
418
419def stability_diagram():
420 """
421 Create stability diagram: q(a) vs q(0) with stability boundaries.
422 """
423 print("\n" + "=" * 70)
424 print("Tokamak Stability Diagram")
425 print("=" * 70)
426
427 tokamak = TokamakEquilibrium(R0=1.0, a=0.3, B0=2.0, I_p=1.0e6)
428 r = np.linspace(0, tokamak.a, 200)
429
430 # Scan over different peaking parameters
431 alphas = np.linspace(0.5, 4.0, 30)
432 q0_values = []
433 qa_values = []
434 stable_flags = []
435
436 for alpha in alphas:
437 j_func = parabolic_current(alpha=alpha)
438 q = tokamak.safety_factor(r, lambda x: j_func(x, tokamak.a))
439
440 q0_values.append(q[0])
441 qa_values.append(q[-1])
442
443 stable, _ = tokamak.kruskal_shafranov_criterion(q[-1])
444 stable_flags.append(stable)
445
446 q0_values = np.array(q0_values)
447 qa_values = np.array(qa_values)
448 stable_flags = np.array(stable_flags)
449
450 # Plot
451 fig, ax = plt.subplots(figsize=(10, 8))
452
453 # Color by stability
454 stable_idx = stable_flags == True
455 unstable_idx = stable_flags == False
456
457 ax.scatter(q0_values[stable_idx], qa_values[stable_idx],
458 c='green', s=50, alpha=0.7, label='Stable', zorder=3)
459 ax.scatter(q0_values[unstable_idx], qa_values[unstable_idx],
460 c='red', s=50, alpha=0.7, label='Unstable (kink)', zorder=3)
461
462 # Stability boundaries
463 ax.axhline(y=2, color='red', linestyle='--', linewidth=2, label='q(a)=2 (KS limit)')
464 ax.axhline(y=3, color='orange', linestyle='--', linewidth=2, label='q(a)=3')
465
466 # Optimal region
467 ax.fill_between([0, 10], 2, 3, alpha=0.1, color='orange', label='Marginal')
468 ax.fill_between([0, 10], 3, 10, alpha=0.1, color='green', label='Safe region')
469
470 ax.set_xlabel('q(0) (on-axis safety factor)', fontsize=12)
471 ax.set_ylabel('q(a) (edge safety factor)', fontsize=12)
472 ax.set_title('Tokamak Stability Diagram (Kruskal-Shafranov)', fontsize=14)
473 ax.set_xlim([0.5, 3])
474 ax.set_ylim([1.5, 6])
475 ax.legend(loc='upper left')
476 ax.grid(True, alpha=0.3)
477
478 plt.tight_layout()
479 plt.savefig('/opt/projects/01_Personal/03_Study/examples/Numerical_Simulation/10_projects/tokamak_stability_diagram.png',
480 dpi=150, bbox_inches='tight')
481 plt.close()
482 print("Figure saved: tokamak_stability_diagram.png")
483
484 print("\nKey findings:")
485 print(f" - Stable configurations: q(a) > 2-3")
486 print(f" - Typical range: q(0) ≈ 1-2, q(a) ≈ 3-5")
487 print(f" - Lower q(0) (peaked current) → higher q(a) needed for stability")
488
489
490def main():
491 """Run complete tokamak stability analysis."""
492 compare_current_profiles()
493 stability_diagram()
494
495 print("\n" + "=" * 70)
496 print("Summary: Tokamak MHD Stability")
497 print("=" * 70)
498 print("""
499 Key Stability Criteria:
500
501 1. Kruskal-Shafranov Criterion (External Kink):
502 - q(a) > 2: Stable to m=2 kink mode
503 - q(a) > 3: Safe operating margin
504 - Prevents large-scale disruptions
505
506 2. Rational Surfaces (Tearing Modes):
507 - Occur where q = m/n (m,n integers)
508 - Most dangerous: m/n = 2/1, 3/2
509 - Δ' > 0 → unstable (islands grow)
510 - Hollow current profiles more susceptible
511
512 3. Current Profile Effects:
513 - Parabolic (broad): High q(a), good stability
514 - Peaked: Low q(a), may violate KS criterion
515 - Hollow: Tearing mode instabilities
516 - Bootstrap: Advanced scenarios, needs careful control
517
518 4. Magnetic Shear s = (r/q) dq/dr:
519 - Positive shear (standard): Generally stable
520 - Negative shear (reversed): Can stabilize turbulence
521 - Low/zero shear: Enhanced confinement, but instabilities
522
523 Operational Implications:
524 - ITER design: q(a) ≈ 3, q(0) ≈ 1 (safety margin)
525 - Real-time control needed to maintain q-profile
526 - Advanced scenarios explore low-shear, high-β regimes
527 - Disruption avoidance critical for machine protection
528
529 Limitations of This Analysis:
530 - Cylindrical approximation (real tokamaks are toroidal)
531 - No finite pressure effects (β)
532 - No kinetic effects
533 - Simplified Δ' calculation
534 """)
535
536
537if __name__ == "__main__":
538 main()