Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Upwind operators of Mattsson (2017) #61

Merged
merged 8 commits into from
Mar 25, 2020
Merged

Upwind operators of Mattsson (2017) #61

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
967228d
fix indentation
ranocha Mar 25, 2020
b37b19e
add checks for LowerOffset > 0 for non-periodic operators
ranocha Mar 25, 2020
cce0765
upwind operators of Mattsson(2017), order 2
ranocha Mar 25, 2020
a5f0041
upwind operators of Mattsson(2017), order 3
ranocha Mar 25, 2020
e2b7247
upwind operators of Mattsson(2017), order 4
ranocha Mar 25, 2020
4ffdcfc
upwind operators of Mattsson(2017), order 5
ranocha Mar 25, 2020
b79f9b8
upwind operators of Mattsson(2017), order 6
ranocha Mar 25, 2020
7ae4991
upwind operators of Mattsson(2017), order 7
ranocha Mar 25, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
123 changes: 123 additions & 0 deletions dev/Mattsson2017.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,123 @@
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
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//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 = [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)

# 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)

# 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)
Loading