MethodOfLines.jl
MethodOfLines.jl copied to clipboard
SciML
2D incompressible MHD equations
Hi,
I want to see if this package can solve a system of 2D compressible ideal magnetohydrodynamic equations in the X-Z plane.
Problem Description
The original equations are
$$ \begin{align} \frac{\partial\rho}{\partial t} &= -\mathbf{u}\cdot\nabla\rho - \rho\nabla\cdot\mathbf{u}, \ \frac{\partial p}{\partial t} &= -\mathbf{u}\cdot\nabla p - \gamma p \nabla\cdot\mathbf{u}, \ \frac{\partial \mathbf{u}}{\partial t} &= -\rho\mathbf{u}\cdot\nabla\mathbf{u} -\nabla p - \frac{1}{\mu_0}\nabla\frac{B^2}{2} + \frac{1}{\mu_0}\mathbf{B}\cdot\nabla\mathbf{B} + \nu\nabla^2\mathbf{u}, \ \frac{\partial \mathbf{B}}{\partial t} &= \nabla\times\mathbf{u}\times\mathbf{B} - \nabla\times[\eta\nabla\times\mathbf{B}]. \end{align} $$
Since $\nabla\cdot\mathbf{B}=0$, instead of solving the magnetic field directly, we can solve for the magnetic vector potential $\mathbf{A}$. Let $$\mathbf{B} = \nabla\times\mathbf{A} = \nabla\times(0, A_y, 0) = (-\frac{\partial A_y}{\partial z}, 0, \frac{\partial A_y}{\partial x})$$, the last equation above can be simplified to
$$ \frac{\partial\mathbf{A}}{\partial t} = \mathbf{u}\times(\nabla\times\mathbf{A}) - \eta\nabla\times(\nabla\times\mathbf{A}). $$
The normalized 2D equations can be written as
$$ \begin{align} \frac{\partial\rho}{\partial t} &= -u_x\frac{\partial \rho}{\partial x} - u_z\frac{\partial \rho}{\partial z} - \rho\Big( \frac{\partial u_x}{\partial x} + \frac{\partial u_z}{\partial z} \Big), \ \frac{\partial p}{\partial t} &= -u_x\frac{\partial p}{\partial x} - u_z\frac{\partial p}{\partial z} - \gamma p\Big( \frac{\partial u_x}{\partial x} + \frac{\partial u_z}{\partial z} \Big), \ \frac{\partial u_x}{\partial t} &= -u_x\frac{\partial u_x}{\partial x} - u_z\frac{\partial u_x}{\partial z} - \frac{1}{\rho}\frac{\partial}{\partial x}\Big( p + \frac{B^2}{2} \Big) + \frac{1}{\rho}\Big( B_x\frac{\partial Bx}{\partial x} + B_z\frac{\partial Bx}{\partial z} \Big) + \frac{1}{\rho}\nu_m\Big( \frac{\partial^2 u_x}{\partial x^2} + \frac{\partial^2 u_x}{\partial z^2} \Big), \ \frac{\partial u_z}{\partial t} &= -u_x\frac{\partial u_z}{\partial x} - u_z\frac{\partial u_z}{\partial z} - \frac{1}{\rho}\frac{\partial}{\partial z}\Big( p + \frac{B^2}{2} \Big) + \frac{1}{\rho}\Big( B_x\frac{\partial Bz}{\partial x} + B_z\frac{\partial Bz}{\partial z} \Big) + \frac{1}{\rho}\nu_m\Big( \frac{\partial^2 u_z}{\partial x^2} + \frac{\partial^2 u_z}{\partial z^2} \Big), \ \frac{\partial A_y}{\partial t} &= u_x\frac{\partial A_y}{\partial x} - u_z\frac{\partial A_y}{\partial z} + \eta_m\Big( \frac{\partial A_y}{\partial x^2} + \frac{\partial A_y}{\partial z^2} \Big) \end{align} $$
where $\gamma= 5/3$ is the adiabatic index, $\nu_m, \eta_m$ are some normalized constants, and
$$ \begin{align} B^2 = B_x^2 + B_z^2. \end{align} $$
Solving with MethodOfLines.jl
Based on my understanding of the examples given in the tutorials, in principle we shall be able to solve this. For simplicity, I set $\eta_m = 0$ and $\nu_m = 0$. Here is my attempt:
Solving with MethodOfLines.jl
```julia # 2D magnetic reconnection for GEM challenge solved using MethodOfLines.jl. # # Initial condition: # Harris sheet equilibrium with perturbation # # Configuration: # z # Lz/2 | conducting, Bz = ∂Bz/∂z = ∂By/∂z = 0 # | # | periodic # -Lx/2 | Lx/2 # --------------------------------------------> x # | # | # | # -Lz/2 | conducting, Bz = ∂Bz/∂z = ∂By/∂z = 0 # # Ref: # [Fu1995], section 7.4 and Appendix 2 # [Birn+2001]( https://doi.org/10.1029/1999JA900449)
using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets
const Lx = 25.6 const Lz = 12.8 const nx = 16 const nz = 16 "background field" const B₀ = 1.0 "mass density" const ρ₀ = 1.0 "mass density at infinity" const ρ∞ = 0.2ρ₀ "width of current sheet" const λ = 0.5 "perturbation amplitude of the magnetic flux" const ψ₀ = 0.1 "initial plasma β" const β = 1.0 "Alfven velocity" #const va = √(B₀^2/ρ₀) "pressure normalization parameter" const p₀ = 0.5βB₀^2 "temperature normalization parameter" #const T₀ = 0.5β*va^2
physical parameters in MHD equations
"adiabatic index" const γ = 5/3 const η = 0.0 # η/(vaL₀) const ν = 0.0 # μ/(vaL₀*ρ₀)
@parameters x z t #@parameters η, ν @variables ρ(..) p(..) ux(..) uz(..) Ay(..) Bx(..) Bz(..) Dt = Differential(t) Dx = Differential(x) Dz = Differential(z) Dxx = Differential(x)^2 Dzz = Differential(z)^2
∇²(u) = Dxx(u) + Dzz(u)
x_min = -Lx/2 z_min = -Lz/2 t_min = 0.0 x_max = Lx/2 z_max = Lz/2 t_max = 10.0
dx = Lx / nx dz = Lz / nz
ψ(x,z,t) = ψ₀cos(2πx/Lx)cos(πz/Lz)
ρ0(x,z,t) = ρ₀*sech(z/λ)^2 + ρ∞
p0(x,z,t) = begin b = B₀tanh(z/λ) p₀ + 0.5(B₀^2 - b^2) end
ux0(x,z,t) = 0.0 uz0(x,z,t) = 0.0
Bx0(x,z,t) = B₀tanh(z/λ) + ψ₀(-π/Lz)cos(2πx/Lx)sin(πz/Lz) Bz0(x,z,t) = 0.0 + ψ₀*(-2π/Lx)sin(2πx/Lx)cos(πz/Lz)
Ay0(x,z,t) = B₀λlog(cosh(z)) + ψ(x,z,t)
eq = [ Dt(ρ(x,z,t)) ~ -ux(x,z,t)*Dx(ρ(x,z,t)) - uz(x,z,t)Dz(ρ(x,z,t)) - ρ(x,z,t)(Dx(ux(x,z,t)) + Dz(uz(x,z,t))), Dt(p(x,z,t)) ~ -ux(x,z,t)Dx(p(x,z,t)) - uz(x,z,t)Dz(p(x,z,t)) - γp(x,z,t)(Dx(ux(x,z,t)) + Dz(uz(x,z,t))), Dt(ux(x,z,t)) ~ -ux(x,z,t)*Dx(ux(x,z,t)) - uz(x,z,t)Dz(ux(x,z,t)) + 1/ρ(x,z,t)(Bx(x,z,t)*Dx(Bx(x,z,t)) + Bz(x,z,t)*Dz(Bx(x,z,t)) - Dx(p(x,z,t)) - (Bx(x,z,t)*Dx(Bx(x,z,t)) + Bz(x,z,t)Dx(Bz(x,z,t))) + ν∇²(ux(x,z,t))), Dt(uz(x,z,t)) ~ -ux(x,z,t)*Dx(uz(x,z,t)) - uz(x,z,t)Dz(uz(x,z,t)) + 1/ρ(x,z,t)(Bx(x,z,t)*Dx(Bz(x,z,t)) + Bz(x,z,t)*Dz(Bz(x,z,t)) - Dz(p(x,z,t)) - (Bx(x,z,t)*Dz(Bx(x,z,t)) + Bz(x,z,t)Dz(Bz(x,z,t))) + ν∇²(uz(x,z,t))), Dt(Ay(x,z,t)) ~ -ux(x,z,t)*Dx(Ay(x,z,t)) - uz(x,z,t)Dz(Ay(x,z,t)) + η∇²(Ay(x,z,t)), Bx(x,z,t) ~ -Dz(Ay(x,z,t)), Bz(x,z,t) ~ Dx(Ay(x,z,t)) ]
domains = [x ∈ Interval(x_min, x_max), z ∈ Interval(z_min, z_max), t ∈ Interval(t_min, t_max)]
BCs: periodic in x, Neumann in z
ICs: set from functions
bcs = [ρ(x,z,0) ~ ρ0(x,z,0), ρ(x_min,z,t) ~ ρ(x_max,z,t), Dz(ρ(x,z_min,t)) ~ 0.0, Dz(ρ(x,z_max,t)) ~ 0.0,
p(x,z,0) ~ p0(x,z,0),
p(x_min,z,t) ~ p(x_max,z,t),
Dz(p(x,z_min,t)) ~ 0.0,
Dz(p(x,z_max,t)) ~ 0.0,
ux(x,z,0) ~ ux0(x,z,0),
ux(x_min,z,t) ~ ux(x_max,z,t),
Dz(ux(x,z_min,t)) ~ 0.0,
Dz(ux(x,z_max,t)) ~ 0.0,
uz(x,z,0) ~ uz0(x,z,0),
uz(x_min,z,t) ~ uz(x_max,z,t),
Dz(uz(x,z_min,t)) ~ 0.0,
Dz(uz(x,z_max,t)) ~ 0.0,
Ay(x,z,0) ~ Ay0(x,z,0),
Ay(x_min,z,t) ~ Ay(x_max,z,t),
Dz(Ay(x,z_min,t)) ~ 0.0,
Dz(Ay(x,z_max,t)) ~ 0.0,
Bx(x,z,0) ~ Bx0(x,z,0),
Bx(x_min,z,t) ~ Bx(x_max,z,t),
Dz(Bx(x,z_min,t)) ~ 0.0,
Dz(Bx(x,z_max,t)) ~ 0.0,
Bz(x,z,0) ~ Bz0(x,z,0),
Bz(x_min,z,t) ~ Bz(x_max,z,t),
Dz(Bz(x,z_min,t)) ~ 0.0,
Dz(Bz(x,z_max,t)) ~ 0.0,
]
@named pdesys = PDESystem(eq, bcs, domains, [x,z,t], [ρ(x,z,t), p(x,z,t), ux(x,z,t), uz(x,z,t), Ay(x,z,t), Bx(x,z,t), Bz(x,z,t)])
Discretization
order = 2
discretization = MOLFiniteDifference([x=>dx, z=>dz], t, approx_order=order, grid_align=center_align)
Convert the PDE problem into an ODE problem
println("Discretization:") @time prob = discretize(pdesys, discretization)
println("Solve:") #@time sol = solve(prob, Tsit5(), saveat=0.1) @time sol = solve(prob, RK4(), dt=0.05, saveat=0.1)
Extracting results
grid = get_discrete(pdesys, discretization) discrete_x = grid[x] discrete_z = grid[z] discrete_t = sol[t]
@time solBx = map(d -> sol[d][end], grid[Bx(x, z, t)]) solBz = map(d -> sol[d][end], grid[Bz(x, z, t)]) solρ = map(d -> sol[d][end], grid[ρ(x, z, t)])
</p>
</details>
For plotting, I use PyPlot
<details><summary>Plotting script</summary>
<p>
```julia
using PyPlot
@static if matplotlib.__version__ < "3.5"
matplotlib.rc("pcolor", shading="nearest") # newer version default "auto"
end
matplotlib.rc("font", size=14)
matplotlib.rc("xtick", labelsize=10)
matplotlib.rc("ytick", labelsize=10)
function plot_snapshot(xrange, zrange, bx, bz)
# meshgrid: note the array ordering difference between Julia and Python!
X = [i for _ in zrange, i in xrange]
Z = [j for j in zrange, _ in xrange]
fig, ax = subplots(1,1, figsize=(12,8), constrained_layout=true)
im = ax.pcolormesh(xrange, zrange, bz', cmap=matplotlib.cm.RdBu_r)
ax.streamplot(X, Z, bx', bz', color="k")
ax.set_xlabel("x")
ax.set_ylabel("z")
ax.set_title("Bz")
fig.colorbar(im; ax)
return
end
figure()
pcolormesh(discrete_x, discrete_z, solρ', cmap=matplotlib.cm.RdBu_r, shading="nearest")
xlabel("x")
ylabel("z")
colorbar()
plot_snapshot(discrete_x, discrete_z, solBx, solBz)
Testing
I hope I don't make mistakes in expressing the system of PDEs, but the test result is not quite what I expect: it quickly develops some numerical instabilities. As a comparison, here is my hand-written script for solving the PDEs with RK4 in time (fixed timestep) and central differencing in space:
Hand-written finite difference code
using PyPlot
@static if matplotlib.__version__ ≥ "3.3"
matplotlib.rc("image", cmap="turbo") # set default colormap
end
@static if matplotlib.__version__ < "3.5"
matplotlib.rc("pcolor", shading="nearest") # newer version default "auto"
end
matplotlib.rc("font", size=14)
matplotlib.rc("xtick", labelsize=10)
matplotlib.rc("ytick", labelsize=10)
Base.@kwdef struct Parameter
Lx::Float64 = 25.6
Lz::Float64 = 12.8
nx::Int = 18 #34
nz::Int = 18 #34
nt::Int = 600
"background field"
B₀::Float64 = 1.0
"mass density"
ρ₀::Float64 = 1.0
"mass density at infinity"
ρ∞::Float64 = 0.2*ρ₀
"width of current sheet"
λ::Float64 = 0.5
"perturbation amplitude of the magnetic flux"
ψ₀::Float64 = 0.1
"initial plasma β"
β::Float64 = 1.0
"Alfven velocity"
va::Float64 = √(B₀^2/ρ₀)
"pressure normalization parameter"
p₀::Float64 = 0.5*β*B₀^2
"temperature normalization parameter"
T₀::Float64 = 0.5*β*va^2
# physical parameters in MHD equations
"adiabatic index"
γ::Float64 = 5/3
η::Float64 = 0.0 # η/(va*L₀)
ν::Float64 = 0.0 # μ/(va*L₀*ρ₀)
"output cadence"
nplot::Int = 600
# array indices for different variables
ρ_::Int = 1
p_::Int = 2
ux_::Int = 3
uz_::Int = 4
ay_::Int = 5
dx::Float64 = Lx/nx
dz::Float64 = Lz/nz
dt::Float64 = 0.05 # giving a dt < min(dx,dz)/(√(1.0+0.5*γ*β)*va)
inv2dx::Float64 = nx/(2*Lx)
inv2dz::Float64 = nz/(2*Lz)
invdx²::Float64 = (nx/Lx)^2
invdz²::Float64 = (nz/Lz)^2
end
struct Variable
state::Array{Float64,3}
statetmp::Array{Float64,3}
bx::Array{Float64,2}
bz::Array{Float64,2}
"total pressure"
pt::Array{Float64,2}
"maximum bz magnitudes"
bzm::Vector{Float64}
# intermediate arrays for rk4
rrho::Array{Float64,2} # 1/rho
f1::Array{Float64,3}
f2::Array{Float64,3}
f3::Array{Float64,3}
f4::Array{Float64,3}
function Variable(nx::Int, nz::Int, nt::Int)
state = zeros(nx, nz, 5)
statetmp = zeros(nx, nz, 5)
bx = zeros(nx, nz)
bz = zeros(nx, nz)
pt = zeros(nx, nz) # pt = p + b^2/2
bzm = zeros(nt)
rrho = zeros(nx, nz)
f1 = zeros(nx, nz, 5)
f2 = zeros(nx, nz, 5)
f3 = zeros(nx, nz, 5)
f4 = zeros(nx, nz, 5)
new(state, statetmp, bx, bz, pt, bzm, rrho, f1, f2, f3, f4)
end
end
function solve!(param::Parameter, var::Variable)
(;nt, dt, nplot, ρ_, p_) = param
(;state, bzm) = var
t = 0.0
set_initial_condition!(param, var)
fig, cs = save_snapshot(param, var)
for it = 1:nt
bzm[it] = get_bzmax(param, var)
if mod(it-1, nplot) == 0
println(it, ", max(Bz) = ", bzm[it])
save_snapshot!(var, it, fig, cs)
#sleep(2.0)
end
t += dt
update!(param, var)
ρmin = @views minimum(state[2:end-1,2:end-1,ρ_])
pmin = @views minimum(state[2:end-1,2:end-1,p_])
if ρmin < 0
index = @views argmin(state[:,:,ρ_])
@info index, state[index[1],index[2],ρ_]
error("Negative density at step $it")
end
if pmin < 0
index = @views argmin(state[:,:,p_])
@info index, state[index[1],index[2],p_]
error("Negative pressure at step $it")
end
end
println("Finished at step $nt, t = $t")
return
end
"""
set_initial_condition(param::Parameter, var::Variable)
Set initial condition as a perturbation to the Harris current sheet equilibrium.
"""
function set_initial_condition!(param::Parameter, var::Variable)
(;nx, nz, Lx, Lz, B₀, ρ₀, ρ∞, p₀, λ, ψ₀, ρ_, p_, ay_) = param
(; state) = var
x = range(-Lx/2, Lx/2, length=nx)
z = range(-Lz/2, Lz/2, length=nz)
state .= 0.0
# Harris current sheet
for k in eachindex(z)
ρ = ρ₀*sech(z[k]/λ)^2 + ρ∞ # (p₀ + 0.5*(B₀^2 - b^2)) / T₀
b = B₀*tanh(z[k]/λ)
state[2:end-1,k,ay_] .= B₀*λ*log(cosh(z[k]))
state[2:end-1,k,ρ_] .= ρ
state[2:end-1,k,p_] .= p₀ + 0.5*(B₀^2 - b^2)
end
# Perturbation in B, or flux function
for k in eachindex(z), i in eachindex(x)
#δBx = ψ₀*(-π/Lz)*cos(2πx[i]/Lx)*sin(πz[k]/Lz)
#δBz = ψ₀*(-2π/Lx)*sin(2πx[i]/Lx)*cos(πz[k]/Lz)
state[i,k,ay_] += ψ₀*cos(2π*x[i]/Lx)*cos(π*z[k]/Lz)
end
# Neumann B.C. in z
state[:,1,:] = state[:,2,:]
state[:,end,:] = state[:,end-1,:]
# periodic B.C. in x
state[1,:,:] = state[end-1,:,:]
state[end,:,:] = state[2,:,:]
return
end
"""
update!(param::Parameter, var::Variable)
One step update with 1st order in time and RK4 in space.
"""
function update!(param::Parameter, var::Variable)
(;dt) = param
(;state, statetmp, f1, f2, f3, f4) = var
rhs!(param, var, f1, state)
@. statetmp = state + 0.5*dt*f1
rhs!(param, var, f2, statetmp)
@. statetmp = state + 0.5*dt*f2
rhs!(param, var, f3, statetmp)
@. statetmp = state + dt*f3
rhs!(param, var, f4, statetmp)
@. state += dt*(f1 + 2.0*f2 + 2.0*f3 + f4)/6.0
return
end
"Compute for rk4 the right hand side of mhd equations."
function rhs!(param::Parameter, var::Variable, varout::Array{Float64,3}, varin::Array{Float64,3})
(;nx, nz, inv2dx, inv2dz, invdx², invdz², γ, ν, η, ρ_, p_, ux_, uz_, ay_) = param
(;rrho, bx, bz, pt) = var
# calculate Bx, Bz
calcb!(param, var, varin)
for i = 2:nx-1, j = 2:nz-1
rrho[i,j] = 1.0 / varin[i,j,ρ_]
pt[i,j] = varin[i,j,p_] + 0.5*(bx[i,j]^2 + bz[i,j]^2)
end
set_BC!(param, rrho)
set_BC!(param, pt)
for j = 2:nz-1
jm = j - 1
jp = j + 1
for i = 2:nx-1
varout[i,j,ρ_] =
-varin[i,j,ux_]*inv2dx*(varin[i+1,j , ρ_] - varin[i-1,j , ρ_]) -
varin[i,j,uz_]*inv2dz*(varin[i ,jp, ρ_] - varin[i ,jm, ρ_]) -
varin[i,j,ρ_]*(inv2dx*(varin[i+1,j ,ux_] - varin[i-1,j ,ux_]) +
inv2dz*(varin[i ,jp,uz_] - varin[i ,jm,uz_]))
varout[i,j,p_] =
-varin[i,j,ux_]*inv2dx*(varin[i+1,j ,p_] - varin[i-1,j ,p_]) -
varin[i,j,uz_]*inv2dz*(varin[i ,jp,p_] - varin[i ,jm,p_]) -
γ*varin[i,j,p_]*(inv2dx*(varin[i+1,j ,ux_] - varin[i-1,j ,ux_]) +
inv2dz*(varin[i ,jp,uz_] - varin[i ,jm,uz_]))
varout[i,j,ux_] =
-varin[i,j,ux_]*inv2dx*(varin[i+1,j ,ux_] - varin[i-1,j ,ux_]) -
varin[i,j,uz_]*inv2dz*(varin[i ,jp,ux_] - varin[i ,jm,ux_]) +
rrho[i,j]*( (bx[i,j]*inv2dx*(bx[i+1,j ] - bx[i-1,j ]) +
bz[i,j]*inv2dz*(bx[i ,jp] - bx[i ,jm]) -
inv2dx*(pt[i+1,j ] - pt[i-1,j ])) +
ν*(invdx²*(varin[i+1,j ,ux_] + varin[i-1,j ,ux_] - 2.0*varin[i,j,ux_]) +
invdz²*(varin[i ,jp,ux_] + varin[i ,jm,ux_] - 2.0*varin[i,j,ux_])) )
varout[i,j,uz_] =
-varin[i,j,ux_]*inv2dx*(varin[i+1,j ,uz_] - varin[i-1,j ,uz_]) -
varin[i,j,uz_]*inv2dz*(varin[i ,jp,uz_] - varin[i ,jm,uz_]) +
rrho[i,j]*( (bx[i,j]*inv2dx*(bz[i+1,j ] - bz[i-1,j ]) +
bz[i,j]*inv2dz*(bz[i ,jp] - bz[i ,jm]) -
inv2dz*(pt[i ,jp] - pt[i ,jm])) +
ν*(invdx²*(varin[i+1,j ,uz_] + varin[i-1,j ,uz_] - 2.0*varin[i,j,uz_]) +
invdz²*(varin[i ,jp,uz_] + varin[i ,jm,uz_] - 2.0*varin[i,j,uz_])) )
varout[i,j,ay_] =
-varin[i,j,ux_]*inv2dx*(varin[i+1,j ,ay_] - varin[i-1,j ,ay_]) -
varin[i,j,uz_]*inv2dz*(varin[i ,jp,ay_] - varin[i ,jm,ay_]) +
η*(invdx²*(varin[i+1,j ,ay_] + varin[i-1,j ,ay_] - 2.0*varin[i,j,ay_]) +
invdz²*(varin[i ,jp,ay_] + varin[i ,jm,ay_] - 2.0*varin[i,j,ay_]))
end
end
set_BC!(param, varout)
return
end
"Calculate Bx, Bz."
function calcb!(param::Parameter, var::Variable, varin::Array{Float64,3})
(;nx, nz, inv2dx, inv2dz, ay_) = param
(;bx, bz) = var
# calculate Bx, Bz
for i = 2:nx-1, j = 2:nz-1
jp = j + 1
jm = j - 1
bx[i,j] = -inv2dz*(varin[i,jp,ay_] - varin[i,jm,ay_])
bz[i,j] = inv2dx*(varin[i+1,j,ay_] - varin[i-1,j,ay_])
end
set_BC!(param, bx)
set_BC!(param, bz)
return
end
function set_BC!(param::Parameter, var::Array{Float64,2})
(;nx, nz) = param
# x
var[1,2:nz-1] = var[end-1,2:nz-1]
var[end,2:nz-1] = var[2,2:nz-1]
# z
var[2:nx-1,1] = var[2:nx-1,2]
var[2:nx-1,end] = var[2:nx-1,end-1]
end
function set_BC!(param::Parameter, var::Array{Float64,3})
(;nx, nz, ρ_, p_, ux_, uz_, ay_) = param
var[1,2:nz-1,ρ_] = var[end-1,2:nz-1,ρ_]
var[1,2:nz-1,p_] = var[end-1,2:nz-1,p_]
var[1,2:nz-1,ux_] = var[end-1,2:nz-1,ux_]
var[1,2:nz-1,uz_] = var[end-1,2:nz-1,uz_]
var[1,2:nz-1,ay_] = var[end-1,2:nz-1,ay_]
var[end,2:nz-1,ρ_] = var[2,2:nz-1,ρ_]
var[end,2:nz-1,p_] = var[2,2:nz-1,p_]
var[end,2:nz-1,ux_] = var[2,2:nz-1,ux_]
var[end,2:nz-1,uz_] = var[2,2:nz-1,uz_]
var[end,2:nz-1,ay_] = var[2,2:nz-1,ay_]
# z boundary
var[2:nx-1,1,ρ_] = var[2:nx-1,2,ρ_]
var[2:nx-1,1,p_] = var[2:nx-1,2,p_]
var[2:nx-1,1,ux_] = var[2:nx-1,2,ux_]
var[2:nx-1,1,uz_] = var[2:nx-1,2,uz_]
var[2:nx-1,1,ay_] = var[2:nx-1,2,ay_]
var[2:nx-1,end,ρ_] = var[2:nx-1,end-1,ρ_]
var[2:nx-1,end,p_] = var[2:nx-1,end-1,p_]
var[2:nx-1,end,ux_] = var[2:nx-1,end-1,ux_]
var[2:nx-1,end,uz_] = var[2:nx-1,end-1,uz_]
var[2:nx-1,end,ay_] = var[2:nx-1,end-1,ay_]
end
"Calculate Bz max magnitude."
function get_bzmax(param::Parameter, var::Variable)
# Calculate B
calcb!(param, var, var.state)
bzm = maximum(abs, var.bz; init=-100.0)
end
"Save snapshots."
function save_snapshot(param::Parameter, var::Variable)
(;nx, nz, nt, dx, dz, dt) = param
xl = dx*(nx-1)/2
zl = dz*(nz-1)
xrange = -xl:dx:xl
zrange = 0:dz:zl
t = range(0, dt*(nt-1), step=dt)
# meshgrid: note the array ordering difference between Julia and Python!
X = [i for _ in zrange, i in xrange]
Z = [j for j in zrange, _ in xrange]
fig, axs = subplots(2,2, figsize=(12,8), constrained_layout=true)
ρmin, ρmax = 0.0, 1.0
umin, umax = -0.35, 0.35
bmin, bmax = -0.07, 0.07
c1 = @views axs[1,1].pcolormesh(xrange, zrange, var.state[:,:,1]'; vmin=ρmin, vmax=ρmax)
c2 = @views axs[1,2].pcolormesh(xrange, zrange, var.state[:,:,3]'; vmin=umin, vmax=umax,
cmap=matplotlib.cm.RdBu_r)
c3 = axs[2,1].pcolormesh(xrange, zrange, var.bz'; vmin=bmin, vmax=bmax,
cmap=matplotlib.cm.RdBu_r)
l1 = axs[2,2].plot(t, zero(var.bzm))
axs[2,2].set_xlim(0, dt*nt)
axs[2,2].set_ylim(-3.8, -3.2)
for ax in axs[1:3]
ax.set_xlabel("x")
ax.set_ylabel("z")
end
im_ratio = length(zrange)/length(xrange)
fraction = 0.046 * im_ratio
ticks = (range(ρmin, ρmax, length=7), range(umin, umax, length=7),
range(bmin, bmax, length=7))
cb1 = colorbar(c1; ax=axs[1,1], ticks=ticks[1], fraction, pad=0.02)
cb2 = colorbar(c2; ax=axs[1,2], ticks=ticks[2], fraction, pad=0.02)
cb3 = colorbar(c3; ax=axs[2,1], ticks=ticks[3], fraction, pad=0.02)
titles = (L"\rho", "Bz", "Ux", "log(max(Bz))")
for (ax, title) in zip(axs, titles)
ax.set_title(title)
end
return fig, (c1, c2, c3, l1)
end
"Save snapshots by overwriting `fig` and `axs`."
function save_snapshot!(var::Variable, it::Int, fig, cs)
fig.suptitle("2D MHD tearing mode, it = $it")
cs[1].set_array(var.state[:,:,1]')
cs[2].set_array(var.state[:,:,3]')
cs[3].set_array(var.bz')
cs[4][1].set_ydata(log.(var.bzm))
savefig("$(lpad(it, 4, '0')).png")
return
end
function plot_snapshot(param::Parameter, var::Variable)
(;nx, nz, nt, dx, dz, dt) = param
(;state, bx, bz, bzm) = var
xl = dx*(nx-1)/2
zl = dz*(nz-1)
xrange = -xl:dx:xl
zrange = 0:dz:zl
t = range(0, dt*(nt-1), step=dt)
# meshgrid: note the array ordering difference between Julia and Python!
X = [i for _ in zrange, i in xrange]
Z = [j for j in zrange, _ in xrange]
fig, axs = subplots(2,2, figsize=(12,8), constrained_layout=true)
im1 = @views axs[1,1].pcolormesh(xrange, zrange, state[:,:,1]')
im2 = @views axs[1,2].pcolormesh(xrange, zrange, state[:,:,3]', cmap=matplotlib.cm.RdBu_r)
im3 = axs[2,1].pcolormesh(xrange, zrange, bz', cmap=matplotlib.cm.RdBu_r)
axs[2,1].streamplot(X, Z, bx', bz', color="k")
axs[2,2].plot(t, log.(bzm))
for ax in axs[1:3]
ax.set_xlabel("x")
ax.set_ylabel("z")
end
titles = (L"\rho", "Bz", "Ux", "log(max(Bz))")
for (ax, title) in zip(axs, titles)
ax.set_title(title)
end
fig.colorbar(im1, ax=axs[1,1])
fig.colorbar(im2, ax=axs[1,2])
fig.colorbar(im3, ax=axs[2,1])
return
end
##### Main
param = Parameter()
var = Variable(param.nx, param.nz, param.nt)
set_initial_condition!(param, var)
calcb!(param, var, var.state)
plot_snapshot(param, var)
solve!(param, var)
plot_snapshot(param, var)
With my hand-written script, the initial condition looks like
and the solutions at t=30 are
With the first script using this package, I get rapidly increasing densities, e.g. at t=1.8 which is a hint for instability
Troubleshooting
Currently I am uncertain where the problem is. Could you please take a look and offer me some guidance? Thanks!
I am currently working on WENO schemes, which are good for discretizing systems such as this one. The upwind scheme (which we currently use) has exhibited instability for similar problems. Soon we should be able to solve your problem more effectively.
Glad to hear that! However, the comparison script I posted does not use WENO schemes; I believe it is also using the central difference schemes for the 1st and 2nd order derivative terms. So maybe I made some mistakes in expressing the system? That one will also get into numerical issues but at a later time, where some of the variables (i.e. density, pressure) become negative.
Ah, that is a bug then. May be related to #130, the workaround was to rearrange the form of the coefficients, though I am not sure this is applicable here.
Thanks for your excellently written issue, this will be very helpful in debugging.
Please try this with the WENO scheme @henry2004y. Note that this scheme is sensitive to solver choice, for best results use a strong stability preserving solver like SSPRK33() with an appropriately small fixed dt.
Thanks for the reminder. I will try later today.
Currently it does not look promising. I've been waiting for around an hour in the discretization step for a grid size of 16*16. Let's see when it will be finished.
You can also try FBDF() or QBDF, they seem quite good with advective problems
With
discretization = MOLFiniteDifference([x=>dx, z=>dz], t,
advection_scheme=FBDF(),#WENOScheme(),
approx_order=2,
grid_align=center_align)
in MethodOfLines v0.7.2, I get
Discretization:
ERROR: LoadError: ArgumentError: Only `UpwindScheme()` and `WENOScheme()` are supported advection schemes.
Are FBDF() and QBDF() available in master?
I mean as the solver, not the advection scheme - WENO it must be noted has issues with non periodic bcs, or less than 2 BCs per boundary - so try this too - a neumann0 condition often affects the solution little
I mean as the solver, not the advection scheme - WENO it must be noted has issues with non periodic bcs, or less than 2 BCs per boundary - so try this too - a neumann0 condition often affects the solution little
Ah, sorry. Need to refresh my memory a bit :sweat: