GRMHD in Flat Spacetime
The general relativistic MHD (GRMHD) solver in FiniteVolumeMethod.jl implements the Valencia formulation of the relativistic MHD equations on a curved spacetime. In flat (Minkowski) spacetime the geometric source terms vanish and the GRMHD equations reduce exactly to the special relativistic MHD (SRMHD) equations.
This tutorial verifies that equivalence by solving the Balsara 1 shock tube with both SRMHD and GRMHD solvers and comparing the density profiles.
Problem Setup
The Balsara 1 initial data are:
\[(\rho, v_x, v_y, v_z, P, B_x, B_y, B_z) = \begin{cases} (1,\, 0,\, 0,\, 0,\, 1,\, 0.5,\, 1,\, 0) & x < 0.5,\\ (0.125,\, 0,\, 0,\, 0,\, 0.1,\, 0.5,\, -1,\, 0) & x \geq 0.5. \end{cases}\]
using FiniteVolumeMethod
using StaticArrays
gamma = 5.0 / 3.0
eos = IdealGasEOS(gamma)IdealGasEOS{Float64}(1.6666666666666667)Create both the SRMHD and GRMHD conservation laws:
law_sr = SRMHDEquations{1}(eos)
law_gr = GRMHDEquations{1}(eos, MinkowskiMetric{1}())
wL = SVector(1.0, 0.0, 0.0, 0.0, 1.0, 0.5, 1.0, 0.0)
wR = SVector(0.125, 0.0, 0.0, 0.0, 0.1, 0.5, -1.0, 0.0)
N = 400
mesh = StructuredMesh1D(0.0, 1.0, N)
ic(x) = x < 0.5 ? wL : wRic (generic function with 1 method)Solving with SRMHD
prob_sr = HyperbolicProblem(
law_sr, mesh, HLLSolver(), CellCenteredMUSCL(MinmodLimiter()),
TransmissiveBC(), TransmissiveBC(), ic;
final_time = 0.4, cfl = 0.4
)
x_sr, U_sr, t_sr = solve_hyperbolic(prob_sr)400-element Vector{Float64}:
0.00125
0.00375
0.00625
0.00875
0.01125
⋮
0.98875
0.9912500000000001
0.99375
0.99625
0.99875Solving with GRMHD (Minkowski)
prob_gr = HyperbolicProblem(
law_gr, mesh, HLLSolver(), CellCenteredMUSCL(MinmodLimiter()),
TransmissiveBC(), TransmissiveBC(), ic;
final_time = 0.4, cfl = 0.4
)
x_gr, U_gr, t_gr = solve_hyperbolic(prob_gr)([0.00125, 0.00375, 0.00625, 0.00875, 0.01125, 0.01375, 0.01625, 0.01875, 0.02125, 0.02375 … 0.9762500000000001, 0.97875, 0.9812500000000001, 0.98375, 0.9862500000000001, 0.98875, 0.9912500000000001, 0.99375, 0.99625, 0.99875], StaticArraysCore.SVector{8, Float64}[[1.0, 1.8975096649034662e-16, 0.0, 0.0, 2.1249999999999996, 0.5, 1.0, 0.0], [1.0, 1.8974995317206407e-16, 0.0, 0.0, 2.1249999999999996, 0.5, 1.0, 0.0], [1.0, 1.897493649679863e-16, 0.0, 0.0, 2.1249999999999996, 0.5, 1.0, 0.0], [1.0, 1.8974780446347618e-16, 0.0, 0.0, 2.1249999999999996, 0.5, 1.0, 0.0], [1.0, 1.897470474040355e-16, 0.0, 0.0, 2.1249999999999996, 0.5, 1.0, 0.0], [1.0, 1.897457699837429e-16, 0.0, 0.0, 2.1249999999999996, 0.5, 1.0, 0.0], [1.0, 1.8974473522456973e-16, 0.0, 0.0, 2.1249999999999996, 0.5, 1.0, 0.0], [1.0, 1.8974411304271518e-16, 0.0, 0.0, 2.1249999999999996, 0.5, 1.0, 0.0], [1.0, 1.8974252681440586e-16, 0.0, 0.0, 2.1249999999999996, 0.5, 1.0, 0.0], [1.0, 1.8974175331936868e-16, 0.0, 0.0, 2.1249999999999996, 0.5, 1.0, 0.0] … [0.12500000000059514, 6.609108543364055e-12, 2.876804147948251e-12, 0.0, 0.7750000000065275, 0.5, -1.000000000005337, 0.0], [0.12500000000027764, 3.0843531893031072e-12, 1.3427325268892243e-12, 0.0, 0.7750000000030457, 0.5, -1.0000000000024905, 0.0], [0.12500000000012834, 1.4265633553652736e-12, 6.211010706303026e-13, 0.0, 0.7750000000014082, 0.5, -1.0000000000011513, 0.0], [0.12500000000005873, 6.534839892713321e-13, 2.8448299171180675e-13, 0.0, 0.7750000000006448, 0.5, -1.0000000000005271, 0.0], [0.12500000000002653, 2.95697618954179e-13, 1.2856189712873704e-13, 0.0, 0.7750000000002915, 0.5, -1.000000000000238, 0.0], [0.12500000000001182, 1.3212140493325892e-13, 5.73020503493442e-14, 0.0, 0.7750000000001298, 0.5, -1.000000000000106, 0.0], [0.1250000000000051, 5.725267464193994e-14, 2.4919728248789634e-14, 0.0, 0.775000000000056, 0.5, -1.0000000000000457, 0.0], [0.1250000000000021, 2.4077730434246706e-14, 1.0548785638637285e-14, 0.0, 0.7750000000000234, 0.5, -1.0000000000000187, 0.0], [0.1250000000000008, 9.572673982701237e-15, 4.1513801063167406e-15, 0.0, 0.775000000000009, 0.5, -1.000000000000007, 0.0], [0.12500000000000017, 2.755560103738745e-15, 1.1627135318109971e-15, 0.0, 0.7750000000000021, 0.5, -1.0000000000000018, 0.0]], 0.4)Comparing Results
rho_sr = [conserved_to_primitive(law_sr, U_sr[i])[1] for i in eachindex(U_sr)]
vx_sr = [conserved_to_primitive(law_sr, U_sr[i])[2] for i in eachindex(U_sr)]
By_sr = [conserved_to_primitive(law_sr, U_sr[i])[7] for i in eachindex(U_sr)]
rho_gr = [conserved_to_primitive(law_gr, U_gr[i])[1] for i in eachindex(U_gr)]
vx_gr = [conserved_to_primitive(law_gr, U_gr[i])[2] for i in eachindex(U_gr)]
By_gr = [conserved_to_primitive(law_gr, U_gr[i])[7] for i in eachindex(U_gr)]400-element Vector{Float64}:
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
⋮
-1.0000000000024905
-1.0000000000011513
-1.0000000000005271
-1.000000000000238
-1.000000000000106
-1.0000000000000457
-1.0000000000000187
-1.000000000000007
-1.0000000000000018The two solutions should match to high accuracy:
max_rho_diff = maximum(abs.(rho_sr .- rho_gr))trueVisualisation
using CairoMakie
fig = Figure(fontsize = 20, size = (1200, 400))
ax1 = Axis(fig[1, 1], xlabel = "x", ylabel = L"\rho", title = "Density")
lines!(ax1, x_sr, rho_sr, color = :blue, linewidth = 2, label = "SRMHD")
scatter!(ax1, x_gr, rho_gr, color = :red, markersize = 3, label = "GRMHD (Minkowski)")
axislegend(ax1, position = :ct)
ax2 = Axis(fig[1, 2], xlabel = "x", ylabel = L"v_x", title = "Velocity")
lines!(ax2, x_sr, vx_sr, color = :blue, linewidth = 2)
scatter!(ax2, x_gr, vx_gr, color = :red, markersize = 3)
ax3 = Axis(fig[1, 3], xlabel = "x", ylabel = L"B_y", title = "Magnetic Field (y)")
lines!(ax3, x_sr, By_sr, color = :blue, linewidth = 2)
scatter!(ax3, x_gr, By_gr, color = :red, markersize = 3)
resize_to_layout!(fig)
fig
The SRMHD (blue lines) and GRMHD (red dots) solutions overlap perfectly, confirming that the GRMHD solver correctly reduces to SRMHD in flat spacetime. The maximum density difference is round(max_rho_diff, sigdigits = 2).
trueJust 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_sr = SRMHDEquations{1}(eos)
law_gr = GRMHDEquations{1}(eos, MinkowskiMetric{1}())
wL = SVector(1.0, 0.0, 0.0, 0.0, 1.0, 0.5, 1.0, 0.0)
wR = SVector(0.125, 0.0, 0.0, 0.0, 0.1, 0.5, -1.0, 0.0)
N = 400
mesh = StructuredMesh1D(0.0, 1.0, N)
ic(x) = x < 0.5 ? wL : wR
prob_sr = HyperbolicProblem(
law_sr, mesh, HLLSolver(), CellCenteredMUSCL(MinmodLimiter()),
TransmissiveBC(), TransmissiveBC(), ic;
final_time = 0.4, cfl = 0.4
)
x_sr, U_sr, t_sr = solve_hyperbolic(prob_sr)
prob_gr = HyperbolicProblem(
law_gr, mesh, HLLSolver(), CellCenteredMUSCL(MinmodLimiter()),
TransmissiveBC(), TransmissiveBC(), ic;
final_time = 0.4, cfl = 0.4
)
x_gr, U_gr, t_gr = solve_hyperbolic(prob_gr)
rho_sr = [conserved_to_primitive(law_sr, U_sr[i])[1] for i in eachindex(U_sr)]
vx_sr = [conserved_to_primitive(law_sr, U_sr[i])[2] for i in eachindex(U_sr)]
By_sr = [conserved_to_primitive(law_sr, U_sr[i])[7] for i in eachindex(U_sr)]
rho_gr = [conserved_to_primitive(law_gr, U_gr[i])[1] for i in eachindex(U_gr)]
vx_gr = [conserved_to_primitive(law_gr, U_gr[i])[2] for i in eachindex(U_gr)]
By_gr = [conserved_to_primitive(law_gr, U_gr[i])[7] for i in eachindex(U_gr)]
max_rho_diff = maximum(abs.(rho_sr .- rho_gr))
using CairoMakie
fig = Figure(fontsize = 20, size = (1200, 400))
ax1 = Axis(fig[1, 1], xlabel = "x", ylabel = L"\rho", title = "Density")
lines!(ax1, x_sr, rho_sr, color = :blue, linewidth = 2, label = "SRMHD")
scatter!(ax1, x_gr, rho_gr, color = :red, markersize = 3, label = "GRMHD (Minkowski)")
axislegend(ax1, position = :ct)
ax2 = Axis(fig[1, 2], xlabel = "x", ylabel = L"v_x", title = "Velocity")
lines!(ax2, x_sr, vx_sr, color = :blue, linewidth = 2)
scatter!(ax2, x_gr, vx_gr, color = :red, markersize = 3)
ax3 = Axis(fig[1, 3], xlabel = "x", ylabel = L"B_y", title = "Magnetic Field (y)")
lines!(ax3, x_sr, By_sr, color = :blue, linewidth = 2)
scatter!(ax3, x_gr, By_gr, color = :red, markersize = 3)
resize_to_layout!(fig)
figThis page was generated using Literate.jl.