Merge pull request #17 from xtalax/arbitrary-approx-order

Arbitrary approx order

Author: ChrisRackauckas

.gitignore

@@ -4,3 +4,9 @@
 Manifest.toml
 .vscode
 .vscode/*
+plots
+plots/*
+src/scratch
+src/scratch/*
+test/plots
+test/plots/*

Project.toml

@@ -18,7 +18,7 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
-DiffEqOperators = "4"
+DiffEqOperators = "4.40.0"
 DomainSets = "0.5"
 ModelingToolkit = "8"
 NonlinearSolve = "0.3"

README.md

@@ -9,24 +9,24 @@ Feature requests and issues welcome.
 discretization = MOLFiniteDifference(dxs, 
                                       <your choice of continuous variable, usually time>; 
                                       upwind_order = <currently hard coded to 1>, 
-                                      centered_order = <currently hard coded to 2>, 
+                                      approx_order = Order of derivative approximation, starting from 2 
                                       grid_align = your grid type choice>)
 prob = discretize(pdesys, discretization)
 ```
 Where `dxs` is a vector of pairs of parameters to the grid step in this dimension, i.e. `[x=>0.2, y=>0.1]`
 
 Note that the second argument to `MOLFiniteDifference` is optional, all parameters can be discretized if all boundary conditions are specified
 
-Currently supported grid types: `center_align` and `edge_align`
+Currently supported grid types: `center_align` and `edge_align`. Edge align will give better accuracy with Neumann Boundary conditions.
 
 `center_align`: naive grid, starting from lower boundary, ending on upper boundary with step of `dx`
 
 `edge_align`: offset grid, set halfway between the points that would be generated with center_align, with extra points at either end that are above and below the supremum and infimum by `dx/2`.
 
-Currently boundary conditions defined in terms of derivatives at the boundary are unsupported above 1 discretized dimension. Periodic conditions are also unsupported.
+Currently Periodic boundary conditions are unsupported
 
 ## Coming soon:
-- Arbitrary approximation order
+- Automatic up/downwinding for odd order derivatives
 - Above mentioned unsupported boundary conditions supported
 
 ## Full Example:

src/MOLFiniteDifference.jl

@@ -0,0 +1,16 @@
+struct MOLFiniteDifference{G} <: DiffEqBase.AbstractDiscretization
+    dxs
+    time
+    approx_order::Int
+    upwind_order::Int
+    grid_align::G
+end
+
+# Constructors. If no order is specified, both upwind and centered differences will be 2nd order
+function MOLFiniteDifference(dxs, time=nothing; approx_order = 2, upwind_order = 1, grid_align=CenterAlignedGrid())
+    
+    if approx_order % 2 != 0
+        @warn "Discretization approx_order must be even, rounding up to $(approx_order+1)"
+    end
+    return MOLFiniteDifference{typeof(grid_align)}(dxs, time, approx_order, upwind_order, grid_align)
+end

src/discretization/MOL_discretization.jl renamed to src/MOL_discretization.jl

@@ -1,28 +1,15 @@
 # Method of lines discretization scheme
 
-@enum GridAlign center_align edge_align
-
-struct MOLFiniteDifference{T,T2} <: DiffEqBase.AbstractDiscretization
-    dxs::T
-    time::T2
-    upwind_order::Int
-    centered_order::Int
-    grid_align::GridAlign
-end
-
-# Constructors. If no order is specified, both upwind and centered differences will be 2nd order
-MOLFiniteDifference(dxs, time=nothing; upwind_order = 1, centered_order = 2, grid_align=center_align) =
-    MOLFiniteDifference(dxs, time, upwind_order, centered_order, grid_align)
-
-function SciMLBase.symbolic_discretize(pdesys::PDESystem, discretization::MethodOfLines.MOLFiniteDifference)
+function SciMLBase.symbolic_discretize(pdesys::PDESystem, discretization::MethodOfLines.MOLFiniteDifference{G}) where G
     pdeeqs = pdesys.eqs isa Vector ? pdesys.eqs : [pdesys.eqs]
     bcs = pdesys.bcs
     domain = pdesys.domain
-
+    
     grid_align = discretization.grid_align
     t = discretization.time
     # Get tspan
     tspan = nothing
+    # Check that inputs make sense
     if t !== nothing
         tdomain = pdesys.domain[findfirst(d->isequal(t.val, d.variables), pdesys.domain)]
         @assert tdomain.domain isa DomainSets.Interval
@@ -43,11 +30,12 @@ function SciMLBase.symbolic_discretize(pdesys::PDESystem, discretization::Method
         depvars_lhs = get_depvars(pde.lhs, depvar_ops)
         depvars_rhs = get_depvars(pde.rhs, depvar_ops)
         depvars = collect(depvars_lhs ∪ depvars_rhs)
+
         # Read the independent variables,
         # ignore if the only argument is [t]
         allindvars = Set(filter(xs->!isequal(xs, [t]), map(arguments, depvars)))
-        allnottime = Set(filter(!isempty, map(u->filter(x-> t === nothing || !isequal(x, t.val), arguments(u)), depvars)))
-        if isempty(allnottime)
+        allx̄ = Set(filter(!isempty, map(u->filter(x-> t === nothing || !isequal(x, t.val), arguments(u)), depvars)))
+        if isempty(allx̄)
             push!(alleqs, pde)
             push!(alldepvarsdisc, depvars)
             for bc in bcs
@@ -57,59 +45,38 @@ function SciMLBase.symbolic_discretize(pdesys::PDESystem, discretization::Method
             end
         else
             # make sure there is only one set of independent variables per equation
-            @assert length(allnottime) == 1
-            nottime = first(allnottime)
+            @assert length(allx̄) == 1
+            x̄ = first(allx̄)
             @assert length(allindvars) == 1
             indvars = first(allindvars)
-            nspace = length(nottime)
+            #TODO: Factor this out of the loop, along with BCs
+            # Get the grid
+            s = DiscreteSpace(domain, depvars, indvars, x̄, discretization)
+            
+            # Get the finite difference weights
+            derivweights = DifferentialDiscretizer(pde, pdesys.bcs, s, discretization)
 
-            s = DiscreteSpace(pdesys.domain, depvars, indvars, nottime, grid_align, discretization)
+            # Get the boundary conditions
+            BoundaryHandler!!(u0, bceqs, pdesys.bcs, s, depvar_ops, tspan, derivweights)
 
-            #---- Count Boundary Equations --------------------
-            # Count the number of boundary equations that lie at the spatial boundary on
-            # both the left and right side. This will be used to determine number of
-            # interior equations s.t. we have a balanced system of equations.
+            # Find the indexes on the interior
 
-            # get the depvar boundary terms for given depvar and indvar index.
-            get_depvarbcs(depvar, i) = substitute.((depvar,),get_edgevals(s, i))
-
-            # return the counts of the boundary-conditions that reference the "left" and
-            # "right" edges of the given independent variable. Note that we return the
-            # max of the count for each depvar in the system of equations.
-            @inline function get_bc_counts(i)
-                left = 0
-                right = 0
-                for depvar in s.vars
-                    depvaredges = get_depvarbcs(depvar, i)
-                    counts = [map(x->occursin(x, bc.lhs), depvaredges) for bc in pdesys.bcs]
-                    left = max(left, sum([c[1] for c in counts]))
-                    right = max(right, sum([c[2] for c in counts]))
-                end
-                return [left, right]
-            end
-           
 
-            #TODO: Update this when allowing higher order derivatives to correctly deal with boundary upwinding
-            interior = s.Igrid[[let bcs = get_bc_counts(i)
-                                    (1 + first(bcs)):length(g)-last(bcs)
-                                    end
-                                    for (i,g) in enumerate(s.grid)]...]
-
-            ### PDE EQUATIONS ###
-            # Create a stencil in the required dimension centered around 0
-            # e.g. (-1,0,1) for 2nd order, (-2,-1,0,1,2) for 4th order, etc
-            if discretization.centered_order % 2 != 0
-                throw(ArgumentError("Discretization centered_order must be even, given $(discretization.centered_order)"))
-            end
-            # For all cartesian indices in the interior, generate finite difference rules
-            eqs = vec(map(interior) do II
-                rules = generate_finite_difference_rules(II, s, pde, discretization)
+            Iinterior =  s.Igrid[[let bcs = get_bc_counts(i, s, pdesys.bcs)
+                                      (1 + first(bcs)):length(s.grid[x])-last(bcs)
+                                      end
+                                      for (i,x) in enumerate(s.x̄)]...]
+    
+        
+        
+            # Discretize the equation on the interior
+            pdeeqs = vec(map(Iinterior) do II
+                rules = vcat(generate_finite_difference_rules(II, s, pde, derivweights), valmaps(s, II))
                 substitute(pde.lhs,rules) ~ substitute(pde.rhs,rules)
             end)
 
-            generate_u0_and_bceqs!!(u0, bceqs, pdesys.bcs, s, depvar_ops, tspan, discretization)
-            push!(alleqs,eqs)
-            push!(alldepvarsdisc, reduce(vcat, s.discvars))
+            push!(alleqs,pdeeqs)
+            push!(alldepvarsdisc, reduce(vcat, values(s.discvars)))
         end
     end
 
@@ -131,6 +98,7 @@ function SciMLBase.symbolic_discretize(pdesys::PDESystem, discretization::Method
         return sys, nothing
     else
         # * In the end we have reduced the problem to a system of equations in terms of Dt that can be solved by the `solve` method.
+        #println(alleqs, bceqs)
         sys = ODESystem(vcat(alleqs, unique(bceqs)), t, vec(reduce(vcat, vec(alldepvarsdisc))), ps, defaults=Dict(defaults), name=pdesys.name)
         return sys, tspan
     end

src/MOL_utils.jl

@@ -1,5 +1,7 @@
-# Counts the Differential operators for given variable x. This is used to determine
-# the order of a PDE.
+"""
+Counts the Differential operators for given variable x. This is used to determine
+the order of a PDE.
+"""
 function count_differentials(term, x::Symbolics.Symbolic)
     S = Symbolics
     SU = SymbolicUtils
@@ -15,7 +17,9 @@ function count_differentials(term, x::Symbolics.Symbolic)
     end
 end
 
-# return list of differential orders in the equation
+"""
+return list of differential orders in the equation
+"""
 function differential_order(eq, x::Symbolics.Symbolic)
     S = Symbolics
     SU = SymbolicUtils
@@ -33,7 +37,9 @@ function differential_order(eq, x::Symbolics.Symbolic)
     return filter(!iszero, orders)
 end
 
-# find all the dependent variables given by depvar_ops in an expression
+"""
+find all the dependent variables given by depvar_ops in an expression
+"""
 function get_depvars(eq,depvar_ops)
     S = Symbolics
     SU = SymbolicUtils
@@ -51,4 +57,36 @@ function get_depvars(eq,depvar_ops)
         end
     end
     return depvars
-end
+end
+
+"""
+A function that creates a tuple of CartesianIndices of unit length and `N` dimensions, one pointing along each dimension.
+"""
+function unitindices(N::Int) #create unit CartesianIndex for each dimension
+    null = zeros(Int, N)
+    return map(1:N) do i
+        unit_i = copy(null)
+        unit_i[i] = 1
+        CartesianIndex(Tuple(unit_i))
+    end
+end
+
+@inline function unitindex(N, j)
+    unitindices(N)[j]
+end
+
+function split_additive_terms(eq)
+    # Calling the methods from symbolicutils matches the expressions
+    rhs_arg = istree(eq.rhs) && (SymbolicUtils.operation(eq.rhs) == +) ? SymbolicUtils.arguments(eq.rhs) : [eq.rhs]
+    lhs_arg = istree(eq.lhs) && (SymbolicUtils.operation(eq.lhs) == +) ? SymbolicUtils.arguments(eq.lhs) : [eq.lhs]
+
+    return vcat(lhs_arg,rhs_arg)
+end
+@inline clip(II::CartesianIndex{M}, j, N) where M = II[j] > N ? II - unitindices(M)[j] : II
+subsmatch(expr, rule) = isequal(substitute(expr, rule), expr) ? false : true
+
+    
+
+half_range(x) = -div(x,2):div(x,2)
+
+

src/MethodOfLines.jl

@@ -2,18 +2,28 @@ module MethodOfLines
     using LinearAlgebra
     using SciMLBase
     using DiffEqBase
+    using DiffEqOperators
     using ModelingToolkit
+    using ModelingToolkit: operation, istree, arguments, variable
     using SymbolicUtils, Symbolics
+    using SymbolicUtils: operation, arguments
     import DomainSets
 
-    include("MOL_utils.jl")
     include("discretization/fornberg.jl")
+    
+    include("grid_types.jl")
+    include("MOLFiniteDifference.jl")
+    include("bcs/boundary_types.jl")
+    
     include("discretization/discretize_vars.jl")
-    include("equation_types.jl")
+    include("MOL_utils.jl")
+
+    include("discretization/differential_discretizer.jl")
     include("discretization/generate_finite_difference_rules.jl")
+
     include("bcs/generate_bc_eqs.jl")
 
-    include("discretization/MOL_discretization.jl")
+    include("MOL_discretization.jl")
 
     export MOLFiniteDifference, discretize, symbolic_discretize, grid_align, edge_align, center_align
 end

src/bcs/boundary_types.jl

@@ -0,0 +1,23 @@
+### INITIAL AND BOUNDARY CONDITIONS ###
+
+abstract type AbstractBoundary end
+
+abstract type AbstractTruncatingBoundary <: AbstractBoundary end
+
+abstract type AbstractExtendingBoundary <: AbstractBoundary end
+
+struct LowerBoundary <: AbstractTruncatingBoundary
+end
+
+struct UpperBoundary<: AbstractTruncatingBoundary
+end
+
+struct CompleteBoundary <: AbstractTruncatingBoundary
+end
+
+struct PeriodicBoundary <: AbstractBoundary
+end
+
+struct BoundaryHandler{hasperiodic}
+    boundaries::Dict{Num, AbstractBoundary}
+end