1#!/usr/bin/env python3
2"""
3Orszag-Tang Vortex MHD Test Problem
4
5This is a standard benchmark for MHD codes that demonstrates the formation
6of shocks, current sheets, and complex vortical structures from smooth
7initial conditions.
8
9Initial conditions:
10 v = (-sin(y), sin(x), 0)
11 B = (-sin(y), sin(2x), 0)
12 ρ = γ², P = γ
13 Domain: [0, 2π]² with periodic boundaries
14
15The Orszag-Tang vortex tests:
16- Shock capturing
17- Current sheet formation
18- Conservation properties
19- Resolution of small-scale structures
20
21Author: MHD Course Examples
22License: MIT
23"""
24
25import numpy as np
26import matplotlib.pyplot as plt
27
28
29class OrszagTangVortex:
30 """
31 Orszag-Tang vortex test problem.
32
33 This provides initial conditions and analysis for the standard
34 MHD benchmark test.
35 """
36
37 def __init__(self, nx=256, ny=256, gamma=5/3):
38 """
39 Initialize Orszag-Tang vortex.
40
41 Parameters:
42 nx, ny (int): Grid resolution
43 gamma (float): Adiabatic index
44 """
45 self.nx = nx
46 self.ny = ny
47 self.gamma = gamma
48
49 # Domain: [0, 2π]²
50 self.xmin, self.xmax = 0.0, 2*np.pi
51 self.ymin, self.ymax = 0.0, 2*np.pi
52
53 # Grid
54 x = np.linspace(self.xmin, self.xmax, nx, endpoint=False)
55 y = np.linspace(self.ymin, self.ymax, ny, endpoint=False)
56 self.X, self.Y = np.meshgrid(x, y, indexing='ij')
57
58 def initial_conditions(self):
59 """
60 Set up Orszag-Tang initial conditions.
61
62 Returns:
63 dict: Primitive variables {rho, vx, vy, vz, Bx, By, Bz, P}
64 """
65 X, Y = self.X, self.Y
66
67 # Density (constant)
68 rho = self.gamma**2 * np.ones_like(X)
69
70 # Velocity
71 vx = -np.sin(Y)
72 vy = np.sin(X)
73 vz = np.zeros_like(X)
74
75 # Magnetic field
76 Bx = -np.sin(Y)
77 By = np.sin(2*X)
78 Bz = np.zeros_like(X)
79
80 # Pressure (constant)
81 P = self.gamma * np.ones_like(X)
82
83 return {
84 'rho': rho, 'vx': vx, 'vy': vy, 'vz': vz,
85 'Bx': Bx, 'By': By, 'Bz': Bz, 'P': P
86 }
87
88 def compute_diagnostics(self, rho, vx, vy, vz, Bx, By, Bz, P):
89 """
90 Compute diagnostic quantities.
91
92 Parameters:
93 rho, vx, vy, vz, Bx, By, Bz, P: Primitive variables
94
95 Returns:
96 dict: Diagnostic quantities
97 """
98 # Current density: J = ∇×B
99 dBx_dy = np.gradient(Bx, axis=1)
100 dBy_dx = np.gradient(By, axis=0)
101 Jz = dBy_dx - dBx_dy
102
103 # Vorticity: ω = ∇×v
104 dvx_dy = np.gradient(vx, axis=1)
105 dvy_dx = np.gradient(vy, axis=0)
106 omega_z = dvy_dx - dvx_dy
107
108 # Magnetic pressure
109 B_squared = Bx**2 + By**2 + Bz**2
110 P_mag = 0.5 * B_squared
111
112 # Total pressure
113 P_total = P + P_mag
114
115 # Mach number
116 sound_speed = np.sqrt(self.gamma * P / rho)
117 v_mag = np.sqrt(vx**2 + vy**2 + vz**2)
118 mach = v_mag / sound_speed
119
120 # Plasma beta
121 beta = P / P_mag
122
123 return {
124 'current': Jz,
125 'vorticity': omega_z,
126 'P_mag': P_mag,
127 'P_total': P_total,
128 'mach': mach,
129 'beta': beta,
130 'B_mag': np.sqrt(B_squared)
131 }
132
133 def plot_solution(self, rho, vx, vy, vz, Bx, By, Bz, P, time=0.0):
134 """
135 Plot the solution with standard diagnostics.
136
137 Parameters:
138 rho, vx, vy, vz, Bx, By, Bz, P: Primitive variables
139 time (float): Current time
140 """
141 diag = self.compute_diagnostics(rho, vx, vy, vz, Bx, By, Bz, P)
142
143 fig, axes = plt.subplots(2, 3, figsize=(16, 10))
144 ax = axes.flatten()
145
146 # Density
147 im0 = ax[0].contourf(self.X, self.Y, rho, levels=30, cmap='viridis')
148 ax[0].set_title(f'Density (t={time:.3f})')
149 ax[0].set_xlabel('x')
150 ax[0].set_ylabel('y')
151 plt.colorbar(im0, ax=ax[0])
152
153 # Pressure
154 im1 = ax[1].contourf(self.X, self.Y, P, levels=30, cmap='plasma')
155 ax[1].set_title('Gas Pressure')
156 ax[1].set_xlabel('x')
157 ax[1].set_ylabel('y')
158 plt.colorbar(im1, ax=ax[1])
159
160 # Magnetic field magnitude
161 im2 = ax[2].contourf(self.X, self.Y, diag['B_mag'], levels=30,
162 cmap='inferno')
163 ax[2].set_title('|B|')
164 ax[2].set_xlabel('x')
165 ax[2].set_ylabel('y')
166 plt.colorbar(im2, ax=ax[2])
167
168 # Current density
169 im3 = ax[3].contourf(self.X, self.Y, diag['current'], levels=30,
170 cmap='RdBu_r')
171 ax[3].set_title(r'Current Density $J_z$')
172 ax[3].set_xlabel('x')
173 ax[3].set_ylabel('y')
174 plt.colorbar(im3, ax=ax[3])
175
176 # Vorticity
177 im4 = ax[4].contourf(self.X, self.Y, diag['vorticity'], levels=30,
178 cmap='RdBu_r')
179 ax[4].set_title(r'Vorticity $\omega_z$')
180 ax[4].set_xlabel('x')
181 ax[4].set_ylabel('y')
182 plt.colorbar(im4, ax=ax[4])
183
184 # Mach number
185 im5 = ax[5].contourf(self.X, self.Y, diag['mach'], levels=30,
186 cmap='hot')
187 ax[5].set_title('Mach Number')
188 ax[5].set_xlabel('x')
189 ax[5].set_ylabel('y')
190 plt.colorbar(im5, ax=ax[5])
191
192 plt.tight_layout()
193 return fig
194
195
196def main():
197 """
198 Demonstrate Orszag-Tang vortex initial conditions.
199 """
200 print("=" * 60)
201 print("Orszag-Tang Vortex MHD Benchmark")
202 print("=" * 60)
203
204 # Create problem
205 ot = OrszagTangVortex(nx=256, ny=256, gamma=5/3)
206
207 # Get initial conditions
208 ic = ot.initial_conditions()
209
210 print("\nInitial conditions set up on 256×256 grid")
211 print("Domain: [0, 2π]²")
212 print(f"γ = {ot.gamma:.3f}")
213
214 # Compute initial diagnostics
215 diag = ot.compute_diagnostics(**ic)
216
217 print(f"\nInitial state:")
218 print(f" ρ range: [{np.min(ic['rho']):.3f}, {np.max(ic['rho']):.3f}]")
219 print(f" P range: [{np.min(ic['P']):.3f}, {np.max(ic['P']):.3f}]")
220 print(f" |B| range: [{np.min(diag['B_mag']):.3f}, {np.max(diag['B_mag']):.3f}]")
221 print(f" |v| max: {np.max(np.sqrt(ic['vx']**2 + ic['vy']**2)):.3f}")
222
223 # Plot initial conditions
224 fig = ot.plot_solution(**ic, time=0.0)
225 fig.suptitle('Orszag-Tang Vortex: Initial Conditions', fontsize=14,
226 weight='bold', y=0.995)
227
228 plt.savefig('/tmp/orszag_tang_vortex.png', dpi=150, bbox_inches='tight')
229 print("\nPlot saved to /tmp/orszag_tang_vortex.png")
230 print("\nTo run full simulation, use with MHD solver (mhd_2d_ct.py)")
231
232 plt.show()
233
234
235if __name__ == "__main__":
236 main()