some assumptions

Author: xtalax

src/bcs/generate_bc_eqs.jl

@@ -13,7 +13,7 @@ struct UpperBoundary<: AbstractTruncatingBoundary
 end
 
 struct CompleteBoundary <: AbstractTruncatingBoundary
-end#
+end
 
 struct PeriodicBoundary <: AbstractBoundary
 end
@@ -37,6 +37,7 @@ function BoundaryHandler!!(u0, bceqs, bcs, s, depvar_ops, tspan, derivweights)
     end
 
     # 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])
 
@@ -61,7 +62,6 @@ function BoundaryHandler!!(u0, bceqs, bcs, s, depvar_ops, tspan, derivweights)
 
                 u_, x_ = (nothing, nothing)
                 terms = vcat(lhs_arg,rhs_arg)
-                # * Assume that the term of the condition is applied additively and has no multiplier/divisor/power etc.
                 boundary = nothing
                 # Check whether the bc is on the lower boundary, or periodic
                 for term in terms, r in lower_boundary_rules
@@ -71,7 +71,7 @@ function BoundaryHandler!!(u0, bceqs, bcs, s, depvar_ops, tspan, derivweights)
                         for term_ in setdiff(terms, term)
                             for r in upper_boundary_rules
                                 if r(term_) !== nothing
-                                    # boundary = :periodic
+                                    # boundary = PeriodicBoundary()
                                     #TODO: Add handling for perioodic boundary conditions here
                                 end
                             end
@@ -104,41 +104,49 @@ function generate_bc_rules(II, dim, s, bc, u_, x_, boundary::AbstractTruncatingB
     # 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]
-
+    # * 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_))
             # 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
         end
     end
 
     fd_rules = generate_finite_difference_rules(II, s, bc, derivweights)
-    varrules = axiesvals(s, II)
+    varrules = axiesvals(s, x_, II)
+
+    # valrules should be caught by depvarbcmaps and varrules if the above assumption holds
+    #valr = valrules(s, II)
 
     return vcat(depvarderivbcmaps, depvarbcmaps, fd_rules, varrules)
 end
 
 function generate_bc_rules(II, dim, s, bc, u_, x_, boundary::AbstractTruncatingBoundary, derivweights, G::EdgeAlignedGrid) 
 
-    boundaryoffset(::LowerBoundary) = 1/2
-    boundaryoffset(::UpperBoundary) = -1/2
+    offset(::LowerBoundary) = 1/2
+    offset(::UpperBoundary) = -1/2
     ufunc(v, I, x) = s.discvars[v][I]
 
     # 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
         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_]]
 
             depvarbcmaps = [u_ => half_offset_centered_difference(derivweights.interpmap[x_], II, s, offset(boundary), (s.x2i[x_],x_), u, ufunc)]
+            break
         end
     end
 
     fd_rules = generate_finite_difference_rules(II, s, bc, derivweights)
-    varrules = axiesvals(s, II)
-    valr = valrules(s, II)
+    varrules = axiesvals(s, x_, II)
+    
+    # valrules should be caught by depvarbcmaps and varrules if the above assumption holds
+    #valr = valrules(s, II)
 
     return vcat(depvarderivbcmaps, depvarbcmaps, fd_rules, varrules)
 end

src/discretization/discretize_vars.jl

@@ -42,7 +42,7 @@ params(s::DiscreteSpace{N,M}) where {N,M}= s.params
 grid_idxs(s::DiscreteSpace) = CartesianIndices(((axes(g)[1] for g in s.grid)...,))
 edge_idxs(s::DiscreteSpace{N}) where {N} = reduce(vcat, [[vcat([Colon() for j = 1:i-1], 1, [Colon() for j = i+1:N]), vcat([Colon() for j = 1:i-1], length(s.axies[i]), [Colon() for j = i+1:N])] for i = 1:N])
 
-axiesvals(s::DiscreteSpace{N}, I) where {N} = [x => s.axies[x][I[j]] for (j,x) in enumerate(s.nottime)]
+axiesvals(s::DiscreteSpace{N}, x_, I) where {N} = [x => (x == x_ ? s.axies[x][I[j]] : s.grid[x][I[j]]) for (j,x) in enumerate(s.nottime)]
 gridvals(s::DiscreteSpace{N}, I) where {N} = [x => s.grid[x][I[j]] for (j,x) in enumerate(s.nottime)]
 
 ## Boundary methods ##

src/discretization/generate_finite_difference_rules.jl

@@ -148,7 +148,7 @@ There are of course more specific schemes that are used to improve stability/spe
 ##Planned special cases include:
     - Up/Downwind schemes to be used for odd ordered derivatives multiplied by a coefficient, downwinding when the coefficient is positive, and upwinding when the coefficient is negative.
 
-Please submit an issue if you know of any special cases that are not implemented, with links to papers and/or code that demonstrates the special case.
+Please submit an issue if you know of any special cases which impact stability or accuracy that are not implemented, with links to papers and/or code that demonstrates the special case.
 """
 function generate_finite_difference_rules(II, s, pde, derivweights)
     # 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)]

test/runtests.jl

@@ -10,7 +10,7 @@ const is_TRAVIS = haskey(ENV,"TRAVIS")
 
     if GROUP == "All" || GROUP == "Component Tests"
        # @time @safetestset "Test for regression against original code" begin include("regression_test.jl") end
-        #@time @safetestset "MOLFiniteDifference Utils" begin include("utils_test.jl") end
+        @time @safetestset "MOLFiniteDifference Utils" begin include("utils_test.jl") end
     end
 
     if GROUP == "All" || GROUP == "Integration Tests"