# Solving Grad-Shafranov Equation with ModelingToolkit

I am trying to solve the Grad-Shafranov partial differential equation and I am running into the following error

``````using ModelingToolkit

R0 = 2.0
C = 1.0
A = 1.0
#ψ(r,z) = (1/2)*(A*R0^2 +a0*r^2)*z^2 + (1/8)*(C−a0)(r^2−R0^2)^2

@parameters r z
@variables ψ(..)

Dr = Differential(r)
Drr = Differential(r)^2
Dzz = Differential(z)^2

eq  = Drr(ψ(r,z)) - Dr(ψ(r,z))/r + Dzz(ψ(r,z)) ~ A*R0^2 + C*r^2
bcs = [ψ(0,z) ~ 0,
ψ(4,z) ~ 0,
ψ(r,-2) ~ 0,
ψ(r,2) ~ 0]

domains = [r ∈ IntervalDomain(0.0,4.0),
z ∈ IntervalDomain(-2.0,2.0)]

pde_system = PDESystem(eq,bcs,domains,[r,z],[ψ])
``````
``````MethodError: no method matching ^(::Differential, ::Int64)
Closest candidates are:
^(::Union{AbstractChar, AbstractString}, ::Integer) at strings/basic.jl:718
^(::Union{VectorizationBase.AbstractSIMDVector{W, T}, VectorizationBase.VecUnroll{var"#s2", W, T, V} where {var"#s2", V<:VectorizationBase.AbstractSIMDVector{W, T}}}, ::Union{Int16, Int32, Int64, Int8, UInt16, UInt32, UInt64, UInt8}) where {W, T<:Union{Float32, Float64}} at /home/lstagner/.julia/packages/VectorizationBase/qmYqb/src/special/misc.jl:2
^(::Union{VectorizationBase.AbstractSIMDVector{W, T}, VectorizationBase.VecUnroll{var"#s2", W, T, V} where {var"#s2", V<:VectorizationBase.AbstractSIMDVector{W, T}}}, ::Union{Int16, Int32, Int64, Int8, UInt16, UInt32, UInt64, UInt8}) where {W, T} at /home/lstagner/.julia/packages/VectorizationBase/qmYqb/src/special/misc.jl:1
...

Stacktrace:
[1] macro expansion
@ ./none:0 [inlined]
[2] literal_pow(f::typeof(^), x::Differential, #unused#::Val{2})
@ Base ./none:0
[3] top-level scope
@ In[4]:11
[4] eval
@ ./boot.jl:360 [inlined]
[5] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
``````

Also, I don’t know how to progress past making the `PDESystem` All the tutorials were about solving the problem with PINNS but I want to solve it the regular way.

Your example seems to work for me using Julia v1.6.1 and ModelingToolkit v5.16.0.

I was using v1.6 so I updated to v1.6.1 but it looks like my ModelingToolkit is at v3.20.1 which is strange since I just added it. Turns out one of my dev package I haven’t used in a while was holding the whole environment back. After I freed the dev package I was able to update and the error went away.

However, after I create the PDESystem I still don’t know what to do next as all the examples I found were for using PINNS. I just want to solve it the normal way.

See the DiffEqOperators.jl documentation for a finite difference way, but note that is still in development.

The examples listed in the DiffEqOperators are all functions of time and use the Method of Lines (MOL) Operator with the continuous variable being time. After discretizing the next step is to put it into an ode solver but I cant do that for the Grad-Shafranov Equation because it is independent of time. A central differencing discretization is more appropriate. I just don’t know how to proceed.

I’ve given up on the Grad-Shafranov for now and I’m now trying to solve the 2D Poisson equation from the NeuralPDE docs using finite differencing.

``````using ModelingToolkit, DiffEqOperators
@parameters x y
@variables u(..)
Dxx = Differential(x)^2
Dyy = Differential(y)^2

# 2D PDE
eq  = Dxx(u(x,y)) + Dyy(u(x,y)) ~ -sin(pi*x)*sin(pi*y)

# Boundary conditions
bcs = [u(0,y) ~ 0.f0, u(1,y) ~ -sin(pi*1)*sin(pi*y),
u(x,0) ~ 0.f0, u(x,1) ~ -sin(pi*x)*sin(pi*1)]
# Space and time domains
domains = [x ∈ IntervalDomain(0.0,1.0),
y ∈ IntervalDomain(0.0,1.0)]

# Discretization
dx = 0.05
discretization = CenteredDifference{2}(2,2,dx) #<<<This is wrong

pde_system = PDESystem(eq,bcs,domains,[x,y],[u])
prob = discretize(pde_system,discretization)

# Once I have the discretized problem how do I solve it?
``````

DiffEqOperators `MOLFiniteDifference` only supports time-dependent problems right now. Time independent is an open issue: MOLFiniteDifference: NonlinearProblem · Issue #282 · SciML/DiffEqOperators.jl · GitHub

So I tried something stupid but it seems to work. I basically added a dummy time variable to the system.

``````using ModelingToolkit, DiffEqOperators, OrdinaryDiffEq

@parameters x y t
@variables u(..)

Dxx = Differential(x)^2
Dyy = Differential(y)^2

# 2D PDE
eq  = Dxx(u(t,x,y)) + Dyy(u(t,x,y)) ~ -sin(pi*x)*sin(pi*y)

# Boundary conditions
bcs = [u(0,x,y) ~ -sin(pi*x)*sin(pi*y),
u(t,0,y) ~ 0.f0, u(t,1,y) ~ -sin(pi*1)*sin(pi*y),
u(t,x,0) ~ 0.f0, u(t,x,1) ~ -sin(pi*x)*sin(pi*1)]

# Space and time domains
domains = [t ∈ IntervalDomain(0.0,1.0),
x ∈ IntervalDomain(0.0,1.0),
y ∈ IntervalDomain(0.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 = discretize(pde_system,discretization)
sol = solve(prob,Rodas5())[1]

9-element Vector{Float64}:
0.009277825067485276
0.012769830684160145
0.024289661368330728
0.015784378786771757
0.030023672601675946
0.012769830684142692
0.024289661368289295
0.00927782506745599
0.017647471974907916
``````

However the size of the output array doesn’t make sense excluding the boundaries it should have 81 elements not 9 (or 121 including the boundaries). Also, if I try to decrease the step size to `dx=0.05; dy=0.05` the `discretize` call hangs.

Yeah I’m not quite sure how it would interpret that. If you did:

``````eq  = Dt(u(t,x,y)) ~ -Dxx(u(t,x,y)) - Dyy(u(t,x,y)) -sin(pi*x)*sin(pi*y)
``````

then the solution at t->infty (i.e. SteadyStateProblem) would be the solution to the original elliptic PDE. That’s a hack that should work for now, but I think your original form just needs the issue to be solved.

Yeah, it doesn’t like that. It complains that the problem is unstable

``````┌ Warning: Instability detected. Aborting
└ @ SciMLBase /home/lstagner/.julia/packages/SciMLBase/XuLdB/src/integrator_interface.jl:351
``````

I may have had the sign flipped?

``````eq  = Dt(u(t,x,y)) ~ Dxx(u(t,x,y)) + Dyy(u(t,x,y)) -sin(pi*x)*sin(pi*y)
``````

Almost

``````eq  = Dt(u(t,x,y)) ~ Dxx(u(t,x,y)) + Dyy(u(t,x,y)) + sin(pi*x)*sin(pi*y)
``````

Anyway its working now. Not the most efficient around 8 secs (edit: 2 sec when actually solving the steady state equation correctly. Also the discretize function doesn’t scale well with decreased spatial step)

``````using ModelingToolkit, DiffEqOperators, DifferentialEquations, SteadyStateDiffEq, PyPlot
@parameters x y t
@variables u(..)
Dxx = Differential(x)^2
Dyy = Differential(y)^2
Dt = Differential(t)

# 2D PDE
eq  = Dt(u(t,x,y)) ~ Dxx(u(t,x,y)) + Dyy(u(t,x,y)) + sin(pi*x)*sin(pi*y)

# Boundary conditions
bcs = [u(0,x,y) ~ 0.0,
u(t,0,y) ~ 0.0, u(t,1,y) ~ 0.0,
u(t,x,0) ~ 0.0, u(t,x,1) ~ 0.0]

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

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

# Discretization
dx = 0.05
dy = 0.05
discretization = MOLFiniteDifference([x=>dx,y=>dy],t)
#prob = discretize(pde_system,discretization)
#sol = solve(prob,DynamicSS(Tsit5())); # slow
sol = solve(prob,SSRootfind());

x = 0:dx:1
y = 0:dy:1
M = zeros(length(x),length(y))
M[2:end-1,2:end-1] = sol
contourf(x,y,M)
colorbar()
``````

2 Likes

Yeah, not the best way to do it but it works, and we’ll work on the nonlinear system style soon.

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)
@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)
``````

Trying with a different set of boundary conditions produces a different problem

``````using ModelingToolkit, DiffEqOperators, DifferentialEquations, SteadyStateDiffEq, PyPlot

α = -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.0, x, y) ~ 0.0, # initial guess
u(t, 0.0, y) ~ 0.0,
u(t, 1.0, y) ~ 0.0,
u(t, x,-1.0) ~ 0.0,
u(t, x, 1.0) ~ 0.0,
Dy(u(t, 0.0, y)) ~ 0.0,
Dy(u(t, 1.0, y)) ~ 0.0,
Dx(u(t, x, -1.0)) ~ 0.0,
Dx(u(t, x,  1.0)) ~ 0.0]

# Space and time domains
domains = [t ∈ IntervalDomain(0.0,Inf),
x ∈ IntervalDomain(0.0,1.0),
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)
@time sol = solve(prob,SSRootfind());

x = 0:dx:1
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()
``````
``````MethodError: no method matching -(::Float64, ::Term{Vector{Num}, Nothing})
Closest candidates are:
-(::VectorizationBase.CartesianVIndex, ::Any) at /home/lstagner/.julia/packages/VectorizationBase/BPiSF/src/cartesianvindex.jl:58
-(::PyCall.PyObject, ::Any) at /home/lstagner/.julia/packages/PyCall/BD546/src/pyoperators.jl:13
-(::Num, ::Term) at /home/lstagner/.julia/packages/SymbolicUtils/9iQGH/src/methods.jl:49
...

Stacktrace:
[1] initialize_system_structure(sys::ODESystem)
@ ModelingToolkit.SystemStructures ~/.julia/packages/ModelingToolkit/DOKSJ/src/systems/systemstructure.jl:112
[2] alias_elimination(sys::ODESystem)
@ ModelingToolkit ~/.julia/packages/ModelingToolkit/DOKSJ/src/systems/alias_elimination.jl:6
[3] structural_simplify(sys::ODESystem)
@ ModelingToolkit ~/.julia/packages/ModelingToolkit/DOKSJ/src/systems/abstractsystem.jl:642
[4] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{Vector{Pair{Num, Float64}}, Num})
@ DiffEqOperators ~/.julia/packages/DiffEqOperators/DMNmH/src/MOLFiniteDifference/MOL_discretization.jl:226
[5] top-level scope
@ In[24]:50
[6] eval
@ ./boot.jl:360 [inlined]
[7] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
``````

Being interested in the dynamics of the problem, I’ve stripped off the steady state part and done:

``````using ModelingToolkit, DiffEqOperators
using DifferentialEquations
using Plots
#
@parameters x y t
@variables u(..)
Dxx = Differential(x)^2
Dyy = Differential(y)^2
Dt = Differential(t)

# 2D PDE
eq  = Dt(u(t,x,y)) ~ Dxx(u(t,x,y)) + Dyy(u(t,x,y)) + sin(pi*x)*sin(pi*y)

# Boundary conditions
bcs = [u(0,x,y) ~ 0.0,
u(t,0,y) ~ 0.0, u(t,1,y) ~ 0.0,
u(t,x,0) ~ 0.0, u(t,x,1) ~ 0.0]

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

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

# Discretization
dx = 0.05
dy = 0.05
discretization = MOLFiniteDifference([x=>dx,y=>dy],t)
prob = discretize(pde_system,discretization)
@time sol = solve(prob,DynamicSS(Tsit5())); # slow
``````

[What does solver `DynamicSS(Tsit5())` imply?]

This runs in ca. 4.2 seconds on my low power i9 processor.

Is there some documentation on how to plot the result? When I do `plot(sol)`, the result is:

Suppose I want to plot temporal snap-shots, e.g., for every t in `0:0.1:1`, suppose I want to plot a surface plot of u(x,y) in order to generate an animation… Is there an example that would explain how to do this?

1 Like

Not yet, all of the PDE tools are still under heavy development and YMMV. At least for the next year or so, they will not be plug-and-play and will be a fun ride

1 Like

OK – the following works for now:

``````x = 0.05:0.05:0.95
y = x
n = length(x)
#
anim = Animation()
for t in 0:0.01:1
Sol = reshape(sol(t),n,n)
plot(x,y,Sol; st=:surface,xlim=(0,1),ylim=(0,1),zlim=(0,0.05))
frame(anim)
end
#
gif(anim,"anim.gif";fps=15)
``````

[(i) The boundary doesn’t seem to be part of `sol`. (ii) keyword `fps` doesn’t seem to have any effect on the playback.]

A snapshot of the animation looks like:

Anyway, I realize that the solution layout, etc. may change during the next year or so

Seems like `solve(prob,Tsit5())` also works (instead of `solve(prob,DynamicSS(Tsit5()))` I’d guess that the `DynamicSS()` construct says that the solution should approach a steady state, or something. And that perhaps that construct should not be used if the solution converges to some stationary oscillation, or non-steady state solutions.

Yes it’s a solver for a SteadyStateProblem, along with rootfinding methods and pseudo-transient.