MOLFiniteDifference : Adding support for higher derivative orders in bcs

src/MOLFiniteDifference/MOL_discretization.jl

@@ -155,23 +155,29 @@ function SciMLBase.symbolic_discretize(pdesys::ModelingToolkit.PDESystem,discret
             # depvarderivbcmaps will dictate what to replace the Differential terms with in the bcs
             if nspace == 1
                 # 1D system
-                left_weights(j) = DiffEqOperators.calculate_weights(discretization.upwind_order, grid[j][1], grid[j][1:2])
-                right_weights(j) = DiffEqOperators.calculate_weights(discretization.upwind_order, grid[j][end], grid[j][end-1:end])
-                central_neighbor_idxs(II,j) = [II-CartesianIndex((1:nspace.==j)...),II,II+CartesianIndex((1:nspace.==j)...)]
-                left_idxs = central_neighbor_idxs(CartesianIndex(2),1)[1:2]
-                right_idxs(j) = central_neighbor_idxs(CartesianIndex(length(grid[j])-1),1)[end-1:end]
+                bc_der_orders = filter(!iszero,sort(unique([count_differentials(bc.lhs, nottime[1]) for bc in pdesys.bcs])))
+                left_weights(d_order,j) = DiffEqOperators.calculate_weights(d_order, grid[j][1], grid[j][1:1+d_order])
+                right_weights(d_order,j) = DiffEqOperators.calculate_weights(d_order, grid[j][end], grid[j][end-d_order:end])
+                # central_neighbor_idxs(II,j) = [II-CartesianIndex((1:nspace.==j)...),II,II+CartesianIndex((1:nspace.==j)...)]
+                left_idxs(d_order) = CartesianIndices(grid[1])[1:1+d_order]
+                right_idxs(d_order,j) = CartesianIndices(grid[j])[end-d_order:end]
                 # Constructs symbolic spatially discretized terms of the form e.g. au₂ - bu₁
-                derivars = [[dot(left_weights(1),depvar[left_idxs]), dot(right_weights(1),depvar[right_idxs(1)])]
-                for depvar in depvarsdisc]
+                derivars = [[[dot(left_weights(d,1),depvar[left_idxs(d)]), dot(right_weights(d,1),depvar[right_idxs(d,1)])] for d in bc_der_orders]
+                                for depvar in depvarsdisc]
                 # Create list of all the symbolic Differential terms evaluated at boundary e.g. Differential(x)(u(t,0))
-                subderivar(depvar,s) = substitute.((Differential(s)(depvar),),edgevals)
+                subderivar(depvar,d,s) = substitute.(((Differential(s)^d)(depvar),),edgevals)
                 # Create map of symbolic Differential terms with symbolic spatially discretized terms
-                depvarderivbcmaps = reduce(vcat,[subderivar(depvar, s) .=> derivars[i]
-                                                for (i, depvar) in enumerate(depvars) for s in nottime])
+                depvarderivbcmaps = []
+                for (k,s) in enumerate(nottime)
+                    rs = (subderivar(depvar,bc_der_orders[j],s) .=> derivars[i][j] for j in 1:length(bc_der_orders), (i,depvar) in enumerate(depvars))
+                    for r in rs
+                        push!(depvarderivbcmaps, r)
+                    end
+                end
 
                 if grid_align == edge_align
                     # Constructs symbolic spatially discretized terms of the form e.g. (u₁ + u₂) / 2 
-                    bcvars = [[dot(ones(2)/2,depvar[left_idxs]), dot(ones(2)/2,depvar[right_idxs(1)])]
+                    bcvars = [[dot(ones(2)/2,depvar[left_idxs(1)]), dot(ones(2)/2,depvar[right_idxs(1,1)])]
                             for depvar in depvarsdisc]
                     # replace u(t,0) with (u₁ + u₂) / 2, etc
                     depvarbcmaps = reduce(vcat,[subvar(depvar) .=> bcvars[i]
@@ -183,7 +189,10 @@ function SciMLBase.symbolic_discretize(pdesys::ModelingToolkit.PDESystem,discret
                 # TODO: Fix Neumann and Robin on higher dimension
                 depvarderivbcmaps = []
             end
-
+            # All unique order of derivates in BCs
+            bc_der_orders = filter(!iszero,sort(unique([count_differentials(bc.lhs, nottime[1]) for bc in pdesys.bcs])))
+            # no. of different orders in BCs 
+            n = length(bc_der_orders)
             # Generate initial conditions and bc equations
             for bc in pdesys.bcs
                 bcdepvar = first(get_depvars(bc.lhs, depvar_ops))
@@ -202,7 +211,10 @@ function SciMLBase.symbolic_discretize(pdesys::ModelingToolkit.PDESystem,discret
                             bcargs = arguments(first(depvarbcmaps[i]))
                             # Replace Differential terms in the bc lhs with the symbolic spatially discretized terms
                             # TODO: Fix Neumann and Robin on higher dimension
-                            lhs = nspace == 1 ? substitute(bc.lhs,depvarderivbcmaps[i]) : bc.lhs
+                            # Update: Fixed for 1D systems
+                            j = findfirst(isequal(count_differentials(bc.lhs, nottime[1])),bc_der_orders)
+                            k = i%2 == 0 ? 2 : 1
+                            lhs = nspace == 1 ? (j isa Nothing ? bc.lhs : substitute(bc.lhs,depvarderivbcmaps[j + n*Int(floor((i-1)/2))][k])) : bc.lhs
 
                             # Replace symbol in the bc lhs with the spatial discretized term
                             lhs = substitute(lhs,depvarbcmaps[i])

src/derivative_operators/convolutions.jl

@@ -117,7 +117,7 @@ function convolve_interior!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::
     for i in (1+A.boundary_point_count) : (length(x_temp)-A.boundary_point_count - A.offside)
         xtempi = zero(T)
         cur_stencil = stencil
-        cur_coeff   = typeof(coeff)   <: AbstractVector ? coeff[i] : coeff isa Number ? coeff : true
+        cur_coeff   = coeff[i]
         cur_stencil = cur_coeff >= 0 ? cur_stencil : A.derivative_order % 2 == 0 ? reverse(cur_stencil) : -1*reverse(cur_stencil)
         for idx in 1:A.stencil_length
             x_idx = cur_coeff < 0 ? x[i - A.stencil_length + 1 + idx + A.offside] : x[i + idx - A.offside]

test/MOL/MOL_1D_Linear_Diffusion.jl

@@ -59,6 +59,37 @@ using ModelingToolkit: Differential
             @test all(isapprox.(u_approx, exact, atol=0.01))
         end
     end
+
+    @parameters t x
+
+    # BCs with 2nd order Neumann
+    bcs = [u(0,x) ~ cos(x),
+           Dxx(u(t,0)) ~ -exp(-t),
+           Dxx(u(t,Float64(π))) ~ exp(-t)]
+
+    pdesys = PDESystem(eq,bcs,domains,[t,x],[u(t,x)])
+
+    for disc in [discretization, discretization_edge, discretization_centered, discretization_centered_order4]
+        # Convert the PDE problem into an ODE problem
+        prob = discretize(pdesys,disc)
+
+        # Solve ODE problem
+        sol = solve(prob,Tsit5(),saveat=0.1)
+
+        if disc.grid_align == center_align
+            x = dx[2:end-1]
+        else
+            x = (dx[1:end-1]+dx[2:end])/2
+        end
+        t = sol.t
+
+        # Test against exact solution
+        for i in 1:length(sol)
+            exact = u_exact(x, t[i])
+            u_approx = sol.u[i]
+            @test all(isapprox.(u_approx, exact, atol=0.01))
+        end
+    end
 end
 
 
@@ -496,7 +527,7 @@ end
     end
 end
 
-@testset "Test 10: linear diffusion, two variables, mixed BCs" begin
+@testset "Test 10: linear diffusion, two variables, mixed & higher order BCs" begin
     # Method of Manufactured Solutions
     u_exact = (x,t) -> exp.(-t) * cos.(x)
     v_exact = (x,t) -> exp.(-t) * sin.(x)
@@ -507,6 +538,7 @@ end
     Dt = Differential(t)
     Dx = Differential(x)
     Dxx = Dx^2
+    Dxxx = Dx^3
 
     # 1D PDE and boundary conditions
     eqs = [Dt(u(t,x)) ~ Dxx(u(t,x)),
@@ -545,6 +577,33 @@ end
         @test all(isapprox.(u_exact(x_sol, t_sol[i]), sol.u[i][1:l-2], atol=0.01))
         @test all(isapprox.(v_exact(x_sol, t_sol[i]), sol.u[i][l-1:end], atol=0.01))
     end
+
+    @parameters t x
+
+    # BCs with 3rd order Neumann
+    bcs = [u(0,x) ~ cos(x),
+           v(0,x) ~ sin(x),
+           u(t,0) ~ exp(-t),
+           Dxxx(u(t,1)) ~ exp(-t) * sin(1),
+           Dxxx(v(t,0)) ~ -exp(-t),
+           v(t,1) ~ exp(-t) * sin(1)]
+
+    pdesys = PDESystem(eqs,bcs,domains,[t,x],[u(t,x),v(t,x)])
+
+    prob = discretize(pdesys,discretization)
+
+    # Solve ODE problem
+    sol = solve(prob,Tsit5(),saveat=0.1)
+
+    x_sol = dx[2:end-1]
+    t_sol = sol.t
+
+    # Test against exact solution
+    for i in 1:length(sol)
+        @test all(isapprox.(u_exact(x_sol, t_sol[i]), sol.u[i][1:l-2], atol=0.01))
+        @test all(isapprox.(v_exact(x_sol, t_sol[i]), sol.u[i][l-1:end], atol=0.01))
+    end
+
 end
 
 @testset "Test 11: linear diffusion, two variables, mixed BCs, with parameters" begin