Building off of the story in https://github.com/JuliaDiffEq/DifferentialEquation…s.jl/issues/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:
```julia
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:
- [x] Add domains and VarDomainPairing to ModelingToolkit
- [x] Add ConstrainedEquation to ModelingToolkit
- [x] Implement PDESystem in ModelingToolkit
- [x] Implement LinearProblem and NonlinearProblem in DiffEqBase
- [ ] Implement some basic DEAlgorithms for it. Solve dispatches on DEAlgorithms that just do `\`, GMRES (`LinSolveGMRES()`), `nlsolve`
- [x] Implement PDEProblem, PDETimeseriesSolution, PDENoTimeseriesSolution
- [x] Create `discretize` stub function in DiffEqBase
- [x] 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.
