1#!/usr/bin/env python3
2"""
3Safety Factor (q) Profile Stability Analysis
4
5This module analyzes the stability of tokamak q-profiles with respect to
6various MHD instabilities.
7
8Stability criteria checked:
91. q(0) > 1: Sawtooth stability (internal kink mode)
102. q(a) > 2: External kink mode stability
113. Suydam criterion: Local interchange stability
124. Rational surfaces q = m/n: Locations of resonant modes
135. Magnetic shear: s = (r/q)(dq/dr)
14
15Author: MHD Course Examples
16License: MIT
17"""
18
19import numpy as np
20import matplotlib.pyplot as plt
21from scipy.interpolate import interp1d
22
23
24class QProfileStability:
25 """
26 Safety factor profile stability analyzer.
27
28 Attributes:
29 r (ndarray): Minor radius grid
30 a (float): Minor radius
31 q (ndarray): Safety factor profile
32 """
33
34 def __init__(self, r, q):
35 """
36 Initialize with q-profile.
37
38 Parameters:
39 r (ndarray): Minor radius array
40 q (ndarray): Safety factor array
41 """
42 self.r = r
43 self.a = r[-1]
44 self.q = q
45
46 # Create interpolator
47 self.q_interp = interp1d(r, q, kind='cubic', fill_value='extrapolate')
48
49 def check_sawtooth_stability(self):
50 """
51 Check sawtooth stability: q(0) > 1
52
53 Returns:
54 dict: Stability info
55 """
56 q0 = self.q[0]
57 stable = q0 > 1.0
58
59 return {
60 'stable': stable,
61 'q0': q0,
62 'criterion': 'q(0) > 1',
63 'message': f'q(0) = {q0:.3f}' + (' ✓' if stable else ' ✗ UNSTABLE')
64 }
65
66 def check_kink_stability(self):
67 """
68 Check external kink stability: q(a) > 2
69
70 Returns:
71 dict: Stability info
72 """
73 qa = self.q[-1]
74 stable = qa > 2.0
75
76 return {
77 'stable': stable,
78 'qa': qa,
79 'criterion': 'q(a) > 2',
80 'message': f'q(a) = {qa:.3f}' + (' ✓' if stable else ' ✗ UNSTABLE')
81 }
82
83 def compute_shear(self):
84 """
85 Compute magnetic shear s = (r/q)(dq/dr).
86
87 Returns:
88 ndarray: Shear profile
89 """
90 # Compute dq/dr using centered differences
91 dq_dr = np.gradient(self.q, self.r)
92
93 # Shear
94 shear = (self.r / self.q) * dq_dr
95
96 # Handle r=0
97 shear[0] = shear[1]
98
99 return shear
100
101 def suydam_criterion(self, pressure_gradient=None):
102 """
103 Check Suydam criterion for local interchange stability.
104
105 Criterion: (r q'²) / (4 q²) + 2μ₀ p'/B² > 0
106
107 Simplified version without pressure: q' > 0 (positive shear)
108
109 Parameters:
110 pressure_gradient (ndarray): dp/dr if available
111
112 Returns:
113 dict: Stability info
114 """
115 dq_dr = np.gradient(self.q, self.r)
116
117 # Simplified Suydam (no pressure)
118 suydam_parameter = dq_dr
119
120 if pressure_gradient is not None:
121 # Full Suydam with pressure (simplified)
122 # Assuming cylindrical approximation
123 term1 = (self.r * dq_dr**2) / (4 * self.q**2)
124 term2 = pressure_gradient # Simplified pressure term
125 suydam_parameter = term1 + term2
126
127 # Check if positive everywhere
128 stable = np.all(suydam_parameter > 0)
129
130 # Find unstable regions
131 unstable_regions = self.r[suydam_parameter < 0]
132
133 return {
134 'stable': stable,
135 'criterion': "Suydam: r q'²/(4q²) + 2μ₀p'/B² > 0",
136 'parameter': suydam_parameter,
137 'unstable_r': unstable_regions,
138 'message': 'Suydam criterion satisfied' if stable else
139 f'Suydam unstable at {len(unstable_regions)} points'
140 }
141
142 def find_rational_surfaces(self, m_max=5, n_max=3):
143 """
144 Find locations of rational surfaces q = m/n.
145
146 Parameters:
147 m_max (int): Maximum poloidal mode number
148 n_max (int): Maximum toroidal mode number
149
150 Returns:
151 list: Rational surface locations
152 """
153 rational_surfaces = []
154
155 for m in range(1, m_max + 1):
156 for n in range(1, n_max + 1):
157 q_rational = m / n
158
159 # Check if this q value exists in the profile
160 if self.q[0] <= q_rational <= self.q[-1]:
161 # Find radial location
162 try:
163 # Interpolate to find r where q = m/n
164 idx = np.where(np.diff(np.sign(self.q - q_rational)))[0]
165 if len(idx) > 0:
166 r_rational = self.r[idx[0]]
167 rational_surfaces.append({
168 'm': m, 'n': n,
169 'q': q_rational,
170 'r': r_rational,
171 'rho': r_rational / self.a
172 })
173 except:
174 pass
175
176 return rational_surfaces
177
178 def compute_tearing_mode_parameter(self):
179 """
180 Compute simplified tearing mode parameter Δ'.
181
182 Δ' > 0 indicates tearing instability.
183 This is a very simplified estimate.
184
185 Returns:
186 ndarray: Tearing parameter (approximate)
187 """
188 # Simplified: Δ' ≈ -1/q² × d²q/dr²
189 dq_dr = np.gradient(self.q, self.r)
190 d2q_dr2 = np.gradient(dq_dr, self.r)
191
192 Delta_prime = -d2q_dr2 / self.q**2
193
194 return Delta_prime
195
196 def plot_stability_analysis(self):
197 """
198 Plot comprehensive stability analysis.
199 """
200 fig = plt.figure(figsize=(14, 10))
201 gs = fig.add_gridspec(3, 2)
202
203 ax1 = fig.add_subplot(gs[0, :])
204 ax2 = fig.add_subplot(gs[1, 0])
205 ax3 = fig.add_subplot(gs[1, 1])
206 ax4 = fig.add_subplot(gs[2, 0])
207 ax5 = fig.add_subplot(gs[2, 1])
208
209 rho = self.r / self.a
210 shear = self.compute_shear()
211 suydam = self.suydam_criterion()
212 rational = self.find_rational_surfaces()
213
214 # q-profile with rational surfaces
215 ax1.plot(rho, self.q, 'b-', linewidth=2.5, label='q(r)')
216 ax1.axhline(y=1, color='r', linestyle='--', alpha=0.6,
217 label='q=1 (sawtooth)')
218 ax1.axhline(y=2, color='orange', linestyle='--', alpha=0.6,
219 label='q=2 (external kink)')
220
221 # Mark rational surfaces
222 for rs in rational:
223 ax1.plot(rs['rho'], rs['q'], 'go', markersize=8)
224 ax1.text(rs['rho'], rs['q'], f" {rs['m']}/{rs['n']}",
225 fontsize=9, verticalalignment='center')
226
227 ax1.set_xlabel(r'$\rho = r/a$', fontsize=11)
228 ax1.set_ylabel('Safety Factor q', fontsize=11)
229 ax1.set_title('Safety Factor Profile with Rational Surfaces', fontsize=13)
230 ax1.grid(True, alpha=0.3)
231 ax1.legend(fontsize=10)
232
233 # Add stability status
234 sawtooth = self.check_sawtooth_stability()
235 kink = self.check_kink_stability()
236 status_text = f"{sawtooth['message']}\n{kink['message']}"
237 ax1.text(0.02, 0.98, status_text, transform=ax1.transAxes,
238 verticalalignment='top', fontsize=10,
239 bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.8))
240
241 # Magnetic shear
242 ax2.plot(rho, shear, 'purple', linewidth=2)
243 ax2.axhline(y=0, color='k', linestyle='--', alpha=0.5)
244 ax2.fill_between(rho, 0, shear, where=(shear>0),
245 alpha=0.3, color='green', label='Positive shear')
246 ax2.fill_between(rho, shear, 0, where=(shear<0),
247 alpha=0.3, color='red', label='Negative shear')
248 ax2.set_xlabel(r'$\rho = r/a$', fontsize=11)
249 ax2.set_ylabel(r'Shear $s = (r/q)(dq/dr)$', fontsize=11)
250 ax2.set_title('Magnetic Shear Profile', fontsize=13)
251 ax2.grid(True, alpha=0.3)
252 ax2.legend(fontsize=10)
253
254 # Suydam parameter
255 ax3.plot(rho, suydam['parameter'], 'brown', linewidth=2)
256 ax3.axhline(y=0, color='k', linestyle='--', alpha=0.5,
257 label='Stability boundary')
258 ax3.fill_between(rho, 0, np.maximum(suydam['parameter'], 0),
259 alpha=0.3, color='green', label='Stable')
260 if not suydam['stable']:
261 ax3.fill_between(rho, np.minimum(suydam['parameter'], 0), 0,
262 alpha=0.3, color='red', label='Unstable')
263 ax3.set_xlabel(r'$\rho = r/a$', fontsize=11)
264 ax3.set_ylabel('Suydam Parameter', fontsize=11)
265 ax3.set_title('Suydam Stability Criterion', fontsize=13)
266 ax3.grid(True, alpha=0.3)
267 ax3.legend(fontsize=10)
268
269 # Tearing mode parameter
270 Delta_prime = self.compute_tearing_mode_parameter()
271 ax4.plot(rho, Delta_prime, 'teal', linewidth=2)
272 ax4.axhline(y=0, color='k', linestyle='--', alpha=0.5,
273 label="Δ'=0 (marginal)")
274 ax4.set_xlabel(r'$\rho = r/a$', fontsize=11)
275 ax4.set_ylabel(r"Tearing Parameter $\Delta'$", fontsize=11)
276 ax4.set_title('Tearing Mode Stability (simplified)', fontsize=13)
277 ax4.grid(True, alpha=0.3)
278 ax4.legend(fontsize=10)
279
280 # Stability diagram
281 criteria = ['Sawtooth\n(q₀>1)', 'Kink\n(qₐ>2)', 'Suydam']
282 status = [
283 1 if sawtooth['stable'] else 0,
284 1 if kink['stable'] else 0,
285 1 if suydam['stable'] else 0
286 ]
287 colors_bar = ['green' if s else 'red' for s in status]
288
289 bars = ax5.barh(criteria, status, color=colors_bar, alpha=0.7)
290 ax5.set_xlim([0, 1])
291 ax5.set_xlabel('Status', fontsize=11)
292 ax5.set_title('Stability Summary', fontsize=13)
293 ax5.set_xticks([0, 1])
294 ax5.set_xticklabels(['UNSTABLE', 'STABLE'])
295 ax5.grid(True, alpha=0.3, axis='x')
296
297 # Add checkmarks/crosses
298 for i, (bar, s) in enumerate(zip(bars, status)):
299 symbol = '✓' if s else '✗'
300 color = 'green' if s else 'red'
301 ax5.text(0.5, i, symbol, fontsize=20, color=color,
302 ha='center', va='center', weight='bold')
303
304 plt.tight_layout()
305 return fig
306
307 def generate_stability_report(self):
308 """
309 Generate text stability report.
310
311 Returns:
312 str: Stability report
313 """
314 report = []
315 report.append("=" * 60)
316 report.append("TOKAMAK Q-PROFILE STABILITY ANALYSIS REPORT")
317 report.append("=" * 60)
318
319 # Sawtooth
320 sawtooth = self.check_sawtooth_stability()
321 report.append(f"\n1. Sawtooth Stability ({sawtooth['criterion']})")
322 report.append(f" {sawtooth['message']}")
323
324 # Kink
325 kink = self.check_kink_stability()
326 report.append(f"\n2. External Kink Stability ({kink['criterion']})")
327 report.append(f" {kink['message']}")
328
329 # Suydam
330 suydam = self.suydam_criterion()
331 report.append(f"\n3. Suydam Criterion")
332 report.append(f" {suydam['message']}")
333
334 # Rational surfaces
335 rational = self.find_rational_surfaces()
336 report.append(f"\n4. Rational Surfaces (q = m/n)")
337 report.append(f" Found {len(rational)} rational surfaces:")
338 for rs in rational[:10]: # Limit to first 10
339 report.append(f" - q = {rs['m']}/{rs['n']} = {rs['q']:.3f} at ρ = {rs['rho']:.3f}")
340
341 # Overall assessment
342 report.append(f"\n" + "=" * 60)
343 all_stable = sawtooth['stable'] and kink['stable'] and suydam['stable']
344 if all_stable:
345 report.append("OVERALL: ALL STABILITY CRITERIA SATISFIED ✓")
346 else:
347 report.append("OVERALL: STABILITY ISSUES DETECTED ✗")
348 report.append("=" * 60)
349
350 return "\n".join(report)
351
352
353def main():
354 """
355 Main function demonstrating q-profile stability analysis.
356 """
357 print("=" * 60)
358 print("Q-Profile Stability Analysis")
359 print("=" * 60)
360
361 # Create sample q-profiles
362 r = np.linspace(0, 1.0, 200)
363
364 # Profile 1: Stable (monotonic, q0>1, qa>2)
365 q_stable = 1.2 + 2.0 * r**2
366
367 # Profile 2: Unstable (q0<1, sawtooth unstable)
368 q_unstable_sawtooth = 0.8 + 2.5 * r**2
369
370 # Profile 3: Unstable kink (qa<2)
371 q_unstable_kink = 1.1 + 0.7 * r**2
372
373 # Analyze stable profile
374 print("\nAnalyzing STABLE q-profile:")
375 analyzer = QProfileStability(r, q_stable)
376 report = analyzer.generate_stability_report()
377 print(report)
378
379 fig1 = analyzer.plot_stability_analysis()
380 fig1.suptitle('Stable Q-Profile Analysis', fontsize=14, weight='bold', y=0.995)
381
382 # Analyze unstable profile
383 print("\n\nAnalyzing UNSTABLE q-profile (sawtooth):")
384 analyzer2 = QProfileStability(r, q_unstable_sawtooth)
385 report2 = analyzer2.generate_stability_report()
386 print(report2)
387
388 fig2 = analyzer2.plot_stability_analysis()
389 fig2.suptitle('Unstable Q-Profile Analysis (Sawtooth)', fontsize=14, weight='bold', y=0.995)
390
391 plt.savefig('/tmp/q_profile_stability.png', dpi=150, bbox_inches='tight')
392 print("\nPlots saved to /tmp/q_profile_stability.png")
393
394 plt.show()
395
396
397if __name__ == "__main__":
398 main()