diff --git a/src/SBP_coefficients/Mattsson2017.jl b/src/SBP_coefficients/Mattsson2017.jl index 9e6b2059..9397b2f4 100644 --- a/src/SBP_coefficients/Mattsson2017.jl +++ b/src/SBP_coefficients/Mattsson2017.jl @@ -35,7 +35,65 @@ end function first_derivative_coefficients(source::Mattsson2017, order::Int, T=Float64, mode=FastMode()) - if order == 2 + if order == 1 + left_boundary_plus = ( + DerivativeCoefficientRow{T,1,2}(SVector(T(-2), # -2 for continuously coupled = same + T(2), )), # 2 for continuously coupled = same + ) + right_boundary_plus = ( + DerivativeCoefficientRow{T,1,2}(SVector(T(0), + T(0), )), + ) + upper_coef_plus = SVector( T(1), ) + central_coef_plus = T(-1) + lower_coef_plus = SVector{0,T}() + left_weights = SVector( T(1//2), ) # 1//2 for continuously coupled = same + right_weights = left_weights + left_boundary_derivatives = Tuple{}() + right_boundary_derivatives = left_boundary_derivatives + + left_boundary_minus = ( + DerivativeCoefficientRow{T,1,2}(SVector(T(0), + T(0), )), + ) + right_boundary_minus = ( + DerivativeCoefficientRow{T,1,2}(SVector(T(2), # 2 for continuously coupled = same + T(-2), )), # -2 for continuously coupled = same + ) + upper_coef_minus = .- lower_coef_plus + central_coef_minus = .- central_coef_plus + lower_coef_minus = .- upper_coef_plus + + left_boundary_central = (left_boundary_plus .+ left_boundary_minus) ./ 2 + right_boundary_central = (right_boundary_plus .+ right_boundary_minus) ./ 2 + upper_coef_central = widening_plus(upper_coef_plus, upper_coef_minus) / 2 + central_coef_central = (central_coef_plus + central_coef_minus) / 2 + lower_coef_central = widening_plus(lower_coef_plus, lower_coef_minus) / 2 + + if source.kind === :plus + left_boundary = left_boundary_plus + right_boundary = right_boundary_plus + upper_coef = upper_coef_plus + central_coef = central_coef_plus + lower_coef = lower_coef_plus + elseif source.kind === :minus + left_boundary = left_boundary_minus + right_boundary = right_boundary_minus + upper_coef = upper_coef_minus + central_coef = central_coef_minus + lower_coef = lower_coef_minus + elseif source.kind === :central + left_boundary = left_boundary_central + right_boundary = right_boundary_central + upper_coef = upper_coef_central + central_coef = central_coef_central + lower_coef = lower_coef_central + end + DerivativeCoefficients(left_boundary, right_boundary, + left_boundary_derivatives, right_boundary_derivatives, + lower_coef, central_coef, upper_coef, + left_weights, right_weights, mode, 1, order, source) + elseif order == 2 left_boundary_plus = ( DerivativeCoefficientRow{T,1,3}(SVector(T(-3), T(5), diff --git a/test/upwind_operators_test.jl b/test/upwind_operators_test.jl index 732e53d3..ddd5ae10 100644 --- a/test/upwind_operators_test.jl +++ b/test/upwind_operators_test.jl @@ -2,52 +2,67 @@ using Test using LinearAlgebra using SummationByPartsOperators +source_list = (Mattsson2017,) + # check construction of interior part of upwind operators @testset "Check interior parts" begin N = 21 - xmin = 0. + xmin = 0.0 xmax = Float64(N + 1) interior = 10:N-10 + for Source in source_list + for acc_order in 1:9 + Dp_bounded, Dm_bounded, Dc_bounded = try + Dp_bounded = derivative_operator(Source(:plus ), 1, acc_order, xmin, xmax, N) + Dm_bounded = derivative_operator(Source(:minus ), 1, acc_order, xmin, xmax, N) + Dc_bounded = derivative_operator(Source(:central), 1, acc_order, xmin, xmax, N) + Dp_bounded, Dm_bounded, Dc_bounded + catch + Dp_bounded = nothing + Dm_bounded = nothing + Dc_bounded = nothing + Dp_bounded, Dm_bounded, Dc_bounded + end + + Dp_bounded === nothing && continue - for acc_order in 2:9 - Dp_bounded = derivative_operator(Mattsson2017(:plus ), 1, acc_order, xmin, xmax, N) - Dm_bounded = derivative_operator(Mattsson2017(:minus ), 1, acc_order, xmin, xmax, N) - Dc_bounded = derivative_operator(Mattsson2017(:central), 1, acc_order, xmin, xmax, N) - for compact in (true, false) - show(IOContext(devnull, :compact=>compact), Dp_bounded) - show(IOContext(devnull, :compact=>compact), Dm_bounded) - show(IOContext(devnull, :compact=>compact), Dc_bounded) - summary(IOContext(devnull, :compact=>compact), Dp_bounded) - summary(IOContext(devnull, :compact=>compact), Dm_bounded) - summary(IOContext(devnull, :compact=>compact), Dc_bounded) - end - M = mass_matrix(Dp_bounded) - @test M == mass_matrix(Dm_bounded) - @test M == mass_matrix(Dc_bounded) - Dp_periodic = periodic_derivative_operator(1, acc_order, xmin, xmax, N-1, -(acc_order - 1) ÷ 2) - Dm_periodic = periodic_derivative_operator(1, acc_order, xmin, xmax, N-1, -acc_order + (acc_order - 1) ÷ 2) - Dp = Matrix(Dp_bounded) - Dm = Matrix(Dm_bounded) - @test Dp[interior,interior] ≈ Matrix(Dp_periodic)[interior,interior] - @test Dm[interior,interior] ≈ Matrix(Dm_periodic)[interior,interior] - res = M * Dp + Dm' * M - res[1,1] += 1 - res[end,end] -= 1 - @test norm(res) < N * eps() - x = grid(Dp_bounded) - for D in (Dp_bounded, Dm_bounded, Dc_bounded) - @test norm(D * x.^0) < N * eps() - for k in 1:acc_order÷2 - @test D * x.^k ≈ k .* x.^(k-1) + for compact in (true, false) + show(IOContext(devnull, :compact => compact), Dp_bounded) + show(IOContext(devnull, :compact => compact), Dm_bounded) + show(IOContext(devnull, :compact => compact), Dc_bounded) + summary(IOContext(devnull, :compact => compact), Dp_bounded) + summary(IOContext(devnull, :compact => compact), Dm_bounded) + summary(IOContext(devnull, :compact => compact), Dc_bounded) end - for k in acc_order÷2+1:acc_order - @test (D * x.^k)[interior] ≈ (k .* x.^(k-1))[interior] + + M = mass_matrix(Dp_bounded) + @test M == mass_matrix(Dm_bounded) + @test M == mass_matrix(Dc_bounded) + Dp_periodic = periodic_derivative_operator(1, acc_order, xmin, xmax, N-1, -(acc_order - 1) ÷ 2) + Dm_periodic = periodic_derivative_operator(1, acc_order, xmin, xmax, N-1, -acc_order + (acc_order - 1) ÷ 2) + Dp = Matrix(Dp_bounded) + Dm = Matrix(Dm_bounded) + @test Dp[interior,interior] ≈ Matrix(Dp_periodic)[interior,interior] + @test Dm[interior,interior] ≈ Matrix(Dm_periodic)[interior,interior] + res = M * Dp + Dm' * M + res[1,1] += 1 + res[end,end] -= 1 + @test norm(res) < N * eps() + x = grid(Dp_bounded) + for D in (Dp_bounded, Dm_bounded, Dc_bounded) + @test norm(D * x.^0) < N * eps() + for k in 1:acc_order÷2 + @test D * x.^k ≈ k .* x.^(k-1) + end + for k in acc_order÷2+1:acc_order + @test (D * x.^k)[interior] ≈ (k .* x.^(k-1))[interior] + end end + diss = M * (Dp - Dm) + @test diss ≈ diss' + @test maximum(eigvals(Symmetric(diss))) < N * eps() end - diss = M * (Dp - Dm) - @test diss ≈ diss' - @test maximum(eigvals(Symmetric(diss))) < N * eps() end end