1#!/usr/bin/env python3
2"""
32D Magnetic Reconnection Simulation
4
5Simplified 2D resistive MHD simulation of magnetic reconnection in a Harris current sheet.
6This demonstrates the fundamental process by which magnetic field lines break and reconnect,
7releasing stored magnetic energy as kinetic and thermal energy.
8
9Key results:
10- Harris current sheet evolves through resistive reconnection
11- Formation of X-point and magnetic islands
12- Current density concentrates at reconnection sites
13- Magnetic energy converts to kinetic/thermal energy
14
15Physics:
16- Ideal MHD with resistivity: ∂B/∂t = ∇×(v×B) - ∇×(ηJ)
17- Momentum: ρ∂v/∂t = J×B - ∇p
18- Incompressible flow: ∇·v = 0
19
20Author: Claude
21Date: 2026-02-14
22"""
23
24import numpy as np
25import matplotlib.pyplot as plt
26from matplotlib.patches import FancyArrowPatch
27from typing import Tuple
28
29
30class MagneticReconnection2D:
31 """
32 2D Magnetic Reconnection Solver using simplified resistive MHD.
33
34 Uses the vector potential formulation where B = ∇×A (with A = A_z ẑ in 2D).
35 Evolution equations:
36 - ∂A/∂t = -E_z = (v×B)_z - η J_z
37 - ∂ω/∂t = ∇×(J×B/ρ) - ν∇²ω (vorticity equation)
38 where ω = ∇²ψ is the vorticity and v = ∇×ψ ŷ
39 """
40
41 def __init__(self, nx: int = 128, ny: int = 128,
42 Lx: float = 20.0, Ly: float = 10.0,
43 eta: float = 0.01, nu: float = 0.001):
44 """
45 Initialize the reconnection simulation.
46
47 Args:
48 nx, ny: Grid resolution
49 Lx, Ly: Domain size
50 eta: Magnetic resistivity (controls reconnection rate)
51 nu: Kinematic viscosity (for numerical stability)
52 """
53 self.nx, self.ny = nx, ny
54 self.Lx, self.Ly = Lx, Ly
55 self.eta = eta
56 self.nu = nu
57
58 # Grid spacing
59 self.dx = Lx / (nx - 1)
60 self.dy = Ly / (ny - 1)
61
62 # Coordinate arrays
63 self.x = np.linspace(-Lx/2, Lx/2, nx)
64 self.y = np.linspace(-Ly/2, Ly/2, ny)
65 self.X, self.Y = np.meshgrid(self.x, self.y, indexing='ij')
66
67 # Field arrays
68 self.A = np.zeros((nx, ny)) # Vector potential (magnetic flux function)
69 self.psi = np.zeros((nx, ny)) # Stream function (velocity potential)
70 self.Bx = np.zeros((nx, ny)) # Magnetic field components
71 self.By = np.zeros((nx, ny))
72 self.Jz = np.zeros((nx, ny)) # Current density (out-of-plane)
73 self.vx = np.zeros((nx, ny)) # Velocity components
74 self.vy = np.zeros((nx, ny))
75
76 # Energy tracking
77 self.time = 0.0
78 self.magnetic_energy = []
79 self.kinetic_energy = []
80 self.time_history = []
81
82 def harris_current_sheet(self, B0: float = 1.0, width: float = 1.0,
83 perturbation: float = 0.1):
84 """
85 Initialize Harris current sheet equilibrium with perturbation.
86
87 The Harris sheet is a classic 1D equilibrium with:
88 - B_x(y) = B0 * tanh(y/width)
89 - J_z = -B0/width * sech²(y/width)
90
91 We add a small perturbation to trigger reconnection.
92
93 Args:
94 B0: Asymptotic magnetic field strength
95 width: Current sheet thickness
96 perturbation: Amplitude of initial perturbation
97 """
98 # Harris sheet: A = -B0 * width * ln(cosh(y/width))
99 for i in range(self.nx):
100 for j in range(self.ny):
101 y = self.Y[i, j]
102 x = self.X[i, j]
103
104 # Base Harris profile
105 self.A[i, j] = -B0 * width * np.log(np.cosh(y / width))
106
107 # Add perturbation to seed reconnection
108 # Symmetric tearing mode: δA ∝ cos(kx) * sech(y/width)
109 k = 2.0 * np.pi / self.Lx
110 self.A[i, j] += perturbation * B0 * np.cos(k * x) / np.cosh(y / width)
111
112 # Initialize stream function (small or zero)
113 self.psi = np.zeros_like(self.A)
114
115 # Compute derived fields
116 self._update_derived_fields()
117
118 def _update_derived_fields(self):
119 """Compute B, J, v from potentials A and psi."""
120 # Magnetic field: B = ∇×A = (∂A/∂y, -∂A/∂x, 0)
121 self.By, self.Bx = np.gradient(self.A, self.dx, self.dy)
122 self.Bx = -self.Bx
123
124 # Current density: J_z = -∇²A (Ampere's law in 2D)
125 self.Jz = -self._laplacian(self.A)
126
127 # Velocity: v = ∇×psi ŷ = (-∂psi/∂y, ∂psi/∂x, 0)
128 self.vy, self.vx = np.gradient(self.psi, self.dx, self.dy)
129 self.vx = -self.vx
130
131 def _laplacian(self, field: np.ndarray) -> np.ndarray:
132 """Compute 2D Laplacian using centered finite differences."""
133 laplacian = np.zeros_like(field)
134
135 laplacian[1:-1, 1:-1] = (
136 (field[2:, 1:-1] - 2*field[1:-1, 1:-1] + field[:-2, 1:-1]) / self.dx**2 +
137 (field[1:-1, 2:] - 2*field[1:-1, 1:-1] + field[1:-1, :-2]) / self.dy**2
138 )
139
140 return laplacian
141
142 def _poisson_solve(self, rhs: np.ndarray, max_iter: int = 1000,
143 tol: float = 1e-5) -> np.ndarray:
144 """
145 Solve Poisson equation ∇²φ = rhs using Jacobi iteration.
146
147 Used for: ∇²ψ = ω (vorticity → stream function)
148 """
149 phi = np.zeros_like(rhs)
150 dx2 = self.dx**2
151 dy2 = self.dy**2
152 denom = 2.0 * (dx2 + dy2)
153
154 for iteration in range(max_iter):
155 phi_old = phi.copy()
156
157 phi[1:-1, 1:-1] = (
158 dx2 * (phi_old[1:-1, 2:] + phi_old[1:-1, :-2]) +
159 dy2 * (phi_old[2:, 1:-1] + phi_old[:-2, 1:-1]) -
160 dx2 * dy2 * rhs[1:-1, 1:-1]
161 ) / denom
162
163 # Boundary conditions: φ = 0 at boundaries
164 phi[0, :] = phi[-1, :] = phi[:, 0] = phi[:, -1] = 0
165
166 # Check convergence
167 error = np.max(np.abs(phi - phi_old))
168 if error < tol:
169 break
170
171 return phi
172
173 def step(self, dt: float):
174 """
175 Advance the simulation by one time step using forward Euler.
176
177 Evolution equations:
178 1. Induction: ∂A/∂t = (v×B)_z - η J_z = v_x B_y - v_y B_x - η J_z
179 2. Vorticity: ∂ω/∂t = (∇×(J×B))_z - ν∇²ω
180 """
181 # Compute vorticity
182 omega = -self._laplacian(self.psi)
183
184 # 1. Advance vector potential A
185 # ∂A/∂t = E_z = v×B - ηJ
186 vxB = self.vx * self.By - self.vy * self.Bx
187 dA_dt = vxB - self.eta * self.Jz
188
189 self.A += dt * dA_dt
190
191 # Boundary conditions for A (conducting walls)
192 self.A[0, :] = self.A[-1, :] = self.A[:, 0] = self.A[:, -1] = 0
193
194 # 2. Advance vorticity ω
195 # J×B force term
196 JxB_x = self.Jz * self.By
197 JxB_y = -self.Jz * self.Bx
198
199 # Curl of J×B (simplified, assuming constant density ρ=1)
200 dJxB_y_dx = np.gradient(JxB_y, self.dx, axis=0)
201 dJxB_x_dy = np.gradient(JxB_x, self.dy, axis=1)
202 curl_JxB = dJxB_y_dx - dJxB_x_dy
203
204 # Viscous term
205 visc_omega = self.nu * self._laplacian(omega)
206
207 domega_dt = curl_JxB + visc_omega
208 omega += dt * domega_dt
209
210 # 3. Recover stream function from vorticity: ∇²ψ = ω
211 self.psi = self._poisson_solve(omega)
212
213 # 4. Update all derived fields
214 self._update_derived_fields()
215
216 # Update time
217 self.time += dt
218
219 # Track energies
220 self._compute_energies()
221
222 def _compute_energies(self):
223 """Compute magnetic and kinetic energies."""
224 # Magnetic energy: ∫ B²/2 dV
225 B_squared = self.Bx**2 + self.By**2
226 E_mag = 0.5 * np.sum(B_squared) * self.dx * self.dy
227
228 # Kinetic energy: ∫ ρv²/2 dV (ρ=1)
229 v_squared = self.vx**2 + self.vy**2
230 E_kin = 0.5 * np.sum(v_squared) * self.dx * self.dy
231
232 self.magnetic_energy.append(E_mag)
233 self.kinetic_energy.append(E_kin)
234 self.time_history.append(self.time)
235
236 def find_xpoint(self) -> Tuple[float, float]:
237 """
238 Find the X-point location (magnetic null point where B=0).
239
240 Returns:
241 (x, y) coordinates of X-point
242 """
243 # Find minimum of |B|
244 B_magnitude = np.sqrt(self.Bx**2 + self.By**2)
245
246 # Search in central region
247 nx_mid = self.nx // 2
248 ny_mid = self.ny // 2
249 search_region = B_magnitude[nx_mid-20:nx_mid+20, ny_mid-20:ny_mid+20]
250
251 i_min, j_min = np.unravel_index(search_region.argmin(), search_region.shape)
252 i_min += nx_mid - 20
253 j_min += ny_mid - 20
254
255 return self.X[i_min, j_min], self.Y[i_min, j_min]
256
257
258def visualize_reconnection(sim: MagneticReconnection2D,
259 save_prefix: str = "reconnection"):
260 """
261 Create comprehensive visualization of magnetic reconnection.
262
263 Shows:
264 1. Current density with magnetic field lines
265 2. Velocity field
266 3. Energy evolution
267 """
268 fig = plt.figure(figsize=(16, 10))
269
270 # 1. Current density and magnetic field lines
271 ax1 = plt.subplot(2, 3, 1)
272
273 # Current density color map
274 levels = np.linspace(-1.5, 1.5, 30)
275 cf = ax1.contourf(sim.X, sim.Y, sim.Jz, levels=levels, cmap='RdBu_r', extend='both')
276 plt.colorbar(cf, ax=ax1, label='Current density $J_z$')
277
278 # Magnetic field lines (contours of A)
279 A_levels = np.linspace(sim.A.min(), sim.A.max(), 20)
280 ax1.contour(sim.X, sim.Y, sim.A, levels=A_levels, colors='black',
281 linewidths=0.8, alpha=0.6)
282
283 # Mark X-point
284 try:
285 x_xpt, y_xpt = sim.find_xpoint()
286 ax1.plot(x_xpt, y_xpt, 'g*', markersize=15, label='X-point')
287 ax1.legend()
288 except:
289 pass
290
291 ax1.set_xlabel('x')
292 ax1.set_ylabel('y')
293 ax1.set_title(f'Current Density and B-field Lines (t={sim.time:.2f})')
294 ax1.set_aspect('equal')
295
296 # 2. Magnetic field vectors
297 ax2 = plt.subplot(2, 3, 2)
298
299 # Subsample for quiver plot
300 skip = 4
301 ax2.quiver(sim.X[::skip, ::skip], sim.Y[::skip, ::skip],
302 sim.Bx[::skip, ::skip], sim.By[::skip, ::skip],
303 alpha=0.7)
304 ax2.set_xlabel('x')
305 ax2.set_ylabel('y')
306 ax2.set_title('Magnetic Field Vectors')
307 ax2.set_aspect('equal')
308
309 # 3. Velocity field
310 ax3 = plt.subplot(2, 3, 3)
311
312 v_magnitude = np.sqrt(sim.vx**2 + sim.vy**2)
313 cf3 = ax3.contourf(sim.X, sim.Y, v_magnitude, levels=20, cmap='plasma')
314 plt.colorbar(cf3, ax=ax3, label='|v|')
315
316 # Velocity streamlines
317 ax3.streamplot(sim.x, sim.y, sim.vx.T, sim.vy.T,
318 color='white', linewidth=0.5, density=1.5, arrowsize=0.8)
319 ax3.set_xlabel('x')
320 ax3.set_ylabel('y')
321 ax3.set_title('Velocity Field')
322 ax3.set_aspect('equal')
323
324 # 4. Energy evolution
325 ax4 = plt.subplot(2, 3, 4)
326
327 times = np.array(sim.time_history)
328 E_mag = np.array(sim.magnetic_energy)
329 E_kin = np.array(sim.kinetic_energy)
330 E_total = E_mag + E_kin
331
332 ax4.plot(times, E_mag, 'b-', label='Magnetic', linewidth=2)
333 ax4.plot(times, E_kin, 'r-', label='Kinetic', linewidth=2)
334 ax4.plot(times, E_total, 'k--', label='Total', linewidth=1.5)
335 ax4.set_xlabel('Time')
336 ax4.set_ylabel('Energy')
337 ax4.set_title('Energy Evolution')
338 ax4.legend()
339 ax4.grid(True, alpha=0.3)
340
341 # 5. Current density profile at y=0
342 ax5 = plt.subplot(2, 3, 5)
343
344 j_mid = sim.ny // 2
345 ax5.plot(sim.x, sim.Jz[:, j_mid], 'b-', linewidth=2)
346 ax5.set_xlabel('x')
347 ax5.set_ylabel('$J_z(x, y=0)$')
348 ax5.set_title('Current Sheet Profile')
349 ax5.grid(True, alpha=0.3)
350 ax5.axhline(y=0, color='k', linestyle='--', alpha=0.5)
351
352 # 6. Magnetic field profile at x=0
353 ax6 = plt.subplot(2, 3, 6)
354
355 i_mid = sim.nx // 2
356 ax6.plot(sim.y, sim.Bx[i_mid, :], 'b-', label='$B_x$', linewidth=2)
357 ax6.plot(sim.y, sim.By[i_mid, :], 'r-', label='$B_y$', linewidth=2)
358 ax6.set_xlabel('y')
359 ax6.set_ylabel('B(x=0, y)')
360 ax6.set_title('Magnetic Field Profile')
361 ax6.legend()
362 ax6.grid(True, alpha=0.3)
363 ax6.axhline(y=0, color='k', linestyle='--', alpha=0.5)
364
365 plt.tight_layout()
366 plt.savefig(f'/opt/projects/01_Personal/03_Study/examples/Numerical_Simulation/10_projects/{save_prefix}.png',
367 dpi=150, bbox_inches='tight')
368 plt.close()
369 print(f"Figure saved: {save_prefix}.png")
370
371
372def run_reconnection_simulation():
373 """Run a complete magnetic reconnection simulation."""
374 print("=" * 70)
375 print("2D Magnetic Reconnection Simulation")
376 print("=" * 70)
377
378 # Initialize simulation
379 print("\nInitializing Harris current sheet...")
380 sim = MagneticReconnection2D(nx=128, ny=64, Lx=20.0, Ly=10.0,
381 eta=0.01, nu=0.001)
382 sim.harris_current_sheet(B0=1.0, width=1.0, perturbation=0.1)
383
384 print(f" Grid: {sim.nx} × {sim.ny}")
385 print(f" Domain: [{-sim.Lx/2:.1f}, {sim.Lx/2:.1f}] × [{-sim.Ly/2:.1f}, {sim.Ly/2:.1f}]")
386 print(f" Resistivity η = {sim.eta}")
387 print(f" Initial magnetic energy: {sim.magnetic_energy[0]:.4f}")
388
389 # Time stepping
390 dt = 0.01
391 t_final = 50.0
392 n_steps = int(t_final / dt)
393 output_interval = 500
394
395 print(f"\nTime integration:")
396 print(f" dt = {dt}, t_final = {t_final}")
397 print(f" Total steps: {n_steps}")
398
399 # Initial state
400 print("\nSaving initial state...")
401 visualize_reconnection(sim, save_prefix="reconnection_t000")
402
403 # Time evolution
404 print("\nEvolving system...")
405 for step in range(n_steps):
406 sim.step(dt)
407
408 if (step + 1) % output_interval == 0:
409 print(f" Step {step+1}/{n_steps}, t={sim.time:.2f}, "
410 f"E_mag={sim.magnetic_energy[-1]:.4f}, "
411 f"E_kin={sim.kinetic_energy[-1]:.4f}")
412
413 # Save snapshot
414 visualize_reconnection(sim,
415 save_prefix=f"reconnection_t{int(sim.time):03d}")
416
417 # Final state
418 print("\nFinal state:")
419 print(f" Time: {sim.time:.2f}")
420 print(f" Magnetic energy: {sim.magnetic_energy[-1]:.4f} "
421 f"(change: {sim.magnetic_energy[-1]-sim.magnetic_energy[0]:.4f})")
422 print(f" Kinetic energy: {sim.kinetic_energy[-1]:.4f}")
423 print(f" Total energy: {sim.magnetic_energy[-1]+sim.kinetic_energy[-1]:.4f}")
424
425 # Find X-point
426 try:
427 x_xpt, y_xpt = sim.find_xpoint()
428 print(f" X-point location: ({x_xpt:.2f}, {y_xpt:.2f})")
429 except:
430 print(" X-point not clearly identified")
431
432 print("\n" + "=" * 70)
433 print("Simulation complete!")
434 print("=" * 70)
435 print("""
436 Physical Interpretation:
437
438 1. Initial Harris sheet: Anti-parallel magnetic field with thin current layer
439 2. Perturbation triggers instability (tearing mode)
440 3. X-point forms where field lines break and reconnect
441 4. Magnetic energy converts to kinetic energy and heat
442 5. Current sheet fragments into magnetic islands (plasmoids)
443
444 Key Parameters:
445 - Resistivity η: Controls reconnection rate (higher η → faster reconnection)
446 - Perturbation: Seeds the instability
447 - Aspect ratio: Affects island formation
448
449 Applications:
450 - Solar flares and coronal mass ejections
451 - Magnetospheric substorms
452 - Tokamak disruptions
453 - Astrophysical jets
454 """)
455
456
457if __name__ == "__main__":
458 run_reconnection_simulation()