GitHub

Orszag-Tang Vortex

The Orszag-Tang vortex is an iconic 2D MHD test problem that demonstrates the transition to MHD turbulence. Starting from smooth initial conditions, the flow develops a complex pattern of interacting shocks and current sheets.

Problem Setup

The domain is $[0, 1]^2$ with periodic boundary conditions and $\gamma = 5/3$. The initial condition is a superposition of velocity and magnetic field vortices:

\[\rho = \gamma^2, \quad P = \gamma, \quad v_x = -\sin(2\pi y), \quad v_y = \sin(2\pi x),\]

\[B_x = -\frac{\sin(2\pi y)}{\sqrt{4\pi}}, \quad B_y = \frac{\sin(4\pi x)}{\sqrt{4\pi}}.\]

using FiniteVolumeMethod
using StaticArrays

gamma = 5.0 / 3.0
eos = IdealGasEOS(gamma)
law = IdealMHDEquations{2}(eos)

N = 50
mesh = StructuredMesh2D(0.0, 1.0, 0.0, 1.0, N, N)

function ot_ic(x, y)
    rho = gamma^2
    P = gamma
    vx = -sin(2 * pi * y)
    vy = sin(2 * pi * x)
    vz = 0.0
    Bx = -sin(2 * pi * y) / sqrt(4 * pi)
    By = sin(4 * pi * x) / sqrt(4 * pi)
    Bz = 0.0
    return SVector(rho, vx, vy, vz, P, Bx, By, Bz)
end
ot_ic (generic function with 1 method)

Vector Potential Initialisation

For MHD with constrained transport, we initialise the magnetic field using a vector potential $A_z(x,y)$ to guarantee $\nabla\cdot\vb B = 0$ to machine precision. The vector potential satisfying $B_x = \partial A_z/\partial y$ and $B_y = -\partial A_z/\partial x$ is:

function Az_ot(x, y)
    return cos(2 * pi * y) / (2 * pi * sqrt(4 * pi)) +
        cos(4 * pi * x) / (4 * pi * sqrt(4 * pi))
end
Az_ot (generic function with 1 method)

Solving

We use the HLLD Riemann solver with MUSCL reconstruction and periodic BCs.

prob = HyperbolicProblem2D(
    law, mesh, HLLDSolver(), CellCenteredMUSCL(MinmodLimiter()),
    PeriodicHyperbolicBC(), PeriodicHyperbolicBC(),
    PeriodicHyperbolicBC(), PeriodicHyperbolicBC(),
    ot_ic; final_time = 0.5, cfl = 0.4
)
HyperbolicProblem2D: IdealMHDEquations with HLLDSolver on 50×50 cells, t ∈ (0.0, 0.5)

The vector_potential keyword tells solve_hyperbolic to initialise the face-centred magnetic field from $A_z$ via Stokes' theorem:

coords, U, t_final, ct = solve_hyperbolic(prob; vector_potential = Az_ot)
50×50 Matrix{Tuple{Float64, Float64}}:
 (0.01, 0.01)  (0.01, 0.03)  (0.01, 0.05)  …  (0.01, 0.97)  (0.01, 0.99)
 (0.03, 0.01)  (0.03, 0.03)  (0.03, 0.05)     (0.03, 0.97)  (0.03, 0.99)
 (0.05, 0.01)  (0.05, 0.03)  (0.05, 0.05)     (0.05, 0.97)  (0.05, 0.99)
 (0.07, 0.01)  (0.07, 0.03)  (0.07, 0.05)     (0.07, 0.97)  (0.07, 0.99)
 (0.09, 0.01)  (0.09, 0.03)  (0.09, 0.05)     (0.09, 0.97)  (0.09, 0.99)
 ⋮                                         ⋱                
 (0.91, 0.01)  (0.91, 0.03)  (0.91, 0.05)  …  (0.91, 0.97)  (0.91, 0.99)
 (0.93, 0.01)  (0.93, 0.03)  (0.93, 0.05)     (0.93, 0.97)  (0.93, 0.99)
 (0.95, 0.01)  (0.95, 0.03)  (0.95, 0.05)     (0.95, 0.97)  (0.95, 0.99)
 (0.97, 0.01)  (0.97, 0.03)  (0.97, 0.05)     (0.97, 0.97)  (0.97, 0.99)
 (0.99, 0.01)  (0.99, 0.03)  (0.99, 0.05)     (0.99, 0.97)  (0.99, 0.99)

Checking $\nabla\cdot\vb B$

The constrained transport algorithm should keep $|\nabla\cdot\vb B|$ at machine precision:

divB_max = max_divB(ct, mesh)
true

Visualisation

using CairoMakie

nx, ny = N, N
xc = [coords[i, 1][1] for i in 1:nx]
yc = [coords[1, j][2] for j in 1:ny]
rho = [conserved_to_primitive(law, U[i, j])[1] for i in 1:nx, j in 1:ny]
P_vals = [conserved_to_primitive(law, U[i, j])[5] for i in 1:nx, j in 1:ny]

fig = Figure(fontsize = 24, size = (1100, 500))
ax1 = Axis(
    fig[1, 1], xlabel = "x", ylabel = "y",
    title = "Density at t = $(round(t_final, digits = 2))", aspect = DataAspect()
)
hm1 = heatmap!(ax1, xc, yc, rho, colormap = :viridis)
Colorbar(fig[1, 2], hm1)

ax2 = Axis(
    fig[1, 3], xlabel = "x", ylabel = "y",
    title = "Pressure", aspect = DataAspect()
)
hm2 = heatmap!(ax2, xc, yc, P_vals, colormap = :magma)
Colorbar(fig[1, 4], hm2)
resize_to_layout!(fig)
fig
Example block output

The density and pressure fields show the characteristic pattern of interacting shocks and current sheets. The maximum $|\nabla\cdot\vb B|$ is round(divB_max, sigdigits = 2), confirming that constrained transport maintains the divergence-free constraint to machine precision.

Just the code

An uncommented version of this example is given below. You can view the source code for this file here.

using FiniteVolumeMethod
using StaticArrays

gamma = 5.0 / 3.0
eos = IdealGasEOS(gamma)
law = IdealMHDEquations{2}(eos)

N = 50
mesh = StructuredMesh2D(0.0, 1.0, 0.0, 1.0, N, N)

function ot_ic(x, y)
    rho = gamma^2
    P = gamma
    vx = -sin(2 * pi * y)
    vy = sin(2 * pi * x)
    vz = 0.0
    Bx = -sin(2 * pi * y) / sqrt(4 * pi)
    By = sin(4 * pi * x) / sqrt(4 * pi)
    Bz = 0.0
    return SVector(rho, vx, vy, vz, P, Bx, By, Bz)
end

function Az_ot(x, y)
    return cos(2 * pi * y) / (2 * pi * sqrt(4 * pi)) +
        cos(4 * pi * x) / (4 * pi * sqrt(4 * pi))
end

prob = HyperbolicProblem2D(
    law, mesh, HLLDSolver(), CellCenteredMUSCL(MinmodLimiter()),
    PeriodicHyperbolicBC(), PeriodicHyperbolicBC(),
    PeriodicHyperbolicBC(), PeriodicHyperbolicBC(),
    ot_ic; final_time = 0.5, cfl = 0.4
)

coords, U, t_final, ct = solve_hyperbolic(prob; vector_potential = Az_ot)

divB_max = max_divB(ct, mesh)


using CairoMakie

nx, ny = N, N
xc = [coords[i, 1][1] for i in 1:nx]
yc = [coords[1, j][2] for j in 1:ny]
rho = [conserved_to_primitive(law, U[i, j])[1] for i in 1:nx, j in 1:ny]
P_vals = [conserved_to_primitive(law, U[i, j])[5] for i in 1:nx, j in 1:ny]

fig = Figure(fontsize = 24, size = (1100, 500))
ax1 = Axis(
    fig[1, 1], xlabel = "x", ylabel = "y",
    title = "Density at t = $(round(t_final, digits = 2))", aspect = DataAspect()
)
hm1 = heatmap!(ax1, xc, yc, rho, colormap = :viridis)
Colorbar(fig[1, 2], hm1)

ax2 = Axis(
    fig[1, 3], xlabel = "x", ylabel = "y",
    title = "Pressure", aspect = DataAspect()
)
hm2 = heatmap!(ax2, xc, yc, P_vals, colormap = :magma)
Colorbar(fig[1, 4], hm2)
resize_to_layout!(fig)
fig

This page was generated using Literate.jl.