fixed and simplified half offset

Author: xtalax

src/MOL_utils.jl

@@ -74,12 +74,23 @@ function unitindices(N::Int) #create unit CartesianIndex for each dimension
 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]
+    lhs_arg = istree(eq.lhs) && (SymbolicUtils.operation(eq.lhs) == +) ? SymbolicUtils.arguments(eq.lhs) : [eq.lhs]
 
     return vcat(lhs_arg,rhs_arg)
 end
 
+function clip(I::CartesianIndex, s::DiscreteSpace{N}, j::Int) where N
+    # Clip the index by 1 at each end of the dimension
+    I1 = unitindices(N)[j]
+    return I[j] > length(s, j) ? I-I1 : I+I1
+end
+    
+subsmatch(expr, rule) = isequal(substitute(expr, rule), expr) ? false : true
+
+    
+
 half_range(x) = -div(x,2):div(x,2)
 
 

src/bcs/generate_bc_eqs.jl

@@ -1,16 +1,12 @@
-function _boundary_rules(s, orders, val)
-    args = copy(s.params)
+function _boundary_rules(s, orders, x, val)
+    #TODO: Change this around to detect boundary variables anywhere in the eq
+   args = copy(s.params)
 
-    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]
+    args = substitute.(args, (x=>val,))
+
+    rules = [operation(u)(args...) => (operation(u)(args...), x) for u in ū]
 
-    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]))
+    return vcat(rules, vec([(Differential(x)^d)(operation(u)(args...)) => (operation(u)(args...), x) for d in orders[x], u in ū]))
 end
 
 function generate_boundary_matching_rules(s, orders)
@@ -19,9 +15,9 @@ function generate_boundary_matching_rules(s, orders)
     upperboundary(x) = last(s.axies[x])
 
     # Rules to match boundary conditions on the lower boundaries
-    lower = reduce(vcat, [_boundary_rules(s, orders, lowerboundary(x)) for x in s.vars])
+    lower = reduce(vcat, [_boundary_rules(s, orders, x, lowerboundary(x)) for x in s.x̄])
 
-    upper = reduce(vcat, [_boundary_rules(s, orders, upperboundary(x)) for x in s.vars])
+    upper = reduce(vcat, [_boundary_rules(s, orders, x, upperboundary(x)) for x in s.x̄])
 
     return (lower, upper)
 end
@@ -35,9 +31,9 @@ function BoundaryHandler!!(u0, bceqs, bcs, s::DiscreteSpace, depvar_ops, tspan,
     t=s.time
 
     if t === nothing
-        initmaps = s.vars
+        initmaps = ū
     else
-        initmaps = substitute.(s.vars,[t=>tspan[1]])
+        initmaps = substitute.(ū,[t=>tspan[1]])
     end
 
     # Create some rules to match which bundary/variable a bc concerns
@@ -56,13 +52,13 @@ function BoundaryHandler!!(u0, bceqs, bcs, s::DiscreteSpace, depvar_ops, tspan,
         # * 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 any(u -> isequal(operation(u), operation(bcdepvar)), ū)
             if t !== nothing && operation(bc.lhs) isa Sym && !any(x -> isequal(x, t.val), arguments(bc.lhs))
                 # initial condition
                 # * 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))))
+                    push!(u0,vec(s.discvars[ū[initindex]] .=> substitute.((bc.rhs,),gridvals(s))))
                 end
             else
                 # Split out additive terms
@@ -72,23 +68,26 @@ function BoundaryHandler!!(u0, bceqs, bcs, s::DiscreteSpace, depvar_ops, tspan,
                 boundary = nothing
                 # Check whether the bc is on the lower boundary, or periodic
                 for term in terms, r in lower_boundary_rules
-                    if r(term) !== nothing
-                        u_, x_ = (term, r(term))
+                    #Check if the rule changes the expression
+                    if subsmatch(term, r)
+                        # Get the matched variables from the rule
+                        u_, x_ = r.second
+                        # Mark the boundary                        
                         boundary = LowerBoundary()
-                        for term_ in setdiff(terms, term)
-                            for r in upper_boundary_rules
-                                if r(term_) !== nothing
-                                    # boundary = PeriodicBoundary()
-                                    #TODO: Add handling for perioodic boundary conditions here
-                                end
-                            end
-                        end
+                        # for term_ in setdiff(terms, term)
+                        #     for r_ in upper_boundary_rules
+                        #         if subsmatch(term_, r_) !== nothing
+                        #             boundary = PeriodicBoundary()
+                        #             #TODO: Add handling for perioodic boundary conditions here
+                        #         end
+                        #     end
+                        # end
                         break
                     end
                 end
                 for term in terms, r in upper_boundary_rules
-                    if r(term) !== nothing
-                        u_, x_ = (term, r(term))
+                    if subsmatch(term, r)
+                        u_, x_ = r.second 
                         boundary = UpperBoundary()
                         break
                     end
@@ -115,10 +114,10 @@ function generate_bc_rules(II, derivweights, s::DiscreteSpace{N,M,G}, bc, u_, x_
     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
+    for u in ū
         if isequal(operation(u), operation(u_))
             # What to replace derivatives at the boundary with
-            depvarderivbcmaps = [(Differential(x)^d)(u_) => central_difference(derivweights.map[Differential(x_)^d], II, s, (s.x2i[x_], x_), u, ufunc) for d in derivweights.orders[x_]]
+            depvarderivbcmaps = [(Differential(x_)^d)(u_) => central_difference(derivweights.map[Differential(x_)^d], II, s, (s.x2i[x_], x_), u, ufunc) for d in derivweights.orders[x_]]
             # ? Does this need to be done for all variables at the boundary?
             depvarbcmaps = [u_ => s.discvars[u][II]]
             break
@@ -136,8 +135,6 @@ end
 
 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 = []
@@ -146,11 +143,11 @@ function generate_bc_rules(II, derivweights, s::DiscreteSpace{N,M,G}, bc, u_, x_
     # 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
-    for u in s.vars
+    for u in ū
         if isequal(operation(u), operation(u_))
-            depvarderivbcmaps = [(Differential(x)^d)(u_) => half_offset_centered_difference(derivweights.halfoffsetmap[Differential(x_)^d], II, s, offset(boundary), (s.x2i[x_],x_), u, ufunc) for d in derivweights.orders[x_]]
+            depvarderivbcmaps = [(Differential(x_)^d)(u_) => half_offset_centered_difference(derivweights.halfoffsetmap[Differential(x_)^d], II, s, (s.x2i[x_],x_), u, ufunc) for d in derivweights.orders[x_]]
 
-            depvarbcmaps = [u_ => half_offset_centered_difference(derivweights.interpmap[x_], II, s, offset(boundary), (s.x2i[x_],x_), u, ufunc)]
+            depvarbcmaps = [u_ => half_offset_centered_difference(derivweights.interpmap[x_], II, s, (s.x2i[x_],x_), u, ufunc)]
             break
         end
     end

src/bcs/scratch.jl

@@ -0,0 +1,62 @@
+
+using DiffEqOperators
+
+"""
+Inner function for `get_half_offset_weights_and_stencil` (see below).
+"""
+function _get_weights_and_stencil(D, I)
+    # k => i of itap - 
+    # offset is important due to boundary proximity
+    # The low boundary coeffs has a heirarchy of coefficients following: number of indices from boundary -> which half offset point does it correspond to -> weights
+    if I <= (D.boundary_point_count)
+        println("low boundary")
+        weights = D.low_boundary_coefs[I]    
+        offset = 1 - I
+        Itap = [I + (i+offset) for i in 0:(D.boundary_stencil_length-1)]
+    elseif I > (20 - D.boundary_point_count)
+        println("high boundary")
+
+        weights = D.high_boundary_coefs[20-I+1]
+        offset = 20 - I
+        Itap = [I + (i+offset) for i in (-D.boundary_stencil_length+1):1:0]
+    else
+        println("interior")
+
+        weights = D.stencil_coefs
+        Itap = [I + i for i in (1-div(D.stencil_length,2)):(div(D.stencil_length,2))]
+    end    
+
+    return (weights, Itap)
+end
+
+
+function clip(I, n)
+    return I>(n/2) ? I-1 : I+1
+    
+end
+
+for I in [2,19, 3, 18, 10]
+    
+    D =  CompleteHalfCenteredDifference(1, 6, 1.0)
+    
+     
+    
+    outerweights, outerstencil = _get_weights_and_stencil(D, clip(I,20))
+    
+    # Get the correct weights and stencils for this II
+    
+    @show I
+    @show outerweights
+    @show outerstencil
+    
+    
+    for II in outerstencil
+        println()
+        println("II: ", II)
+        println()
+        weights, stencil = _get_weights_and_stencil(D, II)
+        println(weights)
+        println(stencil)
+    end
+
+end

src/discretization/MOL_discretization.jl

@@ -34,8 +34,8 @@ function SciMLBase.symbolic_discretize(pdesys::PDESystem, discretization::Method
         # 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
@@ -45,13 +45,13 @@ 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)
 
             # Get the grid
-            s = DiscreteSpace(domain, depvars, indvars, nottime, discretization)
+            s = DiscreteSpace(domain, depvars, indvars, x̄, discretization)
 
             # Get the finite difference weights
             derivweights = DifferentialDiscretizer(pde, s, discretization)
@@ -64,12 +64,12 @@ function SciMLBase.symbolic_discretize(pdesys::PDESystem, discretization::Method
 
             # Discretize the equation on the interior
             pdeeqs = vec(map(interior) do II
-                rules = vcat(generate_finite_difference_rules(II, s, pde, derivweights), valrules(s, II))
+                rules = vcat(generate_finite_difference_rules(II, s, pde, derivweights), valmaps(s, II))
                 substitute(pde.lhs,rules) ~ substitute(pde.rhs,rules)
             end)
 
             push!(alleqs,pdeeqs)
-            push!(alldepvarsdisc, reduce(vcat, s.discvars))
+            push!(alldepvarsdisc, reduce(vcat, values(s.discvars)))
         end
     end
 

src/discretization/differential_discretizer.jl

@@ -12,14 +12,14 @@ function DifferentialDiscretizer(pde, s, discretization)
     # TODO: Include bcs in this calculation
     d_orders(x) = reverse(sort(collect(union(differential_order(pde.rhs, x), differential_order(pde.lhs, x)))))
 
-    # central_deriv_rules = [(Differential(s)^2)(u) => central_deriv(2,II,j,k) for (j,s) in enumerate(s.nottime), (k,u) in enumerate(s.vars)]
+    # central_deriv_rules = [(Differential(s)^2)(u) => central_deriv(2,II,j,k) for (j,s) in enumerate(s.x̄), (k,u) in enumerate(ū)]
     differentialmap = Array{Pair{Num,DiffEqOperators.DerivativeOperator},1}()
     nonlinlap = []
     interp = []
     orders = []
     # Hardcoded to centered difference, generate weights for each differential
     # TODO: Add handling for upwinding
-    for x in s.nottime
+    for x in s.x̄
         push!(orders, x => d_orders(x))
         # TODO: Only generate weights for derivatives that are actually used and avoid redundant calculations
         rs = [(Differential(x)^d) => CompleteCenteredDifference(d, approx_order, s.dxs[x] ) for d in last(orders).second]
@@ -60,51 +60,39 @@ function central_difference(D, II, s, jx, u, ufunc)
     return dot(weights, ufunc(u, Itap, x))
 end
 
+
 """
-Inner function for `get_half_offset_weights_and_stencil` (see below).
+Get the weights and stencil for the inner half offset centered difference for the nonlinear laplacian for a given index and differentiating variable.
+
+Does not discretize so that the weights can be used in a replacement rule
+TODO: consider refactoring this to harmonize with centered difference
 """
-function _get_weights_and_stencil(D, II, I1, s, k, j, x)
-    # k => i of itap - 
+function get_half_offset_weights_and_stencil(D, II, s, jx)
+    j, x = jx
+    # unit index in direction of the derivative
+    I1 = unitindices(nparams(s))[j] 
     # offset is important due to boundary proximity
-    # The low boundary coeffs has a heirarchy of coefficients following: number of indices from boundary -> which half offset point does it correspond to -> weights
-    if II[j] <= (D.boundary_point_count-1)
-        weights = D.low_boundary_coefs[II[j]][k]    
+
+    if II[j] <= D.boundary_point_count
+        weights = D.low_boundary_coefs[II[j]]    
         offset = 1 - II[j]
         Itap = [II + (i+offset)*I1 for i in 0:(D.boundary_stencil_length-1)]
-    elseif II[j] > (length(s, x) - D.boundary_point_count - 1)
-        weights = D.high_boundary_coefs[length(s,x)-II[j]+1][k]
+    elseif II[j] > (length(s, x) - D.boundary_point_count)
+        weights = D.high_boundary_coefs[length(s, x)-II[j]+1]
         offset = length(s, x) - II[j]
         Itap = [II + (i+offset)*I1 for i in (-D.boundary_stencil_length+1):1:0]
     else
         weights = D.stencil_coefs
-        Itap = [II + (i+1)*I1 for i in half_range(D.stencil_length)]
+        Itap = [II + i*I1 for i in (1-div(D.stencil_length,2)):(div(D.stencil_length,2))]
     end    
 
     return (weights, Itap)
 end
 
-"""
-Get the weights and stencil for the inner half offset centered difference for the nonlinear laplacian for a given index and differentiating variable.
-
-Does not discretize so that the weights can be used in a replacement rule
-TODO: consider refactoring this to harmonize with centered difference
-"""
-function get_half_offset_weights_and_stencil(D, II, s, offset, jx)
-    j, x = jx
-    I1 = unitindices(nparams(s))[j]
-    # Shift the current index to the correct offset
-    II_prime = II + Int(offset-0.5)*I1
-    @assert all(i-> 0<i<=length(s,x), II_prime) "Index out of bounds"
-    return _get_weights_and_stencil(D, II_prime, I1, s, offset, j, x)
-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, jx, u, ufunc)
-    j, x = jx
-    I1 = unitindices(nparams(s))[j]
-    # Shift the current index to the correct offset
-    II_prime = II + Int(offset-0.5)*I1
-    # Get the weights and stencil
-    (weights, Itap) = _get_weights_and_stencil(D, II_prime, I1, s, offset, j, x)
+function half_offset_centered_difference(D, II, s, jx, u, ufunc)
+
+    (weights, Itap) = get_half_offset_weights_and_stencil(D, II, s, jx)
     return dot(weights, ufunc(u, Itap, x))
 end