From 4ffdcfc3334bb832476b96bb8dd70a45151ed92d Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Wed, 25 Mar 2020 13:17:28 +0300 Subject: [PATCH] upwind operators of Mattsson(2017), order 5 --- dev/Mattsson2017.jl | 22 +++++ src/SBP_coefficients/Mattsson2017.jl | 138 +++++++++++++++++++++++++++ test/upwind_operators.jl | 2 +- 3 files changed, 161 insertions(+), 1 deletion(-) diff --git a/dev/Mattsson2017.jl b/dev/Mattsson2017.jl index 36b2ad49..ab0470b5 100644 --- a/dev/Mattsson2017.jl +++ b/dev/Mattsson2017.jl @@ -50,3 +50,25 @@ m = [49//144, 61//48, 41//48, 149//144] M = Diagonal(vcat(m, ones(Int, size(Qp,1)-2*length(m)), m[end:-1:1])) Dp = M \ (Qp + B / 2) Dm = M \ (Qm + B / 2) + +# order 5 +Qp = [ + -1//120 941//1440 -47//360 -7//480 0 0 0 0 0 0 0 + -869//1440 -11//120 25//32 -43//360 1//30 0 0 0 0 0 0 + 29//360 -17//32 -29//120 1309//1440 -1//4 1//30 0 0 0 0 0 + 1//32 -11//360 -661//1440 -13//40 1 -1//4 1//30 0 0 0 0 + 0 0 1//20 -1//2 -1//3 1 -1//4 1//30 0 0 0 + 0 0 0 1//20 -1//2 -1//3 1 -1//4 1//30 0 0 + 0 0 0 0 1//20 -1//2 -1//3 1 -1//4 1//30 0 + 0 0 0 0 0 1//20 -1//2 -13//40 1309//1440 -43//360 -7//480 + 0 0 0 0 0 0 1//20 -661//1440 -29//120 25//32 -47//360 + 0 0 0 0 0 0 0 -11//360 -17//32 -11//120 941//1440 + 0 0 0 0 0 0 0 1//32 29//360 -869//1440 -1//120 +] +Qm = Matrix(-Qp') +B = zero(Qp); B[1,1] = -1; B[end,end] = 1 +# m = [251//720, 299//240, 41//48, 149//144] +m = [251//720, 299//240, 211//240, 739//720] +M = Diagonal(vcat(m, ones(Int, size(Qp,1)-2*length(m)), m[end:-1:1])) +Dp = M \ (Qp + B / 2) +Dm = M \ (Qm + B / 2) diff --git a/src/SBP_coefficients/Mattsson2017.jl b/src/SBP_coefficients/Mattsson2017.jl index 77d73550..22b01b52 100644 --- a/src/SBP_coefficients/Mattsson2017.jl +++ b/src/SBP_coefficients/Mattsson2017.jl @@ -288,6 +288,144 @@ function first_derivative_coefficients(source::Mattsson2017, order::Int, T=Float 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, parallel, 1, order, source) + elseif order == 5 + left_boundary_plus = ( + DerivativeCoefficientRow{T,1,4}(SVector(T(-366//251), + T(941//502), + T(-94//251), + T(-21//502), )), + DerivativeCoefficientRow{T,1,5}(SVector(T(-869//1794), + T(-22//299), + T(375//598), + T(-86//897), + T(8//299), )), + DerivativeCoefficientRow{T,1,6}(SVector(T(58//633), + T(-255//422), + T(-58//211), + T(1309//1266), + T(-60//211), + T(8//211), )), + DerivativeCoefficientRow{T,1,7}(SVector(T(45//1478), + T(-22//739), + T(-661//1478), + T(-234//739), + T(720//739), + T(-180//739), + T(24//739), )), + ) + right_boundary_plus = ( + DerivativeCoefficientRow{T,1,4}(SVector(T(354//251), + T(-869//502), + T(58//251), + T(45//502), )), + DerivativeCoefficientRow{T,1,4}(SVector(T(941//1794), + T(-22//299), + T(-255//598), + T(-22//897), )), + DerivativeCoefficientRow{T,1,5}(SVector(T(-94//633), + T(375//422), + T(-58//211), + T(-661//1266), + T(12//211), )), + DerivativeCoefficientRow{T,1,6}(SVector(T(-21//1478), + T(-86//739), + T(1309//1478), + T(-234//739), + T(-360//739), + T(36//739), )), + ) + upper_coef_plus = SVector( T(1), + T(-1//4), + T(1//30), ) + central_coef_plus = T(-1//3) + lower_coef_plus = SVector( T(-1//2), + T(1//20), ) + left_weights = SVector( T(251//720), + T(299//240), + T(211//240), + T(739//720), ) + right_weights = left_weights + left_boundary_derivatives = Tuple{}() + right_boundary_derivatives = left_boundary_derivatives + + left_boundary_minus = ( + DerivativeCoefficientRow{T,1,4}(SVector(T(-354//251), + T(869//502), + T(-58//251), + T(-45//502), )), + DerivativeCoefficientRow{T,1,4}(SVector(T(-941//1794), + T(22//299), + T(255//598), + T(22//897), )), + DerivativeCoefficientRow{T,1,5}(SVector(T(94//633), + T(-375//422), + T(58//211), + T(661//1266), + T(-12//211), )), + DerivativeCoefficientRow{T,1,6}(SVector(T(21//1478), + T(86//739), + T(-1309//1478), + T(234//739), + T(360//739), + T(-36//739), )), + ) + right_boundary_minus = ( + DerivativeCoefficientRow{T,1,4}(SVector(T(366//251), + T(-941//502), + T(94//251), + T(21//502), )), + DerivativeCoefficientRow{T,1,5}(SVector(T(869//1794), + T(22//299), + T(-375//598), + T(86//897), + T(-8//299), )), + DerivativeCoefficientRow{T,1,6}(SVector(T(-58//633), + T(255//422), + T(58//211), + T(-1309//1266), + T(60//211), + T(-8//211), )), + DerivativeCoefficientRow{T,1,7}(SVector(T(-45//1478), + T(22//739), + T(661//1478), + T(234//739), + T(-720//739), + T(180//739), + T(-24//739), )), + ) + 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 diff --git a/test/upwind_operators.jl b/test/upwind_operators.jl index a498c238..54c89f73 100644 --- a/test/upwind_operators.jl +++ b/test/upwind_operators.jl @@ -10,7 +10,7 @@ using SummationByPartsOperators interior = 10:N-10 - for acc_order in 2:4 + for acc_order in 2:5 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)