Generalizing MOLFiniteDifference to N-order PDEs

Project.toml

@@ -18,6 +18,7 @@ SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 SuiteSparse = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"
+SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
 
 [compat]
 BandedMatrices = "0.15.11, 0.16"

src/MOLFiniteDifference/MOL_discretization.jl

@@ -68,7 +68,9 @@ function SciMLBase.symbolic_discretize(pdesys::ModelingToolkit.PDESystem,discret
     edges = reduce(vcat,[[vcat([Colon() for j in 1:i-1],1,[Colon() for j in i+1:nspace]),
                           vcat([Colon() for j in 1:i-1],length(space[i]),[Colon() for j in i+1:nspace])] for i in 1:nspace])
 
-    edgevals = reduce(vcat,[[nottime[i]=>first(space[i]),nottime[i]=>last(space[i])] for i in 1:length(space)])
+    #edgeindices = [indices[e...] for e in edges]
+    get_edgevals(i) = [nottime[i]=>first(space[i]),nottime[i]=>last(space[i])]
+    edgevals = reduce(vcat,[get_edgevals(i) for i in 1:length(space)])
     edgevars = [[d[e...] for e in edges] for d in depvarsdisc]
 
     bclocs = map(e->substitute.(pdesys.indvars,e),edgevals) # location of the boundary conditions e.g. (t,0.0,y)
@@ -141,31 +143,73 @@ function SciMLBase.symbolic_discretize(pdesys::ModelingToolkit.PDESystem,discret
     u0 = reduce(vcat,u0)
     bceqs = reduce(vcat,bceqs)
 
+    #---- 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(depvar, i) = substitute.((depvar,),get_edgevals(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.
+    get_bc_counts(i) =
+        begin
+            left = 0
+            right = 0
+            for depvar in depvars
+                depvaredges = get_depvarbcs(depvar, i)
+                counts = [map(x->occursin(x.val, bc.lhs), depvaredges) for bc in pdesys.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
+    #--------------------------------------------------
+
     ### PDE EQUATIONS ###
-    interior = grid_indices[[2:length(g)-1 for g in grid]...]
+    interior = grid_indices[[let bcs = get_bc_counts(i)
+                             (1 + first(bcs)):length(g)-last(bcs)
+                             end
+                             for (i,g) in enumerate(grid)]...]
     eqs = vec(map(Base.product(interior,pdeeqs)) do p
         II,eq = p
 
         # 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
-        order = discretization.centered_order
-        stencil(j) = CartesianIndices(Tuple(map(x -> -x:x, (1:nspace.==j) * (order÷2))))
+        if discretization.centered_order % 2 != 0
+            throw(ArgumentError("Discretization centered_order must be even, given $(discretization.centered_order)"))
+        end
+        approx_order = discretization.centered_order
+        stencil(j) = CartesianIndices(Tuple(map(x -> -x:x, (1:nspace.==j) * (approx_order÷2))))
 
         # TODO: Generalize central difference handling to allow higher even order derivatives
         # The central neighbour indices should add the stencil to II, unless II is too close
         # to an edge in which case we need to shift away from the edge
 
         # Calculate buffers
         I1 = oneunit(first(grid_indices))
-        Imin = first(grid_indices) + I1 * (order÷2)
-        Imax = last(grid_indices) - I1 * (order÷2)
+        Imin = first(grid_indices) + I1 * (approx_order÷2)
+        Imax = last(grid_indices) - I1 * (approx_order÷2)
         # Use max and min to apply buffers
         central_neighbor_idxs(II,j) = stencil(j) .+ max(Imin,min(II,Imax))
         central_neighbor_space(II,j) = vec(grid[j][map(i->i[j],central_neighbor_idxs(II,j))])
-        central_weights(II,j) = DiffEqOperators.calculate_weights(2, grid[j][II[j]], central_neighbor_space(II,j))
-        central_deriv(II,j,k) = dot(central_weights(II,j),depvarsdisc[k][central_neighbor_idxs(II,j)])
+        central_weights(d_order, II,j) = DiffEqOperators.calculate_weights(d_order, grid[j][II[j]], central_neighbor_space(II,j))
+        central_deriv(d_order, II,j,k) = dot(central_weights(d_order, II,j),depvarsdisc[k][central_neighbor_idxs(II,j)])
+
+        # get a sorted list derivative order such that highest order is first. This is useful when substituting rules
+        # starting from highest to lowest order.
+        d_orders(s) = reverse(sort(collect(union(differential_order(eq.rhs, s.val), differential_order(eq.lhs, s.val)))))
 
-        central_deriv_rules = [(Differential(s)^2)(u) => central_deriv(II,j,k) for (j,s) in enumerate(nottime), (k,u) in enumerate(depvars)]
+        # central_deriv_rules = [(Differential(s)^2)(u) => central_deriv(2,II,j,k) for (j,s) in enumerate(nottime), (k,u) in enumerate(depvars)]
+        central_deriv_rules = Array{Pair{Num,Num},1}()
+        for (j,s) in enumerate(nottime)
+            rs = [(Differential(s)^d)(u) => central_deriv(d,II,j,k) for d in d_orders(s), (k,u) in enumerate(depvars)]
+            for r in rs
+              push!(central_deriv_rules, r)
+            end
+        end
         valrules = vcat([depvars[k] => depvarsdisc[k][II] for k in 1:length(depvars)],
                         [nottime[j] => grid[j][II[j]] for j in 1:nspace])
 

src/utils.jl

@@ -1,3 +1,4 @@
+using SymbolicUtils
 
 """
 A function that creates a tuple of CartesianIndices of unit length and `N` dimensions, one pointing along each dimension.
@@ -21,6 +22,41 @@ function cartesian_to_linear(I::CartesianIndex, s)  #Not sure if there is a buil
     return out
 end
 
+# Counts the Differential operators for given variable x. This is used to determine
+# the order of a PDE.
+function count_differentials(term, x::Symbolics.Symbolic)
+    S = Symbolics
+    SU = SymbolicUtils
+    if !S.istree(term)
+        return 0
+    else
+        op = SU.operation(term)
+        count_children = sum(map(arg -> count_differentials(arg, x), SU.arguments(term)))
+        if op isa Differential && op.x === x
+            return 1 + count_children
+        end
+        return count_children
+    end
+end
+
+# return list of differential orders in the equation
+function differential_order(eq, x::Symbolics.Symbolic)
+    S = Symbolics
+    SU = SymbolicUtils
+    orders = Set{Int}()
+    if S.istree(eq)
+        op = SU.operation(eq)
+        if op isa Differential
+            push!(orders, count_differentials(eq, x))
+        else
+            for o in map(ch -> differential_order(ch, x), SU.arguments(eq))
+                union!(orders, o)
+            end
+        end
+    end
+    return filter(!iszero, orders)
+end
+
 add_dims(A::AbstractArray, n::Int; dims::Int = 1) = cat(ndims(A) + n, A, dims = dims)
 
 ""

test/MOL/MOL_1D_HigherOrder.jl

@@ -0,0 +1,122 @@
+# 1D diffusion problem
+
+# Packages and inclusions
+using ModelingToolkit,DiffEqOperators,LinearAlgebra,Test,OrdinaryDiffEq
+using ModelingToolkit: Differential
+
+# Beam Equation
+@test_broken begin
+    @parameters x, t
+    @variables u(..)
+    Dt = Differential(t)
+    Dtt = Differential(t)^2
+    Dx = Differential(x)
+    Dxx = Differential(x)^2
+    Dx3 = Differential(x)^3
+    Dx4 = Differential(x)^4
+
+    g = -9.81
+    EI = 1
+    mu = 1
+    L = 10.0
+    dx = 0.4
+
+    eq = Dtt(u(t,x)) ~ -mu*EI*Dx4(u(t,x)) + mu*g
+
+    bcs = [u(0, x) ~ 0,
+           u(t,0) ~ 0,
+           Dx(u(t,0)) ~ 0,
+           Dxx(u(t, L)) ~ 0,
+           Dx3(u(t, L)) ~ 0]
+
+    # Space and time domains
+    domains = [t ∈ IntervalDomain(0.0,1.0),
+               x ∈ IntervalDomain(0.0,L)]
+
+    pdesys = PDESystem(eq,bcs,domains,[t,x],[u(t,x)])
+    discretization = MOLFiniteDifference([x=>dx],t, centered_order=4)
+    prob = discretize(pdesys,discretization)
+end
+
+# Beam Equation with Velocity
+@test_broken begin
+    @parameters x, t
+    @variables u(..), v(..)
+    Dt = Differential(t)
+    Dtt = Differential(t)^2
+    Dx = Differential(x)
+    Dxx = Differential(x)^2
+    Dx3 = Differential(x)^3
+    Dx4 = Differential(x)^4
+
+    g = -9.81
+    EI = 1
+    mu = 1
+    L = 10.0
+    dx = 0.4
+
+    eqs = [v(t, x) ~ Dt(u(t,x)),
+           Dt(v(t,x)) ~ -mu*EI*Dx4(u(t,x)) + mu*g]
+
+    bcs = [u(0, x) ~ 0,
+           v(0, x) ~ 0,
+           u(t,0) ~ 0,
+           v(t,0) ~ 0,
+           Dxx(u(t, L)) ~ 0,
+           Dx3(u(t, L)) ~ 0]
+
+    # Space and time domains
+    domains = [t ∈ IntervalDomain(0.0,1.0),
+               x ∈ IntervalDomain(0.0,L)]
+
+    pdesys = PDESystem(eqs,bcs,domains,[t,x],[u(t,x),v(t,x)])
+    discretization = MOLFiniteDifference([x=>dx],t, centered_order=4)
+    prob = discretize(pdesys,discretization)
+end
+
+@testset "Kuramoto–Sivashinsky equation" begin
+    @parameters x, t
+    @variables u(..)
+    Dt = Differential(t)
+    Dx = Differential(x)
+    Dx2 = Differential(x)^2
+    Dx3 = Differential(x)^3
+    Dx4 = Differential(x)^4
+
+    α = 1
+    β = 4
+    γ = 1
+    eq = Dt(u(x,t)) ~ -u(x,t)*Dx(u(x,t)) - α*Dx2(u(x,t)) - β*Dx3(u(x,t)) - γ*Dx4(u(x,t))
+
+    u_analytic(x,t;z = -x/2+t) = 11 + 15*tanh(z) -15*tanh(z)^2 - 15*tanh(z)^3
+    du(x,t;z = -x/2+t) = 15/2*(tanh(z) + 1)*(3*tanh(z) - 1)*sech(z)^2
+
+    bcs = [u(x,0) ~ u_analytic(x,0),
+           u(-10,t) ~ u_analytic(-10,t),
+           u(10,t) ~ u_analytic(10,t),
+           Dx(u(-10,t)) ~ du(-10,t),
+           Dx(u(10,t)) ~ du(10,t)]
+
+    # Space and time domains
+    domains = [x ∈ IntervalDomain(-10.0,10.0),
+               t ∈ IntervalDomain(0.0,1.0)]
+    # Discretization
+    dx = 0.4; dt = 0.2
+
+    discretization = MOLFiniteDifference([x=>dx],t;centered_order=4,grid_align=center_align)
+    pdesys = PDESystem(eq,bcs,domains,[x,t],[u(x,t)])
+    prob = discretize(pdesys,discretization)
+
+    sol = solve(prob,Tsit5(),saveat=0.1,dt=dt)
+
+    @test sol.retcode == :Success
+
+    xs = domains[1].domain.lower+dx+dx:dx:domains[1].domain.upper-dx-dx
+    ts = sol.t
+
+    u_predict = sol.u
+    u_real = [[u_analytic(x, t) for x in xs] for t in ts]
+    u_diff = u_real - u_predict
+    @test_broken u_diff[:] ≈ zeros(length(u_diff)) atol=0.01;
+    #plot(xs, u_diff)
+end

test/MOL/MOL_1D_Linear_Diffusion.jl

@@ -478,4 +478,41 @@ end
     @test prob.p == [0.5,2]
     # Make sure it can be solved
     sol = solve(prob,Tsit5())
-end
+end
+
+@testset "Test 12: Test Invalid Centered Order" begin
+    # Method of Manufactured Solutions
+    u_exact = (x,t) -> exp.(-t) * cos.(x)
+
+    # Parameters, variables, and derivatives
+    @parameters t x
+    @variables u(..)
+    Dt = Differential(t)
+    Dxx = Differential(x)^2
+
+    # 1D PDE and boundary conditions
+    eq  = Dt(u(t,x)) ~ Dxx(u(t,x))
+    bcs = [u(0,x) ~ cos(x),
+           u(t,0) ~ exp(-t),
+           u(t,Float64(π)) ~ -exp(-t)]
+
+    # Space and time domains
+    domains = [t ∈ IntervalDomain(0.0,1.0),
+               x ∈ IntervalDomain(0.0,Float64(π))]
+
+    # PDE system
+    pdesys = PDESystem(eq,bcs,domains,[t,x],[u(t,x)])
+
+    # Method of lines discretization
+    dx = range(0.0,Float64(π),length=30)
+
+    # Explicitly specify and invalid order of centered difference
+    for order in 1:6
+        discretization = MOLFiniteDifference([x=>dx],t;centered_order=order)
+        if order % 2 != 0
+            @test_throws ArgumentError discretize(pdesys,discretization)
+        else
+            discretize(pdesys,discretization)
+        end
+    end
+end

test/Misc/utils.jl

@@ -1,9 +1,83 @@
 using Test, LinearAlgebra
 using DiffEqOperators
+using ModelingToolkit
 
 @testset "utility functions" begin
     @test DiffEqOperators.unit_indices(2) == (CartesianIndex(1,0), CartesianIndex(0,1))
     @test DiffEqOperators.add_dims(zeros(2,2), ndims(zeros(2,2)) + 2) == [6. 6.; 0. 0.; 0. 0.]
     @test DiffEqOperators.perpindex(collect(1:5), 3) == [1, 2, 4, 5]
     @test DiffEqOperators.perpsize(zeros(2,2,3,2), 3) == (2, 2, 2)
 end
+
+@testset "count differentials 1D" begin
+    @parameters t x
+    @variables u(..)
+    Dt = Differential(t)
+
+    Dx = Differential(x)
+    eq  = Dt(u(t,x)) ~ -Dx(u(t,x))
+    @test first(DiffEqOperators.differential_order(eq.rhs, x.val)) == 1
+    @test isempty(DiffEqOperators.differential_order(eq.rhs, t.val))
+    @test first(DiffEqOperators.differential_order(eq.lhs, t.val)) == 1
+    @test isempty(DiffEqOperators.differential_order(eq.lhs, x.val))
+
+    Dxx = Differential(x)^2
+    eq  = Dt(u(t,x)) ~ Dxx(u(t,x))
+    @test first(DiffEqOperators.differential_order(eq.rhs, x.val)) == 2
+    @test isempty(DiffEqOperators.differential_order(eq.rhs, t.val))
+    @test first(DiffEqOperators.differential_order(eq.lhs, t.val)) == 1
+    @test isempty(DiffEqOperators.differential_order(eq.lhs, x.val))
+
+    Dxxxx = Differential(x)^4
+    eq  = Dt(u(t,x)) ~ -Dxxxx(u(t,x))
+    @test first(DiffEqOperators.differential_order(eq.rhs, x.val)) == 4
+    @test isempty(DiffEqOperators.differential_order(eq.rhs, t.val))
+    @test first(DiffEqOperators.differential_order(eq.lhs, t.val)) == 1
+    @test isempty(DiffEqOperators.differential_order(eq.lhs, x.val))
+end
+
+@testset "count differentials 2D" begin
+    @parameters t x y
+    @variables u(..)
+    Dxx = Differential(x)^2
+    Dyy = Differential(y)^2
+    Dt = Differential(t)
+
+    eq  = Dt(u(t,x,y)) ~ Dxx(u(t,x,y)) + Dyy(u(t,x,y))
+    @test first(DiffEqOperators.differential_order(eq.rhs, x.val)) == 2
+    @test first(DiffEqOperators.differential_order(eq.rhs, y.val)) == 2
+    @test isempty(DiffEqOperators.differential_order(eq.rhs, t.val))
+    @test first(DiffEqOperators.differential_order(eq.lhs, t.val)) == 1
+    @test isempty(DiffEqOperators.differential_order(eq.lhs, x.val))
+    @test isempty(DiffEqOperators.differential_order(eq.lhs, y.val))
+end
+
+@testset "count with mixed terms" begin
+    @parameters t x y
+    @variables u(..)
+    Dxx = Differential(x)^2
+    Dyy = Differential(y)^2
+    Dx = Differential(x)
+    Dy = Differential(y)
+    Dt = Differential(t)
+
+    eq  = Dt(u(t,x,y)) ~ Dxx(u(t,x,y)) + Dyy(u(t,x,y)) + Dx(Dy(u(t,x,y)))
+    @test DiffEqOperators.differential_order(eq.rhs, x.val) == Set([2, 1])
+    @test DiffEqOperators.differential_order(eq.rhs, y.val) == Set([2, 1])
+end
+
+@testset "Kuramoto–Sivashinsky equation" begin
+    @parameters x, t
+    @variables u(..)
+    Dt = Differential(t)
+    Dx = Differential(x)
+    Dx2 = Differential(x)^2
+    Dx3 = Differential(x)^3
+    Dx4 = Differential(x)^4
+
+    α = 1
+    β = 4
+    γ = 1
+    eq = Dt(u(x,t)) + u(x,t)*Dx(u(x,t)) + α*Dx2(u(x,t)) + β*Dx3(u(x,t)) + γ*Dx4(u(x,t)) ~ 0
+    @test DiffEqOperators.differential_order(eq.lhs, x.val) == Set([4, 3, 2, 1])
+end