Merge pull request #179 from JuliaDiffEq/MOL

MOL discretization start

Author: ChrisRackauckas

Project.toml

@@ -11,12 +11,14 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LazyArrays =  "5078a376-72f3-5289-bfd5-ec5146d43c02"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NNlib = "872c559c-99b0-510c-b3b7-b6c96a88d5cd"
+ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 SuiteSparse = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"
 
 [compat]
-DiffEqBase = ">= 5.19.0"
+DiffEqBase = "6"
+ModelingToolkit = "0.8.0"
 julia = "1"
 
 [extras]

src/DiffEqOperators.jl

@@ -5,6 +5,7 @@ using DiffEqBase, StaticArrays, LinearAlgebra
 import LinearAlgebra: mul!, ldiv!, lmul!, rmul!, axpy!, opnorm, factorize, I
 import DiffEqBase: AbstractDiffEqLinearOperator, update_coefficients!, is_constant
 using SparseArrays, ForwardDiff, BandedMatrices, NNlib, LazyArrays, BlockBandedMatrices
+using ModelingToolkit
 
 abstract type AbstractDerivativeOperator{T} <: AbstractDiffEqLinearOperator{T} end
 abstract type AbstractDiffEqCompositeOperator{T} <: AbstractDiffEqLinearOperator{T} end
@@ -33,10 +34,11 @@ include("derivative_operators/concretization.jl")
 include("derivative_operators/ghost_derivative_operator.jl")
 include("derivative_operators/derivative_operator_functions.jl")
 
-
 ### Composite Operators
 include("composite_operators.jl")
 
+include("MOL_discretization.jl")
+
 # The (u,p,t) and (du,u,p,t) interface
 for T in [DiffEqScaledOperator, DiffEqOperatorCombination, DiffEqOperatorComposition]
   (L::T)(u,p,t) = (update_coefficients!(L,u,p,t); L * u)
@@ -50,6 +52,6 @@ export AbstractDerivativeOperator, DerivativeOperator,
 export DirichletBC, Dirichlet0BC, NeumannBC, Neumann0BC, RobinBC, GeneralBC, MultiDimBC, PeriodicBC,
        MultiDimDirectionalBC, ComposedMultiDimBC,
        compose, decompose, perpsize
-
 export GhostDerivativeOperator
+export MOLFiniteDifference
 end # module

src/MOL_discretization.jl

@@ -0,0 +1,22 @@
+struct MOLFiniteDifference{T} <: DiffEqBase.AbstractDiscretization
+  dxs::T
+  order::Int
+end
+MOLFiniteDifference(args...;order=2) = MOLFiniteDifference(args,order)
+
+function DiffEqBase.discretize(pdesys::PDESystem,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 = DirichletBC(0.0,0.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

test/MOLtest.jl

@@ -0,0 +1,20 @@
+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)

test/runtests.jl

@@ -17,4 +17,5 @@ import Base: isapprox
 @time @safetestset "Differentiation Dimension" begin include("differentiation_dimension.jl") end
 @time @safetestset "2D and 3D fast multiplication" begin include("2D_3D_fast_multiplication.jl") end
 @time @safetestset "Higher Dimensional Concretization" begin include("concretization.jl") end
-@time @safetestset "Upwind Operator Interface" begin include("upwind_operators_interface.jl") end
+@time @safetestset "Upwind Operator Interface" begin include("upwind_operators_interface.jl") end
+@time @safetestset "MOLFiniteDifference Interface" begin include("MOLtest.jl") end