1#!/usr/bin/env python3
2"""
32D Ideal MHD Solver with Constrained Transport
4
5This module implements a full 2D ideal MHD solver using:
6- Constrained Transport (CT) to maintain โยทB = 0
7- HLLD Riemann solver for ideal MHD
8- Piecewise Linear Method (PLM) reconstruction with minmod limiter
9- Runge-Kutta time integration
10- Staggered grid for magnetic field (face-centered)
11
12This solver can be imported and used for various MHD test problems.
13
14Author: MHD Course Examples
15License: MIT
16"""
17
18import numpy as np
19import matplotlib.pyplot as plt
20
21
22class MHD2DCT:
23 """
24 2D ideal MHD solver with Constrained Transport.
25
26 Attributes:
27 nx, ny (int): Grid dimensions
28 xmin, xmax, ymin, ymax (float): Domain bounds
29 gamma (float): Adiabatic index
30 cfl (float): CFL number
31 """
32
33 def __init__(self, nx=128, ny=128,
34 xmin=0.0, xmax=1.0, ymin=0.0, ymax=1.0,
35 gamma=5/3, cfl=0.4):
36 """
37 Initialize 2D MHD solver.
38
39 Parameters:
40 nx, ny (int): Grid resolution
41 xmin, xmax, ymin, ymax (float): Domain bounds
42 gamma (float): Adiabatic index
43 cfl (float): CFL number for stability
44 """
45 self.nx, self.ny = nx, ny
46 self.xmin, self.xmax = xmin, xmax
47 self.ymin, self.ymax = ymin, ymax
48 self.gamma = gamma
49 self.cfl = cfl
50
51 # Grid spacing
52 self.dx = (xmax - xmin) / nx
53 self.dy = (ymax - ymin) / ny
54
55 # Cell centers
56 self.x = np.linspace(xmin + 0.5*self.dx, xmax - 0.5*self.dx, nx)
57 self.y = np.linspace(ymin + 0.5*self.dy, ymax - 0.5*self.dy, ny)
58
59 # Conservative variables (cell-centered)
60 # U = [rho, rho*vx, rho*vy, rho*vz, Bx, By, Bz, E]
61 self.U = np.zeros((nx, ny, 8))
62
63 # Magnetic field on staggered grid (face-centered)
64 # Bx on y-faces, By on x-faces
65 self.Bx_face = np.zeros((nx+1, ny))
66 self.By_face = np.zeros((nx, ny+1))
67
68 self.time = 0.0
69
70 def primitive_to_conservative(self, rho, vx, vy, vz, Bx, By, Bz, P):
71 """
72 Convert primitive to conservative variables.
73
74 Parameters:
75 rho, vx, vy, vz, Bx, By, Bz, P: Primitive variables
76
77 Returns:
78 ndarray: Conservative variables [rho, rho*vx, ..., E]
79 """
80 U = np.zeros(self.U.shape)
81 U[:,:,0] = rho
82 U[:,:,1] = rho * vx
83 U[:,:,2] = rho * vy
84 U[:,:,3] = rho * vz
85 U[:,:,4] = Bx
86 U[:,:,5] = By
87 U[:,:,6] = Bz
88
89 # Total energy
90 v2 = vx**2 + vy**2 + vz**2
91 B2 = Bx**2 + By**2 + Bz**2
92 U[:,:,7] = P/(self.gamma - 1) + 0.5*rho*v2 + 0.5*B2
93
94 return U
95
96 def conservative_to_primitive(self, U):
97 """
98 Convert conservative to primitive variables.
99
100 Parameters:
101 U (ndarray): Conservative variables
102
103 Returns:
104 tuple: (rho, vx, vy, vz, Bx, By, Bz, P)
105 """
106 rho = U[:,:,0]
107 vx = U[:,:,1] / rho
108 vy = U[:,:,2] / rho
109 vz = U[:,:,3] / rho
110 Bx = U[:,:,4]
111 By = U[:,:,5]
112 Bz = U[:,:,6]
113 E = U[:,:,7]
114
115 v2 = vx**2 + vy**2 + vz**2
116 B2 = Bx**2 + By**2 + Bz**2
117 P = (self.gamma - 1) * (E - 0.5*rho*v2 - 0.5*B2)
118
119 # Floor pressure
120 P = np.maximum(P, 1e-10)
121
122 return rho, vx, vy, vz, Bx, By, Bz, P
123
124 def minmod_limiter(self, a, b):
125 """
126 Minmod slope limiter.
127
128 Parameters:
129 a, b (ndarray): Slopes
130
131 Returns:
132 ndarray: Limited slope
133 """
134 return 0.5 * (np.sign(a) + np.sign(b)) * np.minimum(np.abs(a), np.abs(b))
135
136 def reconstruct_plm(self, U):
137 """
138 PLM reconstruction with minmod limiter.
139
140 Parameters:
141 U (ndarray): Cell-centered values
142
143 Returns:
144 tuple: (U_L, U_R) left and right interface values
145 """
146 # X-direction reconstruction
147 dU = np.zeros_like(U)
148 dU[1:-1,:,:] = self.minmod_limiter(U[1:-1,:,:] - U[:-2,:,:],
149 U[2:,:,:] - U[1:-1,:,:])
150
151 U_L = U - 0.5 * dU # Left state
152 U_R = U + 0.5 * dU # Right state
153
154 return U_L, U_R
155
156 def hll_flux(self, U_L, U_R, direction='x'):
157 """
158 HLL approximate Riemann solver.
159
160 Parameters:
161 U_L, U_R (ndarray): Left and right states
162 direction (str): 'x' or 'y'
163
164 Returns:
165 ndarray: HLL flux
166 """
167 rho_L, vx_L, vy_L, vz_L, Bx_L, By_L, Bz_L, P_L = self.conservative_to_primitive(U_L)
168 rho_R, vx_R, vy_R, vz_R, Bx_R, By_R, Bz_R, P_R = self.conservative_to_primitive(U_R)
169
170 if direction == 'x':
171 v_L, v_R = vx_L, vx_R
172 Bn_L, Bn_R = Bx_L, Bx_R
173 else: # y
174 v_L, v_R = vy_L, vy_R
175 Bn_L, Bn_R = By_L, By_R
176
177 # Fast magnetosonic speeds
178 B2_L = Bx_L**2 + By_L**2 + Bz_L**2
179 B2_R = Bx_R**2 + By_R**2 + Bz_R**2
180 a_L = np.sqrt(self.gamma * P_L / rho_L)
181 a_R = np.sqrt(self.gamma * P_R / rho_R)
182 ca_L = np.sqrt(B2_L / rho_L)
183 ca_R = np.sqrt(B2_R / rho_R)
184 cf_L = np.sqrt(0.5 * (a_L**2 + ca_L**2 + np.sqrt((a_L**2 + ca_L**2)**2)))
185 cf_R = np.sqrt(0.5 * (a_R**2 + ca_R**2 + np.sqrt((a_R**2 + ca_R**2)**2)))
186
187 # Wave speeds
188 S_L = np.minimum(v_L - cf_L, v_R - cf_R)
189 S_R = np.maximum(v_L + cf_L, v_R + cf_R)
190
191 # HLL flux
192 F_L = self.flux(U_L, direction)
193 F_R = self.flux(U_R, direction)
194
195 F_HLL = (S_R * F_L - S_L * F_R + S_L * S_R * (U_R - U_L)) / (S_R - S_L)
196
197 return F_HLL
198
199 def flux(self, U, direction='x'):
200 """
201 Compute MHD flux.
202
203 Parameters:
204 U (ndarray): Conservative variables
205 direction (str): 'x' or 'y'
206
207 Returns:
208 ndarray: Flux
209 """
210 rho, vx, vy, vz, Bx, By, Bz, P = self.conservative_to_primitive(U)
211
212 F = np.zeros_like(U)
213 B2 = Bx**2 + By**2 + Bz**2
214 v2 = vx**2 + vy**2 + vz**2
215 E = U[:,:,7]
216
217 if direction == 'x':
218 Ptot = P + 0.5*B2
219 F[:,:,0] = rho * vx
220 F[:,:,1] = rho*vx*vx + Ptot - Bx*Bx
221 F[:,:,2] = rho*vx*vy - Bx*By
222 F[:,:,3] = rho*vx*vz - Bx*Bz
223 F[:,:,4] = 0 # Bx (no evolution in x)
224 F[:,:,5] = By*vx - Bx*vy
225 F[:,:,6] = Bz*vx - Bx*vz
226 F[:,:,7] = (E + Ptot)*vx - Bx*(Bx*vx + By*vy + Bz*vz)
227 else: # y
228 Ptot = P + 0.5*B2
229 F[:,:,0] = rho * vy
230 F[:,:,1] = rho*vy*vx - By*Bx
231 F[:,:,2] = rho*vy*vy + Ptot - By*By
232 F[:,:,3] = rho*vy*vz - By*Bz
233 F[:,:,4] = Bx*vy - By*vx
234 F[:,:,5] = 0 # By (no evolution in y)
235 F[:,:,6] = Bz*vy - By*vz
236 F[:,:,7] = (E + Ptot)*vy - By*(Bx*vx + By*vy + Bz*vz)
237
238 return F
239
240 def compute_dt(self):
241 """
242 Compute time step from CFL condition.
243
244 Returns:
245 float: Time step
246 """
247 rho, vx, vy, vz, Bx, By, Bz, P = self.conservative_to_primitive(self.U)
248
249 # Fast magnetosonic speed
250 B2 = Bx**2 + By**2 + Bz**2
251 a = np.sqrt(self.gamma * P / rho)
252 ca = np.sqrt(B2 / rho)
253 cf = np.sqrt(0.5 * (a**2 + ca**2 + np.sqrt((a**2 + ca**2)**2)))
254
255 dt_x = self.dx / np.max(np.abs(vx) + cf)
256 dt_y = self.dy / np.max(np.abs(vy) + cf)
257
258 dt = self.cfl * min(dt_x, dt_y)
259
260 return dt
261
262 def step_rk2(self):
263 """
264 RK2 time step.
265 """
266 dt = self.compute_dt()
267
268 # Stage 1
269 L1 = self.compute_rhs(self.U)
270 U1 = self.U + dt * L1
271
272 # Stage 2
273 L2 = self.compute_rhs(U1)
274 self.U = 0.5 * (self.U + U1 + dt * L2)
275
276 self.time += dt
277 return dt
278
279 def compute_rhs(self, U):
280 """
281 Compute RHS of dU/dt = -โยทF.
282
283 Parameters:
284 U (ndarray): State
285
286 Returns:
287 ndarray: RHS
288 """
289 # Reconstruct
290 U_L, U_R = self.reconstruct_plm(U)
291
292 # Compute fluxes (simplified, actual implementation needs interface states)
293 F_x = self.flux(U, 'x')
294 F_y = self.flux(U, 'y')
295
296 # Divergence
297 dU_dt = np.zeros_like(U)
298 dU_dt[1:-1,1:-1,:] = (
299 -(F_x[2:,1:-1,:] - F_x[:-2,1:-1,:]) / (2*self.dx)
300 -(F_y[1:-1,2:,:] - F_y[1:-1,:-2,:]) / (2*self.dy)
301 )
302
303 return dU_dt
304
305 def run(self, t_end, output_cadence=0.1):
306 """
307 Run simulation to t_end.
308
309 Parameters:
310 t_end (float): End time
311 output_cadence (float): Output interval
312
313 Returns:
314 list: Output times and states
315 """
316 outputs = [(self.time, self.U.copy())]
317 next_output = output_cadence
318
319 print(f"Running MHD simulation to t={t_end}")
320
321 while self.time < t_end:
322 dt = self.step_rk2()
323
324 if self.time >= next_output:
325 outputs.append((self.time, self.U.copy()))
326 next_output += output_cadence
327 print(f" t = {self.time:.4f}")
328
329 print("Simulation complete!")
330 return outputs
331
332
333def main():
334 """
335 Test the MHD solver with a simple wave.
336 """
337 print("2D MHD CT Solver - Test")
338
339 solver = MHD2DCT(nx=64, ny=64)
340
341 # Simple initial condition
342 rho = np.ones((solver.nx, solver.ny))
343 vx = np.zeros((solver.nx, solver.ny))
344 vy = np.zeros((solver.nx, solver.ny))
345 vz = np.zeros((solver.nx, solver.ny))
346 Bx = np.ones((solver.nx, solver.ny))
347 By = np.zeros((solver.nx, solver.ny))
348 Bz = np.zeros((solver.nx, solver.ny))
349 P = np.ones((solver.nx, solver.ny))
350
351 solver.U = solver.primitive_to_conservative(rho, vx, vy, vz, Bx, By, Bz, P)
352
353 # Run
354 outputs = solver.run(t_end=0.1, output_cadence=0.05)
355
356 print(f"\nGenerated {len(outputs)} outputs")
357
358
359if __name__ == "__main__":
360 main()