1#!/usr/bin/env python3
2"""
3Alpha-Omega Mean-Field Dynamo Model
4
5This module implements a 1D α-Ω mean-field dynamo model, which is used to
6understand magnetic field generation in rotating systems like stars and planets.
7
8The governing equations are:
9 ∂A/∂t = αB + η∂²A/∂r² (poloidal field evolution)
10 ∂B/∂t = S(∂A/∂r) + η∂²B/∂r² (toroidal field evolution)
11
12where:
13 A = poloidal field potential
14 B = toroidal field
15 α = alpha effect (helicity parameter)
16 η = magnetic diffusivity
17 S = shear parameter (differential rotation, Ω-effect)
18
19The dynamo number D = αS/η² determines whether the dynamo is self-sustaining.
20
21Author: MHD Course Examples
22License: MIT
23"""
24
25import numpy as np
26import matplotlib.pyplot as plt
27from matplotlib.animation import FuncAnimation
28from scipy.sparse import diags
29from scipy.sparse.linalg import spsolve
30
31
32class AlphaOmegaDynamo:
33 """
34 1D α-Ω mean-field dynamo model solver.
35
36 Attributes:
37 nr (int): Number of radial grid points
38 r (ndarray): Radial grid
39 dr (float): Grid spacing
40 alpha (float): Alpha effect strength
41 eta (float): Magnetic diffusivity
42 shear (float): Shear parameter S
43 dynamo_number (float): D = αS/η²
44 """
45
46 def __init__(self, nr=100, r_min=0.5, r_max=1.0, alpha=1.0,
47 eta=0.1, shear=10.0):
48 """
49 Initialize the dynamo model.
50
51 Parameters:
52 nr (int): Number of radial grid points
53 r_min (float): Inner radius
54 r_max (float): Outer radius
55 alpha (float): Alpha effect strength
56 eta (float): Magnetic diffusivity
57 shear (float): Shear parameter
58 """
59 self.nr = nr
60 self.r = np.linspace(r_min, r_max, nr)
61 self.dr = self.r[1] - self.r[0]
62 self.alpha = alpha
63 self.eta = eta
64 self.shear = shear
65 self.dynamo_number = (alpha * shear) / (eta**2)
66
67 # Initialize fields
68 self.A = np.zeros(nr) # Poloidal potential
69 self.B = np.zeros(nr) # Toroidal field
70
71 # Time history for butterfly diagram
72 self.time_history = []
73 self.B_history = []
74
75 def diffusion_operator(self):
76 """
77 Construct the finite difference diffusion operator ∂²/∂r².
78
79 Returns:
80 scipy.sparse matrix: Diffusion operator
81 """
82 # Second derivative with periodic boundary conditions
83 diag_main = -2.0 * np.ones(self.nr) / self.dr**2
84 diag_off = np.ones(self.nr - 1) / self.dr**2
85
86 # Periodic boundaries
87 data = [diag_main, diag_off, diag_off]
88 offsets = [0, 1, -1]
89 D = diags(data, offsets, shape=(self.nr, self.nr), format='csr')
90
91 # Wrap around for periodic BC
92 D[0, -1] = 1.0 / self.dr**2
93 D[-1, 0] = 1.0 / self.dr**2
94
95 return D
96
97 def gradient(self, f):
98 """
99 Compute ∂f/∂r using centered differences.
100
101 Parameters:
102 f (ndarray): Field to differentiate
103
104 Returns:
105 ndarray: Gradient of f
106 """
107 df = np.zeros_like(f)
108 df[1:-1] = (f[2:] - f[:-2]) / (2 * self.dr)
109 # Periodic boundaries
110 df[0] = (f[1] - f[-1]) / (2 * self.dr)
111 df[-1] = (f[0] - f[-2]) / (2 * self.dr)
112 return df
113
114 def step_explicit(self, dt):
115 """
116 Time step using explicit Euler method.
117
118 Parameters:
119 dt (float): Time step
120 """
121 # Compute spatial derivatives
122 d2A_dr2 = self.eta * self.diffusion_operator().dot(self.A)
123 d2B_dr2 = self.eta * self.diffusion_operator().dot(self.B)
124 dA_dr = self.gradient(self.A)
125
126 # Update equations
127 dA_dt = self.alpha * self.B + d2A_dr2
128 dB_dt = self.shear * dA_dr + d2B_dr2
129
130 # Euler step
131 self.A += dt * dA_dt
132 self.B += dt * dB_dt
133
134 def step_implicit(self, dt):
135 """
136 Time step using implicit method (Crank-Nicolson).
137
138 Parameters:
139 dt (float): Time step
140 """
141 I = np.eye(self.nr)
142 D = self.diffusion_operator().toarray()
143
144 # Crank-Nicolson matrices
145 LHS_A = I - 0.5 * dt * self.eta * D
146 RHS_A = I + 0.5 * dt * self.eta * D
147
148 LHS_B = I - 0.5 * dt * self.eta * D
149 RHS_B = I + 0.5 * dt * self.eta * D
150
151 # Source terms
152 source_A = dt * self.alpha * self.B
153 source_B = dt * self.shear * self.gradient(self.A)
154
155 # Solve linear systems
156 self.A = np.linalg.solve(LHS_A, RHS_A @ self.A + source_A)
157 self.B = np.linalg.solve(LHS_B, RHS_B @ self.B + source_B)
158
159 def set_initial_condition(self, mode='random'):
160 """
161 Set initial magnetic field perturbation.
162
163 Parameters:
164 mode (str): 'random' or 'sinusoidal'
165 """
166 if mode == 'random':
167 self.A = 0.01 * np.random.randn(self.nr)
168 self.B = 0.01 * np.random.randn(self.nr)
169 elif mode == 'sinusoidal':
170 self.A = 0.01 * np.sin(2 * np.pi * (self.r - self.r[0]) /
171 (self.r[-1] - self.r[0]))
172 self.B = 0.01 * np.cos(2 * np.pi * (self.r - self.r[0]) /
173 (self.r[-1] - self.r[0]))
174
175 def run_simulation(self, t_end=10.0, dt=0.01, save_interval=10):
176 """
177 Run the dynamo simulation.
178
179 Parameters:
180 t_end (float): End time
181 dt (float): Time step
182 save_interval (int): Save data every N steps
183 """
184 n_steps = int(t_end / dt)
185 t = 0.0
186
187 print(f"Running α-Ω dynamo simulation...")
188 print(f"Dynamo number D = {self.dynamo_number:.2f}")
189 print(f"Critical D ~ 10 for oscillatory dynamo")
190
191 for step in range(n_steps):
192 self.step_implicit(dt)
193 t += dt
194
195 # Save history for butterfly diagram
196 if step % save_interval == 0:
197 self.time_history.append(t)
198 self.B_history.append(self.B.copy())
199
200 if step % 1000 == 0:
201 B_max = np.max(np.abs(self.B))
202 print(f"Step {step}/{n_steps}, t={t:.2f}, max|B|={B_max:.4f}")
203
204 self.B_history = np.array(self.B_history)
205 print("Simulation complete!")
206
207 def plot_butterfly_diagram(self):
208 """
209 Plot the butterfly diagram (time-radius diagram of toroidal field).
210 """
211 fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 8))
212
213 # Butterfly diagram
214 T, R = np.meshgrid(self.time_history, self.r)
215 im = ax1.contourf(T, R, self.B_history.T, levels=20, cmap='RdBu_r')
216 ax1.set_xlabel('Time')
217 ax1.set_ylabel('Radius r')
218 ax1.set_title(f'Butterfly Diagram (Toroidal Field B)\nD = {self.dynamo_number:.2f}')
219 plt.colorbar(im, ax=ax1, label='B')
220
221 # Dynamo wave propagation
222 ax1.text(0.02, 0.98,
223 'Dynamo wave propagates\ndue to α-Ω coupling',
224 transform=ax1.transAxes,
225 verticalalignment='top',
226 bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.8))
227
228 # Time series at mid-radius
229 mid_idx = self.nr // 2
230 ax2.plot(self.time_history, self.B_history[:, mid_idx], 'b-', linewidth=2)
231 ax2.set_xlabel('Time')
232 ax2.set_ylabel('B')
233 ax2.set_title(f'Toroidal Field at r = {self.r[mid_idx]:.2f}')
234 ax2.grid(True, alpha=0.3)
235 ax2.axhline(y=0, color='k', linestyle='--', alpha=0.5)
236
237 plt.tight_layout()
238 return fig
239
240 def plot_field_profiles(self):
241 """
242 Plot current field profiles.
243 """
244 fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))
245
246 ax1.plot(self.r, self.A, 'b-', linewidth=2, label='A (poloidal)')
247 ax1.set_xlabel('Radius r')
248 ax1.set_ylabel('A')
249 ax1.set_title('Poloidal Field Potential')
250 ax1.grid(True, alpha=0.3)
251 ax1.legend()
252
253 ax2.plot(self.r, self.B, 'r-', linewidth=2, label='B (toroidal)')
254 ax2.set_xlabel('Radius r')
255 ax2.set_ylabel('B')
256 ax2.set_title('Toroidal Field')
257 ax2.grid(True, alpha=0.3)
258 ax2.legend()
259
260 plt.tight_layout()
261 return fig
262
263
264def main():
265 """
266 Main function demonstrating the α-Ω dynamo model.
267 """
268 # Create dynamo instance
269 dynamo = AlphaOmegaDynamo(
270 nr=100,
271 r_min=0.5,
272 r_max=1.0,
273 alpha=1.0,
274 eta=0.05,
275 shear=15.0
276 )
277
278 # Set initial condition
279 dynamo.set_initial_condition(mode='random')
280
281 # Run simulation
282 dynamo.run_simulation(t_end=20.0, dt=0.005, save_interval=20)
283
284 # Plot results
285 dynamo.plot_butterfly_diagram()
286 dynamo.plot_field_profiles()
287
288 plt.savefig('/tmp/alpha_omega_dynamo.png', dpi=150, bbox_inches='tight')
289 print("\nPlot saved to /tmp/alpha_omega_dynamo.png")
290
291 plt.show()
292
293
294if __name__ == "__main__":
295 main()