Mercurial > repos > public > sbplib_julia

comparison test/SbpOperators/volumeops/derivatives/dissipation_test.jl @ 1221:b3b4d29b46c3 refactor/grids

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge default
author Jonatan Werpers <jonatan@werpers.com>
date Fri, 10 Feb 2023 08:36:56 +0100
parents 0905c3c674e9
children 7d52c4835d15
comparison
equal deleted inserted replaced
1220:93bba649aea2 1221:b3b4d29b46c3
1 using Test
2
3 using Sbplib.SbpOperators
4 using Sbplib.Grids
5 using Sbplib.LazyTensors
6
7 using Sbplib.SbpOperators: Stencil
8
9 using Sbplib.SbpOperators: dissipation_interior_weights
10 using Sbplib.SbpOperators: dissipation_interior_stencil, dissipation_transpose_interior_stencil
11 using Sbplib.SbpOperators: midpoint, midpoint_transpose
12 using Sbplib.SbpOperators: dissipation_lower_closure_size, dissipation_upper_closure_size
13 using Sbplib.SbpOperators: dissipation_lower_closure_stencils,dissipation_upper_closure_stencils
14 using Sbplib.SbpOperators: dissipation_transpose_lower_closure_stencils, dissipation_transpose_upper_closure_stencils
15
16 """
17 monomial(x,k)
18
19 Evaluates ``x^k/k!` with the convetion that it is ``0`` for all ``k<0``.
20 Has the property that ``d/dx monomial(x,k) = monomial(x,k-1)``
21 """
22 function monomial(x,k)
23 if k < 0
24 return zero(x)
25 end
26 x^k/factorial(k)
27 end
28
29 @testset "undivided_skewed04" begin
30 g = EquidistantGrid(20, 0., 11.)
31 D,Dᵀ = undivided_skewed04(g, 1)
32
33 @test D isa LazyTensor{Float64,1,1}
34 @test Dᵀ isa LazyTensor{Float64,1,1}
35
36 @testset "Accuracy conditions" begin
37 N = 20
38 g = EquidistantGrid(N, 0//1,2//1)
39 h = only(spacing(g))
40 @testset "D_$p" for p ∈ [1,2,3,4]
41 D,Dᵀ = undivided_skewed04(g, p)
42
43 @testset "x^$k" for k ∈ 0:p
44 v = evalOn(g, x->monomial(x,k))
45 vₚₓ = evalOn(g, x->monomial(x,k-p))
46
47 @test D*v == h^p * vₚₓ
48 end
49 end
50 end
51
52 @testset "transpose equality" begin
53 function get_matrix(D)
54 N = only(range_size(D))
55 M = only(domain_size(D))
56
57 Dmat = zeros(N,M)
58 e = zeros(M)
59 for i ∈ 1:M
60 if i > 1
61 e[i-1] = 0.
62 end
63 e[i] = 1.
64 Dmat[:,i] = D*e
65 end
66
67 return Dmat
68 end
69
70 g = EquidistantGrid(11, 0., 1.)
71 @testset "D_$p" for p ∈ [1,2,3,4]
72 D,Dᵀ = undivided_skewed04(g, p)
73
74 D̄ = get_matrix(D)
75 D̄ᵀ = get_matrix(Dᵀ)
76
77 @test D̄ == D̄ᵀ'
78 end
79 end
80 end
81
82 @testset "dissipation_interior_weights" begin
83 @test dissipation_interior_weights(1) == (-1, 1)
84 @test dissipation_interior_weights(2) == (1,-2, 1)
85 @test dissipation_interior_weights(3) == (-1, 3,-3, 1)
86 @test dissipation_interior_weights(4) == (1, -4, 6, -4, 1)
87 end
88
89 @testset "dissipation_interior_stencil" begin
90 @test dissipation_interior_stencil(dissipation_interior_weights(1)) == Stencil(-1, 1, center=2)
91 @test dissipation_interior_stencil(dissipation_interior_weights(2)) == Stencil( 1,-2, 1, center=2)
92 @test dissipation_interior_stencil(dissipation_interior_weights(3)) == Stencil(-1, 3,-3, 1, center=3)
93 @test dissipation_interior_stencil(dissipation_interior_weights(4)) == Stencil( 1,-4, 6,-4, 1, center=3)
94 end
95
96 @testset "dissipation_transpose_interior_stencil" begin
97 @test dissipation_transpose_interior_stencil(dissipation_interior_weights(1)) == Stencil(1,-1, center=1)
98 @test dissipation_transpose_interior_stencil(dissipation_interior_weights(2)) == Stencil(1,-2, 1, center=2)
99 @test dissipation_transpose_interior_stencil(dissipation_interior_weights(3)) == Stencil(1,-3, 3,-1, center=2)
100 @test dissipation_transpose_interior_stencil(dissipation_interior_weights(4)) == Stencil(1,-4, 6,-4, 1, center=3)
101 end
102
103 @testset "midpoint" begin
104 @test midpoint((1,1)) == 2
105 @test midpoint((1,1,1)) == 2
106 @test midpoint((1,1,1,1)) == 3
107 @test midpoint((1,1,1,1,1)) == 3
108 end
109
110 @testset "midpoint_transpose" begin
111 @test midpoint_transpose((1,1)) == 1
112 @test midpoint_transpose((1,1,1)) == 2
113 @test midpoint_transpose((1,1,1,1)) == 2
114 @test midpoint_transpose((1,1,1,1,1)) == 3
115 end
116
117 @testset "dissipation_lower_closure_size" begin
118 @test dissipation_lower_closure_size((1,1)) == 1
119 @test dissipation_lower_closure_size((1,1,1)) == 1
120 @test dissipation_lower_closure_size((1,1,1,1)) == 2
121 @test dissipation_lower_closure_size((1,1,1,1,1)) == 2
122 end
123
124 @testset "dissipation_upper_closure_size" begin
125 @test dissipation_upper_closure_size((1,1)) == 0
126 @test dissipation_upper_closure_size((1,1,1)) == 1
127 @test dissipation_upper_closure_size((1,1,1,1)) == 1
128 @test dissipation_upper_closure_size((1,1,1,1,1)) == 2
129 end
130
131 @testset "dissipation_lower_closure_stencils" begin
132 cases = (
133 (-1,1) => (
134 Stencil(-1, 1, center=1),
135 ),
136 (1,-2,1) => (
137 Stencil( 1,-2, 1, center=1),
138 ),
139 (-1,3,-3,1) => (
140 Stencil(-1,3,-3,1, center=1),
141 Stencil(-1,3,-3,1, center=2),
142 ),
143 (1, -4, 6, -4, 1) => (
144 Stencil(1, -4, 6, -4, 1, center=1),
145 Stencil(1, -4, 6, -4, 1, center=2),
146 )
147 )
148 @testset "interior_weights = $w" for (w, closure_stencils) ∈ cases
149 @test dissipation_lower_closure_stencils(w) == closure_stencils
150 end
151 end
152
153 @testset "dissipation_upper_closure_stencils" begin
154 cases = (
155 (-1,1) => (),
156 (1,-2,1) => (
157 Stencil( 1,-2, 1, center=3),
158 ),
159 (-1,3,-3,1) => (
160 Stencil(-1,3,-3,1, center=4),
161 ),
162 (1, -4, 6, -4, 1) => (
163 Stencil(1, -4, 6, -4, 1, center=4),
164 Stencil(1, -4, 6, -4, 1, center=5),
165 )
166 )
167 @testset "interior_weights = $w" for (w, closure_stencils) ∈ cases
168 @test dissipation_upper_closure_stencils(w) == closure_stencils
169 end
170 end
171
172
173 @testset "dissipation_transpose_lower_closure_stencils" begin
174 cases = (
175 (-1,1) => (
176 Stencil(-1,-1, 0, center=1),
177 Stencil( 1, 1,-1, center=2),
178 ),
179 (1,-2,1) => (
180 Stencil( 1, 1, 0, 0, center=1),
181 Stencil(-2,-2, 1, 0, center=2),
182 Stencil( 1, 1,-2, 1, center=3),
183 ),
184 (-1,3,-3,1) => (
185 Stencil(-1,-1,-1, 0, 0, 0, center=1),
186 Stencil( 3, 3, 3,-1, 0, 0, center=2),
187 Stencil(-3,-3,-3, 3,-1, 0, center=3),
188 Stencil( 1, 1, 1,-3, 3,-1, center=4),
189 ),
190 )
191 @testset "interior_weights = $w" for (w, closure_stencils) ∈ cases
192 @test dissipation_transpose_lower_closure_stencils(w) == closure_stencils
193 end
194 end
195
196 @testset "dissipation_transpose_upper_closure_stencils" begin
197 cases = (
198 (-1,1) => (
199 Stencil( 1,-1, center = 1),
200 Stencil( 0, 1, center = 2),
201 ),
202 (1,-2,1) => (
203 Stencil( 1,-2, 1, 1, center=2),
204 Stencil( 0, 1,-2,-2, center=3),
205 Stencil( 0, 0, 1, 1, center=4),
206 ),
207 (-1,3,-3,1) => (
208 Stencil( 1,-3, 3,-1,-1, center=2),
209 Stencil( 0, 1,-3, 3, 3, center=3),
210 Stencil( 0, 0, 1,-3,-3, center=4),
211 Stencil( 0, 0, 0, 1, 1, center=5),
212 ),
213 )
214 @testset "interior_weights = $w" for (w, closure_stencils) ∈ cases
215 @test dissipation_transpose_upper_closure_stencils(w) == closure_stencils
216 end
217 end