1#!/usr/bin/env python3
2"""
3MHD Eigenvalue Solver
4
5Solves the linearized MHD force operator eigenvalue problem:
6 F·ξ = -ω² ρ ξ
7
8where F is the force operator, ξ is the displacement, ω is the frequency,
9and ρ is the mass density.
10
11Discretizes the operator on a radial grid and solves the generalized
12eigenvalue problem to find growth rates for unstable modes.
13
14Author: Claude
15Date: 2026-02-14
16"""
17
18import numpy as np
19import matplotlib.pyplot as plt
20from scipy.sparse import diags
21from scipy.sparse.linalg import eigs
22from scipy.integrate import simpson
23
24
25class MHDEigenvalueSolver:
26 """
27 Solver for MHD stability eigenvalue problem.
28
29 Attributes
30 ----------
31 r : ndarray
32 Radial grid
33 n_points : int
34 Number of grid points
35 """
36
37 def __init__(self, r_max=0.01, n_points=100):
38 """
39 Initialize solver.
40
41 Parameters
42 ----------
43 r_max : float
44 Maximum radius
45 n_points : int
46 Number of radial points
47 """
48 # Avoid r=0 for numerical stability
49 self.r = np.linspace(r_max/n_points, r_max, n_points)
50 self.r_max = r_max
51 self.n_points = n_points
52 self.dr = self.r[1] - self.r[0]
53
54 def equilibrium(self, B_z, I_total, rho_0=1e-3):
55 """
56 Set up equilibrium profiles.
57
58 Parameters
59 ----------
60 B_z : float
61 Axial field (T)
62 I_total : float
63 Total current (A)
64 rho_0 : float
65 Central density (kg/m³)
66
67 Returns
68 -------
69 dict
70 Equilibrium quantities
71 """
72 mu_0 = 4 * np.pi * 1e-7
73
74 # Parabolic current density
75 J_z = 2 * I_total / (np.pi * self.r_max**2) * \
76 (1 - (self.r/self.r_max)**2)
77
78 # Azimuthal field
79 B_theta = np.zeros_like(self.r)
80 for i, ri in enumerate(self.r):
81 # Approximate enclosed current
82 r_frac = ri / self.r_max
83 I_enc = I_total * (2*r_frac**2 - r_frac**4)
84 B_theta[i] = mu_0 * I_enc / (2 * np.pi * ri)
85
86 # Pressure from force balance
87 p = np.zeros_like(self.r)
88 p_edge = 100.0
89 p[-1] = p_edge
90
91 for i in range(len(self.r)-2, -1, -1):
92 p[i] = p[i+1] - J_z[i] * B_theta[i] * self.dr
93
94 # Density profile
95 rho = rho_0 * (1 - (self.r/self.r_max)**2)
96
97 return {
98 'B_z': B_z,
99 'B_theta': B_theta,
100 'J_z': J_z,
101 'p': p,
102 'rho': rho
103 }
104
105 def construct_force_operator(self, eq, m):
106 """
107 Construct discretized force operator F.
108
109 F = -∇(δp) + (1/μ₀)(∇×δB)×B + (1/μ₀)(∇×B)×δB
110
111 Simplified for cylindrical geometry with mode number m.
112
113 Parameters
114 ----------
115 eq : dict
116 Equilibrium quantities
117 m : int
118 Poloidal mode number
119
120 Returns
121 -------
122 F : ndarray
123 Force operator matrix
124 M : ndarray
125 Mass matrix (diagonal)
126 """
127 mu_0 = 4 * np.pi * 1e-7
128 n = self.n_points
129
130 # Simplified operator for radial force balance
131 # F_r ~ d²ξ/dr² + (1/r)dξ/dr - (m²/r²)ξ + (pressure and field terms)
132
133 # Construct tridiagonal matrix for Laplacian
134 main_diag = np.zeros(n)
135 upper_diag = np.zeros(n-1)
136 lower_diag = np.zeros(n-1)
137
138 B_theta = eq['B_theta']
139 B_z = eq['B_z']
140 p = eq['p']
141 rho = eq['rho']
142
143 for i in range(1, n-1):
144 r_i = self.r[i]
145
146 # Laplacian coefficients
147 main_diag[i] = -2/self.dr**2 - m**2/r_i**2
148 upper_diag[i] = 1/self.dr**2 + 1/(2*r_i*self.dr)
149 if i > 0:
150 lower_diag[i-1] = 1/self.dr**2 - 1/(2*r_i*self.dr)
151
152 # Magnetic tension term
153 B_tot_sq = B_z**2 + B_theta[i]**2
154 main_diag[i] += B_tot_sq / (mu_0 * r_i**2)
155
156 # Pressure gradient term (stabilizing)
157 if i < n-1 and i > 0:
158 dp_dr = (p[i+1] - p[i-1]) / (2*self.dr)
159 main_diag[i] -= dp_dr / r_i
160
161 # Boundary conditions: ξ = 0 at boundaries
162 main_diag[0] = 1.0
163 main_diag[-1] = 1.0
164 if n > 1:
165 upper_diag[0] = 0.0
166 lower_diag[-1] = 0.0
167
168 # Construct sparse matrix
169 F = diags([lower_diag, main_diag, upper_diag],
170 offsets=[-1, 0, 1],
171 shape=(n, n)).toarray()
172
173 # Mass matrix (diagonal, proportional to density)
174 M = diags([rho], offsets=[0], shape=(n, n)).toarray()
175
176 return F, M
177
178 def solve_eigenvalue_problem(self, F, M, n_eigs=6):
179 """
180 Solve generalized eigenvalue problem: F·ξ = λ M·ξ
181
182 where λ = -ω².
183
184 Parameters
185 ----------
186 F : ndarray
187 Force operator
188 M : ndarray
189 Mass matrix
190 n_eigs : int
191 Number of eigenvalues to compute
192
193 Returns
194 -------
195 eigenvalues : ndarray
196 Eigenvalues (λ = -ω²)
197 eigenvectors : ndarray
198 Eigenvectors (displacement patterns)
199 """
200 # Solve standard eigenvalue problem: M^{-1}·F·ξ = λ·ξ
201 # Add small regularization to M
202 M_reg = M + np.eye(M.shape[0]) * 1e-10
203
204 # Use numpy for small problems
205 eigenvalues, eigenvectors = np.linalg.eig(
206 np.linalg.solve(M_reg, F)
207 )
208
209 # Sort by eigenvalue (most unstable first)
210 idx = np.argsort(eigenvalues.real)[::-1]
211 eigenvalues = eigenvalues[idx]
212 eigenvectors = eigenvectors[:, idx]
213
214 return eigenvalues[:n_eigs], eigenvectors[:, :n_eigs]
215
216 def compute_growth_rates(self, eigenvalues):
217 """
218 Compute growth rates from eigenvalues.
219
220 For F·ξ = -ω²·M·ξ, if λ = -ω², then:
221 - λ > 0: ω² < 0, oscillatory (stable)
222 - λ < 0: ω² > 0, ω = ±sqrt(-λ), exponential growth
223
224 Parameters
225 ----------
226 eigenvalues : ndarray
227 Eigenvalues
228
229 Returns
230 -------
231 gamma : ndarray
232 Growth rates (real part of ω)
233 """
234 gamma = np.zeros(len(eigenvalues))
235
236 for i, lam in enumerate(eigenvalues):
237 if lam.real < 0:
238 # Unstable: ω = ±i·sqrt(-λ) has real part
239 gamma[i] = np.sqrt(-lam.real)
240 else:
241 # Stable or oscillatory
242 gamma[i] = 0.0
243
244 return gamma
245
246
247def plot_eigenfunctions(solver, eigenvectors, eigenvalues, n_plot=4):
248 """
249 Plot eigenfunctions for unstable modes.
250
251 Parameters
252 ----------
253 solver : MHDEigenvalueSolver
254 Solver instance
255 eigenvectors : ndarray
256 Eigenvectors
257 eigenvalues : ndarray
258 Eigenvalues
259 n_plot : int
260 Number of modes to plot
261 """
262 fig, axes = plt.subplots(2, 2, figsize=(12, 10))
263 axes = axes.flatten()
264
265 r_cm = solver.r * 100 # Convert to cm
266
267 for i in range(min(n_plot, len(eigenvalues))):
268 ax = axes[i]
269
270 # Eigenvector
271 xi = eigenvectors[:, i].real
272 # Normalize
273 xi = xi / np.max(np.abs(xi))
274
275 ax.plot(r_cm, xi, 'b-', linewidth=2)
276 ax.axhline(y=0, color='k', linestyle='--', alpha=0.3)
277 ax.set_xlabel('Radius (cm)', fontsize=11)
278 ax.set_ylabel('Normalized ξ_r', fontsize=11)
279
280 # Growth rate
281 gamma = np.sqrt(-eigenvalues[i].real) if eigenvalues[i].real < 0 else 0
282 stability = "UNSTABLE" if gamma > 0 else "STABLE"
283 color = "red" if gamma > 0 else "green"
284
285 ax.set_title(f'Mode {i+1}: λ={eigenvalues[i].real:.2e}, γ={gamma:.2e}\n{stability}',
286 fontsize=11, fontweight='bold', color=color)
287 ax.grid(True, alpha=0.3)
288
289 plt.suptitle('MHD Eigenfunctions (Radial Displacement)',
290 fontsize=14, fontweight='bold')
291 plt.tight_layout()
292 plt.savefig('eigenfunctions.png', dpi=150, bbox_inches='tight')
293 plt.show()
294
295
296def plot_spectrum(eigenvalues, growth_rates):
297 """
298 Plot eigenvalue spectrum and growth rates.
299
300 Parameters
301 ----------
302 eigenvalues : ndarray
303 Eigenvalues
304 growth_rates : ndarray
305 Growth rates
306 """
307 fig, axes = plt.subplots(1, 2, figsize=(12, 5))
308
309 # Eigenvalue spectrum
310 ax1 = axes[0]
311 ax1.scatter(eigenvalues.real, eigenvalues.imag,
312 c=growth_rates, cmap='coolwarm',
313 s=100, edgecolors='black', linewidths=1.5)
314 ax1.axhline(y=0, color='k', linestyle='-', alpha=0.3)
315 ax1.axvline(x=0, color='k', linestyle='-', alpha=0.3)
316 ax1.set_xlabel('Re(λ)', fontsize=12)
317 ax1.set_ylabel('Im(λ)', fontsize=12)
318 ax1.set_title('Eigenvalue Spectrum', fontsize=13, fontweight='bold')
319 ax1.grid(True, alpha=0.3)
320
321 # Add stability boundary
322 ax1.axvline(x=0, color='red', linestyle='--', linewidth=2,
323 label='Stability boundary')
324 ax1.legend(fontsize=10)
325
326 # Growth rates
327 ax2 = axes[1]
328 mode_numbers = np.arange(1, len(growth_rates) + 1)
329 colors = ['red' if g > 1e-6 else 'green' for g in growth_rates]
330
331 ax2.bar(mode_numbers, growth_rates, color=colors,
332 alpha=0.7, edgecolor='black')
333 ax2.set_xlabel('Mode number', fontsize=12)
334 ax2.set_ylabel('Growth rate γ (s⁻¹)', fontsize=12)
335 ax2.set_title('Growth Rate Spectrum', fontsize=13, fontweight='bold')
336 ax2.grid(True, alpha=0.3, axis='y')
337
338 # Add threshold line
339 if np.max(growth_rates) > 0:
340 ax2.axhline(y=0, color='k', linestyle='-', linewidth=1)
341
342 plt.tight_layout()
343 plt.savefig('eigenvalue_spectrum.png', dpi=150, bbox_inches='tight')
344 plt.show()
345
346
347def main():
348 """Main execution function."""
349 # Initialize solver
350 r_max = 0.01 # 1 cm
351 solver = MHDEigenvalueSolver(r_max=r_max, n_points=100)
352
353 # Equilibrium parameters
354 B_z = 0.3 # T
355 I_total = 50e3 # 50 kA
356 rho_0 = 1e-3 # kg/m³
357
358 print("MHD Eigenvalue Solver")
359 print("=" * 60)
360 print(f"Configuration: Z-pinch")
361 print(f" Radius: {r_max*100:.1f} cm")
362 print(f" Axial field: {B_z:.2f} T")
363 print(f" Current: {I_total/1e3:.1f} kA")
364 print(f" Grid points: {solver.n_points}")
365 print()
366
367 # Compute equilibrium
368 eq = solver.equilibrium(B_z, I_total, rho_0)
369
370 print("Equilibrium computed:")
371 print(f" Peak B_θ: {np.max(eq['B_theta'])*1e3:.2f} mT")
372 print(f" Peak pressure: {np.max(eq['p'])/1e3:.2f} kPa")
373 print(f" Peak density: {np.max(eq['rho'])*1e3:.2f} g/m³")
374 print()
375
376 # Solve eigenvalue problem for m=1 mode
377 m = 1
378 print(f"Solving eigenvalue problem for m={m} mode...")
379
380 F, M = solver.construct_force_operator(eq, m)
381 eigenvalues, eigenvectors = solver.solve_eigenvalue_problem(F, M, n_eigs=8)
382 growth_rates = solver.compute_growth_rates(eigenvalues)
383
384 print(f"\nEigenvalue results:")
385 print(f" {'Mode':<6} {'Re(λ)':<12} {'Im(λ)':<12} {'γ (s⁻¹)':<12} {'Status'}")
386 print(" " + "-" * 60)
387
388 for i in range(len(eigenvalues)):
389 lam = eigenvalues[i]
390 gamma = growth_rates[i]
391 status = "UNSTABLE" if gamma > 1e-6 else "STABLE"
392
393 print(f" {i+1:<6} {lam.real:<12.3e} {lam.imag:<12.3e} "
394 f"{gamma:<12.3e} {status}")
395
396 # Count unstable modes
397 n_unstable = np.sum(growth_rates > 1e-6)
398 print(f"\nNumber of unstable modes: {n_unstable}")
399
400 if n_unstable > 0:
401 max_gamma = np.max(growth_rates)
402 max_idx = np.argmax(growth_rates)
403 print(f"Most unstable mode: #{max_idx+1} with γ = {max_gamma:.3e} s⁻¹")
404 print(f"Growth time: {1/max_gamma:.3e} s")
405
406 # Plot results
407 plot_eigenfunctions(solver, eigenvectors, eigenvalues, n_plot=4)
408 print("\nEigenfunction plot saved as 'eigenfunctions.png'")
409
410 plot_spectrum(eigenvalues, growth_rates)
411 print("Eigenvalue spectrum plot saved as 'eigenvalue_spectrum.png'")
412
413
414if __name__ == '__main__':
415 main()