2D incompressible MHD equations

[henry2004y]

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!

[xtalax]

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.

[henry2004y]

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.

[xtalax]

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.

[henry2004y]

Is the problem solved? I haven't tested with the WENO solver.

[xtalax]

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.

[henry2004y]

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.

[xtalax]

You can also try FBDF() or QBDF, they seem quite good with advective problems

[henry2004y]

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?

[xtalax]

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

[henry2004y]

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 😓