From 967228de33e9efddab0de3ad2b488dd5453ca779 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Wed, 25 Mar 2020 08:21:24 +0300 Subject: [PATCH 1/8] fix indentation --- .../MattssonNordstr\303\266m2004.jl" | 974 +++++++++--------- 1 file changed, 487 insertions(+), 487 deletions(-) diff --git "a/src/SBP_coefficients/MattssonNordstr\303\266m2004.jl" "b/src/SBP_coefficients/MattssonNordstr\303\266m2004.jl" index 161b9328..1b7aa475 100644 --- "a/src/SBP_coefficients/MattssonNordstr\303\266m2004.jl" +++ "b/src/SBP_coefficients/MattssonNordstr\303\266m2004.jl" @@ -1,12 +1,12 @@ """ -MattssonNordström2004 + MattssonNordström2004 Coefficients of the SBP operators given in -Mattsson, Nordström (2004) -Summation by parts operators for finite difference approximations of second - derivatives. -Journal of Computational Physics 199, pp.503-540. + Mattsson, Nordström (2004) + Summation by parts operators for finite difference approximations of second + derivatives. + Journal of Computational Physics 199, pp.503-540. """ struct MattssonNordström2004 <: SourceOfCoefficients end @@ -20,501 +20,501 @@ end function first_derivative_coefficients(source::MattssonNordström2004, order::Int, T=Float64, parallel=Val{:serial}()) -if order == 2 - left_boundary = ( - DerivativeCoefficientRow{T,1,2}(SVector(-one(T), one(T))), - ) - right_boundary = .- left_boundary - upper_coef = SVector(T(1//2)) - central_coef = zero(T) - lower_coef = -upper_coef - left_weights = SVector(T(1//2)) - right_weights = left_weights - left_boundary_derivatives = Tuple{}() - right_boundary_derivatives = left_boundary_derivatives + if order == 2 + left_boundary = ( + DerivativeCoefficientRow{T,1,2}(SVector(-one(T), one(T))), + ) + right_boundary = .- left_boundary + upper_coef = SVector(T(1//2)) + central_coef = zero(T) + lower_coef = -upper_coef + left_weights = SVector(T(1//2)) + right_weights = left_weights + left_boundary_derivatives = Tuple{}() + right_boundary_derivatives = left_boundary_derivatives - 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 == 4 - left_boundary = ( - # q1 - DerivativeCoefficientRow{T,1,4}(SVector(T(-24//17), - T(59//34), - T(-4//17), - T(-3//34) )), - # q2 - DerivativeCoefficientRow{T,1,3}(SVector(T(-1//2), - T(0), - T(1//2))), - # q3 - DerivativeCoefficientRow{T,1,5}(SVector(T(4//43), - T(-59//86), - T(0), - T(59//86), - T(-4//43))), - # q4 - DerivativeCoefficientRow{T,1,6}(SVector(T(3//98), - T(0), - T(-59//98), - T(0), - T(32//49), - T(-4//49))), - ) - right_boundary = .- left_boundary - upper_coef = SVector(T(2//3), T(-1//12)) - central_coef = zero(T) - lower_coef = -upper_coef - left_weights = SVector( T(17//48), - T(59//48), - T(43//48), - T(49//48) ) - right_weights = left_weights - left_boundary_derivatives = Tuple{}() - right_boundary_derivatives = left_boundary_derivatives + 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 == 4 + left_boundary = ( + # q1 + DerivativeCoefficientRow{T,1,4}(SVector(T(-24//17), + T(59//34), + T(-4//17), + T(-3//34) )), + # q2 + DerivativeCoefficientRow{T,1,3}(SVector(T(-1//2), + T(0), + T(1//2))), + # q3 + DerivativeCoefficientRow{T,1,5}(SVector(T(4//43), + T(-59//86), + T(0), + T(59//86), + T(-4//43))), + # q4 + DerivativeCoefficientRow{T,1,6}(SVector(T(3//98), + T(0), + T(-59//98), + T(0), + T(32//49), + T(-4//49))), + ) + right_boundary = .- left_boundary + upper_coef = SVector(T(2//3), T(-1//12)) + central_coef = zero(T) + lower_coef = -upper_coef + left_weights = SVector( T(17//48), + T(59//48), + T(43//48), + T(49//48) ) + right_weights = left_weights + left_boundary_derivatives = Tuple{}() + right_boundary_derivatives = left_boundary_derivatives - 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 == 6 - left_boundary = ( - # q1 - DerivativeCoefficientRow{T,1,6}(SVector(T(-21600//13649), - T(104009//54596), - T(30443//81894), - T(-33311//27298), - T(16863//27298), - T(-15025//163788) )), - # q2 - DerivativeCoefficientRow{T,1,6}(SVector(T(-104009//240260), - T(0), - T(-311//72078), - T(20229//24026), - T(-24337//48052), - T(36661//360390) )), - # q3 - DerivativeCoefficientRow{T,1,6}(SVector(T(-30443//162660), - T(311//32532), - T(0), - T(-11155//16266), - T(41287//32532), - T(-21999//54220) )), - # q4 #TODO - DerivativeCoefficientRow{T,1,7}(SVector(T(33311//107180), - T(-20229//21436), - T(485//1398), - T(0), - T(4147//21436), - T(25427//321540), - T(72//5359) )), - # q5 - DerivativeCoefficientRow{T,1,8}(SVector(T(-16863//78770), - T(24337//31508), - T(-41287//47262), - T(-4147//15754), - T(0), - T(342523//472620), - T(-1296//7877), - T(144//7877) )), - # q6 - DerivativeCoefficientRow{T,1,9}(SVector(T(15025//525612), - T(-36661//262806), - T(21999//87602), - T(-25427//262806), - T(-342523//525612), - T(0), - T(32400//43801), - T(-6480//43801), - T(720//43801) )) - ) - right_boundary = .- left_boundary - upper_coef = SVector(T(3//4), T(-3//20), T(1//60)) - central_coef = zero(T) - lower_coef = -upper_coef - left_weights = SVector( T(13649//43200), - T(12013//8640), - T(2711//4320), - T(5359//4320), - T(7877//8640), - T(43801//43200) ) - right_weights = left_weights - left_boundary_derivatives = Tuple{}() - right_boundary_derivatives = left_boundary_derivatives + 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 == 6 + left_boundary = ( + # q1 + DerivativeCoefficientRow{T,1,6}(SVector(T(-21600//13649), + T(104009//54596), + T(30443//81894), + T(-33311//27298), + T(16863//27298), + T(-15025//163788) )), + # q2 + DerivativeCoefficientRow{T,1,6}(SVector(T(-104009//240260), + T(0), + T(-311//72078), + T(20229//24026), + T(-24337//48052), + T(36661//360390) )), + # q3 + DerivativeCoefficientRow{T,1,6}(SVector(T(-30443//162660), + T(311//32532), + T(0), + T(-11155//16266), + T(41287//32532), + T(-21999//54220) )), + # q4 #TODO + DerivativeCoefficientRow{T,1,7}(SVector(T(33311//107180), + T(-20229//21436), + T(485//1398), + T(0), + T(4147//21436), + T(25427//321540), + T(72//5359) )), + # q5 + DerivativeCoefficientRow{T,1,8}(SVector(T(-16863//78770), + T(24337//31508), + T(-41287//47262), + T(-4147//15754), + T(0), + T(342523//472620), + T(-1296//7877), + T(144//7877) )), + # q6 + DerivativeCoefficientRow{T,1,9}(SVector(T(15025//525612), + T(-36661//262806), + T(21999//87602), + T(-25427//262806), + T(-342523//525612), + T(0), + T(32400//43801), + T(-6480//43801), + T(720//43801) )) + ) + right_boundary = .- left_boundary + upper_coef = SVector(T(3//4), T(-3//20), T(1//60)) + central_coef = zero(T) + lower_coef = -upper_coef + left_weights = SVector( T(13649//43200), + T(12013//8640), + T(2711//4320), + T(5359//4320), + T(7877//8640), + T(43801//43200) ) + right_weights = left_weights + left_boundary_derivatives = Tuple{}() + right_boundary_derivatives = left_boundary_derivatives - 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 == 8 - left_boundary = ( - # q1 - DerivativeCoefficientRow{T,1,8}(SVector(T(-2540160//1498139), - T(5544277//5992556), - T(198794991//29962780), - T(-256916579//17977668), - T(20708767//1498139), - T(-41004357//5992556), - T(27390659//17977668), - T(-2323531//29962780) )), - # q2 - DerivativeCoefficientRow{T,1,8}(SVector(T(-5544277//31004596), - T(0), - T(-85002381//22146140), - T(49607267//4429228), - T(-165990199//13287684), - T(7655859//1107307), - T(-7568311//4429228), - T(48319961//465068940) )), - # q3 - DerivativeCoefficientRow{T,1,8}(SVector(T(-66264997//8719620), - T(9444709//415220), - T(0), - T(-20335981//249132), - T(32320879//249132), - T(-35518713//415220), - T(2502774//103805), - T(-3177073//1743924) )), - # q4 - DerivativeCoefficientRow{T,1,8}(SVector(T(256916579//109619916), - T(-49607267//5219996), - T(61007943//5219996), - T(0), - T(-68748371//5219996), - T(65088123//5219996), - T(-66558305//15659988), - T(3870214//9134993) )), - # q5 - DerivativeCoefficientRow{T,1,9}(SVector(T(-20708767//2096689), - T(165990199//3594324), - T(-96962637//1198108), - T(68748371//1198108), - T(0), - T(-27294549//1198108), - T(14054993//1198108), - T(-42678199//25160268), - T(-2592//299527) )), - # q6 - DerivativeCoefficientRow{T,1,10}(SVector(T(13668119//8660148), - T(-850651//103097), - T(35518713//2061940), - T(-21696041//1237164), - T(9098183//1237164), - T(0), - T(-231661//412388), - T(7120007//43300740), - T(3072//103097), - T(-288//103097) )), - # q7 - DerivativeCoefficientRow{T,1,11}(SVector(T(-27390659//56287644), - T(7568311//2680364), - T(-22524966//3350455), - T(66558305//8041092), - T(-14054993//2680364), - T(2084949//2680364), - T(0), - T(70710683//93812740), - T(-145152//670091), - T(27648//670091), - T(-2592//670091) )), - # q8 - DerivativeCoefficientRow{T,1,12}(SVector(T(2323531//102554780), - T(-48319961//307664340), - T(9531219//20510956), - T(-3870214//5127739), - T(2246221//3238572), - T(-21360021//102554780), - T(-70710683//102554780), - T(0), - T(4064256//5127739), - T(-1016064//5127739), - T(193536//5127739), - T(-18144//5127739) )), - ) - right_boundary = .- left_boundary - upper_coef = SVector(T(4//5), T(-1//5), T(4//105), T(-1//280)) - central_coef = zero(T) - lower_coef = -upper_coef - left_weights = SVector( T(1498139//5080320), - T(1107307//725760), - T(20761//80640), - T(1304999//725760), - T(299527//725760), - T(103097//80640), - T(670091//725760), - T(5127739//5080320) ) - right_weights = left_weights - left_boundary_derivatives = Tuple{}() - right_boundary_derivatives = left_boundary_derivatives + 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 == 8 + left_boundary = ( + # q1 + DerivativeCoefficientRow{T,1,8}(SVector(T(-2540160//1498139), + T(5544277//5992556), + T(198794991//29962780), + T(-256916579//17977668), + T(20708767//1498139), + T(-41004357//5992556), + T(27390659//17977668), + T(-2323531//29962780) )), + # q2 + DerivativeCoefficientRow{T,1,8}(SVector(T(-5544277//31004596), + T(0), + T(-85002381//22146140), + T(49607267//4429228), + T(-165990199//13287684), + T(7655859//1107307), + T(-7568311//4429228), + T(48319961//465068940) )), + # q3 + DerivativeCoefficientRow{T,1,8}(SVector(T(-66264997//8719620), + T(9444709//415220), + T(0), + T(-20335981//249132), + T(32320879//249132), + T(-35518713//415220), + T(2502774//103805), + T(-3177073//1743924) )), + # q4 + DerivativeCoefficientRow{T,1,8}(SVector(T(256916579//109619916), + T(-49607267//5219996), + T(61007943//5219996), + T(0), + T(-68748371//5219996), + T(65088123//5219996), + T(-66558305//15659988), + T(3870214//9134993) )), + # q5 + DerivativeCoefficientRow{T,1,9}(SVector(T(-20708767//2096689), + T(165990199//3594324), + T(-96962637//1198108), + T(68748371//1198108), + T(0), + T(-27294549//1198108), + T(14054993//1198108), + T(-42678199//25160268), + T(-2592//299527) )), + # q6 + DerivativeCoefficientRow{T,1,10}(SVector(T(13668119//8660148), + T(-850651//103097), + T(35518713//2061940), + T(-21696041//1237164), + T(9098183//1237164), + T(0), + T(-231661//412388), + T(7120007//43300740), + T(3072//103097), + T(-288//103097) )), + # q7 + DerivativeCoefficientRow{T,1,11}(SVector(T(-27390659//56287644), + T(7568311//2680364), + T(-22524966//3350455), + T(66558305//8041092), + T(-14054993//2680364), + T(2084949//2680364), + T(0), + T(70710683//93812740), + T(-145152//670091), + T(27648//670091), + T(-2592//670091) )), + # q8 + DerivativeCoefficientRow{T,1,12}(SVector(T(2323531//102554780), + T(-48319961//307664340), + T(9531219//20510956), + T(-3870214//5127739), + T(2246221//3238572), + T(-21360021//102554780), + T(-70710683//102554780), + T(0), + T(4064256//5127739), + T(-1016064//5127739), + T(193536//5127739), + T(-18144//5127739) )), + ) + right_boundary = .- left_boundary + upper_coef = SVector(T(4//5), T(-1//5), T(4//105), T(-1//280)) + central_coef = zero(T) + lower_coef = -upper_coef + left_weights = SVector( T(1498139//5080320), + T(1107307//725760), + T(20761//80640), + T(1304999//725760), + T(299527//725760), + T(103097//80640), + T(670091//725760), + T(5127739//5080320) ) + right_weights = left_weights + left_boundary_derivatives = Tuple{}() + right_boundary_derivatives = left_boundary_derivatives - 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) -else - throw(ArgumentError("Order $order not implemented/derived.")) -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) + else + throw(ArgumentError("Order $order not implemented/derived.")) + end end function second_derivative_coefficients(source::MattssonNordström2004, order::Int, T=Float64, parallel=Val{:serial}()) -if order == 2 - left_boundary = ( - # d1 - DerivativeCoefficientRow{T,1,3}(SVector(T(1), - T(-2), - T(1) )), - ) - right_boundary = left_boundary - upper_coef = SVector(T(1)) - central_coef = T(-2) - lower_coef = upper_coef - left_weights = SVector(T(1//2)) - right_weights = left_weights - left_boundary_derivatives = ( - DerivativeCoefficientRow{T,1,3}(SVector(T(-3//2), - T(2), - T(-1//2) )), - ) - right_boundary_derivatives = (-left_boundary_derivatives[1], ) + if order == 2 + left_boundary = ( + # d1 + DerivativeCoefficientRow{T,1,3}(SVector(T(1), + T(-2), + T(1) )), + ) + right_boundary = left_boundary + upper_coef = SVector(T(1)) + central_coef = T(-2) + lower_coef = upper_coef + left_weights = SVector(T(1//2)) + right_weights = left_weights + left_boundary_derivatives = ( + DerivativeCoefficientRow{T,1,3}(SVector(T(-3//2), + T(2), + T(-1//2) )), + ) + right_boundary_derivatives = (-left_boundary_derivatives[1], ) - DerivativeCoefficients(left_boundary, right_boundary, - left_boundary_derivatives, right_boundary_derivatives, - lower_coef, central_coef, upper_coef, - left_weights, right_weights, parallel, 2, order, source) -elseif order == 4 - left_boundary = ( - #TODO: - # d1 - DerivativeCoefficientRow{T,1,4}(SVector(T(2), - T(-5), - T(4), - T(-1) )), - # d2 - DerivativeCoefficientRow{T,1,3}(SVector(T(1), - T(-2), - T(1) )), - # d3 - DerivativeCoefficientRow{T,1,5}(SVector(T(-4//43), - T(59//43), - T(-110//43), - T(59//43), - T(-4//43) )), - # d4 - DerivativeCoefficientRow{T,1,6}(SVector(T(-1//49), - T(0), - T(59//49), - T(-118//49), - T(64//49), - T(-4//49) )) - ) - right_boundary = left_boundary - upper_coef = SVector(T(4//3), T(-1//12)) - central_coef = T(-5//2) - lower_coef = upper_coef - left_weights = SVector(T(17//48), T(59//48), T(43//48), T(49//48)) - right_weights = left_weights - left_boundary_derivatives = ( - DerivativeCoefficientRow{T,1,4}(SVector(T(-11/6), - T(3), - T(-3//2), - T(1//3) )), - ) - right_boundary_derivatives = (-left_boundary_derivatives[1], ) + DerivativeCoefficients(left_boundary, right_boundary, + left_boundary_derivatives, right_boundary_derivatives, + lower_coef, central_coef, upper_coef, + left_weights, right_weights, parallel, 2, order, source) + elseif order == 4 + left_boundary = ( + #TODO: + # d1 + DerivativeCoefficientRow{T,1,4}(SVector(T(2), + T(-5), + T(4), + T(-1) )), + # d2 + DerivativeCoefficientRow{T,1,3}(SVector(T(1), + T(-2), + T(1) )), + # d3 + DerivativeCoefficientRow{T,1,5}(SVector(T(-4//43), + T(59//43), + T(-110//43), + T(59//43), + T(-4//43) )), + # d4 + DerivativeCoefficientRow{T,1,6}(SVector(T(-1//49), + T(0), + T(59//49), + T(-118//49), + T(64//49), + T(-4//49) )) + ) + right_boundary = left_boundary + upper_coef = SVector(T(4//3), T(-1//12)) + central_coef = T(-5//2) + lower_coef = upper_coef + left_weights = SVector(T(17//48), T(59//48), T(43//48), T(49//48)) + right_weights = left_weights + left_boundary_derivatives = ( + DerivativeCoefficientRow{T,1,4}(SVector(T(-11/6), + T(3), + T(-3//2), + T(1//3) )), + ) + right_boundary_derivatives = (-left_boundary_derivatives[1], ) - DerivativeCoefficients(left_boundary, right_boundary, - left_boundary_derivatives, right_boundary_derivatives, - lower_coef, central_coef, upper_coef, - left_weights, right_weights, parallel, 2, order, source) -elseif order == 6 - left_boundary = ( - # d1 - DerivativeCoefficientRow{T,1,6}(SVector(T(114170//40947), - T(-438107//54596), - T(336409//40947), - T(-276997//81894), - T(3747//13649), - T(21035//163788) )), - # d2 - DerivativeCoefficientRow{T,1,6}(SVector(T(6173//5860), - T(-2066//879), - T(3283//1758), - T(-303//293), - T(2111//3516), - T(-601//4395) )), - # d3 - DerivativeCoefficientRow{T,1,6}(SVector(T(-52391//81330), - T(134603//32532), - T(-21982//2711), - T(112915//16266), - T(-46969//16266), - T(30409//54220) )), - # d4 - DerivativeCoefficientRow{T,1,7}(SVector(T(68603//321540), - T(-12423//10718), - T(112915//32154), - T(-75934//16077), - T(53369//21436), - T(-54899//160770), - T(48//5359) )), - # d5 - DerivativeCoefficientRow{T,1,8}(SVector(T(-7053//39385), - T(86551//94524), - T(-46969//23631), - T(53369//15754), - T(-87904//23631), - T(820271//472620), - T(-1296//7877), - T(96//7877) )), - # d6 - DerivativeCoefficientRow{T,1,9}(SVector(T(21035//525612), - T(-24641//131403), - T(30409//87602), - T(-54899//131403), - T(820271//525612), - T(-117600//43801), - T(64800//43801), - T(-6480//43801), - T(480//43801) )), - ) - right_boundary = left_boundary - upper_coef = SVector(T(3//2), T(-3//20), T(1//90)) - central_coef = T(-49//18) - lower_coef = upper_coef - left_weights = SVector( T(13649//43200), - T(12013//8640), - T(2711//4320), - T(5359//4320), - T(7877//8640), - T(43801//43200) ) - right_weights = left_weights - left_boundary_derivatives = ( - DerivativeCoefficientRow{T,1,5}(SVector(T(-25//12), - T(4), - T(-3), - T(4//3), - T(-1//4) )), - ) - right_boundary_derivatives = (-left_boundary_derivatives[1], ) + DerivativeCoefficients(left_boundary, right_boundary, + left_boundary_derivatives, right_boundary_derivatives, + lower_coef, central_coef, upper_coef, + left_weights, right_weights, parallel, 2, order, source) + elseif order == 6 + left_boundary = ( + # d1 + DerivativeCoefficientRow{T,1,6}(SVector(T(114170//40947), + T(-438107//54596), + T(336409//40947), + T(-276997//81894), + T(3747//13649), + T(21035//163788) )), + # d2 + DerivativeCoefficientRow{T,1,6}(SVector(T(6173//5860), + T(-2066//879), + T(3283//1758), + T(-303//293), + T(2111//3516), + T(-601//4395) )), + # d3 + DerivativeCoefficientRow{T,1,6}(SVector(T(-52391//81330), + T(134603//32532), + T(-21982//2711), + T(112915//16266), + T(-46969//16266), + T(30409//54220) )), + # d4 + DerivativeCoefficientRow{T,1,7}(SVector(T(68603//321540), + T(-12423//10718), + T(112915//32154), + T(-75934//16077), + T(53369//21436), + T(-54899//160770), + T(48//5359) )), + # d5 + DerivativeCoefficientRow{T,1,8}(SVector(T(-7053//39385), + T(86551//94524), + T(-46969//23631), + T(53369//15754), + T(-87904//23631), + T(820271//472620), + T(-1296//7877), + T(96//7877) )), + # d6 + DerivativeCoefficientRow{T,1,9}(SVector(T(21035//525612), + T(-24641//131403), + T(30409//87602), + T(-54899//131403), + T(820271//525612), + T(-117600//43801), + T(64800//43801), + T(-6480//43801), + T(480//43801) )), + ) + right_boundary = left_boundary + upper_coef = SVector(T(3//2), T(-3//20), T(1//90)) + central_coef = T(-49//18) + lower_coef = upper_coef + left_weights = SVector( T(13649//43200), + T(12013//8640), + T(2711//4320), + T(5359//4320), + T(7877//8640), + T(43801//43200) ) + right_weights = left_weights + left_boundary_derivatives = ( + DerivativeCoefficientRow{T,1,5}(SVector(T(-25//12), + T(4), + T(-3), + T(4//3), + T(-1//4) )), + ) + right_boundary_derivatives = (-left_boundary_derivatives[1], ) - DerivativeCoefficients(left_boundary, right_boundary, - left_boundary_derivatives, right_boundary_derivatives, - lower_coef, central_coef, upper_coef, - left_weights, right_weights, parallel, 2, order, source) -elseif order == 8 - left_boundary = ( - # d1 - DerivativeCoefficientRow{T,1,7}(SVector(T(4870382994799//1358976868290), - T(-893640087518//75498714905), - T(926594825119//60398971924), - T(-1315109406200//135897686829), - T(39126983272//15099742981), - T(12344491342//75498714905), - T(-451560522577//2717953736580) )), - # d2 - DerivativeCoefficientRow{T,1,8}(SVector(T(333806012194//390619153855), - T(-154646272029//111605472530), - T(1168338040//33481641759), - T(82699112501//133926567036), - T(-171562838//11160547253), - T(-28244698346//167408208795), - T(11904122576//167408208795), - T(-2598164715//312495323084) )), - # d3 - DerivativeCoefficientRow{T,1,8}(SVector(T(7838984095//52731029988), - T(1168338040//5649753213), - T(-88747895//144865467), - T(423587231//627750357), - T(-43205598281//22599012852), - T(4876378562//1883251071), - T(-5124426509//3766502142), - T(10496900965//39548272491) )), - # d4 - DerivativeCoefficientRow{T,1,8}(SVector(T(-94978241528//828644350023), - T(82699112501//157837019052), - T(1270761693//13153084921), - T(-167389605005//118377764289), - T(48242560214//39459254763), - T(-31673996013//52612339684), - T(43556319241//118377764289), - T(-44430275135//552429566682) )), - # d5 - DerivativeCoefficientRow{T,1,9}(SVector(T(1455067816//21132528431), - T(-171562838//3018932633), - T(-43205598281//36227191596), - T(48242560214//9056797899), - T(-52276055645//6037865266), - T(57521587238//9056797899), - T(-80321706377//36227191596), - T(8078087158//21132528431), - T(-1296//299527) )), - # d6 - DerivativeCoefficientRow{T,1,10}(SVector(T(10881504334//327321118845), - T(-28244698346//140280479505), - T(4876378562//9352031967), - T(-10557998671//12469375956), - T(57521587238//28056095901), - T(-278531401019//93520319670), - T(73790130002//46760159835), - T(-137529995233//785570685228), - T(2048//103097), - T(-144//103097) )), - # d7 - DerivativeCoefficientRow{T,1,11}(SVector(T(-135555328849//8509847458140), - T(11904122576//101307707835), - T(-5124426509//13507694378), - T(43556319241//60784624701), - T(-80321706377//81046166268), - T(73790130002//33769235945), - T(-950494905688//303923123505), - T(239073018673//141830790969), - T(-145152//670091), - T(18432//670091), - T(-1296//670091) )), - # d8 - DerivativeCoefficientRow{T,1,12}(SVector(T(0), - T(-2598164715//206729925524), - T(10496900965//155047444143), - T(-44430275135//310094888286), - T(425162482//2720130599), - T(-137529995233//620189776572), - T(239073018673//155047444143), - T(-144648000000//51682481381), - T(8128512//5127739), - T(-1016064//5127739), - T(129024//5127739), - T(-9072//5127739) )), + DerivativeCoefficients(left_boundary, right_boundary, + left_boundary_derivatives, right_boundary_derivatives, + lower_coef, central_coef, upper_coef, + left_weights, right_weights, parallel, 2, order, source) + elseif order == 8 + left_boundary = ( + # d1 + DerivativeCoefficientRow{T,1,7}(SVector(T(4870382994799//1358976868290), + T(-893640087518//75498714905), + T(926594825119//60398971924), + T(-1315109406200//135897686829), + T(39126983272//15099742981), + T(12344491342//75498714905), + T(-451560522577//2717953736580) )), + # d2 + DerivativeCoefficientRow{T,1,8}(SVector(T(333806012194//390619153855), + T(-154646272029//111605472530), + T(1168338040//33481641759), + T(82699112501//133926567036), + T(-171562838//11160547253), + T(-28244698346//167408208795), + T(11904122576//167408208795), + T(-2598164715//312495323084) )), + # d3 + DerivativeCoefficientRow{T,1,8}(SVector(T(7838984095//52731029988), + T(1168338040//5649753213), + T(-88747895//144865467), + T(423587231//627750357), + T(-43205598281//22599012852), + T(4876378562//1883251071), + T(-5124426509//3766502142), + T(10496900965//39548272491) )), + # d4 + DerivativeCoefficientRow{T,1,8}(SVector(T(-94978241528//828644350023), + T(82699112501//157837019052), + T(1270761693//13153084921), + T(-167389605005//118377764289), + T(48242560214//39459254763), + T(-31673996013//52612339684), + T(43556319241//118377764289), + T(-44430275135//552429566682) )), + # d5 + DerivativeCoefficientRow{T,1,9}(SVector(T(1455067816//21132528431), + T(-171562838//3018932633), + T(-43205598281//36227191596), + T(48242560214//9056797899), + T(-52276055645//6037865266), + T(57521587238//9056797899), + T(-80321706377//36227191596), + T(8078087158//21132528431), + T(-1296//299527) )), + # d6 + DerivativeCoefficientRow{T,1,10}(SVector(T(10881504334//327321118845), + T(-28244698346//140280479505), + T(4876378562//9352031967), + T(-10557998671//12469375956), + T(57521587238//28056095901), + T(-278531401019//93520319670), + T(73790130002//46760159835), + T(-137529995233//785570685228), + T(2048//103097), + T(-144//103097) )), + # d7 + DerivativeCoefficientRow{T,1,11}(SVector(T(-135555328849//8509847458140), + T(11904122576//101307707835), + T(-5124426509//13507694378), + T(43556319241//60784624701), + T(-80321706377//81046166268), + T(73790130002//33769235945), + T(-950494905688//303923123505), + T(239073018673//141830790969), + T(-145152//670091), + T(18432//670091), + T(-1296//670091) )), + # d8 + DerivativeCoefficientRow{T,1,12}(SVector(T(0), + T(-2598164715//206729925524), + T(10496900965//155047444143), + T(-44430275135//310094888286), + T(425162482//2720130599), + T(-137529995233//620189776572), + T(239073018673//155047444143), + T(-144648000000//51682481381), + T(8128512//5127739), + T(-1016064//5127739), + T(129024//5127739), + T(-9072//5127739) )), - ) - right_boundary = left_boundary - upper_coef = SVector(T(8//5), T(-1//5), T(8//315), T(-1//560)) - central_coef = T(-205//72) - lower_coef = upper_coef - left_weights = SVector( T(1498139//5080320), - T(1107307//725760), - T(20761//80640), - T(1304999//725760), - T(299527//725760), - T(103097//80640), - T(670091//725760), - T(5127739//5080320) ) - right_weights = left_weights - left_boundary_derivatives = ( - DerivativeCoefficientRow{T,1,7}(SVector(T(-4723//2100), - T(839//175), - T(-157//35), - T(278//105), - T(-103//140), - T(-1//175), - T(6//175) )), - ) - right_boundary_derivatives = (-left_boundary_derivatives[1], ) + ) + right_boundary = left_boundary + upper_coef = SVector(T(8//5), T(-1//5), T(8//315), T(-1//560)) + central_coef = T(-205//72) + lower_coef = upper_coef + left_weights = SVector( T(1498139//5080320), + T(1107307//725760), + T(20761//80640), + T(1304999//725760), + T(299527//725760), + T(103097//80640), + T(670091//725760), + T(5127739//5080320) ) + right_weights = left_weights + left_boundary_derivatives = ( + DerivativeCoefficientRow{T,1,7}(SVector(T(-4723//2100), + T(839//175), + T(-157//35), + T(278//105), + T(-103//140), + T(-1//175), + T(6//175) )), + ) + right_boundary_derivatives = (-left_boundary_derivatives[1], ) - DerivativeCoefficients(left_boundary, right_boundary, - left_boundary_derivatives, right_boundary_derivatives, - lower_coef, central_coef, upper_coef, - left_weights, right_weights, parallel, 2, order, source) + DerivativeCoefficients(left_boundary, right_boundary, + left_boundary_derivatives, right_boundary_derivatives, + lower_coef, central_coef, upper_coef, + left_weights, right_weights, parallel, 2, order, source) else throw(ArgumentError("Order $order not implemented/derived.")) end From b37b19ee266d8c674338631c2584c17952a02eae Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Wed, 25 Mar 2020 09:16:34 +0300 Subject: [PATCH 2/8] add checks for LowerOffset > 0 for non-periodic operators --- src/SBP_operators.jl | 24 ++++++++++++++++-------- 1 file changed, 16 insertions(+), 8 deletions(-) diff --git a/src/SBP_operators.jl b/src/SBP_operators.jl index da8cef51..8133827a 100644 --- a/src/SBP_operators.jl +++ b/src/SBP_operators.jl @@ -175,11 +175,15 @@ end @generated function convolve_interior_coefficients!(dest::AbstractVector, lower_coef::SVector{LowerOffset}, central_coef, upper_coef::SVector{UpperOffset}, u::AbstractVector, α, β, left_boundary_width, right_boundary_width, parallel) where {LowerOffset, UpperOffset} - ex = :( lower_coef[$LowerOffset]*u[i-$LowerOffset] ) - for j in LowerOffset-1:-1:1 - ex = :( $ex + lower_coef[$j]*u[i-$j] ) + if LowerOffset > 0 + ex = :( lower_coef[$LowerOffset]*u[i-$LowerOffset] ) + for j in LowerOffset-1:-1:1 + ex = :( $ex + lower_coef[$j]*u[i-$j] ) + end + ex = :( $ex + central_coef*u[i] ) + else + ex = :( central_coef*u[i] ) end - ex = :( $ex + central_coef*u[i] ) for j in 1:UpperOffset ex = :( $ex + upper_coef[$j]*u[i+$j] ) end @@ -206,11 +210,15 @@ function convolve_interior_coefficients!(dest::AbstractVector{T}, lower_coef::SV end @generated function convolve_interior_coefficients!(dest::AbstractVector, lower_coef::SVector{LowerOffset}, central_coef, upper_coef::SVector{UpperOffset}, u::AbstractVector, α, left_boundary_width, right_boundary_width, parallel) where {LowerOffset, UpperOffset} - ex = :( lower_coef[$LowerOffset]*u[i-$LowerOffset] ) - for j in LowerOffset-1:-1:1 - ex = :( $ex + lower_coef[$j]*u[i-$j] ) + if LowerOffset > 0 + ex = :( lower_coef[$LowerOffset]*u[i-$LowerOffset] ) + for j in LowerOffset-1:-1:1 + ex = :( $ex + lower_coef[$j]*u[i-$j] ) + end + ex = :( $ex + central_coef*u[i] ) + else + ex = :( central_coef*u[i] ) end - ex = :( $ex + central_coef*u[i] ) for j in 1:UpperOffset ex = :( $ex + upper_coef[$j]*u[i+$j] ) end From cce07656f0b3835e64f31698bf6ef49b8038818f Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Wed, 25 Mar 2020 10:12:49 +0300 Subject: [PATCH 3/8] upwind operators of Mattsson(2017), order 2 --- dev/Mattsson2017.jl | 18 ++++ src/SBP_coefficients/Mattsson2017.jl | 139 +++++++++++++++++++++++++++ src/SBP_operators.jl | 34 ++++++- src/SummationByPartsOperators.jl | 6 +- src/general_operators.jl | 2 +- test/periodic_operators_test.jl | 2 +- test/runtests.jl | 1 + test/upwind_operators.jl | 30 ++++++ 8 files changed, 226 insertions(+), 6 deletions(-) create mode 100644 dev/Mattsson2017.jl create mode 100644 src/SBP_coefficients/Mattsson2017.jl create mode 100644 test/upwind_operators.jl diff --git a/dev/Mattsson2017.jl b/dev/Mattsson2017.jl new file mode 100644 index 00000000..5c151f2e --- /dev/null +++ b/dev/Mattsson2017.jl @@ -0,0 +1,18 @@ +using LinearAlgebra + +# order 2 +Qp = [ + -1//4 5//4 -1//2 0 0 0 + -1//4 -5//4 2 -1//2 0 0 + 0 0 -3//2 2 -1//2 0 + 0 0 0 -3//2 2 -1//2 + 0 0 0 0 -5//4 5//4 + 0 0 0 0 -1//4 -1//4 +] +Qm = Matrix(-Qp') +B = zero(Qp); B[1,1] = -1; B[end,end] = 1 +M = Diagonal([1//4, 5//4, 1, 1, 5//4, 1//4]) +Dp = M \ (Qp + B / 2) +Dm = M \ (Qm + B / 2) + +# order 3 diff --git a/src/SBP_coefficients/Mattsson2017.jl b/src/SBP_coefficients/Mattsson2017.jl new file mode 100644 index 00000000..0b64eba5 --- /dev/null +++ b/src/SBP_coefficients/Mattsson2017.jl @@ -0,0 +1,139 @@ + +""" + Mattsson2017 + +Coefficients of the upwind SBP operators given in + Mattsson (2017) + Diagonal-norm upwind SBP operators. + Journal of Computational Physics 335, pp.283-310. +""" +struct Mattsson2017 <: SourceOfCoefficients + kind::Symbol + + function Mattsson2017(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::Mattsson2017) +print(io, + " Upwind coefficients (", source.kind, ") of \n", + " Mattsson (2017) \n", + " Diagonal-norm upwind SBP operators. \n", + " Journal of Computational Physics 335, pp.283-310. \n") +end + + +function first_derivative_coefficients(source::Mattsson2017, order::Int, T=Float64, parallel=Val{:serial}()) + if order == 2 + left_boundary_plus = ( + DerivativeCoefficientRow{T,1,3}(SVector(T(-3), + T(5), + T(-2), )), + DerivativeCoefficientRow{T,1,4}(SVector(T(-1//5), + T(-1), + T(8//5), + T(-2//5), )), + ) + right_boundary_plus = ( + DerivativeCoefficientRow{T,1,2}(SVector(T(1), + T(-1), )), + DerivativeCoefficientRow{T,1,2}(SVector(T(1), + T(-1), )), + ) + upper_coef_plus = SVector( T(2), + T(-1//2), ) + central_coef_plus = T(-3//2) + lower_coef_plus = SVector{0,T}() + left_weights = SVector( T(1//4), + T(5//4), ) + right_weights = left_weights + left_boundary_derivatives = Tuple{}() + right_boundary_derivatives = left_boundary_derivatives + + left_boundary_minus = ( + DerivativeCoefficientRow{T,1,2}(SVector(T(-1), + T(1), )), + DerivativeCoefficientRow{T,1,2}(SVector(T(-1), + T(1), )), + ) + right_boundary_minus = ( + DerivativeCoefficientRow{T,1,3}(SVector(T(3), + T(-5), + T(2), )), + DerivativeCoefficientRow{T,1,4}(SVector(T(1//5), + T(1), + T(-8//5), + T(2//5), )), + ) + 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, parallel, 1, order, source) + else + throw(ArgumentError("Order $order not implemented/derived.")) + end +end + + +# function second_derivative_coefficients(source::Mattsson2017, order::Int, T=Float64, parallel=Val{:serial}()) +# if order == 2 +# left_boundary = ( +# # d1 +# DerivativeCoefficientRow{T,1,3}(SVector(T(1), +# T(-2), +# T(1) )), +# ) +# right_boundary = left_boundary +# upper_coef = SVector(T(1)) +# central_coef = T(-2) +# lower_coef = upper_coef +# left_weights = SVector(T(1//2)) +# right_weights = left_weights +# left_boundary_derivatives = ( +# DerivativeCoefficientRow{T,1,3}(SVector(T(-3//2), +# T(2), +# T(-1//2) )), +# ) +# right_boundary_derivatives = (-left_boundary_derivatives[1], ) + +# DerivativeCoefficients(left_boundary, right_boundary, +# left_boundary_derivatives, right_boundary_derivatives, +# lower_coef, central_coef, upper_coef, +# left_weights, right_weights, parallel, 2, order, source) +# else +# throw(ArgumentError("Order $order not implemented/derived.")) +# end +# end diff --git a/src/SBP_operators.jl b/src/SBP_operators.jl index 8133827a..7ae8903a 100644 --- a/src/SBP_operators.jl +++ b/src/SBP_operators.jl @@ -130,10 +130,42 @@ function offset(::DerivativeCoefficientRow{T,Start,Length}) where {T,Start,Lengt Start end -function -(coef_row::DerivativeCoefficientRow{T,Start,Length}) where {T,Start,Length} +function Base.:-(coef_row::DerivativeCoefficientRow{T,Start,Length}) where {T,Start,Length} DerivativeCoefficientRow{T,Start,Length}(-coef_row.coef) end +function widening_plus(row1::SVector{Length1,T}, row2::SVector{Length2,T}) where {T,Length1,Length2} + Length = max(Length1, Length2) + row = SVector{Length,T}(ntuple(i -> zero(T), Length)) + for i in 1:Length1 + row = Base.setindex(row, row[i] + row1[i], i) + end + for i in 1:Length2 + row = Base.setindex(row, row[i] + row2[i], i) + end + row +end + +function Base.:+(coef_row1::DerivativeCoefficientRow{T,Start1,Length1}, coef_row2::DerivativeCoefficientRow{T,Start2,Length2}) where {T,Start1,Length1,Start2,Length2} + End = max(Start1+Length1, Start2+Length2) + Start = min(Start1, Start2) + Length = End - Start + row = SVector{Length,T}(ntuple(i -> zero(T), Length)) + for i in 1:Length1 + j = i-Start1+Start + row = Base.setindex(row, row[j] + coef_row1.coef[i], j) + end + for i in 1:Length2 + j = i-Start2+Start + row = Base.setindex(row, row[j] + coef_row2.coef[i], j) + end + DerivativeCoefficientRow{T,Start,Length}(row) +end + +function Base.:/(coef_row::DerivativeCoefficientRow{T,Start,Length}, α::Number) where {T,Start,Length} + DerivativeCoefficientRow{T,Start,Length}(coef_row.coef / α) +end + @inline function convolve_left_row(coef_row::DerivativeCoefficientRow{T,Start,Length}, u) where {T,Start,Length} @inbounds tmp = coef_row.coef[1]*u[Start] @inbounds for i in 2:Length diff --git a/src/SummationByPartsOperators.jl b/src/SummationByPartsOperators.jl index 5fa1dece..e9635861 100644 --- a/src/SummationByPartsOperators.jl +++ b/src/SummationByPartsOperators.jl @@ -14,9 +14,7 @@ using StaticArrays @reexport using DiffEqBase using DiffEqCallbacks -import Base: *, - import LinearAlgebra: mul! -export mul! @reexport using PolynomialBases import PolynomialBases: integrate, evaluate_coefficients, evaluate_coefficients!, compute_coefficients, compute_coefficients! @@ -59,6 +57,7 @@ include("SBP_coefficients/Mattsson2012.jl") include("SBP_coefficients/Mattsson2014.jl") include("SBP_coefficients/MattssonAlmquistCarpenter2014Extended.jl") include("SBP_coefficients/MattssonAlmquistCarpenter2014Optimal.jl") +include("SBP_coefficients/Mattsson2017.jl") include("conservation_laws/general_laws.jl") include("conservation_laws/burgers.jl") @@ -91,7 +90,8 @@ export periodic_central_derivative_operator, periodic_derivative_operator, deriv export Fornberg1998, Holoborodko2008, BeljaddLeFlochMishraParés2017 export MattssonNordström2004, MattssonSvärdNordström2004, MattssonSvärdShoeybi2008, Mattsson2012, Mattsson2014, - MattssonAlmquistCarpenter2014Extended, MattssonAlmquistCarpenter2014Optimal + MattssonAlmquistCarpenter2014Extended, MattssonAlmquistCarpenter2014Optimal, + Mattsson2017 export Tadmor1989, MadayTadmor1989, Tadmor1993, TadmorWaagan2012Standard, TadmorWaagan2012Convergent diff --git a/src/general_operators.jl b/src/general_operators.jl index 36ce608b..df9a613a 100644 --- a/src/general_operators.jl +++ b/src/general_operators.jl @@ -31,7 +31,7 @@ Base.@propagate_inbounds function mul!(dest, D::AbstractDerivativeOperator, u) mul!(dest, D, u, one(eltype(dest))) end -@noinline function *(D::AbstractDerivativeOperator, u) +@noinline function Base.:*(D::AbstractDerivativeOperator, u) @boundscheck begin @argcheck size(D,1) == size(D,2) == length(u) DimensionMismatch end diff --git a/test/periodic_operators_test.jl b/test/periodic_operators_test.jl index 4c813958..64dfa711 100644 --- a/test/periodic_operators_test.jl +++ b/test/periodic_operators_test.jl @@ -880,7 +880,7 @@ let N = 5 D = periodic_derivative_operator(1, 2, 0//1, (N-1)//1, N, -2) @test Matrix(D) ≈ [ 3//2 0//1 1//2 -2//1 - -2//1 3//2 0//1 1//2 + -2//1 3//2 0//1 1//2 1//2 -2//1 3//2 0//1 0//1 1//2 -2//1 3//2] end diff --git a/test/runtests.jl b/test/runtests.jl index 9356dcff..84c4cdcc 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -8,6 +8,7 @@ using Test @time @testset "Fourier Operators" begin include("fourier_operators_test.jl") end @time @testset "Legendre Operators" begin include("legendre_operators_test.jl") end @time @testset "Sum of Operators" begin include("sum_of_operators_test.jl") end +@time @testset "Upwind Operators" begin include("upwind_operators.jl") end @time @testset "Conservation Laws" begin include("conservation_laws/burgers_test.jl") include("conservation_laws/cubic_test.jl") diff --git a/test/upwind_operators.jl b/test/upwind_operators.jl new file mode 100644 index 00000000..dfb59348 --- /dev/null +++ b/test/upwind_operators.jl @@ -0,0 +1,30 @@ +using Test +using LinearAlgebra +using SummationByPartsOperators + +# check construction of interior part of upwind operators +@testset "Check interior parts" begin + N = 21 + xmin = 0//1 + xmax = (N + 1) // 1 + interior = 10:N-10 + + acc_order = 2 + 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) + println(devnull, Dp_bounded) + 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, -(acc_order - 1) ÷ 2) + Dm_periodic = periodic_derivative_operator(1, acc_order, xmin, xmax, N, -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) < 10N * eps() +end From a5f004199c62d75a318f9bbf7f927819f4aa1e66 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Wed, 25 Mar 2020 10:39:40 +0300 Subject: [PATCH 4/8] upwind operators of Mattsson(2017), order 3 --- dev/Mattsson2017.jl | 28 ++++++++++ src/SBP_coefficients/Mattsson2017.jl | 78 +++++++++++++++++++++++++++- test/upwind_operators.jl | 48 ++++++++++------- 3 files changed, 135 insertions(+), 19 deletions(-) diff --git a/dev/Mattsson2017.jl b/dev/Mattsson2017.jl index 5c151f2e..a4e51292 100644 --- a/dev/Mattsson2017.jl +++ b/dev/Mattsson2017.jl @@ -16,3 +16,31 @@ Dp = M \ (Qp + B / 2) Dm = M \ (Qm + B / 2) # order 3 +Qp = [ + -1//12 3//4 -1//6 0 0 0 + -5//12 -5//12 1 -1//6 0 0 + 0 -1//3 -1//2 1 -1//6 0 + 0 0 -1//3 -1//2 1 -1//6 + 0 0 0 -1//3 -5//12 3//4 + 0 0 0 0 -5//12 -1//12 +] +Qm = Matrix(-Qp') +B = zero(Qp); B[1,1] = -1; B[end,end] = 1 +M = Diagonal([5//12, 13//12, 1, 1, 13//12, 5//12]) +Dp = M \ (Qp + B / 2) +Dm = M \ (Qm + B / 2) + +# order 4 +Qp = [ + -1//12 3//4 -1//6 0 0 0 + -5//12 -5//12 1 -1//6 0 0 + 0 -1//3 -1//2 1 -1//6 0 + 0 0 -1//3 -1//2 1 -1//6 + 0 0 0 -1//3 -5//12 3//4 + 0 0 0 0 -5//12 -1//12 +] +Qm = Matrix(-Qp') +B = zero(Qp); B[1,1] = -1; B[end,end] = 1 +M = Diagonal([5//12, 13//12, 1, 1, 13//12, 5//12]) +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 0b64eba5..9113582e 100644 --- a/src/SBP_coefficients/Mattsson2017.jl +++ b/src/SBP_coefficients/Mattsson2017.jl @@ -102,7 +102,83 @@ function first_derivative_coefficients(source::Mattsson2017, order::Int, T=Float left_boundary_derivatives, right_boundary_derivatives, lower_coef, central_coef, upper_coef, left_weights, right_weights, parallel, 1, order, source) - else + elseif order == 3 + left_boundary_plus = ( + DerivativeCoefficientRow{T,1,3}(SVector(T(-7//5), + T(9//5), + T(-2//5), )), + DerivativeCoefficientRow{T,1,4}(SVector(T(-5//13), + T(-5//13), + T(12//13), + T(-2//13), )), + ) + right_boundary_plus = ( + DerivativeCoefficientRow{T,1,2}(SVector(T(1), + T(-1), )), + DerivativeCoefficientRow{T,1,3}(SVector(T(9//13), + T(-5//13), + T(-4//13), )), + ) + upper_coef_plus = SVector( T(1), + T(-1//6), ) + central_coef_plus = T(-1//2) + lower_coef_plus = SVector( T(-1//3), ) + left_weights = SVector( T(5//12), + T(13//12), ) + right_weights = left_weights + left_boundary_derivatives = Tuple{}() + right_boundary_derivatives = left_boundary_derivatives + + left_boundary_minus = ( + DerivativeCoefficientRow{T,1,2}(SVector(T(-1), + T(1), )), + DerivativeCoefficientRow{T,1,3}(SVector(T(-9//13), + T(5//13), + T(4//13), )), + ) + right_boundary_minus = ( + DerivativeCoefficientRow{T,1,3}(SVector(T(7//5), + T(-9//5), + T(2//5), )), + DerivativeCoefficientRow{T,1,4}(SVector(T(5//13), + T(5//13), + T(-12//13), + T(2//13), )), + ) + 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, parallel, 1, order, source) + else throw(ArgumentError("Order $order not implemented/derived.")) end end diff --git a/test/upwind_operators.jl b/test/upwind_operators.jl index dfb59348..46b7fb9e 100644 --- a/test/upwind_operators.jl +++ b/test/upwind_operators.jl @@ -9,22 +9,34 @@ using SummationByPartsOperators xmax = (N + 1) // 1 interior = 10:N-10 - acc_order = 2 - 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) - println(devnull, Dp_bounded) - 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, -(acc_order - 1) ÷ 2) - Dm_periodic = periodic_derivative_operator(1, acc_order, xmin, xmax, N, -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) < 10N * eps() + + for acc_order in 2:3 + 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) + println(devnull, Dp_bounded) + 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, -(acc_order - 1) ÷ 2) + Dm_periodic = periodic_derivative_operator(1, acc_order, xmin, xmax, N, -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) < 10N * eps() + x = grid(Dp_bounded) + for D in (Dp_bounded, Dm_bounded, Dc_bounded) + @test norm(D * x.^0) < 10N * 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 + end end From e2b7247405cdba09a60b2539c0f6a3f6f7ba418f Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Wed, 25 Mar 2020 11:18:20 +0300 Subject: [PATCH 5/8] upwind operators of Mattsson(2017), order 4 --- dev/Mattsson2017.jl | 20 +- src/SBP_coefficients/Mattsson2017.jl | 266 +++++++++++++++++++++++++++ test/upwind_operators.jl | 2 +- 3 files changed, 280 insertions(+), 8 deletions(-) diff --git a/dev/Mattsson2017.jl b/dev/Mattsson2017.jl index a4e51292..36b2ad49 100644 --- a/dev/Mattsson2017.jl +++ b/dev/Mattsson2017.jl @@ -32,15 +32,21 @@ Dm = M \ (Qm + B / 2) # order 4 Qp = [ - -1//12 3//4 -1//6 0 0 0 - -5//12 -5//12 1 -1//6 0 0 - 0 -1//3 -1//2 1 -1//6 0 - 0 0 -1//3 -1//2 1 -1//6 - 0 0 0 -1//3 -5//12 3//4 - 0 0 0 0 -5//12 -1//12 + -1//48 205//288 -29//144 1//96 0 0 0 0 0 0 0 + -169//288 -11//48 33//32 -43//144 1//12 0 0 0 0 0 0 + 11//144 -13//32 -29//48 389//288 -1//2 1//12 0 0 0 0 0 + 1//32 -11//144 -65//288 -13//16 3//2 -1//2 1//12 0 0 0 0 + 0 0 0 -1//4 -5//6 3//2 -1//2 1//12 0 0 0 + 0 0 0 0 -1//4 -5//6 3//2 -1//2 1//12 0 0 + 0 0 0 0 0 -1//4 -5//6 3//2 -1//2 1//12 0 + 0 0 0 0 0 0 -1//4 -13//16 389//288 -43//144 1//96 + 0 0 0 0 0 0 0 -65//288 -29//48 33//32 -29//144 + 0 0 0 0 0 0 0 -11//144 -13//32 -11//48 205//288 + 0 0 0 0 0 0 0 1//32 11//144 -169//288 -1//48 ] Qm = Matrix(-Qp') B = zero(Qp); B[1,1] = -1; B[end,end] = 1 -M = Diagonal([5//12, 13//12, 1, 1, 13//12, 5//12]) +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) diff --git a/src/SBP_coefficients/Mattsson2017.jl b/src/SBP_coefficients/Mattsson2017.jl index 9113582e..77d73550 100644 --- a/src/SBP_coefficients/Mattsson2017.jl +++ b/src/SBP_coefficients/Mattsson2017.jl @@ -178,6 +178,272 @@ function first_derivative_coefficients(source::Mattsson2017, order::Int, T=Float left_boundary_derivatives, right_boundary_derivatives, lower_coef, central_coef, upper_coef, left_weights, right_weights, parallel, 1, order, source) + elseif order == 4 + left_boundary_plus = ( + DerivativeCoefficientRow{T,1,4}(SVector(T(-75//49), + T(205//98), + T(-29//49), + T(3//98), )), + DerivativeCoefficientRow{T,1,5}(SVector(T(-169//366), + T(-11//61), + T(99//122), + T(-43//183), + T(4//61), )), + DerivativeCoefficientRow{T,1,6}(SVector(T(11//123), + T(-39//82), + T(-29//41), + T(389//246), + T(-24//41), + T(4//41), )), + DerivativeCoefficientRow{T,1,7}(SVector(T(9//298), + T(-11//149), + T(-65//298), + T(-117//149), + T(216//149), + T(-72//149), + T(12//149), )), + ) + right_boundary_plus = ( + DerivativeCoefficientRow{T,1,4}(SVector(T(69//49), + T(-169//98), + T(11//49), + T(9//98), )), + DerivativeCoefficientRow{T,1,4}(SVector(T(205//366), + T(-11//61), + T(-39//122), + T(-11//183), )), + DerivativeCoefficientRow{T,1,4}(SVector(T(-29//123), + T(99//82), + T(-29//41), + T(-65//246), )), + DerivativeCoefficientRow{T,1,5}(SVector(T(3//298), + T(-43//149), + T(389//298), + T(-117//149), + T(-36//149), )), + ) + upper_coef_plus = SVector( T(3//2), + T(-1//2), + T(1//12), ) + central_coef_plus = T(-5//6) + lower_coef_plus = SVector( T(-1//4), ) + left_weights = SVector( T(49//144), + T(61//48), + T(41//48), + T(149//144), ) + right_weights = left_weights + left_boundary_derivatives = Tuple{}() + right_boundary_derivatives = left_boundary_derivatives + + left_boundary_minus = ( + DerivativeCoefficientRow{T,1,4}(SVector(T(-69//49), + T(169//98), + T(-11//49), + T(-9//98), )), + DerivativeCoefficientRow{T,1,4}(SVector(T(-205//366), + T(11//61), + T(39//122), + T(11//183), )), + DerivativeCoefficientRow{T,1,4}(SVector(T(29//123), + T(-99//82), + T(29//41), + T(65//246), )), + DerivativeCoefficientRow{T,1,5}(SVector(T(-3//298), + T(43//149), + T(-389//298), + T(117//149), + T(36//149), )), + ) + right_boundary_minus = ( + DerivativeCoefficientRow{T,1,4}(SVector(T(75//49), + T(-205//98), + T(29//49), + T(-3//98), )), + DerivativeCoefficientRow{T,1,5}(SVector(T(169//366), + T(11//61), + T(-99//122), + T(43//183), + T(-4//61), )), + DerivativeCoefficientRow{T,1,6}(SVector(T(-11//123), + T(39//82), + T(29//41), + T(-389//246), + T(24//41), + T(-4//41), )), + DerivativeCoefficientRow{T,1,7}(SVector(T(-9//298), + T(11//149), + T(65//298), + T(117//149), + T(-216//149), + T(72//149), + T(-12//149), )), + ) + 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, parallel, 1, order, source) + # elseif order == 4 + # left_boundary_plus = ( + # DerivativeCoefficientRow{T,1,4}(SVector(T(), + # T(), + # T(), + # T(), )), + # DerivativeCoefficientRow{T,1,5}(SVector(T(), + # T(), + # T(), + # T(), + # T(), )), + # DerivativeCoefficientRow{T,1,6}(SVector(T(), + # T(), + # T(), + # T(), + # T(), + # T(), )), + # DerivativeCoefficientRow{T,1,7}(SVector(T(), + # T(), + # T(), + # T(), + # T(), + # T(), + # T(), )), + # ) + # right_boundary_plus = ( + # DerivativeCoefficientRow{T,1,4}(SVector(T(), + # T(), + # T(), + # T(), )), + # DerivativeCoefficientRow{T,1,4}(SVector(T(), + # T(), + # T(), + # T(), )), + # DerivativeCoefficientRow{T,1,4}(SVector(T(), + # T(), + # T(), + # T(), )), + # DerivativeCoefficientRow{T,1,5}(SVector(T(), + # T(), + # T(), + # T(), + # T(), )), + # ) + # upper_coef_plus = SVector( T(), + # T(), + # T(), ) + # central_coef_plus = T() + # lower_coef_plus = SVector( T(), ) + # left_weights = SVector( T(), + # T(), + # T(), + # T(), ) + # right_weights = left_weights + # left_boundary_derivatives = Tuple{}() + # right_boundary_derivatives = left_boundary_derivatives + + # left_boundary_minus = ( + # DerivativeCoefficientRow{T,1,4}(SVector(T(), + # T(), + # T(), + # T(), )), + # DerivativeCoefficientRow{T,1,4}(SVector(T(), + # T(), + # T(), + # T(), )), + # DerivativeCoefficientRow{T,1,4}(SVector(T(), + # T(), + # T(), + # T(), )), + # DerivativeCoefficientRow{T,1,5}(SVector(T(), + # T(), + # T(), + # T(), + # T(), )), + # ) + # right_boundary_minus = ( + # DerivativeCoefficientRow{T,1,4}(SVector(T(), + # T(), + # T(), + # T(), )), + # DerivativeCoefficientRow{T,1,5}(SVector(T(), + # T(), + # T(), + # T(), + # T(), )), + # DerivativeCoefficientRow{T,1,6}(SVector(T(), + # T(), + # T(), + # T(), + # T(), + # T(), )), + # DerivativeCoefficientRow{T,1,7}(SVector(T(), + # T(), + # T(), + # T(), + # T(), + # T(), + # T(), )), + # ) + # 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, parallel, 1, order, source) else throw(ArgumentError("Order $order not implemented/derived.")) end diff --git a/test/upwind_operators.jl b/test/upwind_operators.jl index 46b7fb9e..a498c238 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:3 + for acc_order in 2:4 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) From 4ffdcfc3334bb832476b96bb8dd70a45151ed92d Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Wed, 25 Mar 2020 13:17:28 +0300 Subject: [PATCH 6/8] 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) From b79f9b872f5ba0683ee9d0e3df07728d2bec7264 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Wed, 25 Mar 2020 14:12:26 +0300 Subject: [PATCH 7/8] upwind operators of Mattsson(2017), order 6 --- dev/Mattsson2017.jl | 24 ++ src/SBP_coefficients/Mattsson2017.jl | 358 +++++++++++++++++++++++++-- test/upwind_operators.jl | 6 +- 3 files changed, 368 insertions(+), 20 deletions(-) diff --git a/dev/Mattsson2017.jl b/dev/Mattsson2017.jl index ab0470b5..66a05f87 100644 --- a/dev/Mattsson2017.jl +++ b/dev/Mattsson2017.jl @@ -72,3 +72,27 @@ 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) + +# order 6 +Qptmp = [ + -265//128688 1146190567//1737288000 -1596619//18384000 -55265831//579096000 26269819//3474576000 2464501//144774000 0 0 0 0 0 0 + -1116490567//1737288000 -8839//214480 190538869//347457600 102705469//694915200 413741//9651600 -191689861//3474576000 0 0 0 0 0 0 + 1096619//18384000 -135385429//347457600 -61067//321720 45137333//57909600 -253641811//694915200 70665929//579096000 -1//60 0 0 0 0 0 + 66965831//579096000 -208765789//694915200 -17623253//57909600 -18269//45960 410905829//347457600 -477953317//1158192000 2//15 -1//60 0 0 0 0 + -49219819//3474576000 293299//9651600 26422771//694915200 -141938309//347457600 -346583//643440 2217185207//1737288000 -1//2 2//15 -1//60 0 0 0 + -2374501//144774000 142906261//3474576000 -3137129//579096000 -29884283//1158192000 -630168407//1737288000 -3559//6128 4//3 -1//2 2//15 -1//60 0 0 + 0 0 0 0 1//30 -2//5 -7//12 4//3 -1//2 2//15 -1//60 0 + 0 0 0 0 0 1//30 -2//5 -7//12 4//3 -1//2 2//15 -1//60 +] +bnd = size(Qptmp, 1)-2 +Qp = zeros(eltype(Qptmp), 2*bnd+2, 2*bnd+2) +Qp[1:size(Qptmp,1), 1:size(Qptmp,2)] = Qptmp +for i in 1:bnd + Qp[end+1-i, end+1-size(Qptmp,1):end] = Qptmp[end:-1:1, i] +end +Qm = Matrix(-Qp') +B = zero(Qp); B[1,1] = -1; B[end,end] = 1 +m = [13613//43200, 12049//8640, 535//864, 1079//864, 7841//8640, 43837//43200] +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 22b01b52..56fac7b0 100644 --- a/src/SBP_coefficients/Mattsson2017.jl +++ b/src/SBP_coefficients/Mattsson2017.jl @@ -449,14 +449,240 @@ function first_derivative_coefficients(source::Mattsson2017, order::Int, T=Float left_boundary_derivatives, right_boundary_derivatives, lower_coef, central_coef, upper_coef, left_weights, right_weights, parallel, 1, order, source) - # elseif order == 4 + elseif order == 6 + left_boundary_plus = ( + DerivativeCoefficientRow{T,1,6}(SVector(T(-58148100//36496453), + T(1146190567//547446795), + T(-14369571//52137790), + T(-55265831//182482265), + T(26269819//1094893590), + T(9858004//182482265), )), + DerivativeCoefficientRow{T,1,6}(SVector(T(-1116490567//2422752675), + T(-954612//32303369), + T(190538869//484550535), + T(102705469//969101070), + T(4964892//161516845), + T(-191689861//4845505350), )), + DerivativeCoefficientRow{T,1,7}(SVector(T(9869571//102452500), + T(-135385429//215150250), + T(-2198412//7171675), + T(45137333//35858375), + T(-253641811//430300500), + T(70665929//358583750), + T(-72//2675), )), + DerivativeCoefficientRow{T,1,8}(SVector(T(66965831//723199750), + T(-208765789//867839700), + T(-17623253//72319975), + T(-657684//2066285), + T(410905829//433919850), + T(-477953317//1446399500), + T(576//5395), + T(-72//5395), )), + DerivativeCoefficientRow{T,1,9}(SVector(T(-49219819//3153258150), + T(3519588//105108605), + T(26422771//630651630), + T(-141938309//315325815), + T(-12476988//21021721), + T(2217185207//1576629075), + T(-4320//7841), + T(1152//7841), + T(-144//7841), )), + DerivativeCoefficientRow{T,1,10}(SVector(T(-9498004//587634985), + T(142906261//3525809910), + T(-3137129//587634985), + T(-29884283//1175269970), + T(-630168407//1762904955), + T(-9609300//16789571), + T(57600//43837), + T(-21600//43837), + T(5760//43837), + T(-720//43837), )), + ) + right_boundary_plus = ( + DerivativeCoefficientRow{T,1,6}(SVector(T(57671100//36496453), + T(-1116490567//547446795), + T(9869571//52137790), + T(66965831//182482265), + T(-49219819//1094893590), + T(-9498004//182482265), )), + DerivativeCoefficientRow{T,1,6}(SVector(T(1146190567//2422752675), + T(-954612//32303369), + T(-135385429//484550535), + T(-208765789//969101070), + T(3519588//161516845), + T(142906261//4845505350), )), + DerivativeCoefficientRow{T,1,6}(SVector(T(-14369571//102452500), + T(190538869//215150250), + T(-2198412//7171675), + T(-17623253//35858375), + T(26422771//430300500), + T(-3137129//358583750), )), + DerivativeCoefficientRow{T,1,6}(SVector(T(-55265831//723199750), + T(102705469//867839700), + T(45137333//72319975), + T(-657684//2066285), + T(-141938309//433919850), + T(-2298791//111261500), )), + DerivativeCoefficientRow{T,1,7}(SVector(T(26269819//3153258150), + T(4964892//105108605), + T(-253641811//630651630), + T(410905829//315325815), + T(-12476988//21021721), + T(-630168407//1576629075), + T(288//7841), )), + DerivativeCoefficientRow{T,1,8}(SVector(T(9858004//587634985), + T(-191689861//3525809910), + T(70665929//587634985), + T(-477953317//1175269970), + T(2217185207//1762904955), + T(-9609300//16789571), + T(-17280//43837), + T(1440//43837), )), + ) + upper_coef_plus = SVector( T(4//3), + T(-1//2), + T(2//15), + T(-1//60), ) + central_coef_plus = T(-7//12) + lower_coef_plus = SVector( T(-2//5), + T(1//30), ) + left_weights = SVector( T(13613//43200), + T(12049//8640), + T(535//864), + T(1079//864), + T(7841//8640), + T(43837//43200), ) + right_weights = left_weights + left_boundary_derivatives = Tuple{}() + right_boundary_derivatives = left_boundary_derivatives + + left_boundary_minus = ( + DerivativeCoefficientRow{T,1,6}(SVector(T(-57671100//36496453), + T(1116490567//547446795), + T(-9869571//52137790), + T(-66965831//182482265), + T(49219819//1094893590), + T(9498004//182482265), )), + DerivativeCoefficientRow{T,1,6}(SVector(T(-1146190567//2422752675), + T(954612//32303369), + T(135385429//484550535), + T(208765789//969101070), + T(-3519588//161516845), + T(-142906261//4845505350), )), + DerivativeCoefficientRow{T,1,6}(SVector(T(14369571//102452500), + T(-190538869//215150250), + T(2198412//7171675), + T(17623253//35858375), + T(-26422771//430300500), + T(3137129//358583750), )), + DerivativeCoefficientRow{T,1,6}(SVector(T(55265831//723199750), + T(-102705469//867839700), + T(-45137333//72319975), + T(657684//2066285), + T(141938309//433919850), + T(2298791//111261500), )), + DerivativeCoefficientRow{T,1,7}(SVector(T(-26269819//3153258150), + T(-4964892//105108605), + T(253641811//630651630), + T(-410905829//315325815), + T(12476988//21021721), + T(630168407//1576629075), + T(-288//7841), )), + DerivativeCoefficientRow{T,1,8}(SVector(T(-9858004//587634985), + T(191689861//3525809910), + T(-70665929//587634985), + T(477953317//1175269970), + T(-2217185207//1762904955), + T(9609300//16789571), + T(17280//43837), + T(-1440//43837), )), + ) + right_boundary_minus = ( + DerivativeCoefficientRow{T,1,6}(SVector(T(58148100//36496453), + T(-1146190567//547446795), + T(14369571//52137790), + T(55265831//182482265), + T(-26269819//1094893590), + T(-9858004//182482265), )), + DerivativeCoefficientRow{T,1,6}(SVector(T(1116490567//2422752675), + T(954612//32303369), + T(-190538869//484550535), + T(-102705469//969101070), + T(-4964892//161516845), + T(191689861//4845505350), )), + DerivativeCoefficientRow{T,1,7}(SVector(T(-9869571//102452500), + T(135385429//215150250), + T(2198412//7171675), + T(-45137333//35858375), + T(253641811//430300500), + T(-70665929//358583750), + T(72//2675), )), + DerivativeCoefficientRow{T,1,8}(SVector(T(-66965831//723199750), + T(208765789//867839700), + T(17623253//72319975), + T(657684//2066285), + T(-410905829//433919850), + T(477953317//1446399500), + T(-576//5395), + T(72//5395), )), + DerivativeCoefficientRow{T,1,9}(SVector(T(49219819//3153258150), + T(-3519588//105108605), + T(-26422771//630651630), + T(141938309//315325815), + T(12476988//21021721), + T(-2217185207//1576629075), + T(4320//7841), + T(-1152//7841), + T(144//7841), )), + DerivativeCoefficientRow{T,1,10}(SVector(T(9498004//587634985), + T(-142906261//3525809910), + T(3137129//587634985), + T(29884283//1175269970), + T(630168407//1762904955), + T(9609300//16789571), + T(-57600//43837), + T(21600//43837), + T(-5760//43837), + T(720//43837), )), + ) + 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, parallel, 1, order, source) + # elseif order == 6 # left_boundary_plus = ( - # DerivativeCoefficientRow{T,1,4}(SVector(T(), + # DerivativeCoefficientRow{T,1,6}(SVector(T(), # T(), # T(), - # T(), )), - # DerivativeCoefficientRow{T,1,5}(SVector(T(), - # T(), # T(), # T(), # T(), )), @@ -473,21 +699,70 @@ function first_derivative_coefficients(source::Mattsson2017, order::Int, T=Float # T(), # T(), # T(), )), + # DerivativeCoefficientRow{T,1,8}(SVector(T(), + # T(), + # T(), + # T(), + # T(), + # T(), + # T(), + # T(), )), + # DerivativeCoefficientRow{T,1,9}(SVector(T(), + # T(), + # T(), + # T(), + # T(), + # T(), + # T(), + # T(), + # T(), )), + # DerivativeCoefficientRow{T,1,10}(SVector(T(), + # T(), + # T(), + # T(), + # T(), + # T(), + # T(), + # T(), + # T(), + # T(), )), # ) # right_boundary_plus = ( - # DerivativeCoefficientRow{T,1,4}(SVector(T(), + # DerivativeCoefficientRow{T,1,6}(SVector(T(), + # T(), + # T(), # T(), # T(), # T(), )), - # DerivativeCoefficientRow{T,1,4}(SVector(T(), + # DerivativeCoefficientRow{T,1,6}(SVector(T(), + # T(), + # T(), + # T(), + # T(), + # T(), )), + # DerivativeCoefficientRow{T,1,6}(SVector(T(), + # T(), + # T(), # T(), # T(), # T(), )), - # DerivativeCoefficientRow{T,1,4}(SVector(T(), + # DerivativeCoefficientRow{T,1,6}(SVector(T(), + # T(), + # T(), # T(), # T(), # T(), )), - # DerivativeCoefficientRow{T,1,5}(SVector(T(), + # DerivativeCoefficientRow{T,1,7}(SVector(T(), + # T(), + # T(), + # T(), + # T(), + # T(), + # T(), )), + # DerivativeCoefficientRow{T,1,8}(SVector(T(), + # T(), + # T(), + # T(), # T(), # T(), # T(), @@ -498,39 +773,61 @@ function first_derivative_coefficients(source::Mattsson2017, order::Int, T=Float # T(), ) # central_coef_plus = T() # lower_coef_plus = SVector( T(), ) + # T(), ) # left_weights = SVector( T(), # T(), # T(), + # T(), + # T(), # T(), ) # right_weights = left_weights # left_boundary_derivatives = Tuple{}() # right_boundary_derivatives = left_boundary_derivatives # left_boundary_minus = ( - # DerivativeCoefficientRow{T,1,4}(SVector(T(), + # DerivativeCoefficientRow{T,1,6}(SVector(T(), + # T(), + # T(), # T(), # T(), # T(), )), - # DerivativeCoefficientRow{T,1,4}(SVector(T(), + # DerivativeCoefficientRow{T,1,6}(SVector(T(), + # T(), + # T(), # T(), # T(), # T(), )), - # DerivativeCoefficientRow{T,1,4}(SVector(T(), + # DerivativeCoefficientRow{T,1,6}(SVector(T(), + # T(), + # T(), # T(), # T(), # T(), )), - # DerivativeCoefficientRow{T,1,5}(SVector(T(), + # DerivativeCoefficientRow{T,1,6}(SVector(T(), + # T(), # T(), # T(), # T(), # T(), )), - # ) - # right_boundary_minus = ( - # DerivativeCoefficientRow{T,1,4}(SVector(T(), + # DerivativeCoefficientRow{T,1,7}(SVector(T(), + # T(), + # T(), + # T(), # T(), # T(), # T(), )), - # DerivativeCoefficientRow{T,1,5}(SVector(T(), + # DerivativeCoefficientRow{T,1,8}(SVector(T(), + # T(), + # T(), + # T(), + # T(), + # T(), + # T(), + # T(), )), + # ) + # right_boundary_minus = ( + # DerivativeCoefficientRow{T,1,6}(SVector(T(), + # T(), # T(), # T(), # T(), @@ -548,6 +845,33 @@ function first_derivative_coefficients(source::Mattsson2017, order::Int, T=Float # T(), # T(), # T(), )), + # DerivativeCoefficientRow{T,1,8}(SVector(T(), + # T(), + # T(), + # T(), + # T(), + # T(), + # T(), + # T(), )), + # DerivativeCoefficientRow{T,1,9}(SVector(T(), + # T(), + # T(), + # T(), + # T(), + # T(), + # T(), + # T(), + # T(), )), + # DerivativeCoefficientRow{T,1,10}(SVector(T(), + # T(), + # T(), + # T(), + # T(), + # T(), + # T(), + # T(), + # T(), + # T(), )), # ) # upper_coef_minus = .- lower_coef_plus # central_coef_minus = .- central_coef_plus diff --git a/test/upwind_operators.jl b/test/upwind_operators.jl index 54c89f73..265cee2c 100644 --- a/test/upwind_operators.jl +++ b/test/upwind_operators.jl @@ -5,12 +5,12 @@ using SummationByPartsOperators # check construction of interior part of upwind operators @testset "Check interior parts" begin N = 21 - xmin = 0//1 - xmax = (N + 1) // 1 + xmin = 0. + xmax = Float64(N + 1) interior = 10:N-10 - for acc_order in 2:5 + for acc_order in 2:6 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) From 7ae4991782cb18989b910f0125af56acb440be51 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Wed, 25 Mar 2020 14:52:32 +0300 Subject: [PATCH 8/8] upwind operators of Mattsson(2017), order 7 --- dev/Mattsson2017.jl | 25 +++ src/SBP_coefficients/Mattsson2017.jl | 239 ++++++++++++++++++++++++++- test/upwind_operators.jl | 2 +- 3 files changed, 264 insertions(+), 2 deletions(-) diff --git a/dev/Mattsson2017.jl b/dev/Mattsson2017.jl index 66a05f87..d07d37fb 100644 --- a/dev/Mattsson2017.jl +++ b/dev/Mattsson2017.jl @@ -96,3 +96,28 @@ m = [13613//43200, 12049//8640, 535//864, 1079//864, 7841//8640, 43837//43200] 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 7 +Qptmp = Rational{Int128}[ + -265//300272 1587945773//2432203200 -1926361//25737600 -84398989//810734400 48781961//4864406400 3429119//202683600 0 0 0 0 0 0 0 + -1570125773//2432203200 -26517//1501360 240029831//486440640 202934303//972881280 118207//13512240 -231357719//4864406400 0 0 0 0 0 0 0 + 1626361//25737600 -206937767//486440640 -61067//750680 49602727//81073440 -43783933//194576256 51815011//810734400 -1//140 0 0 0 0 0 0 + 91418989//810734400 -53314099//194576256 -33094279//81073440 -18269//107240 440626231//486440640 -365711063//1621468800 1//15 -1//140 0 0 0 0 0 + -62551961//4864406400 799//35280 82588241//972881280 -279245719//486440640 -346583//1501360 2312302333//2432203200 -3//10 1//15 -1//140 0 0 0 0 + -3375119//202683600 202087559//4864406400 -11297731//810734400 61008503//1621468800 -1360092253//2432203200 -10677//42896 1 -3//10 1//15 -1//140 0 0 0 + 0 0 0 -1//105 1//10 -3//5 -1//4 1 -3//10 1//15 -1//140 0 0 + 0 0 0 0 -1//105 1//10 -3//5 -1//4 1 -3//10 1//15 -1//140 0 + 0 0 0 0 0 -1//105 1//10 -3//5 -1//4 1 -3//10 1//15 -1//140 +] +bnd = size(Qptmp, 1)-3 +Qp = zeros(eltype(Qptmp), 2*bnd+3, 2*bnd+3) +Qp[1:size(Qptmp,1), 1:size(Qptmp,2)] = Qptmp +for i in 1:bnd + Qp[end+1-i, end+1-size(Qptmp,1):end] = Qptmp[end:-1:1, i] +end +Qm = Matrix(-Qp') +B = zero(Qp); B[1,1] = -1; B[end,end] = 1 +m = [19087//60480, 84199//60480, 18869//30240, 37621//30240, 55031//60480, 61343//60480] +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 56fac7b0..506fdbcd 100644 --- a/src/SBP_coefficients/Mattsson2017.jl +++ b/src/SBP_coefficients/Mattsson2017.jl @@ -655,6 +655,242 @@ 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 == 7 + left_boundary_plus = ( + DerivativeCoefficientRow{T,1,6}(SVector(T(-81216540//51172247), + T(1587945773//767583705), + T(-17337249//73103210), + T(-84398989//255861235), + T(48781961//1535167410), + T(13716476//255861235), )), + DerivativeCoefficientRow{T,1,6}(SVector(T(-1570125773//3386062785), + T(-2863836//225737519), + T(240029831//677212557), + T(202934303//1354425114), + T(1418484//225737519), + T(-231357719//6772125570), )), + DerivativeCoefficientRow{T,1,7}(SVector(T(14637249//144536540), + T(-206937767//303526734), + T(-6595236//50587789), + T(49602727//50587789), + T(-218919665//607053468), + T(51815011//505877890), + T(-216//18869), )), + DerivativeCoefficientRow{T,1,8}(SVector(T(91418989//1008619010), + T(-266570495//1210342812), + T(-33094279//100861901), + T(-1973052//14408843), + T(440626231//605171406), + T(-365711063//2017238020), + T(2016//37621), + T(-216//37621), )), + DerivativeCoefficientRow{T,1,9}(SVector(T(-62551961//4426143330), + T(9588//385217), + T(82588241//885228666), + T(-279245719//442614333), + T(-37430964//147538111), + T(2312302333//2213071665), + T(-18144//55031), + T(4032//55031), + T(-432//55031), )), + DerivativeCoefficientRow{T,1,10}(SVector(T(-13500476//822302915), + T(202087559//4933817490), + T(-11297731//822302915), + T(61008503//1644605830), + T(-1360092253//2466908745), + T(-5765580//23494369), + T(60480//61343), + T(-18144//61343), + T(4032//61343), + T(-432//61343), )), + ) + right_boundary_plus = ( + DerivativeCoefficientRow{T,1,6}(SVector(T(80930340//51172247), + T(-1570125773//767583705), + T(14637249//73103210), + T(91418989//255861235), + T(-62551961//1535167410), + T(-13500476//255861235), )), + DerivativeCoefficientRow{T,1,6}(SVector(T(1587945773//3386062785), + T(-2863836//225737519), + T(-206937767//677212557), + T(-266570495//1354425114), + T(9588//589393), + T(202087559//6772125570), )), + DerivativeCoefficientRow{T,1,6}(SVector(T(-17337249//144536540), + T(240029831//303526734), + T(-6595236//50587789), + T(-33094279//50587789), + T(82588241//607053468), + T(-11297731//505877890), )), + DerivativeCoefficientRow{T,1,7}(SVector(T(-84398989//1008619010), + T(202934303//1210342812), + T(49602727//100861901), + T(-1973052//14408843), + T(-279245719//605171406), + T(61008503//2017238020), + T(-288//37621), )), + DerivativeCoefficientRow{T,1,8}(SVector(T(48781961//4426143330), + T(1418484//147538111), + T(-218919665//885228666), + T(440626231//442614333), + T(-37430964//147538111), + T(-1360092253//2213071665), + T(6048//55031), + T(-576//55031), )), + DerivativeCoefficientRow{T,1,9}(SVector(T(13716476//822302915), + T(-231357719//4933817490), + T(51815011//822302915), + T(-365711063//1644605830), + T(2312302333//2466908745), + T(-5765580//23494369), + T(-36288//61343), + T(6048//61343), + T(-576//61343), )), + ) + upper_coef_plus = SVector( T(1), + T(-3//10), + T(1//15), + T(-1//140), ) + central_coef_plus = T(-1//4) + lower_coef_plus = SVector( T(-3//5), + T(1//10), + T(-1//105), ) + left_weights = SVector( T(19087//60480), + T(84199//60480), + T(18869//30240), + T(37621//30240), + T(55031//60480), + T(61343//60480), ) + right_weights = left_weights + left_boundary_derivatives = Tuple{}() + right_boundary_derivatives = left_boundary_derivatives + + left_boundary_minus = ( + DerivativeCoefficientRow{T,1,6}(SVector(T(-80930340//51172247), + T(1570125773//767583705), + T(-14637249//73103210), + T(-91418989//255861235), + T(62551961//1535167410), + T(13500476//255861235), )), + DerivativeCoefficientRow{T,1,6}(SVector(T(-1587945773//3386062785), + T(2863836//225737519), + T(206937767//677212557), + T(266570495//1354425114), + T(-9588//589393), + T(-202087559//6772125570), )), + DerivativeCoefficientRow{T,1,6}(SVector(T(17337249//144536540), + T(-240029831//303526734), + T(6595236//50587789), + T(33094279//50587789), + T(-82588241//607053468), + T(11297731//505877890), )), + DerivativeCoefficientRow{T,1,7}(SVector(T(84398989//1008619010), + T(-202934303//1210342812), + T(-49602727//100861901), + T(1973052//14408843), + T(279245719//605171406), + T(-61008503//2017238020), + T(288//37621), )), + DerivativeCoefficientRow{T,1,8}(SVector(T(-48781961//4426143330), + T(-1418484//147538111), + T(218919665//885228666), + T(-440626231//442614333), + T(37430964//147538111), + T(1360092253//2213071665), + T(-6048//55031), + T(576//55031), )), + DerivativeCoefficientRow{T,1,9}(SVector(T(-13716476//822302915), + T(231357719//4933817490), + T(-51815011//822302915), + T(365711063//1644605830), + T(-2312302333//2466908745), + T(5765580//23494369), + T(36288//61343), + T(-6048//61343), + T(576//61343), )), + ) + right_boundary_minus = ( + DerivativeCoefficientRow{T,1,6}(SVector(T(81216540//51172247), + T(-1587945773//767583705), + T(17337249//73103210), + T(84398989//255861235), + T(-48781961//1535167410), + T(-13716476//255861235), )), + DerivativeCoefficientRow{T,1,6}(SVector(T(1570125773//3386062785), + T(2863836//225737519), + T(-240029831//677212557), + T(-202934303//1354425114), + T(-1418484//225737519), + T(231357719//6772125570), )), + DerivativeCoefficientRow{T,1,7}(SVector(T(-14637249//144536540), + T(206937767//303526734), + T(6595236//50587789), + T(-49602727//50587789), + T(218919665//607053468), + T(-51815011//505877890), + T(216//18869), )), + DerivativeCoefficientRow{T,1,8}(SVector(T(-91418989//1008619010), + T(266570495//1210342812), + T(33094279//100861901), + T(1973052//14408843), + T(-440626231//605171406), + T(365711063//2017238020), + T(-2016//37621), + T(216//37621), )), + DerivativeCoefficientRow{T,1,9}(SVector(T(62551961//4426143330), + T(-9588//385217), + T(-82588241//885228666), + T(279245719//442614333), + T(37430964//147538111), + T(-2312302333//2213071665), + T(18144//55031), + T(-4032//55031), + T(432//55031), )), + DerivativeCoefficientRow{T,1,10}(SVector(T(13500476//822302915), + T(-202087559//4933817490), + T(11297731//822302915), + T(-61008503//1644605830), + T(1360092253//2466908745), + T(5765580//23494369), + T(-60480//61343), + T(18144//61343), + T(-4032//61343), + T(432//61343), )), + ) + 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 @@ -772,7 +1008,8 @@ function first_derivative_coefficients(source::Mattsson2017, order::Int, T=Float # T(), # T(), ) # central_coef_plus = T() - # lower_coef_plus = SVector( T(), ) + # lower_coef_plus = SVector( T(), + # T(), # T(), ) # left_weights = SVector( T(), # T(), diff --git a/test/upwind_operators.jl b/test/upwind_operators.jl index 265cee2c..640d04e1 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:6 + 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)