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_HigherOrder.jl

@@ -74,28 +74,28 @@ end
     prob = discretize(pdesys,discretization)
 end
 
-@testset "Kuramoto–Sivashinsky equation" begin
+@testset "KdV Single Soliton 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
+    α = 6
+    β = 1
+    eq = Dt(u(x,t)) ~ -α*u(x,t)*Dx(u(x,t)) - β*Dx3(u(x,t)) 
 
+    u_analytic(x,t;z = (x - t)/2) =  1/2*sech(z)^2
+    du(x,t;z = (x-t)/2) = 1/2*tanh(z)*sech(z)^2
+    ddu(x,t; z= (x-t)/2) = 1/4*(2*tanh(z)^2 + sech(z)^2)*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)]
+           Dx(u(10,t)) ~ du(10,t),
+           Dx2(u(-10,t)) ~ ddu(-10,t),
+           Dx2(u(10,t)) ~ ddu(10,t)]
 
     # Space and time domains
     domains = [x ∈ Interval(-10.0,10.0),
@@ -111,12 +111,10 @@ end
 
     @test sol.retcode == :Success
 
-    xs = domains[1].domain.lower+dx+dx:dx:domains[1].domain.upper-dx-dx
+    xs = domains[1].domain.lower+dx+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 all(isapprox.(u_predict, u_real, rtol = 0.03))
+end