MOL discretization start

Author: ChrisRackauckas

Project.toml

@@ -11,6 +11,7 @@ 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"

src/DiffEqOperators.jl

@@ -33,10 +33,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 +51,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 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 = 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