1#!/usr/bin/env python3
2"""
3Growth Rate Scan Over Mode Numbers
4
5Scans over (m,n) mode numbers and computes growth rate for each using
6simplified stability analysis. Creates 2D heatmap of growth rates to
7identify the most dangerous modes.
8
9Author: Claude
10Date: 2026-02-14
11"""
12
13import numpy as np
14import matplotlib.pyplot as plt
15from scipy.integrate import simpson
16
17
18def compute_growth_rate_simple(m, n, eq, r):
19 """
20 Compute growth rate using simplified dispersion relation.
21
22 For MHD instabilities, growth rate scales roughly as:
23 γ² ~ (pressure drive - magnetic stabilization) / inertia
24
25 Parameters
26 ----------
27 m, n : int
28 Mode numbers
29 eq : dict
30 Equilibrium quantities
31 r : ndarray
32 Radial grid
33
34 Returns
35 -------
36 gamma : float
37 Growth rate (s⁻¹)
38 """
39 mu_0 = 4 * np.pi * 1e-7
40
41 B_theta = eq['B_theta']
42 B_z = eq['B_z']
43 p = eq['p']
44 rho = eq['rho']
45
46 # Wave number
47 k = np.sqrt(m**2 + n**2) / r[-1]
48
49 # Pressure gradient (destabilizing)
50 dp_dr = np.gradient(p, r)
51 pressure_drive = -np.mean(dp_dr[dp_dr < 0]) # Average negative gradient
52
53 # Magnetic stabilization
54 B_tot_sq = B_z**2 + np.mean(B_theta**2)
55 magnetic_stabilization = k**2 * B_tot_sq / mu_0
56
57 # Average density for inertia
58 rho_avg = np.mean(rho)
59
60 # Growth rate estimate
61 gamma_sq = (pressure_drive - magnetic_stabilization) / (rho_avg + 1e-10)
62
63 if gamma_sq > 0:
64 gamma = np.sqrt(gamma_sq)
65 else:
66 gamma = 0.0
67
68 # Scale by mode number (higher m,n typically have lower growth)
69 gamma = gamma / (1 + 0.1 * (m + n))
70
71 return gamma
72
73
74def scan_mode_numbers(eq, r, m_max=8, n_max=8):
75 """
76 Scan over mode numbers and compute growth rates.
77
78 Parameters
79 ----------
80 eq : dict
81 Equilibrium
82 r : ndarray
83 Radial grid
84 m_max, n_max : int
85 Maximum mode numbers
86
87 Returns
88 -------
89 gamma_map : ndarray
90 2D array of growth rates
91 m_arr, n_arr : ndarray
92 Mode number arrays
93 """
94 m_arr = np.arange(0, m_max + 1)
95 n_arr = np.arange(1, n_max + 1)
96
97 gamma_map = np.zeros((len(m_arr), len(n_arr)))
98
99 for i, m in enumerate(m_arr):
100 for j, n in enumerate(n_arr):
101 m_eff = max(m, 1) # Avoid m=0
102 gamma_map[i, j] = compute_growth_rate_simple(m_eff, n, eq, r)
103
104 return gamma_map, m_arr, n_arr
105
106
107def create_equilibrium(r):
108 """
109 Create Z-pinch equilibrium.
110
111 Parameters
112 ----------
113 r : ndarray
114 Radial grid
115
116 Returns
117 -------
118 eq : dict
119 Equilibrium quantities
120 """
121 mu_0 = 4 * np.pi * 1e-7
122 r_max = r[-1]
123
124 # Parameters
125 B_z = 0.4 # T
126 I_total = 60e3 # 60 kA
127
128 # Current density
129 J_z = 2 * I_total / (np.pi * r_max**2) * (1 - (r/r_max)**2)
130
131 # Azimuthal field
132 B_theta = np.zeros_like(r)
133 for i, ri in enumerate(r):
134 if ri > 1e-10:
135 r_frac = ri / r_max
136 I_enc = I_total * (2*r_frac**2 - r_frac**4)
137 B_theta[i] = mu_0 * I_enc / (2 * np.pi * ri)
138
139 # Pressure
140 p = np.zeros_like(r)
141 p_edge = 100.0
142 p[-1] = p_edge
143
144 for i in range(len(r)-2, -1, -1):
145 dr = r[i+1] - r[i]
146 p[i] = p[i+1] - J_z[i] * B_theta[i] * dr
147
148 # Density
149 rho_0 = 1e-3 # kg/m³
150 rho = rho_0 * (1 - (r/r_max)**2)
151
152 return {
153 'B_z': B_z,
154 'B_theta': B_theta,
155 'J_z': J_z,
156 'p': p,
157 'rho': rho
158 }
159
160
161def plot_growth_rate_heatmap(gamma_map, m_arr, n_arr):
162 """
163 Plot 2D heatmap of growth rates.
164
165 Parameters
166 ----------
167 gamma_map : ndarray
168 Growth rate map
169 m_arr, n_arr : ndarray
170 Mode numbers
171 """
172 fig, ax = plt.subplots(figsize=(11, 8))
173
174 # Create heatmap
175 im = ax.imshow(gamma_map, cmap='hot', aspect='auto',
176 extent=[n_arr[0]-0.5, n_arr[-1]+0.5,
177 m_arr[-1]+0.5, m_arr[0]-0.5],
178 interpolation='bilinear')
179
180 # Colorbar
181 cbar = plt.colorbar(im, ax=ax)
182 cbar.set_label('Growth rate γ (s⁻¹)', fontsize=13)
183
184 # Find and mark most dangerous mode
185 max_idx = np.unravel_index(np.argmax(gamma_map), gamma_map.shape)
186 m_max = m_arr[max_idx[0]]
187 n_max = n_arr[max_idx[1]]
188 gamma_max = gamma_map[max_idx]
189
190 ax.plot(n_max, m_max, 'c*', markersize=25, markeredgecolor='white',
191 markeredgewidth=2, label=f'Most unstable: ({m_max},{n_max})')
192
193 # Contour lines
194 if np.max(gamma_map) > 0:
195 levels = np.linspace(0, np.max(gamma_map), 6)[1:]
196 cs = ax.contour(n_arr, m_arr, gamma_map, levels=levels,
197 colors='cyan', linewidths=1.5, alpha=0.6)
198 ax.clabel(cs, inline=True, fontsize=9, fmt='%.1e')
199
200 ax.set_xlabel('Toroidal mode number n', fontsize=13)
201 ax.set_ylabel('Poloidal mode number m', fontsize=13)
202 ax.set_title('MHD Growth Rate Map: γ(m,n)',
203 fontsize=14, fontweight='bold')
204 ax.legend(fontsize=11, loc='upper right')
205
206 # Grid
207 ax.set_xticks(n_arr)
208 ax.set_yticks(m_arr)
209 ax.grid(True, color='white', linestyle='--', linewidth=0.5, alpha=0.3)
210
211 plt.tight_layout()
212 plt.savefig('growth_rate_heatmap.png', dpi=150, bbox_inches='tight')
213 plt.show()
214
215
216def plot_growth_vs_mode(gamma_map, m_arr, n_arr):
217 """
218 Plot growth rate vs mode number for different n.
219
220 Parameters
221 ----------
222 gamma_map : ndarray
223 Growth rates
224 m_arr, n_arr : ndarray
225 Mode numbers
226 """
227 fig, axes = plt.subplots(1, 2, figsize=(14, 5))
228
229 # Plot 1: γ vs m for different n
230 ax1 = axes[0]
231 colors = plt.cm.viridis(np.linspace(0, 1, len(n_arr)))
232
233 for j, n in enumerate(n_arr):
234 ax1.plot(m_arr, gamma_map[:, j], 'o-', color=colors[j],
235 linewidth=2, markersize=6, label=f'n={n}')
236
237 ax1.set_xlabel('Poloidal mode number m', fontsize=12)
238 ax1.set_ylabel('Growth rate γ (s⁻¹)', fontsize=12)
239 ax1.set_title('Growth Rate vs m (different n)', fontsize=13, fontweight='bold')
240 ax1.legend(fontsize=9, ncol=2)
241 ax1.grid(True, alpha=0.3)
242 ax1.set_xlim([m_arr[0]-0.5, m_arr[-1]+0.5])
243
244 # Plot 2: γ vs n for different m
245 ax2 = axes[1]
246 colors2 = plt.cm.plasma(np.linspace(0, 1, len(m_arr)))
247
248 for i, m in enumerate(m_arr):
249 if m <= 5: # Plot only first few for clarity
250 ax2.plot(n_arr, gamma_map[i, :], 's-', color=colors2[i],
251 linewidth=2, markersize=6, label=f'm={m}')
252
253 ax2.set_xlabel('Toroidal mode number n', fontsize=12)
254 ax2.set_ylabel('Growth rate γ (s⁻¹)', fontsize=12)
255 ax2.set_title('Growth Rate vs n (different m)', fontsize=13, fontweight='bold')
256 ax2.legend(fontsize=9)
257 ax2.grid(True, alpha=0.3)
258 ax2.set_xlim([n_arr[0]-0.5, n_arr[-1]+0.5])
259
260 plt.tight_layout()
261 plt.savefig('growth_rate_vs_mode.png', dpi=150, bbox_inches='tight')
262 plt.show()
263
264
265def plot_most_dangerous_modes(gamma_map, m_arr, n_arr, top_n=10):
266 """
267 Plot bar chart of most dangerous modes.
268
269 Parameters
270 ----------
271 gamma_map : ndarray
272 Growth rates
273 m_arr, n_arr : ndarray
274 Mode numbers
275 top_n : int
276 Number of top modes to show
277 """
278 # Flatten and sort
279 gamma_flat = gamma_map.flatten()
280 idx_flat = np.argsort(gamma_flat)[::-1][:top_n]
281
282 # Get mode numbers for top modes
283 modes = []
284 gammas = []
285
286 for idx in idx_flat:
287 i, j = np.unravel_index(idx, gamma_map.shape)
288 m = m_arr[i]
289 n = n_arr[j]
290 gamma = gamma_map[i, j]
291
292 modes.append(f'({m},{n})')
293 gammas.append(gamma)
294
295 # Plot
296 fig, ax = plt.subplots(figsize=(11, 6))
297
298 colors = plt.cm.Reds(np.linspace(0.4, 0.9, len(gammas)))
299 bars = ax.barh(modes, gammas, color=colors, edgecolor='black', linewidth=1.5)
300
301 ax.set_xlabel('Growth rate γ (s⁻¹)', fontsize=13)
302 ax.set_ylabel('Mode (m,n)', fontsize=13)
303 ax.set_title(f'Top {top_n} Most Dangerous Modes',
304 fontsize=14, fontweight='bold')
305 ax.grid(True, alpha=0.3, axis='x')
306
307 # Add value labels
308 for bar, val in zip(bars, gammas):
309 width = bar.get_width()
310 ax.text(width, bar.get_y() + bar.get_height()/2,
311 f' {val:.2e}',
312 ha='left', va='center', fontsize=10, fontweight='bold')
313
314 plt.tight_layout()
315 plt.savefig('most_dangerous_modes.png', dpi=150, bbox_inches='tight')
316 plt.show()
317
318
319def main():
320 """Main execution function."""
321 # Create radial grid
322 r_max = 0.01 # 1 cm
323 r = np.linspace(r_max/200, r_max, 200)
324
325 print("Growth Rate Scan Over Mode Numbers")
326 print("=" * 60)
327 print(f"Configuration: Z-pinch")
328 print(f" Radius: {r_max*100:.1f} cm")
329 print()
330
331 # Create equilibrium
332 print("Creating equilibrium...")
333 eq = create_equilibrium(r)
334
335 print(f" B_z: {eq['B_z']:.2f} T")
336 print(f" Peak B_θ: {np.max(eq['B_theta'])*1e3:.2f} mT")
337 print(f" Peak pressure: {np.max(eq['p'])/1e3:.2f} kPa")
338 print()
339
340 # Scan mode numbers
341 m_max = 8
342 n_max = 8
343 print(f"Scanning mode numbers (0 ≤ m ≤ {m_max}, 1 ≤ n ≤ {n_max})...")
344
345 gamma_map, m_arr, n_arr = scan_mode_numbers(eq, r, m_max, n_max)
346
347 # Find most dangerous mode
348 max_idx = np.unravel_index(np.argmax(gamma_map), gamma_map.shape)
349 m_danger = m_arr[max_idx[0]]
350 n_danger = n_arr[max_idx[1]]
351 gamma_max = gamma_map[max_idx]
352
353 print(f"\nMost dangerous mode: (m={m_danger}, n={n_danger})")
354 print(f" Growth rate: {gamma_max:.3e} s⁻¹")
355 if gamma_max > 0:
356 print(f" Growth time: {1/gamma_max:.3e} s")
357 print(f" e-folding time: {1/gamma_max:.3e} s")
358
359 # Count unstable modes
360 n_unstable = np.sum(gamma_map > 1e-6)
361 n_total = gamma_map.size
362 print(f"\nUnstable modes: {n_unstable}/{n_total} ({100*n_unstable/n_total:.1f}%)")
363
364 # Plot results
365 print("\nGenerating plots...")
366 plot_growth_rate_heatmap(gamma_map, m_arr, n_arr)
367 print(" Heatmap saved as 'growth_rate_heatmap.png'")
368
369 plot_growth_vs_mode(gamma_map, m_arr, n_arr)
370 print(" Mode dependence saved as 'growth_rate_vs_mode.png'")
371
372 plot_most_dangerous_modes(gamma_map, m_arr, n_arr, top_n=10)
373 print(" Top modes saved as 'most_dangerous_modes.png'")
374
375
376if __name__ == '__main__':
377 main()