Higher Order Passing

Author: xtalax

src/MOL_discretization.jl

@@ -60,8 +60,15 @@ function SciMLBase.symbolic_discretize(pdesys::PDESystem, discretization::Method
             BoundaryHandler!!(u0, bceqs, pdesys.bcs, s, depvar_ops, tspan, derivweights)
 
             # Find the indexes on the interior
-            Iinterior = interior(s.Igrid, nparams(s))
-
+            
+            
+            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))

src/bcs/generate_bc_eqs.jl

@@ -22,6 +22,31 @@ function generate_boundary_matching_rules(s, orders)
     return (lower, upper)
 end
 
+
+#---- 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.
+
+# get the depvar boundary terms for given depvar and indvar index.
+get_depvarbcs(u, s, i) = substitute.((u,), 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, s, bcs)
+    left = 0
+    right = 0
+    for u in s.ū
+        depvaredges = get_depvarbcs(u, s, i)
+        counts = [map(u -> occursin(u, bc.lhs), depvaredges) for bc in 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
+
+
 """
 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. Important to remember that BCs can have 

src/discretization/discretize_vars.jl

@@ -10,7 +10,7 @@ end
 Return a map of of variables to the gridpoints at the edge of the domain
 """
 @inline function get_edgevals(s, i)
-    return [s.x̄[i] => first(s.axies[x̄[i]]), s.x̄[i] => last(s.axies[x̄[i]])]
+    return [s.x̄[i] => first(s.axies[s.x̄[i]]), s.x̄[i] => last(s.axies[s.x̄[i]])]
 end
 
 @inline function subvar(depvar, edge_vals)
@@ -75,9 +75,9 @@ varmaps(s::DiscreteSpace, II) = [u => s.discvars[u][II] for u in s.ū]
 valmaps(s::DiscreteSpace, II) = vcat(varmaps(s,II), gridvals(s,II))
 
 
-
 interior(Igrid, N) = Igrid[[2:(size(Igrid, i)-1) for i in 1:N]...]
 
+
 map_symbolic_to_discrete(II::CartesianIndex, s::DiscreteSpace{N,M}) where {N,M} = vcat([s.ū[k] => s.discvars[k][II] for k = 1:M], [s.x̄[j] => s.grid[j][II[j]] for j = 1:N])
 
 # TODO: Allow other grids
@@ -170,4 +170,4 @@ edges(s::DiscreteSpace) = s.Iedges
     Ivalid = setdiff(s.Igrid, vcat(vec(Iinterior), vec(Iedges)))
     #TODO: change this when vars have different domains
     return mapreduce(u -> vec(s.discvars[u][Ivalid]) .~ 0., vcat, s.ū)::Vector
-end
+end

test/pde_systems/MOL_1D_HigherOrder.jl

@@ -107,14 +107,24 @@ end
     @named pdesys = PDESystem(eq,bcs,domains,[x,t],[u(x,t)])
     prob = discretize(pdesys,discretization)
 
-    sol = solve(prob,Tsit5(),saveat=0.1,dt=dt)
+    sol = solve(prob,Rosenbrock23(),saveat=0.1,dt=dt)
 
     @test sol.retcode == :Success
 
-    xs = domains[1].domain.lower+dx+dx+dx:dx:domains[1].domain.upper-dx-dx
+    xs = domains[1].domain.lower+dx+dx+dx:dx:domains[1].domain.upper
     ts = sol.t
 
     u_predict = sol.u
     u_real = [[u_analytic(x, t) for x in xs] for t in ts]
+    
+
+       # anim = @animate for (i,T) in enumerate(ts) 
+       #        plot(xs, u_real[i], seriestype = :scatter,label="Analytic solution")
+       #        plot!(xs, sol.u[i], label="Numeric solution")
+       #        plot!(xs, log.(abs.(u_real[i]-sol.u[i])), label="Log Error at t = $(ts[i])")
+       # end
+       # gif(anim, "plots/MOL_Higher_order_1D_KdV_single_soliton.gif", fps = 5)
+
+
     @test all(isapprox.(u_predict, u_real, rtol = 0.03))
 end

test/runtests.jl

@@ -9,14 +9,14 @@ const is_TRAVIS = haskey(ENV,"TRAVIS")
 @time begin
 
     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 "Discretization of space and grid types" begin include("components/DiscreteSpace.jl") end
-        #@time @safetestset "Finite Difference Schemes" begin include("components/finite_diff_schemes.jl") end
+       #@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 "Discretization of space and grid types" begin include("components/DiscreteSpace.jl") end
+        @time @safetestset "Finite Difference Schemes" begin include("components/finite_diff_schemes.jl") end
     end
 
     if GROUP == "All" || GROUP == "Integration Tests"
-        #@time @safetestset "MOLFiniteDifference Interface" begin include("pde_systems/MOLtest.jl") end
+        @time @safetestset "MOLFiniteDifference Interface" begin include("pde_systems/MOLtest.jl") end
         #@time @safetestset "MOLFiniteDifference Interface: Linear Convection" begin include("pde_systems/MOL_1D_Linear_Convection.jl") end
         @time @safetestset "MOLFiniteDifference Interface: 1D Linear Diffusion" begin include("pde_systems/MOL_1D_Linear_Diffusion.jl") end
         @time @safetestset "MOLFiniteDifference Interface: 1D Non-Linear Diffusion" begin include("pde_systems/MOL_1D_NonLinear_Diffusion.jl") end