Automated PDE solving fit for DifferentialEquations.jl

[ChrisRackauckas]

Building off of the story in #260 , DiffEqOperators getting auto finite difference discretizations likely done this summer. Thus the next step is automating the PDE discretization from a high level description. Since ModelingToolkit.jl is well-developed, it makes sense to utilize that to give the high level description. Thus I put a prototype together:

using ModelingToolkit, DiffEqOperators, DiffEqBase, LinearAlgebra

# Define some variables
@parameters t x
@variables u(..)
@derivatives Dt'~t
@derivatives Dxx''~x
eq  = Dt(u(t,x)) ~ Dxx(u(t,x))
bcs = [u(0,x) ~ - x * (x-1) * sin(x),
           u(t,0) ~ 0, u(t,1) ~ 0]

abstract type AbstractDomain{T,N} end
struct IntervalDomain{T} <: AbstractDomain{T,1}
  lower::T
  upper::T
end

struct VarDomainPairing
  variables
  domain::AbstractDomain
end
Base.:∈(variable::ModelingToolkit.Operation,domain::AbstractDomain) = VarDomainPairing(variable,domain)
Base.:∈(variables::NTuple{N,ModelingToolkit.Operation},domain::AbstractDomain) where N = VarDomainPairing(variables,domain)

domains = [t ∈ IntervalDomain(0.0,1.0),
                   x ∈ IntervalDomain(0.0,1.0)]

#########################################

# Note: Variables can be multidimensional
struct ProductDomain{D,T,N} <: AbstractDomain{T,N}
  domains::D
end
⊗(args::AbstractDomain{T}...) where T = ProductDomain{typeof(args),T,length(args)}(args)
@parameters z
z ∈ (IntervalDomain(0.0,1.0) ⊗ IntervalDomain(0.0,1.0))
@register Base.getindex(x,i)
Base.getindex(x::Operation,i::Int64) = Operation(getindex,[x,i])
z[1]^2 + z[2]^2 ~ 1
z[1]^2 + z[2]^2 < 1

struct CircleDomain <: AbstractDomain{Float64,2}
  polar::Bool
  CircleDomain(polar=false) = new(polar)
end
z ∈ CircleDomain()
@parameters y
(x,y) ∈ CircleDomain()
@parameters r θ
(r,θ) ∈ CircleDomain(true)

# Use constrained equations for more detailed BCs
struct ConstrainedEquation
  constraints
  eq
end
ConstrainedEquation([x ~ 0,y < 1/2], u(t,x,y) ~ x + y^2)

#########################################

struct PDESystem <: ModelingToolkit.AbstractSystem
  eq
  bcs
  domain
  indvars
  depvars
end
pdesys = PDESystem(eq,bcs,domains,[t,x],[u])

abstract type AbstractDiscretization end
struct MOLFiniteDifference{T} <: AbstractDiscretization
  dxs::T
  order::Int
end
MOLFiniteDifference(args...;order=2) = MOLFiniteDifference(args,order)
discretization = MOLFiniteDifference(0.1)

struct PDEProblem{P,E,S} <: DiffEqBase.DEProblem
  prob::P
  extrapolation::E
  space::S
end
function Base.show(io::IO, A::PDEProblem)
  println(io,summary(A.prob))
  println(io)
end
Base.summary(prob::PDEProblem) = string(DiffEqBase.TYPE_COLOR, 
                                                                     nameof(typeof(prob)),
                                                                     DiffEqBase.NO_COLOR)

function discretize(pdesys,
                               discretization::MOLFiniteDifference)
  tdomain = pdesys.domain[1].domain
  domain = pdesys.domain[2].domain
  @assert domain isa IntervalDomain
  len = domain.upper - domain.lower
  dx = discretization.dxs[1]
  interior = domain.lower+dx:dx:domain.upper-dx
  X = domain.lower:dx:domain.upper
  L = CenteredDifference(2,2,dx,Int(len/dx)-2)
  Q = DiffEqOperators.DirichletBC([0.0,0.0],[1.0,1.0])
  function f(du,u,p,t)
    mul!(du,L,Array(Q*u))
  end
  u0 = @. - interior * (interior - 1) * sin(interior)
  PDEProblem(ODEProblem(f,u0,(tdomain.lower,tdomain.upper),nothing),Q,X)
end

prob = discretize(pdesys,discretization)

function DiffEqBase.solve(prob::PDEProblem,alg::DiffEqBase.DEAlgorithm,args...;
                                          kwargs...)
    solve(prob.prob,alg,args...;kwargs...)
end

using OrdinaryDiffEq
sol = solve(prob,Tsit5(),saveat=0.1)

using Plots
plot(prob.space,Array(prob.extrapolation*sol[1]))
plot!(prob.space,Array(prob.extrapolation*sol[2]))
plot!(prob.space,Array(prob.extrapolation*sol[3]))
plot!(prob.space,Array(prob.extrapolation*sol[4]))

Description of the code

The equation and boundary conditions are just given by ModelingToolkit.jl operations. The ConstrainedEquation allows for specifying equations only on specific parts. Then each independent variable needs a domain. Domains can be multi-dimensional, and variables can be multi-dimensional, so this lets us smartly define variables in domains (and in manifolds!) and equations using them.

The totality of the PDE is in the PDESystem, just a ModelingToolkit.AbstractSystem object. The workflow for solving a PDE is thus:

1. Define a PDESystem.
2. Discretize to a PDEProblem with an AbstractDiscretization. This step will kick out things like a PDEProblem that wraps an ODEProblem, or SDEProblem, or even LinearProblem and NonlinearProblem.
3. Solve the PDEProblem.

So (1) is all ModelingToolkit.jl. (2) is just a discretize function and some problem types. We will need to add LinearProblem and NonlinearProblem with some solve dispatches for this. After (3), we will want to wrap the solution so it can get a nicer plot recipe.

Thus the steps are:

• Add domains and VarDomainPairing to ModelingToolkit
• Add ConstrainedEquation to ModelingToolkit
• Implement PDESystem in ModelingToolkit
• Implement LinearProblem and NonlinearProblem in DiffEqBase
• Implement some basic DEAlgorithms for it. Solve dispatches on DEAlgorithms that just do \, GMRES (LinSolveGMRES()), nlsolve
• Implement PDEProblem, PDETimeseriesSolution, PDENoTimeseriesSolution
• Create discretize stub function in DiffEqBase
• Implement a smart discretize function in DiffEqOperators.jl for automatic finite difference discretizations. This is where a ton of work will have to go, since it will need to parse the equation and find out what CenteredDifference and UpwindDifference operators to create. This work will continue after this closes.

[ChrisRackauckas]

@YingboMa @HarrisonGrodin solidifying the discussion here.

[xtalax]

Thinking about some ways to automatically handle boundary condition operator generation for polar circle domains down in DiffEqOperators.jl - the boundary condition on theta can simply be periodic, but the condition on r is a bit more difficult, The higher index end can be arbitrarily chosen, but it is the lower index end i'm concerned with here.

An important question is what the lower index end of r should correspond to: r=0 or r=dr.
In the first case it would be necessary to ensure somehow that the lower index points of r for all theta contained the same value.
Would it be sufficient to place in the ghost nodes along this boundary the value at r[2] at the diametrically opposed value of theta? In this case it would be necessary to have a uniformly sampled theta with an even number of points.
In the second case, where r[1] = dr[1], the value at r=0 could be interpolated from the rest of r, and this value placed in the ghost nodes along this boundary - this may be preferable because then theta could be non uniform.

[ChrisRackauckas]

That's a question for DiffEqOperators. It has the information to do the discretization, and it has to find out how it wants to do that. However, this is all about making sure there's general enough information such that it can make the correct choice, and I believe that is the case.

[tpdsantos]

I'm anxiously waiting that this gets finished. I practically only solve PDEs in my work, so I usually discretize them using some Collocation Method, transforming it into a DAE and solve it with the IDA solver. The problem is that the only linear solver that does not give me IC errors is the Dense solver, which is not appropriate at all for my problems. When this is done will I be able to solve my problems more efficiently?

[ChrisRackauckas]

Hopefully it'll set them up in a very efficient manner haha. That's the goal

[tpdsantos]

That would be amazing! I actually wanted to use Modia to do process simulations but unfortunately they cannot handle, yet, sparse jacobians and event handling, things that I absolutely need. Looking forward to see this development!

[ChrisRackauckas]

SciML/DiffEqBase.jl#342
SciML/DiffEqOperators.jl#179
SciML/ModelingToolkit.jl#180

Makes this work:

using ModelingToolkit, DiffEqOperators, DiffEqBase, LinearAlgebra

# Define some variables
@parameters t x
@variables u(..)
@derivatives Dt'~t
@derivatives Dxx''~x
eq  = Dt(u(t,x)) ~ Dxx(u(t,x))
bcs = [u(0,x) ~ - x * (x-1) * sin(x),
           u(t,0) ~ 0, u(t,1) ~ 0]

domains = [t ∈ IntervalDomain(0.0,1.0),
           x ∈ IntervalDomain(0.0,1.0)]

pdesys = PDESystem(eq,bcs,domains,[t,x],[u])
discretization = MOLFiniteDifference(0.1)
prob = discretize(pdesys,discretization) # This gives an ODEProblem since it's time-dependent

using OrdinaryDiffEq
sol = solve(prob,Tsit5(),saveat=0.1)

using Plots
plot(prob.space,Array(prob.extrapolation*sol[1]))
plot!(prob.space,Array(prob.extrapolation*sol[2]))
plot!(prob.space,Array(prob.extrapolation*sol[3]))
plot!(prob.space,Array(prob.extrapolation*sol[4]))
savefig("pde_solve.png")

[ChrisRackauckas]

Implement a smart discretize function in DiffEqOperators.jl for automatic finite difference discretizations. This is where a ton of work will have to go, since it will need to parse the equation and find out what CenteredDifference and UpwindDifference operators to create. This work will continue after this closes.

This is now SciML/DiffEqOperators.jl#180

[SayedAIrfan]

Hi Chris
I am trying to run the above equation to learn about the PDE solution using NN. I am running the code in JupyterLab I am getting the following error. Can you kindly help to understand the error and how to resolve it.

**MethodError: no method matching one(::Type{Array{Float64,1}})
Closest candidates are:
one(!Matched::Type{Missing}) at missing.jl:103
one(!Matched::BitArray{2}) at bitarray.jl:422
one(!Matched::Missing) at missing.jl:100
...

Stacktrace:
[1] DirichletBC(::Array{Float64,1}, ::Array{Float64,1}) at C:\Users\Administrator.julia\packages\DiffEqOperators\QiB3a\src\derivative_operators\BC_operators.jl:140
[2] discretize(::PDESystem, ::MOLFiniteDifference{Tuple{Float64}}) at .\In[2]:101
[3] top-level scope at In[2]:108**

Best Regards
Irfan