# Solving Grad-Shafranov Equation with ModelingToolkit

**URL:** <https://discourse.julialang.org/t/solving-grad-shafranov-equation-with-modelingtoolkit/60513>\
**Category:** General Usage\
**Created:** [May 4, 2021, 11:36am UTC](https://discourse.julialang.org/t/solving-grad-shafranov-equation-with-modelingtoolkit/60513 "2021-05-04T11:36:01Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![lstagner](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lstagner/32/448_2.png) [@lstagner](https://discourse.julialang.org/u/lstagner)\
**Post date:** [May 4, 2021, 11:36am UTC](https://discourse.julialang.org/t/solving-grad-shafranov-equation-with-modelingtoolkit/60513/1 "2021-05-04T11:36:01Z")

</div>

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

```julia
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],[ψ])

```

```julia
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.

---

<div class="post-metadata">

**Author:** ![mihalybaci](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mihalybaci/32/13528_2.png) [@mihalybaci](https://discourse.julialang.org/u/mihalybaci)\
**Post date:** [May 4, 2021, 1:13pm UTC](https://discourse.julialang.org/t/solving-grad-shafranov-equation-with-modelingtoolkit/60513/2 "2021-05-04T13:13:44Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![lstagner](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lstagner/32/448_2.png) [@lstagner](https://discourse.julialang.org/u/lstagner)\
**Post date:** [May 4, 2021, 5:35pm UTC](https://discourse.julialang.org/t/solving-grad-shafranov-equation-with-modelingtoolkit/60513/3 "2021-05-04T17:35:35Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 9, 2021, 10:50am UTC](https://discourse.julialang.org/t/solving-grad-shafranov-equation-with-modelingtoolkit/60513/4 "2021-05-09T10:50:58Z")

</div>

> [@lstagner](#):
>
> 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.

---

<div class="post-metadata">

**Author:** ![lstagner](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lstagner/32/448_2.png) [@lstagner](https://discourse.julialang.org/u/lstagner)\
**Post date:** [May 11, 2021, 8:41am UTC](https://discourse.julialang.org/t/solving-grad-shafranov-equation-with-modelingtoolkit/60513/5 "2021-05-11T08:41:45Z")

</div>

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.

```julia
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?

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 13, 2021, 1:22am UTC](https://discourse.julialang.org/t/solving-grad-shafranov-equation-with-modelingtoolkit/60513/6 "2021-05-13T01:22:55Z")

</div>

DiffEqOperators `MOLFiniteDifference` only supports time-dependent problems right now. Time independent is an open issue: [MOLFiniteDifference: NonlinearProblem · Issue #282 · SciML/DiffEqOperators.jl · GitHub](https://github.com/SciML/DiffEqOperators.jl/issues/282)

---

<div class="post-metadata">

**Author:** ![lstagner](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lstagner/32/448_2.png) [@lstagner](https://discourse.julialang.org/u/lstagner)\
**Post date:** [May 13, 2021, 3:09am UTC](https://discourse.julialang.org/t/solving-grad-shafranov-equation-with-modelingtoolkit/60513/7 "2021-05-13T03:09:39Z")

</div>

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

```julia
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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 13, 2021, 5:46am UTC](https://discourse.julialang.org/t/solving-grad-shafranov-equation-with-modelingtoolkit/60513/8 "2021-05-13T05:46:30Z")

</div>

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

```julia
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.

---

<div class="post-metadata">

**Author:** ![lstagner](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lstagner/32/448_2.png) [@lstagner](https://discourse.julialang.org/u/lstagner)\
**Post date:** [May 13, 2021, 6:14am UTC](https://discourse.julialang.org/t/solving-grad-shafranov-equation-with-modelingtoolkit/60513/9 "2021-05-13T06:14:30Z")

</div>

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

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

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 13, 2021, 6:34am UTC](https://discourse.julialang.org/t/solving-grad-shafranov-equation-with-modelingtoolkit/60513/10 "2021-05-13T06:34:28Z")

</div>

I may have had the sign flipped?

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

```

---

<div class="post-metadata">

**Author:** ![lstagner](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lstagner/32/448_2.png) [@lstagner](https://discourse.julialang.org/u/lstagner)\
**Post date:** [May 13, 2021, 7:01am UTC](https://discourse.julialang.org/t/solving-grad-shafranov-equation-with-modelingtoolkit/60513/11 "2021-05-13T07:01:54Z")

</div>

Almost

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

```julia
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](https://global.discourse-cdn.com/julialang/original/3X/e/2/e25faadffc0ebc78039f9f5f695845c375357a41.png)

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 13, 2021, 8:03am UTC](https://discourse.julialang.org/t/solving-grad-shafranov-equation-with-modelingtoolkit/60513/12 "2021-05-13T08:03:20Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![lstagner](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lstagner/32/448_2.png) [@lstagner](https://discourse.julialang.org/u/lstagner)\
**Post date:** [May 13, 2021, 8:14am UTC](https://discourse.julialang.org/t/solving-grad-shafranov-equation-with-modelingtoolkit/60513/13 "2021-05-13T08:14:27Z")

</div>

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

```julia
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()

```

```julia
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

```

---

<div class="post-metadata">

**Author:** ![lstagner](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lstagner/32/448_2.png) [@lstagner](https://discourse.julialang.org/u/lstagner)\
**Post date:** [May 13, 2021, 8:55am UTC](https://discourse.julialang.org/t/solving-grad-shafranov-equation-with-modelingtoolkit/60513/14 "2021-05-13T08:55:53Z")

</div>

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

```julia
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()

```

```julia
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

```

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [May 13, 2021, 9:41am UTC](https://discourse.julialang.org/t/solving-grad-shafranov-equation-with-modelingtoolkit/60513/15 "2021-05-13T09:41:17Z")

</div>

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

```julia
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](https://global.discourse-cdn.com/julialang/original/3X/5/a/5a837d5670f0b5b3e6a27ca5cea0a47134c3ca39.png)

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?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 13, 2021, 12:05pm UTC](https://discourse.julialang.org/t/solving-grad-shafranov-equation-with-modelingtoolkit/60513/16 "2021-05-13T12:05:35Z")

</div>

> [@BLI](#):
>
> [What does solver `DynamicSS(Tsit5())` imply?]

[https://diffeq.sciml.ai/stable/solvers/steady\_state\_solve/](https://diffeq.sciml.ai/stable/solvers/steady_state_solve/)

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 13, 2021, 12:06pm UTC](https://discourse.julialang.org/t/solving-grad-shafranov-equation-with-modelingtoolkit/60513/17 "2021-05-13T12:06:16Z")

</div>

> [@BLI](#):
>
> Suppose I want to plot temporal snap-shots, e.g., for every t tt 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?

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 😉

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [May 14, 2021, 8:24am UTC](https://discourse.julialang.org/t/solving-grad-shafranov-equation-with-modelingtoolkit/60513/18 "2021-05-14T08:24:34Z")

</div>

OK – the following works for now:

```julia
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](https://global.discourse-cdn.com/julialang/original/3X/b/f/bf9ba4ea614a66ec71ee482c7e823fff6b2c5452.jpeg)

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

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [May 14, 2021, 8:28am UTC](https://discourse.julialang.org/t/solving-grad-shafranov-equation-with-modelingtoolkit/60513/19 "2021-05-14T08:28:08Z")

</div>

> [@ChrisRackauckas](#):
>
> > [@BLI](#):
> >
> > [What does solver `DynamicSS(Tsit5())` imply?]
> 
> [https://diffeq.sciml.ai/stable/solvers/steady\_state\_solve/](https://diffeq.sciml.ai/stable/solvers/steady_state_solve/)

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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 14, 2021, 1:30pm UTC](https://discourse.julialang.org/t/solving-grad-shafranov-equation-with-modelingtoolkit/60513/20 "2021-05-14T13:30:37Z")

</div>

> [@BLI](#):
>
> 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.