Open henry2004y opened 2 years ago
xtalax
commented
2 years ago 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
commented
2 years ago 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
commented
2 years ago 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
commented
2 years ago Is the problem solved? I haven't tested with the WENO solver.
xtalax
commented
1 year ago 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
commented
1 year ago 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
commented
1 year ago You can also try FBDF() or QBDF, they seem quite good with advective problems
henry2004y
commented
1 year ago 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
commented
1 year ago 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
commented
1 year ago 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:
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 # η/(va*L₀) const ν = 0.0 # μ/(va*L₀*ρ₀) @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)]) ```
For plotting, I use PyPlot
Plotting script
```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
```julia 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=30areWith the first script using this package, I get rapidly increasing densities, e.g. at
t=1.8which is a hint for instabilityTroubleshooting
Currently I am uncertain where the problem is. Could you please take a look and offer me some guidance? Thanks!