1#!/usr/bin/env python3
2"""
3Turbulent Small-Scale Dynamo
4
5This module simulates a small-scale (turbulent) dynamo where magnetic fields
6are amplified by turbulent motions at scales smaller than the energy injection
7scale.
8
9Key features:
10- Kazantsev spectrum: E_B(k) ∝ k^{3/2} during kinematic phase
11- Exponential growth phase followed by saturation
12- Stochastic flow model (Ornstein-Uhlenbeck process)
13- Magnetic energy evolution and spectral analysis
14
15The small-scale dynamo is important in astrophysical contexts such as
16the interstellar medium and early universe.
17
18Author: MHD Course Examples
19License: MIT
20"""
21
22import numpy as np
23import matplotlib.pyplot as plt
24from scipy.integrate import odeint
25
26
27class TurbulentDynamo:
28 """
29 Small-scale turbulent dynamo model.
30
31 This uses a simplified stochastic differential equation approach
32 to model magnetic field amplification by turbulent flows.
33
34 Attributes:
35 n_modes (int): Number of Fourier modes
36 k_modes (ndarray): Wavenumber array
37 Re (float): Reynolds number
38 Rm (float): Magnetic Reynolds number
39 Pm (float): Magnetic Prandtl number (Rm/Re)
40 """
41
42 def __init__(self, n_modes=50, k_min=1.0, k_max=50.0,
43 Re=1000.0, Pm=1.0):
44 """
45 Initialize turbulent dynamo.
46
47 Parameters:
48 n_modes (int): Number of Fourier modes
49 k_min (float): Minimum wavenumber
50 k_max (float): Maximum wavenumber
51 Re (float): Reynolds number
52 Pm (float): Magnetic Prandtl number
53 """
54 self.n_modes = n_modes
55 self.k_modes = np.logspace(np.log10(k_min), np.log10(k_max), n_modes)
56 self.Re = Re
57 self.Pm = Pm
58 self.Rm = Re * Pm
59
60 # Turbulent velocity amplitude (from Kolmogorov)
61 self.u_rms = 1.0
62 self.eta = self.u_rms / self.Rm # Magnetic diffusivity
63 self.nu = self.u_rms / self.Re # Kinematic viscosity
64
65 # Eddy turnover time
66 self.tau_eddy = 1.0
67
68 # Initialize magnetic energy spectrum
69 self.E_B = 1e-6 * np.ones(n_modes) # Seed field
70
71 # Time history
72 self.time_history = []
73 self.energy_history = []
74 self.spectrum_history = []
75
76 def kolmogorov_spectrum(self, k):
77 """
78 Kolmogorov turbulent energy spectrum E_v(k) ∝ k^{-5/3}.
79
80 Parameters:
81 k (ndarray): Wavenumber
82
83 Returns:
84 ndarray: Velocity energy spectrum
85 """
86 k0 = 1.0 # Energy injection scale
87 epsilon = 1.0 # Energy dissipation rate
88 C_K = 1.5 # Kolmogorov constant
89
90 E_v = C_K * epsilon**(2/3) * k**(-5/3)
91
92 # Add exponential cutoff at dissipation scale
93 k_d = (epsilon / self.nu**3)**(1/4) # Kolmogorov scale
94 E_v *= np.exp(-k / k_d)
95
96 return E_v
97
98 def kazantsev_spectrum(self, k):
99 """
100 Kazantsev magnetic energy spectrum E_B(k) ∝ k^{3/2}.
101
102 This is the theoretical prediction for small-scale dynamo
103 in the kinematic regime.
104
105 Parameters:
106 k (ndarray): Wavenumber
107
108 Returns:
109 ndarray: Magnetic energy spectrum
110 """
111 # Kazantsev spectrum
112 k0 = self.k_modes[0]
113 C_K = 0.1
114
115 E_B = C_K * (k / k0)**(3/2)
116
117 # Cutoff at resistive scale
118 k_eta = np.sqrt(self.u_rms / self.eta)
119 E_B *= np.exp(-k / k_eta)
120
121 return E_B
122
123 def growth_rate(self, k):
124 """
125 Compute growth rate γ(k) for each mode.
126
127 In kinematic regime: γ(k) ≈ constant (Lyapunov exponent)
128 With diffusion: γ(k) = γ₀ - η k²
129
130 Parameters:
131 k (ndarray): Wavenumber
132
133 Returns:
134 ndarray: Growth rate
135 """
136 # Lyapunov exponent (stretching rate)
137 gamma_0 = 0.1 * self.u_rms / self.tau_eddy
138
139 # Diffusive damping
140 gamma = gamma_0 - self.eta * k**2
141
142 return gamma
143
144 def rhs(self, E_B, t):
145 """
146 Right-hand side for magnetic energy evolution.
147
148 dE_B/dt = 2γ(k)E_B - saturation term
149
150 Parameters:
151 E_B (ndarray): Magnetic energy spectrum
152 t (float): Time
153
154 Returns:
155 ndarray: Time derivative
156 """
157 gamma = self.growth_rate(self.k_modes)
158
159 # Velocity energy for saturation
160 E_v = self.kolmogorov_spectrum(self.k_modes)
161
162 # Kinematic growth
163 dE_dt = 2 * gamma * E_B
164
165 # Saturation: when E_B ~ E_v, growth slows
166 # Simple saturation model
167 saturation_factor = E_v / (E_v + E_B)
168 dE_dt *= saturation_factor
169
170 # Add stochastic forcing (simplified)
171 forcing = 1e-4 * np.random.randn(self.n_modes)
172 dE_dt += forcing
173
174 return dE_dt
175
176 def run_simulation(self, t_end=50.0, n_steps=500):
177 """
178 Run the turbulent dynamo simulation.
179
180 Parameters:
181 t_end (float): End time
182 n_steps (int): Number of time steps
183 """
184 t_array = np.linspace(0, t_end, n_steps)
185
186 print("Running turbulent dynamo simulation...")
187 print(f"Magnetic Reynolds number Rm = {self.Rm:.1f}")
188 print(f"Magnetic Prandtl number Pm = {self.Pm:.1f}")
189
190 # Storage
191 self.time_history = []
192 self.energy_history = []
193 self.spectrum_history = []
194
195 # Time stepping (using manual Euler due to stochasticity)
196 dt = t_end / n_steps
197 E_B = self.E_B.copy()
198
199 for i, t in enumerate(t_array):
200 # Store data
201 self.time_history.append(t)
202 self.energy_history.append(np.sum(E_B))
203 self.spectrum_history.append(E_B.copy())
204
205 # Euler step
206 dE_dt = self.rhs(E_B, t)
207 E_B += dt * dE_dt
208
209 # Ensure positivity
210 E_B = np.maximum(E_B, 1e-10)
211
212 if i % 100 == 0:
213 total_energy = np.sum(E_B)
214 print(f"Step {i}/{n_steps}, t={t:.2f}, E_mag={total_energy:.4e}")
215
216 self.E_B = E_B
217 self.spectrum_history = np.array(self.spectrum_history)
218 print("Simulation complete!")
219
220 def plot_energy_evolution(self):
221 """
222 Plot magnetic energy evolution.
223 """
224 fig, ax = plt.subplots(figsize=(10, 6))
225
226 # Compute growth rate
227 growth_rate_theory = self.growth_rate(self.k_modes[0])
228
229 ax.semilogy(self.time_history, self.energy_history, 'b-',
230 linewidth=2, label='Total magnetic energy')
231
232 # Theoretical exponential growth
233 if len(self.time_history) > 0:
234 E0 = self.energy_history[0]
235 t_array = np.array(self.time_history)
236 E_theory = E0 * np.exp(2 * 0.1 * t_array / self.tau_eddy)
237 ax.semilogy(t_array, E_theory, 'r--', alpha=0.6,
238 label='Exponential growth (kinematic)')
239
240 ax.set_xlabel('Time', fontsize=12)
241 ax.set_ylabel('Total Magnetic Energy', fontsize=12)
242 ax.set_title('Turbulent Dynamo: Magnetic Energy Growth', fontsize=14)
243 ax.grid(True, alpha=0.3)
244 ax.legend(fontsize=11)
245
246 # Mark phases
247 ax.text(0.5, 0.95, 'Kinematic phase → Saturation',
248 transform=ax.transAxes, fontsize=11,
249 verticalalignment='top',
250 bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.8))
251
252 plt.tight_layout()
253 return fig
254
255 def plot_spectrum_evolution(self):
256 """
257 Plot magnetic energy spectrum evolution.
258 """
259 fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))
260
261 # Spectrum at different times
262 n_times = min(5, len(self.spectrum_history))
263 time_indices = np.linspace(0, len(self.spectrum_history)-1,
264 n_times, dtype=int)
265
266 for idx in time_indices:
267 t = self.time_history[idx]
268 spectrum = self.spectrum_history[idx]
269 ax1.loglog(self.k_modes, spectrum, '-', linewidth=2,
270 label=f't = {t:.1f}')
271
272 # Theoretical spectra
273 kazantsev = self.kazantsev_spectrum(self.k_modes)
274 ax1.loglog(self.k_modes, kazantsev * np.max(self.spectrum_history[-1]) /
275 np.max(kazantsev), 'k--', linewidth=2,
276 alpha=0.6, label=r'Kazantsev $k^{3/2}$')
277
278 ax1.set_xlabel('Wavenumber k', fontsize=12)
279 ax1.set_ylabel(r'$E_B(k)$', fontsize=12)
280 ax1.set_title('Magnetic Energy Spectrum Evolution', fontsize=14)
281 ax1.grid(True, alpha=0.3)
282 ax1.legend(fontsize=10)
283
284 # Final spectrum with comparison
285 E_v = self.kolmogorov_spectrum(self.k_modes)
286 E_v *= np.max(self.E_B) / np.max(E_v) # Normalize for comparison
287
288 ax2.loglog(self.k_modes, self.E_B, 'b-', linewidth=2,
289 label='Magnetic energy')
290 ax2.loglog(self.k_modes, E_v, 'r--', linewidth=2, alpha=0.6,
291 label=r'Velocity energy $k^{-5/3}$')
292 ax2.set_xlabel('Wavenumber k', fontsize=12)
293 ax2.set_ylabel('Energy', fontsize=12)
294 ax2.set_title('Final State: Magnetic vs Kinetic Energy', fontsize=14)
295 ax2.grid(True, alpha=0.3)
296 ax2.legend(fontsize=11)
297
298 plt.tight_layout()
299 return fig
300
301
302def main():
303 """
304 Main function demonstrating turbulent dynamo.
305 """
306 print("=" * 60)
307 print("Turbulent Small-Scale Dynamo Simulation")
308 print("=" * 60)
309
310 # Create dynamo
311 dynamo = TurbulentDynamo(
312 n_modes=50,
313 k_min=1.0,
314 k_max=100.0,
315 Re=1000.0,
316 Pm=1.0
317 )
318
319 # Run simulation
320 dynamo.run_simulation(t_end=50.0, n_steps=1000)
321
322 # Plot results
323 dynamo.plot_energy_evolution()
324 dynamo.plot_spectrum_evolution()
325
326 plt.savefig('/tmp/turbulent_dynamo.png', dpi=150, bbox_inches='tight')
327 print("\nPlots saved to /tmp/turbulent_dynamo.png")
328
329 plt.show()
330
331
332if __name__ == "__main__":
333 main()