From 5ca14b1d4b5e16315c55663e38e25f78fb8e6c3a Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Thu, 24 Nov 2022 09:48:43 +0100 Subject: [PATCH 1/4] WIP: first-order upwind SBP coefficients --- src/SBP_coefficients/wip.jl | 102 +++++++++++++++++++++++++++++++ src/SummationByPartsOperators.jl | 3 +- 2 files changed, 104 insertions(+), 1 deletion(-) create mode 100644 src/SBP_coefficients/wip.jl diff --git a/src/SBP_coefficients/wip.jl b/src/SBP_coefficients/wip.jl new file mode 100644 index 00000000..e83924b7 --- /dev/null +++ b/src/SBP_coefficients/wip.jl @@ -0,0 +1,102 @@ + +""" + WIP(version::Symbol) + +Coefficients of the first-order upwind SBP operators given in +- WIP + +Higher-order coefficients are taken from +- Mattsson (2017) + Diagonal-norm upwind SBP operators. + Journal of Computational Physics 335, pp. 283-310. + +You can choose between the different versions `:central`, `:plus`, and `:minus`. +""" +struct WIP <: SourceOfCoefficients + kind::Symbol + + function WIP(kind::Symbol) + if (kind !== :plus) && (kind !== :minus) && (kind !== :central) + throw(ArgumentError("The only choices are :plus, :minus, and :central, not :$kind.")) + end + new(kind) + end +end + +# function Base.show(io::IO, source::WIP) +# if get(io, :compact, false) +# summary(io, source) +# else +# print(io, +# "Mattsson (2017) \n", +# " Diagonal-norm upwind SBP operators. \n", +# " Journal of Computational Physics 335, pp. 283-310. \n", +# " (upwind coefficients ", source.kind, ")") +# end +# end + + +function first_derivative_coefficients(source::WIP, order::Int, + T=Float64, mode=FastMode()) + if order == 1 + left_boundary_plus = ( + DerivativeCoefficientRow{T,1,2}(SVector(T(-1), # -2 for continuously coupled = same + T(1), )), # 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), ) # 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(1), # 2 for continuously coupled = same + T(-1), )), # -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) + else + return first_derivative_coefficients(Mattsson2017(source.kind), order, T, mode) + end +end diff --git a/src/SummationByPartsOperators.jl b/src/SummationByPartsOperators.jl index e2010dd0..65816db2 100644 --- a/src/SummationByPartsOperators.jl +++ b/src/SummationByPartsOperators.jl @@ -112,6 +112,7 @@ include("SBP_coefficients/Mattsson2017.jl") include("SBP_coefficients/MattssonAlmquistVanDerWeide2018Minimal.jl") include("SBP_coefficients/MattssonAlmquistVanDerWeide2018Accurate.jl") include("SBP_coefficients/DienerDorbandSchnetterTiglio2007.jl") +include("SBP_coefficients/wip.jl") include("conservation_laws/general_laws.jl") include("conservation_laws/burgers.jl") @@ -150,7 +151,7 @@ export Fornberg1998, Holoborodko2008, BeljaddLeFlochMishraParés2017 export MattssonNordström2004, MattssonSvärdNordström2004, MattssonSvärdShoeybi2008, Mattsson2012, Mattsson2014, MattssonAlmquistCarpenter2014Extended, MattssonAlmquistCarpenter2014Optimal, - Mattsson2017, + Mattsson2017, WIP, MattssonAlmquistVanDerWeide2018Minimal, MattssonAlmquistVanDerWeide2018Accurate, DienerDorbandSchnetterTiglio2007 export Tadmor1989, MadayTadmor1989, Tadmor1993, From 03011fbbbf2ca30fb5ffec3a6903ef4f6185a04e Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Thu, 24 Nov 2022 10:41:56 +0100 Subject: [PATCH 2/4] add tests --- test/upwind_operators_test.jl | 80 ++++++++++++++++++++--------------- 1 file changed, 47 insertions(+), 33 deletions(-) diff --git a/test/upwind_operators_test.jl b/test/upwind_operators_test.jl index b8c62990..a893909e 100644 --- a/test/upwind_operators_test.jl +++ b/test/upwind_operators_test.jl @@ -2,48 +2,62 @@ using Test using LinearAlgebra using SummationByPartsOperators +source_list = (Mattsson2017, WIP,) + # 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:7 + 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:7 - 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=>false), Dp_bounded) - show(IOContext(devnull, :compact=>false), Dm_bounded) - show(IOContext(devnull, :compact=>false), 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=>false), Dp_bounded) + show(IOContext(devnull, :compact=>false), Dm_bounded) + show(IOContext(devnull, :compact=>false), 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 From 601e9deb7f7f64f8633cbb78f42ac5e344f2a207 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Fri, 17 Mar 2023 16:44:27 +0100 Subject: [PATCH 3/4] put coefficients to the other upwind SBP coefficients --- src/SBP_coefficients/Mattsson2017.jl | 60 +++++++++++++++- src/SBP_coefficients/wip.jl | 102 --------------------------- src/SummationByPartsOperators.jl | 3 +- test/upwind_operators_test.jl | 4 +- 4 files changed, 62 insertions(+), 107 deletions(-) delete mode 100644 src/SBP_coefficients/wip.jl diff --git a/src/SBP_coefficients/Mattsson2017.jl b/src/SBP_coefficients/Mattsson2017.jl index 7916e0f2..c5a509d6 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(-1), # -2 for continuously coupled = same + T(1), )), # 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), ) # 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(1), # 2 for continuously coupled = same + T(-1), )), # -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/src/SBP_coefficients/wip.jl b/src/SBP_coefficients/wip.jl deleted file mode 100644 index e83924b7..00000000 --- a/src/SBP_coefficients/wip.jl +++ /dev/null @@ -1,102 +0,0 @@ - -""" - WIP(version::Symbol) - -Coefficients of the first-order upwind SBP operators given in -- WIP - -Higher-order coefficients are taken from -- Mattsson (2017) - Diagonal-norm upwind SBP operators. - Journal of Computational Physics 335, pp. 283-310. - -You can choose between the different versions `:central`, `:plus`, and `:minus`. -""" -struct WIP <: SourceOfCoefficients - kind::Symbol - - function WIP(kind::Symbol) - if (kind !== :plus) && (kind !== :minus) && (kind !== :central) - throw(ArgumentError("The only choices are :plus, :minus, and :central, not :$kind.")) - end - new(kind) - end -end - -# function Base.show(io::IO, source::WIP) -# if get(io, :compact, false) -# summary(io, source) -# else -# print(io, -# "Mattsson (2017) \n", -# " Diagonal-norm upwind SBP operators. \n", -# " Journal of Computational Physics 335, pp. 283-310. \n", -# " (upwind coefficients ", source.kind, ")") -# end -# end - - -function first_derivative_coefficients(source::WIP, order::Int, - T=Float64, mode=FastMode()) - if order == 1 - left_boundary_plus = ( - DerivativeCoefficientRow{T,1,2}(SVector(T(-1), # -2 for continuously coupled = same - T(1), )), # 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), ) # 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(1), # 2 for continuously coupled = same - T(-1), )), # -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) - else - return first_derivative_coefficients(Mattsson2017(source.kind), order, T, mode) - end -end diff --git a/src/SummationByPartsOperators.jl b/src/SummationByPartsOperators.jl index f952e9d7..9645c6f7 100644 --- a/src/SummationByPartsOperators.jl +++ b/src/SummationByPartsOperators.jl @@ -115,7 +115,6 @@ include("SBP_coefficients/Mattsson2017.jl") include("SBP_coefficients/MattssonAlmquistVanDerWeide2018Minimal.jl") include("SBP_coefficients/MattssonAlmquistVanDerWeide2018Accurate.jl") include("SBP_coefficients/DienerDorbandSchnetterTiglio2007.jl") -include("SBP_coefficients/wip.jl") include("conservation_laws/general_laws.jl") include("conservation_laws/burgers.jl") @@ -156,7 +155,7 @@ export Fornberg1998, Holoborodko2008, BeljaddLeFlochMishraParés2017 export MattssonNordström2004, MattssonSvärdNordström2004, MattssonSvärdShoeybi2008, Mattsson2012, Mattsson2014, MattssonAlmquistCarpenter2014Extended, MattssonAlmquistCarpenter2014Optimal, - Mattsson2017, WIP, + Mattsson2017, MattssonAlmquistVanDerWeide2018Minimal, MattssonAlmquistVanDerWeide2018Accurate, DienerDorbandSchnetterTiglio2007 export Tadmor1989, MadayTadmor1989, Tadmor1993, diff --git a/test/upwind_operators_test.jl b/test/upwind_operators_test.jl index af519195..ddd5ae10 100644 --- a/test/upwind_operators_test.jl +++ b/test/upwind_operators_test.jl @@ -2,7 +2,7 @@ using Test using LinearAlgebra using SummationByPartsOperators -source_list = (Mattsson2017, WIP,) +source_list = (Mattsson2017,) # check construction of interior part of upwind operators @testset "Check interior parts" begin @@ -12,7 +12,7 @@ source_list = (Mattsson2017, WIP,) interior = 10:N-10 for Source in source_list - for acc_order in 1:7 + 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) From 71dfad64d3e211f437cff52581211f6a37e9b6b6 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Fri, 17 Mar 2023 17:47:19 +0100 Subject: [PATCH 4/4] choose another first-order operator --- src/SBP_coefficients/Mattsson2017.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/SBP_coefficients/Mattsson2017.jl b/src/SBP_coefficients/Mattsson2017.jl index fcf21324..9397b2f4 100644 --- a/src/SBP_coefficients/Mattsson2017.jl +++ b/src/SBP_coefficients/Mattsson2017.jl @@ -37,8 +37,8 @@ function first_derivative_coefficients(source::Mattsson2017, order::Int, T=Float64, mode=FastMode()) if order == 1 left_boundary_plus = ( - DerivativeCoefficientRow{T,1,2}(SVector(T(-1), # -2 for continuously coupled = same - T(1), )), # 2 for continuously coupled = same + 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), @@ -47,7 +47,7 @@ function first_derivative_coefficients(source::Mattsson2017, order::Int, upper_coef_plus = SVector( T(1), ) central_coef_plus = T(-1) lower_coef_plus = SVector{0,T}() - left_weights = SVector( T(1), ) # 1//2 for continuously coupled = same + 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 @@ -57,8 +57,8 @@ function first_derivative_coefficients(source::Mattsson2017, order::Int, T(0), )), ) right_boundary_minus = ( - DerivativeCoefficientRow{T,1,2}(SVector(T(1), # 2 for continuously coupled = same - T(-1), )), # -2 for continuously coupled = same + 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