test fixes

Author: xtalax

src/MOLFiniteDifference.jl

@@ -0,0 +1,15 @@
+struct MOLFiniteDifference{G} <: DiffEqBase.AbstractDiscretization
+    dxs
+    time
+    approx_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, 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, grid_align)
+end

src/MOL_utils.jl

@@ -73,6 +73,13 @@ function unitindices(N::Int) #create unit CartesianIndex for each dimension
     Tuple(out)
 end
 
+function split_additive_terms(eq)
+    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
+
 half_range(x) = -div(x,2):div(x,2)
 
 

src/MethodOfLines.jl

@@ -11,10 +11,17 @@ module MethodOfLines
 
     include("MOL_utils.jl")
     include("discretization/fornberg.jl")
-    include("bcs/generate_bc_eqs.jl")
+    
+    include("grid_types.jl")
+    include("MOLFiniteDifference.jl")
+    include("bcs/boundary_types.jl")
+
     include("discretization/discretize_vars.jl")
+
     include("discretization/differential_discretizer.jl")
-    include("discretization/generate_rules.jl")
+    include("discretization/generate_finite_difference_rules.jl")
+
+    include("bcs/generate_bc_eqs.jl")
 
     include("discretization/MOL_discretization.jl")
 

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

src/bcs/generate_bc_eqs.jl

@@ -1,32 +1,36 @@
-### INITIAL AND BOUNDARY CONDITIONS ###
+function _boundary_rules(s, orders, val)
+    args = copy(s.params)
 
-abstract type AbstractBoundary end
-
-abstract type AbstractTruncatingBoundary <: AbstractBoundary end
-
-abstract type AbstractExtendingBoundary <: AbstractBoundary end
+    if isequal(val, floor(val))
+        args = [substitute.(args, (x=>val,)), substitute.(args, (x=>Int(val),))]
+    else
+        args = [substitute.(args, (x=>val,))]
+    end
+    substitute.(args, (x=>lowerboundary(x),))
+    
+    rules = [@rule operation(u)(arg...) => (u, x) for u in s.vars, arg in args]
 
-struct LowerBoundary <: AbstractTruncatingBoundary
+    return vcat(rules, vec([@rule (Differential(x)^d)(operation(u)(arg...)) => (u, x) for d in orders[x], u in s.vars, arg in args]))
 end
 
-struct UpperBoundary<: AbstractTruncatingBoundary
-end
+function generate_boundary_matching_rules(s, orders)
+    # TODO: Check for bc equations of multiple variables
+    lowerboundary(x) = first(s.axies[x])
+    upperboundary(x) = last(s.axies[x])
 
-struct CompleteBoundary <: AbstractTruncatingBoundary
-end
+    # Rules to match boundary conditions on the lower boundaries
+    lower = reduce(vcat, [_boundary_rules(s, orders, lowerboundary(x)) for x in s.vars])
 
-struct PeriodicBoundary <: AbstractBoundary
-end
+    upper = reduce(vcat, [_boundary_rules(s, orders, upperboundary(x)) for x in s.vars])
 
-struct BoundaryHandler{hasperiodic}
-    boundaries::Dict{Num, AbstractBoundary}
+    return (lower, upper)
 end
 
 """
 Mutates bceqs and u0 by finding relevant equations and discretizing them.
-TODO: return a handler for use with generate_finite_difference_rules and pull out initial condition
+TODO: return a handler for use with generate_finite_difference_rules and pull out initial condition. Important to remember that BCs can have 
 """
-function BoundaryHandler!!(u0, bceqs, bcs, s, depvar_ops, tspan, derivweights)
+function BoundaryHandler!!(u0, bceqs, bcs, s::DiscreteSpace, depvar_ops, tspan, derivweights::DifferentialDiscretizer) 
 
     t=s.time
 
@@ -38,30 +42,33 @@ function BoundaryHandler!!(u0, bceqs, bcs, s, depvar_ops, tspan, derivweights)
 
     # Create some rules to match which bundary/variable a bc concerns
     # * Assume that the term of the condition is applied additively and has no multiplier/divisor/power etc.
-    # ? Is it nessecary to check whether all other args are present?
-    lower_boundary_rules = vec([@rule operation(u)(~~a, lowerboundary(x), ~~b) => IfElse.ifelse(all(y-> y in vcat(~~a, ~~b), setdiff(x, arguments(u))), x, nothing) for x in setdiff(arguments(u), t), u in s.vars])
+    
+    ## BC matching rules, returns the variable and parameter the bc concerns
+
+    lower_boundary_rules, upper_boundary_rules = generate_boundary_matching_rules(s, derivweights.orders)
 
-    upper_boundary_rules = vec([@rule operation(u)(~~a, upperboundary(x), ~~b) => IfElse.ifelse(all(y-> y in vcat(~~a, ~~b), setdiff(x, arguments(u))), x, nothing) for x in setdiff(arguments(u), t), u in s.vars])
+    # indexes for Iedge depending on boundary type
+    idx(::LowerBoundary) = 1
+    idx(::UpperBoundary) = 2
 
     # Generate initial conditions and bc equations
     for bc in bcs
+        # * Assume in the form `u(...) ~ ...` for now
         bcdepvar = first(get_depvars(bc.lhs, depvar_ops))
+        
         if any(u -> isequal(operation(u), operation(bcdepvar)), s.vars)
             if t !== nothing && operation(bc.lhs) isa Sym && !any(x -> isequal(x, t.val), arguments(bc.lhs))
                 # initial condition
-                # * Assume in the form `u(...) ~ ...` for now
                 # * Assume that the initial condition is not in terms of the initial derivative
                 initindex = findfirst(isequal(bc.lhs), initmaps) 
                 if initindex !== nothing
                     push!(u0,vec(s.discvars[s.vars[initindex]] .=> substitute.((bc.rhs,),gridvals(s))))
                 end
             else
                 # Split out additive terms
-                rhs_arg = istree(pde.rhs) && (SymbolicUtils.operation(pde.rhs) == +) ? SymbolicUtils.arguments(pde.rhs) : [pde.rhs]
-                lhs_arg = istree(pde.lhs) && (SymbolicUtils.operation(pde.lhs) == +) ? SymbolicUtils.arguments(pde.lhs) : [pde.lhs]
+                terms = split_additive_terms(bc)
 
                 u_, x_ = (nothing, nothing)
-                terms = vcat(lhs_arg,rhs_arg)
                 boundary = nothing
                 # Check whether the bc is on the lower boundary, or periodic
                 for term in terms, r in lower_boundary_rules
@@ -75,6 +82,7 @@ function BoundaryHandler!!(u0, bceqs, bcs, s, depvar_ops, tspan, derivweights)
                                     #TODO: Add handling for perioodic boundary conditions here
                                 end
                             end
+                        end
                         break
                     end
                 end
@@ -86,24 +94,26 @@ function BoundaryHandler!!(u0, bceqs, bcs, s, depvar_ops, tspan, derivweights)
                     end
                 end
 
-                @assert boundary !== nothing "Boundary condition ${bc} is not on a boundary of the domain, or is not a valid boundary condition"
+                @assert boundary !== nothing "Boundary condition $bc is not on a boundary of the domain, or is not a valid boundary condition"
 
-                push!(bceqs, vec(map(s.Iedge[x_][boundary]) do II
-                    rules = generate_bc_rules(II, s, bc, u_, boundary, derivweights)
-                    rules = vcat(rules, )
+                push!(bceqs, vec(map(s.Iedge[x_][idx(boundary)]) do II
+                    rules = generate_bc_rules(II, derivweights, s, bc, u_, x_, boundary)
 
                     substitute(bc.lhs, rules) ~ substitute(bc.rhs, rules)
                 end))
             end
-        else 
-            throw(ArgumentError("No active variables in boundary condition $bc lhs, please ensure that bcs are "))
+        end
     end
 end
 
-function generate_bc_rules(II, dim, s, bc, u_, x_, boundary::AbstractTruncatingBoundary, derivweights, G::CenterAlignedGrid) 
+function generate_bc_rules(II, derivweights, s::DiscreteSpace{N,M,G}, bc, u_, x_, ::AbstractTruncatingBoundary) where {N, M, G<:CenterAlignedGrid}
     # depvarbcmaps will dictate what to replace the variable terms with in the bcs
     # replace u(t,0) with u₁, etc
     ufunc(v, I, x) = s.discvars[v][I]
+
+    depvarderivbcmaps = []
+    depvarbcmaps = []
+
     # * Assume that the BC is in terms of an explicit expression, not containing references to variables other than u_ at the boundary
     for u in s.vars
         if isequal(operation(u), operation(u_))
@@ -124,12 +134,15 @@ function generate_bc_rules(II, dim, s, bc, u_, x_, boundary::AbstractTruncatingB
     return vcat(depvarderivbcmaps, depvarbcmaps, fd_rules, varrules)
 end
 
-function generate_bc_rules(II, dim, s, bc, u_, x_, boundary::AbstractTruncatingBoundary, derivweights, G::EdgeAlignedGrid) 
+function generate_bc_rules(II, derivweights, s::DiscreteSpace{N,M,G}, bc, u_, x_, boundary::AbstractTruncatingBoundary) where {N, M, G<:EdgeAlignedGrid}
 
     offset(::LowerBoundary) = 1/2
     offset(::UpperBoundary) = -1/2
     ufunc(v, I, x) = s.discvars[v][I]
 
+    depvarderivbcmaps = []
+    depvarbcmaps = []
+
     # depvarbcmaps will dictate what to replace the variable terms with in the bcs
     # replace u(t,0) with u₁, etc
     # * Assume that the BC is in terms of an explicit expression, not containing references to variables other than u_ at the boundary
@@ -144,7 +157,7 @@ function generate_bc_rules(II, dim, s, bc, u_, x_, boundary::AbstractTruncatingB
 
     fd_rules = generate_finite_difference_rules(II, s, bc, derivweights)
     varrules = axiesvals(s, x_, II)
-    
+
     # valrules should be caught by depvarbcmaps and varrules if the above assumption holds
     #valr = valrules(s, II)
 

src/discretization/MOL_discretization.jl

@@ -1,32 +1,6 @@
 # Method of lines discretization scheme
-abstract type AbstractGrid end
 
-struct CenterAlignedGrid <: AbstractGrid 
-end
-
-struct EdgeAlignedGrid <: AbstractGrid 
-end
-
-const center_align=CenterAlignedGrid()
-const edge_align=EdgeAlignedGrid()
-
-struct MOLFiniteDifference{G} <: DiffEqBase.AbstractDiscretization
-    dxs
-    time
-    approx_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; upwind_order = 1, centered_order = 2, grid_align=CenterAlignedGrid())
-    
-    if centered_order % 2 != 0
-        warn("Discretization centered_order must be even, rounding up to $(centered_order+1)")
-    end
-    return MOLFiniteDifference{typeof(grid_align)}(dxs, time, centered_order, grid_align)
-end
-
-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
@@ -56,6 +30,7 @@ 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)))
@@ -82,7 +57,7 @@ function SciMLBase.symbolic_discretize(pdesys::PDESystem, discretization::Method
             derivweights = DifferentialDiscretizer(pde, s, discretization)
 
             # Get the boundary conditions
-            generate_u0_and_bceqs!!(u0, bceqs, pdesys.bcs, s, depvar_ops, tspan, derivweights)
+            BoundaryHandler!!(u0, bceqs, pdesys.bcs, s, depvar_ops, tspan, derivweights)
 
             # Find the indexes on the interior
             interior = Iinterior(s)

src/discretization/differential_discretizer.jl

@@ -2,7 +2,7 @@
 struct DifferentialDiscretizer{T, D1}
     approx_order::Int
     map::D1
-    halfoffsetmap::Dict{Num, DiffEqOperators.DerivativeOperator}
+    halfoffsetmap::Dict
     interpmap::Dict{Num, DiffEqOperators.DerivativeOperator}
     orders::Dict{Num, Vector{Int}}
 end
@@ -26,7 +26,8 @@ function DifferentialDiscretizer(pde, s, discretization)
 
         nonlinlap = vcat(nonlinlap, [Differential(x)^d => CompleteHalfCenteredDifference(d, approx_order, s.dxs[x]) for d in last(orders).second])
         differentialmap = vcat(differentialmap, rs)
-        push!(interp, x => CompleteHalfCenteredDifference(0, approx_order, s.dxs[x])
+        # A 0th order derivative off the grid is an interpolation
+        push!(interp, x => CompleteHalfCenteredDifference(0, approx_order, s.dxs[x]))
     end
 
     return DifferentialDiscretizer{eltype(orders), typeof(Dict(differentialmap))}(approx_order, Dict(differentialmap), Dict(nonlinlap), Dict(interp), Dict(orders))
@@ -98,12 +99,12 @@ function get_half_offset_weights_and_stencil(D, II, s, offset, jx)
 end
 
 # i is the index of the offset, assuming that there is one precalculated set of weights for each offset required for a first order finite difference
-function half_offset_centered_difference(D, II, s, offset, i, jx, u, ufunc)
+function half_offset_centered_difference(D, II, s, offset, jx, u, ufunc)
     j, x = jx
     I1 = unitindices(nparams(s))[j]
     # Shift the current index to the correct offset
-    II_prime = II + offset*I1
+    II_prime = II + Int(offset-0.5)*I1
     # Get the weights and stencil
-    (weights, Itap) = _get_weights_and_stencil(D, II_prime, I1, s, i, j, x)
+    (weights, Itap) = _get_weights_and_stencil(D, II_prime, I1, s, offset, j, x)
     return dot(weights, ufunc(u, Itap, x))
 end