Skip to content

Commit

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

0 comments on commit 4ffdcfc

Please sign in to comment.