Solving Grad-Shafranov Equation with ModelingToolkit

General Usage
13

Yeah but it should let me do my original problem of solving the Grad-Shafranov equation but I’m running into a new error

using ModelingToolkit, DiffEqOperators, DifferentialEquations, SteadyStateDiffEq, PyPlot

δ = 0.33
ϵ = 0.32
κ = 1.7

δ₀ = asin(δ)
N₁ = -(1 + δ₀)^2/(ϵ*κ^2) # [d²x/dy²]_(τ=0)
N₂ =  (1 - δ₀)^2/(ϵ*κ^2) # [d²x/dy²]_(τ=pi)
N₃ = -κ/(ϵ*cos(δ₀)^2)    # [d²y/dx²]_(τ=pi/2)

α = -0.155

@parameters x y t
@variables u(..)
Dxx = Differential(x)^2
Dyy = Differential(y)^2
Dx = Differential(x)
Dy = Differential(y)
Dt = Differential(t)

# Solovev Solution to Grad-Shafranov Δ⋆ψ(x,y) = α + (1 - α)x^2
eq  = Dt(u(t,x,y)) ~ Dxx(u(t,x,y)) - Dx(u(t,x,y))/x + Dyy(u(t,x,y)) - α - (1 - α)*x^2

# Boundary conditions
bcs = [u(0,x,y) ~ 0.0, # initial guess
       u(t,1 + ϵ, 0.0) ~ 0.0, # outer equatorial point
       u(t,1 - ϵ, 0.0) ~ 0.0, # inner equatorial point
       u(t,1 - δ*ϵ, κ*ϵ) ~ 0.0, # high point
       Dx(u(t,1 - δ*ϵ, κ*ϵ)) ~ 0.0, # high point maximum
       Dyy(u(t,1 + ϵ, 0.0)) + N₁*Dx(u(t,1 + ϵ, 0.0)) ~ 0.0, # outer equatorial point curvature
       Dyy(u(t,1 - ϵ, 0.0)) + N₂*Dx(u(t,1 - ϵ, 0.0)) ~ 0.0, # inner equatorial point curvature
       Dxx(u(t,1 - δ*ϵ, κ*ϵ)) + N₃*Dy(u(t,1 - δ*ϵ, κ*ϵ)) ~ 0.0] # High point curvature

# Space and time domains
domains = [t ∈ IntervalDomain(0.0,Inf),
           x ∈ IntervalDomain(0.5,1.5),
           y ∈ IntervalDomain(-1.0,1.0)]

pde_system = PDESystem([eq],bcs,domains,[t,x,y],[u(t,x,y)])

# Discretization
dx = 0.1
dy = 0.1
discretization = MOLFiniteDifference([x=>dx,y=>dy],t)
prob = SteadyStateProblem(discretize(pde_system,discretization))
@time sol = solve(prob,SSRootfind());

x = 0.5:dx:1.5
y = -1:dy:1
M = zeros(length(x),length(y))
M[2:end-1,2:end-1] .= reshape(sol.u,length(x)-2,length(y)-2)
contourf(x,y,M')
colorbar()

ArgumentError: invalid index: nothing of type Nothing

Stacktrace:
  [1] to_index(i::Nothing)
    @ Base ./indices.jl:300
  [2] to_index(A::Vector{Pair{Num, Vector{Num}}}, i::Nothing)
    @ Base ./indices.jl:277
  [3] to_indices
    @ ./indices.jl:333 [inlined]
  [4] to_indices
    @ ./indices.jl:325 [inlined]
  [5] getindex(A::Vector{Pair{Num, Vector{Num}}}, I::Nothing)
    @ Base ./abstractarray.jl:1170
  [6] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{Vector{Pair{Num, Float64}}, Num})
    @ DiffEqOperators ~/.julia/packages/DiffEqOperators/DMNmH/src/MOLFiniteDifference/MOL_discretization.jl:128
  [7] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{Vector{Pair{Num, Float64}}, Num})
    @ DiffEqOperators ~/.julia/packages/DiffEqOperators/DMNmH/src/MOLFiniteDifference/MOL_discretization.jl:225
  [8] top-level scope
    @ In[8]:49
  [9] eval
    @ ./boot.jl:360 [inlined]
 [10] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
    @ Base ./loading.jl:1094