stash

Author: xtalax

src/discretization/MOL_discretization.jl

@@ -70,14 +70,16 @@ function SciMLBase.symbolic_discretize(pdesys::PDESystem, discretization::Method
 
             s = DiscreteSpace(pdesys.domain, depvars, indvars, nottime, grid_align, discretization)
 
+            derivweights = DifferentialDiscretizer(pde, s, discretization)
+
             interior = s.Igrid[[2:(length(grid[x])-1) for x in s.nottime]
 
             ### 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
             # For all cartesian indices in the interior, generate finite difference rules
             pdeeqs = vec(map(interior) do II
-                rules = generate_finite_difference_rules(II, s, pde, discretization)
+                rules = generate_finite_difference_rules(II, s, pde, derivweights)
                 substitute(pde.lhs,rules) ~ substitute(pde.rhs,rules)
             end)
 

src/discretization/discretize_vars.jl

@@ -27,6 +27,7 @@ struct DiscreteSpace{N,M}
     Iaxies
     Igrid
     Iedge
+    x2i
 end
 
 nparams(::DiscreteSpace{N,M}) where {N,M} = N
@@ -52,7 +53,8 @@ map_symbolic_to_discrete(II::CartesianIndex, s::DiscreteSpace{N,M}) where {N,M}
 
 # TODO: Allow other grids
 
-function DiscreteSpace(domain, depvars, indvars, nottime, grid_align, discretization)
+function DiscreteSpace(domain, depvars, indvars, nottime, discretization)
+    grid_align = discretization.grid_align
     t = discretization.time
     nspace = length(nottime)
     # Discretize space
@@ -78,7 +80,7 @@ function DiscreteSpace(domain, depvars, indvars, nottime, grid_align, discretiza
         grid =  map(nottime) do x
             xdomain = domain[findfirst(d -> isequal(x, d.variables), domain)]
             dx = discretization.dxs[findfirst(dxs -> isequal(x, dxs[1].val), discretization.dxs)][2]
-            dx isa Number ? x => ((DomainSets.infimum(xdomain.domain)-dx/2):dx:(DomainSets.supremum(xdomain.domain)+dx/2) : x => dx
+            dx isa Number ? (x => ((DomainSets.infimum(xdomain.domain)-dx/2):dx:(DomainSets.supremum(xdomain.domain)+dx/2))) : x => dx
         end 
         # TODO: allow depvar-specific center/edge choice?
     end
@@ -94,18 +96,18 @@ function DiscreteSpace(domain, depvars, indvars, nottime, grid_align, discretiza
             u => [u for II in s.Igrid]
         else
             sym = nameof(operation(u))
-            u => collect(first(@variables $sym[collect(axes(g.second)[1] for g in grid)...](t)))
+            u => collect(first(@variables $sym[collect(axes(g.second)[1] for g in grid)...](t))) 
         end
     end
 
 
     # Build symbolic maps for boundaries
-    Iedge = reduce(vcat, [[vcat([Colon() for j = 1:i-1], 1, [Colon() for j = i+1:nspace]),
-        vcat([Colon() for j = 1:i-1], length(axies[i].second), [Colon() for j = i+1:nspace])] for i = 1:nspace])
+    Iedge = reduce(vcat, [[Igrid[vcat([Colon() for j = 1:i-1], 1, [Colon() for j = i+1:nspace])],
+        Igrid[vcat([Colon() for j = 1:i-1], length(axies[i].second), [Colon() for j = i+1:nspace])]] for i = 1:nspace])
 
     nottime2dim = [nottime[i] => i for i in 1:nspace]
     dim2nottime = [i => nottime[i] for i in 1:nspace]
-    return DiscreteSpace{nspace,length(depvars)}(depvars, Dict(depvarsdisc), nottime, indvars, Dict(axies), Dict(grid), Dict(dxs), Iaxies, Igrid, Iedge)
+    return DiscreteSpace{nspace,length(depvars)}(depvars, Dict(depvarsdisc), nottime, indvars, Dict(axies), Dict(grid), Dict(dxs), Iaxies, Igrid, Iedge, nottime2dim)
 end
 
 

test/DiscreteSpace.jl

@@ -0,0 +1,69 @@
+using ModelingToolkit, MethodOfLines, DomainSets, Test
+
+@parameters x, y, t 
+@variables u(..) v(..)
+
+indvars = [x, y, t]
+nottime = [x, y]
+depvars = [u, v]
+
+t_min= 0.
+t_max = 2.0
+x_min = 0.
+x_max = 2.
+y_min = 0.
+y_max = 2.
+dx = 0.1; dy = 0.2
+order = 2
+
+domains = [t ∈ Interval(t_min, t_max), x ∈ Interval(x_min, x_max), y ∈ Interval(y_min,y_max)]
+
+@testset "Test 01: discretization of variables, center aligned grid" begin
+    # Test centered order
+    disc = MOLFiniteDifference([x=>dx, y=>dy], t; centered_order=order)
+
+    s = MethodOfLines.DiscreteSpace(domains, depvars, indvars, nottime, disc)
+
+    discx = x_min:dx:x_max
+    discy = y_min:dy:y_max
+
+    @test s.grid[x] == discx 
+    @test s.grid[y] == discy
+
+    @test s.axies[x] == s.grid[x]
+    @test s.axies[y] == s.grid[y]
+
+    @test s.discvars[u] == collect(first(@variables u[axes(discx), axes(discy)]))
+    @test s.discvars[v] == collect(first(@variables v[axes(discx), axes(discy)]))
+
+    @test CartesianIndex(1, 10) ∈ s.Iedge
+    @test all([I ∈ s.Igrid for I in s.Iedge])  
+    
+    @test s.Iaxies == s.Igrid
+    
+end
+
+@testset "Test 02: discretization of variables, edge aligned grid" begin
+    # Test centered order
+    disc = MOLFiniteDifference([x=>dx, y=>dy], t; centered_order=order, grid_align=edge_align)
+
+    s = MethodOfLines.DiscreteSpace(domains, depvars, indvars, nottime, disc)
+
+    discx = (x_min-dx/2): dx : (x_max+dx/2)
+    discy = (y_min-dy/2): dy : (y_max+dy/2)
+
+    @test s.grid[x] == discx
+    @test s.grid[y] == discy
+
+    @test s.axies[x] != s.grid[x]
+    @test s.axies[y] != s.grid[y]
+
+    @test s.discvars[u] == collect(first(@variables u[axes(discx), axes(discy)]))
+    @test s.discvars[v] == collect(first(@variables v[axes(discx), axes(discy)]))
+
+    @test CartesianIndex(1, 10) ∈ s.Iedge
+    @test all([I ∈ s.Igrid for I in s.Iedge])  
+    
+    @test s.Iaxies != s.Igrid
+
+end