Skip to content

Commit

Permalink
upwind operators of Mattsson(2017), order 7
Browse files Browse the repository at this point in the history
  • Loading branch information
@ranocha
ranocha committed Mar 25, 2020
1 parent 8141440 commit 3ec9586
Show file tree
Hide file tree
Showing 3 changed files with 264 additions and 2 deletions.
25 changes: 25 additions & 0 deletions dev/Mattsson2017.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
239 changes: 238 additions & 1 deletion src/SBP_coefficients/Mattsson2017.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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(),
Expand Down
2 changes: 1 addition & 1 deletion test/upwind_operators.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down

0 comments on commit 3ec9586

Please sign in to comment.