diff --git a/src/MOLFiniteDifference/MOL_discretization.jl b/src/MOLFiniteDifference/MOL_discretization.jl index c12a0059c..f7991df69 100644 --- a/src/MOLFiniteDifference/MOL_discretization.jl +++ b/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]) diff --git a/src/derivative_operators/convolutions.jl b/src/derivative_operators/convolutions.jl index e13264067..87f1eed0e 100644 --- a/src/derivative_operators/convolutions.jl +++ b/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] diff --git a/test/MOL/MOL_1D_HigherOrder.jl b/test/MOL/MOL_1D_HigherOrder.jl index 33623c004..1516a7af2 100644 --- a/test/MOL/MOL_1D_HigherOrder.jl +++ b/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 \ No newline at end of file