1#!/usr/bin/env python3
2"""
3Edge Localized Mode (ELM) Crash - Simplified Model
4
5This module implements a simplified model of Edge Localized Modes (ELMs)
6in tokamak H-mode plasmas.
7
8Key physics:
9- ELMs are periodic instabilities at the plasma edge (pedestal region)
10- Driven by peeling-ballooning modes when pressure gradient or current
11 exceeds stability boundary
12- Lead to rapid loss of edge energy and particles
13- ELM cycle: buildup → crash → recovery
14
15Peeling-ballooning stability in (J_edge, α) space where:
16 J_edge = edge current density
17 α = pressure gradient parameter
18
19Author: MHD Course Examples
20License: MIT
21"""
22
23import numpy as np
24import matplotlib.pyplot as plt
25from scipy.integrate import odeint
26
27
28class ELMModel:
29 """
30 Simplified ELM crash model with peeling-ballooning stability.
31
32 Attributes:
33 p_edge (float): Edge pedestal pressure
34 j_edge (float): Edge current density
35 alpha (float): Pressure gradient parameter
36 tau_buildup (float): Pedestal buildup timescale
37 """
38
39 def __init__(self, tau_buildup=0.1, tau_crash=0.001,
40 p_source=1.0, alpha_crit=1.0, j_crit=1.0):
41 """
42 Initialize ELM model.
43
44 Parameters:
45 tau_buildup (float): Buildup time (s)
46 tau_crash (float): Crash time (s)
47 p_source (float): Heating power (normalized)
48 alpha_crit (float): Critical pressure gradient
49 j_crit (float): Critical current density
50 """
51 self.tau_buildup = tau_buildup
52 self.tau_crash = tau_crash
53 self.p_source = p_source
54 self.alpha_crit = alpha_crit
55 self.j_crit = j_crit
56
57 # State
58 self.p_ped = 0.1 # Pedestal height
59 self.w_ped = 0.05 # Pedestal width
60
61 # Derived quantities
62 self.alpha = 0.0
63 self.j_edge = 0.0
64
65 # History
66 self.time_history = []
67 self.p_history = []
68 self.alpha_history = []
69 self.j_history = []
70 self.elm_times = []
71 self.elm_energy = []
72
73 def compute_stability_parameters(self):
74 """
75 Compute peeling-ballooning parameters from pedestal.
76
77 Returns:
78 tuple: (alpha, j_edge)
79 """
80 # Pressure gradient parameter
81 # α ~ ∇p in pedestal region
82 alpha = self.p_ped / self.w_ped
83
84 # Edge current (bootstrap + driven)
85 # Simplified: j_edge ~ sqrt(p_ped)
86 j_edge = np.sqrt(self.p_ped)
87
88 return alpha, j_edge
89
90 def peeling_ballooning_boundary(self, j_array):
91 """
92 Peeling-ballooning stability boundary in (j, α) space.
93
94 Simplified boundary:
95 - Low j: ballooning limit α_max ~ α_crit
96 - High j: peeling limit α_max ~ α_crit - (j/j_crit)²
97
98 Parameters:
99 j_array (ndarray): Current density array
100
101 Returns:
102 ndarray: Critical alpha boundary
103 """
104 # Combined peeling-ballooning boundary
105 j_norm = j_array / self.j_crit
106
107 # Ballooning branch (low j)
108 alpha_ballooning = self.alpha_crit * np.ones_like(j_array)
109
110 # Peeling branch (high j)
111 alpha_peeling = self.alpha_crit * (1 - j_norm**2)
112
113 # Take minimum (most restrictive)
114 alpha_boundary = np.minimum(alpha_ballooning, alpha_peeling)
115
116 return np.maximum(alpha_boundary, 0)
117
118 def is_elm_unstable(self):
119 """
120 Check if ELM is unstable (above stability boundary).
121
122 Returns:
123 bool: True if unstable
124 """
125 alpha, j_edge = self.compute_stability_parameters()
126
127 # Check against boundary
128 alpha_max = self.peeling_ballooning_boundary(np.array([j_edge]))[0]
129
130 unstable = alpha > alpha_max
131
132 return unstable
133
134 def elm_crash(self):
135 """
136 Execute ELM crash: rapid loss of pedestal.
137
138 Returns:
139 float: Energy lost in ELM
140 """
141 # Energy loss proportional to pedestal height
142 energy_loss = self.p_ped * 0.5 # Lose ~50% of pedestal
143
144 # Reduce pedestal height
145 self.p_ped *= 0.3 # Crash to 30% of pre-ELM value
146
147 # Widen pedestal
148 self.w_ped *= 1.5
149
150 return energy_loss
151
152 def evolve_pedestal(self, dt):
153 """
154 Evolve pedestal during inter-ELM period.
155
156 dp_ped/dt = P_source - p_ped/τ_buildup
157
158 Parameters:
159 dt (float): Time step
160 """
161 # Buildup equation
162 dp_dt = self.p_source - self.p_ped / self.tau_buildup
163
164 # Narrow pedestal during buildup
165 dw_dt = -0.1 * self.w_ped / self.tau_buildup
166
167 # Euler step
168 self.p_ped += dt * dp_dt
169 self.w_ped += dt * dw_dt
170
171 # Minimum width
172 self.w_ped = max(self.w_ped, 0.02)
173
174 def run_elm_cycle(self, t_end=2.0, dt=0.0001):
175 """
176 Run ELM cycle simulation.
177
178 Parameters:
179 t_end (float): End time (s)
180 dt (float): Time step (s)
181 """
182 n_steps = int(t_end / dt)
183 t = 0.0
184
185 print("Running ELM cycle simulation...")
186 print(f"Critical α = {self.alpha_crit:.2f}")
187 print(f"Critical j = {self.j_crit:.2f}")
188
189 self.time_history = []
190 self.p_history = []
191 self.alpha_history = []
192 self.j_history = []
193 self.elm_times = []
194 self.elm_energy = []
195
196 for step in range(n_steps):
197 # Store state
198 alpha, j_edge = self.compute_stability_parameters()
199 self.time_history.append(t)
200 self.p_history.append(self.p_ped)
201 self.alpha_history.append(alpha)
202 self.j_history.append(j_edge)
203
204 # Check ELM stability
205 if self.is_elm_unstable():
206 # ELM crash
207 energy = self.elm_crash()
208 self.elm_times.append(t)
209 self.elm_energy.append(energy)
210 print(f"ELM at t={t:.4f} s, ΔW={energy:.3f}")
211 else:
212 # Inter-ELM buildup
213 self.evolve_pedestal(dt)
214
215 t += dt
216
217 print(f"\nSimulation complete!")
218 print(f"Total ELMs: {len(self.elm_times)}")
219 if len(self.elm_times) > 1:
220 elm_freq = 1.0 / np.mean(np.diff(self.elm_times))
221 print(f"ELM frequency: {elm_freq:.1f} Hz")
222
223 def plot_elm_cycle(self):
224 """
225 Plot ELM cycle evolution.
226 """
227 fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(14, 10))
228
229 t = np.array(self.time_history)
230
231 # Pedestal height
232 ax1.plot(t, self.p_history, 'b-', linewidth=1.5)
233 for t_elm in self.elm_times:
234 ax1.axvline(x=t_elm, color='r', alpha=0.3, linewidth=0.8)
235 ax1.set_xlabel('Time (s)', fontsize=11)
236 ax1.set_ylabel('Pedestal Height', fontsize=11)
237 ax1.set_title('Pedestal Evolution (ELM Crashes)', fontsize=13)
238 ax1.grid(True, alpha=0.3)
239
240 # Stability parameters
241 ax2.plot(t, self.alpha_history, 'purple', linewidth=1.5, label='α')
242 ax2.axhline(y=self.alpha_crit, color='r', linestyle='--',
243 linewidth=2, label='α_crit')
244 for t_elm in self.elm_times:
245 ax2.axvline(x=t_elm, color='r', alpha=0.3, linewidth=0.8)
246 ax2.set_xlabel('Time (s)', fontsize=11)
247 ax2.set_ylabel('Pressure Gradient α', fontsize=11)
248 ax2.set_title('Stability Parameter Evolution', fontsize=13)
249 ax2.grid(True, alpha=0.3)
250 ax2.legend(fontsize=10)
251
252 # Peeling-ballooning diagram
253 j_array = np.linspace(0, 1.5 * self.j_crit, 200)
254 alpha_boundary = self.peeling_ballooning_boundary(j_array)
255
256 ax3.plot(j_array, alpha_boundary, 'r-', linewidth=2.5,
257 label='Stability boundary')
258 ax3.fill_between(j_array, 0, alpha_boundary, alpha=0.3,
259 color='green', label='Stable')
260 ax3.fill_between(j_array, alpha_boundary, 2*self.alpha_crit,
261 alpha=0.3, color='red', label='Unstable (ELM)')
262
263 # Plot trajectory
264 ax3.plot(self.j_history, self.alpha_history, 'b-',
265 linewidth=1, alpha=0.5, label='Trajectory')
266
267 # Mark ELM times
268 elm_indices = [np.argmin(np.abs(t - t_elm)) for t_elm in self.elm_times]
269 j_elm = [self.j_history[i] for i in elm_indices]
270 alpha_elm = [self.alpha_history[i] for i in elm_indices]
271 ax3.plot(j_elm, alpha_elm, 'ro', markersize=8, label='ELM crashes')
272
273 ax3.set_xlabel('Edge Current j_edge', fontsize=11)
274 ax3.set_ylabel('Pressure Gradient α', fontsize=11)
275 ax3.set_title('Peeling-Ballooning Stability Diagram', fontsize=13)
276 ax3.set_xlim([0, 1.5 * self.j_crit])
277 ax3.set_ylim([0, 2 * self.alpha_crit])
278 ax3.grid(True, alpha=0.3)
279 ax3.legend(fontsize=10)
280
281 # ELM energy vs pedestal height
282 if len(self.elm_energy) > 0:
283 # Find pedestal height just before each ELM
284 p_before_elm = []
285 for t_elm in self.elm_times:
286 idx = np.argmin(np.abs(t - t_elm))
287 if idx > 0:
288 p_before_elm.append(self.p_history[idx-1])
289
290 ax4.plot(p_before_elm, self.elm_energy, 'go', markersize=8)
291
292 # Fit line
293 if len(p_before_elm) > 1:
294 coeffs = np.polyfit(p_before_elm, self.elm_energy, 1)
295 p_fit = np.linspace(min(p_before_elm), max(p_before_elm), 100)
296 E_fit = np.polyval(coeffs, p_fit)
297 ax4.plot(p_fit, E_fit, 'r--', linewidth=2,
298 label=f'Fit: ΔW = {coeffs[0]:.2f}p + {coeffs[1]:.3f}')
299
300 ax4.set_xlabel('Pedestal Height (pre-ELM)', fontsize=11)
301 ax4.set_ylabel('ELM Energy Loss ΔW', fontsize=11)
302 ax4.set_title('ELM Energy vs Pedestal Height', fontsize=13)
303 ax4.grid(True, alpha=0.3)
304 ax4.legend(fontsize=10)
305 else:
306 ax4.text(0.5, 0.5, 'No ELMs occurred',
307 transform=ax4.transAxes, ha='center', va='center',
308 fontsize=14)
309
310 plt.tight_layout()
311 return fig
312
313
314def main():
315 """
316 Main function demonstrating ELM crash model.
317 """
318 print("=" * 60)
319 print("Edge Localized Mode (ELM) Simplified Model")
320 print("=" * 60)
321
322 # Create ELM model
323 elm = ELMModel(
324 tau_buildup=0.1, # 100 ms buildup
325 tau_crash=0.001, # 1 ms crash
326 p_source=1.5, # Heating power
327 alpha_crit=1.0, # Critical gradient
328 j_crit=0.8 # Critical current
329 )
330
331 # Run simulation
332 elm.run_elm_cycle(t_end=2.0, dt=0.0001)
333
334 # Plot results
335 elm.plot_elm_cycle()
336
337 plt.savefig('/tmp/elm_crash.png', dpi=150, bbox_inches='tight')
338 print("\nPlot saved to /tmp/elm_crash.png")
339
340 plt.show()
341
342
343if __name__ == "__main__":
344 main()