Solving Grad-Shafranov Equation with ModelingToolkit

General Usage
1

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)
   @ Base ./loading.jl:1094

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.

2

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

3

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.

4

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

5

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?
6

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

7

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.

8

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.

9

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
10

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)
11

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
prob = SteadyStateProblem(discretize(pde_system,discretization))
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()

image

2 Likes
12

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

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
14

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)
prob = SteadyStateProblem(discretize(pde_system,discretization))
@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)
   @ Base ./loading.jl:1094
15

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:

image
image1484×998 144 KB

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?

16

https://diffeq.sciml.ai/stable/solvers/steady_state_solve/

1 Like
17

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 :wink:

1 Like
18

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:

image
image1528×1009 115 KB

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

19

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.

20

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