PDE infrastructure

Author: ChrisRackauckas

Project.toml

@@ -1,27 +1,28 @@
-name = "ModelingToolkit"
-uuid = "961ee093-0014-501f-94e3-6117800e7a78"
-authors = ["Chris Rackauckas <accounts@chrisrackauckas.com>"]
-version = "0.7.2"
-
-[deps]
-DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
-DiffRules = "b552c78f-8df3-52c6-915a-8e097449b14b"
-Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
-DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
-GeneralizedGenerated = "6b9d7cbe-bcb9-11e9-073f-15a7a543e2eb"
-LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
-NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
-SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
-SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
-StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
-
-[compat]
-julia = "1"
-
-[extras]
-OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-
-[targets]
-test = ["OrdinaryDiffEq", "Test"]
+name = "ModelingToolkit"
+uuid = "961ee093-0014-501f-94e3-6117800e7a78"
+authors = ["Chris Rackauckas <accounts@chrisrackauckas.com>"]
+version = "0.7.2"
+
+[deps]
+DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
+DiffRules = "b552c78f-8df3-52c6-915a-8e097449b14b"
+Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
+DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+GeneralizedGenerated = "6b9d7cbe-bcb9-11e9-073f-15a7a543e2eb"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
+NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+
+[compat]
+julia = "1"
+
+[extras]
+DiffEqOperators = "9fdde737-9c7f-55bf-ade8-46b3f136cc48"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[targets]
+test = ["DiffEqOperators", "OrdinaryDiffEq", "Test"]

src/ModelingToolkit.jl

@@ -89,6 +89,7 @@ include("function_registration.jl")
 include("simplify.jl")
 include("utils.jl")
 include("direct.jl")
+include("domains.jl")
 include("systems/diffeqs/diffeqsystem.jl")
 include("systems/diffeqs/first_order_transform.jl")
 include("systems/nonlinear/nonlinear_system.jl")

src/domains.jl

@@ -0,0 +1,26 @@
+abstract type AbstractDomain{T,N} 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)
+
+## Specific Domains
+
+struct IntervalDomain{T} <: AbstractDomain{T,1}
+  lower::T
+  upper::T
+end
+
+
+struct ProductDomain{D,T,N} <: AbstractDomain{T,N}
+  domains::D
+end
+⊗(args::AbstractDomain{T}...) where T = ProductDomain{typeof(args),T,length(args)}(args)
+
+struct CircleDomain <: AbstractDomain{Float64,2}
+  polar::Bool
+  CircleDomain(polar=false) = new(polar)
+end

src/equations.jl

@@ -40,3 +40,8 @@ Equation(x() - y(), ModelingToolkit.Constant(0))
 Base.:~(lhs::Expression, rhs::Expression) = Equation(lhs, rhs)
 Base.:~(lhs::Expression, rhs::Number    ) = Equation(lhs, rhs)
 Base.:~(lhs::Number    , rhs::Expression) = Equation(lhs, rhs)
+
+struct ConstrainedEquation
+  constraints
+  eq
+end

src/function_registration.jl

@@ -54,3 +54,6 @@ end
 
 # special cases
 Base.:^(x::Expression,y::T) where T <: Integer = Operation(Base.:^, Expression[x, y])
+
+@register Base.getindex(x,i)
+Base.getindex(x::Operation,i::Int64) = Operation(getindex,[x,i])

src/systems/pde/pdesystem.jl

@@ -0,0 +1,7 @@
+struct PDESystem <: ModelingToolkit.AbstractSystem
+  eq
+  bcs
+  domain
+  indvars
+  depvars
+end

test/constraints.jl

@@ -0,0 +1,7 @@
+using ModelingToolkit, DiffEqOperators, DiffEqBase, LinearAlgebra
+
+# Define some variables
+@parameters t x y
+@variables u(..)
+
+ConstrainedEquation([x ~ 0,y < 1/2], u(t,x,y) ~ x + y^2)

test/domains.jl

@@ -0,0 +1,13 @@
+using ModelingToolkit
+
+@parameters t x
+domains = [t ∈ IntervalDomain(0.0,1.0),
+           x ∈ IntervalDomain(0.0,1.0)]
+
+@parameters z
+z ∈ (IntervalDomain(0.0,1.0) ⊗ IntervalDomain(0.0,1.0))
+
+@parameters y
+(x,y) ∈ CircleDomain()
+@parameters r θ
+(r,θ) ∈ CircleDomain(true)

test/pde.jl

@@ -0,0 +1,28 @@
+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)
+
+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]))
+=#