Solving Grad-Shafranov Equation with ModelingToolkit

lstagner · May 4, 2021, 11:36am

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.

mihalybaci · May 4, 2021, 1:13pm

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

lstagner · May 4, 2021, 5:35pm

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.

ChrisRackauckas · May 9, 2021, 10:50am

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

lstagner · May 11, 2021, 8:41am

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?

ChrisRackauckas · May 13, 2021, 1:22am

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

lstagner · May 13, 2021, 3:09am

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.

ChrisRackauckas · May 13, 2021, 5:46am

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.

lstagner · May 13, 2021, 6:14am

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

ChrisRackauckas · May 13, 2021, 6:34am

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)

lstagner · May 13, 2021, 7:01am

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

ChrisRackauckas · May 13, 2021, 8:03am

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

lstagner · May 13, 2021, 8:14am

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

lstagner · May 13, 2021, 8:55am

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

BLI · May 13, 2021, 9:41am

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?

ChrisRackauckas · May 13, 2021, 12:05pm

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

ChrisRackauckas · May 13, 2021, 12:06pm

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

BLI · May 14, 2021, 8:24am

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

BLI · May 14, 2021, 8:28am

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.

ChrisRackauckas · May 14, 2021, 1:30pm

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

Topic		Replies	Views
MethodError: no method matching -(::SymbolicUtils.BasicSymbolic{Real}, ::Symbolics.Equation) Modelling & Simulations package , modelingtoolkit	2	420	October 16, 2023
MethodError: no method matching *(::ODESystem, ::Float64) Modelling & Simulations modelingtoolkit	1	329	October 29, 2023
Using ModelingToolkit.jl to model a matrix ODE Modelling & Simulations question , diffeq , modelingtoolkit	3	941	September 17, 2022
MethodOfLines.jl not working General Usage package	6	233	July 5, 2024
ModelingToolkit \| Error with uncertain parameters using Measurements.jl Modelling & Simulations modelingtoolkit	7	309	November 2, 2023

Solving Grad-Shafranov Equation with ModelingToolkit

Related topics