Mercurial > repos > public > sbplib_julia

comparison test/SbpOperators/volumeops/derivatives/dissipation_test.jl @ 1395:bdcdbd4ea9cd feature/boundary_conditions

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with default. Comment out broken tests for boundary_conditions at sat
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Wed, 26 Jul 2023 21:35:50 +0200
parents 760a4a1ec4b7
children 43aaf710463e
comparison
equal deleted inserted replaced
1217:ea2e8254820a 1395:bdcdbd4ea9cd
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 = equidistant_grid(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 = equidistant_grid(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 = eval_on(g, x->monomial(x,k))
45 vₚₓ = eval_on(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 = equidistant_grid(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
81 @testset "2D" begin
82 N = 20
83 g = equidistant_grid((N,2N), (0,0), (2,1))
84 h = spacing.(g.grids)
85
86 D,Dᵀ = undivided_skewed04(g, 3, 2)
87
88 v = eval_on(g, x->monomial(x[1],4)*monomial(x[2],3))
89 d³vdy³ = eval_on(g, x->monomial(x[1],4)*monomial(x[2],0))
90
91 @test D*v ≈ h[2]^3*d³vdy³
92 end
93 end
94
95 @testset "dissipation_interior_weights" begin
96 @test dissipation_interior_weights(1) == (-1, 1)
97 @test dissipation_interior_weights(2) == (1,-2, 1)
98 @test dissipation_interior_weights(3) == (-1, 3,-3, 1)
99 @test dissipation_interior_weights(4) == (1, -4, 6, -4, 1)
100 end
101
102 @testset "dissipation_interior_stencil" begin
103 @test dissipation_interior_stencil(dissipation_interior_weights(1)) == Stencil(-1, 1, center=2)
104 @test dissipation_interior_stencil(dissipation_interior_weights(2)) == Stencil( 1,-2, 1, center=2)
105 @test dissipation_interior_stencil(dissipation_interior_weights(3)) == Stencil(-1, 3,-3, 1, center=3)
106 @test dissipation_interior_stencil(dissipation_interior_weights(4)) == Stencil( 1,-4, 6,-4, 1, center=3)
107 end
108
109 @testset "dissipation_transpose_interior_stencil" begin
110 @test dissipation_transpose_interior_stencil(dissipation_interior_weights(1)) == Stencil(1,-1, center=1)
111 @test dissipation_transpose_interior_stencil(dissipation_interior_weights(2)) == Stencil(1,-2, 1, center=2)
112 @test dissipation_transpose_interior_stencil(dissipation_interior_weights(3)) == Stencil(1,-3, 3,-1, center=2)
113 @test dissipation_transpose_interior_stencil(dissipation_interior_weights(4)) == Stencil(1,-4, 6,-4, 1, center=3)
114 end
115
116 @testset "midpoint" begin
117 @test midpoint((1,1)) == 2
118 @test midpoint((1,1,1)) == 2
119 @test midpoint((1,1,1,1)) == 3
120 @test midpoint((1,1,1,1,1)) == 3
121 end
122
123 @testset "midpoint_transpose" begin
124 @test midpoint_transpose((1,1)) == 1
125 @test midpoint_transpose((1,1,1)) == 2
126 @test midpoint_transpose((1,1,1,1)) == 2
127 @test midpoint_transpose((1,1,1,1,1)) == 3
128 end
129
130 @testset "dissipation_lower_closure_size" begin
131 @test dissipation_lower_closure_size((1,1)) == 1
132 @test dissipation_lower_closure_size((1,1,1)) == 1
133 @test dissipation_lower_closure_size((1,1,1,1)) == 2
134 @test dissipation_lower_closure_size((1,1,1,1,1)) == 2
135 end
136
137 @testset "dissipation_upper_closure_size" begin
138 @test dissipation_upper_closure_size((1,1)) == 0
139 @test dissipation_upper_closure_size((1,1,1)) == 1
140 @test dissipation_upper_closure_size((1,1,1,1)) == 1
141 @test dissipation_upper_closure_size((1,1,1,1,1)) == 2
142 end
143
144 @testset "dissipation_lower_closure_stencils" begin
145 cases = (
146 (-1,1) => (
147 Stencil(-1, 1, center=1),
148 ),
149 (1,-2,1) => (
150 Stencil( 1,-2, 1, center=1),
151 ),
152 (-1,3,-3,1) => (
153 Stencil(-1,3,-3,1, center=1),
154 Stencil(-1,3,-3,1, center=2),
155 ),
156 (1, -4, 6, -4, 1) => (
157 Stencil(1, -4, 6, -4, 1, center=1),
158 Stencil(1, -4, 6, -4, 1, center=2),
159 )
160 )
161 @testset "interior_weights = $w" for (w, closure_stencils) ∈ cases
162 @test dissipation_lower_closure_stencils(w) == closure_stencils
163 end
164 end
165
166 @testset "dissipation_upper_closure_stencils" begin
167 cases = (
168 (-1,1) => (),
169 (1,-2,1) => (
170 Stencil( 1,-2, 1, center=3),
171 ),
172 (-1,3,-3,1) => (
173 Stencil(-1,3,-3,1, center=4),
174 ),
175 (1, -4, 6, -4, 1) => (
176 Stencil(1, -4, 6, -4, 1, center=4),
177 Stencil(1, -4, 6, -4, 1, center=5),
178 )
179 )
180 @testset "interior_weights = $w" for (w, closure_stencils) ∈ cases
181 @test dissipation_upper_closure_stencils(w) == closure_stencils
182 end
183 end
184
185
186 @testset "dissipation_transpose_lower_closure_stencils" begin
187 cases = (
188 (-1,1) => (
189 Stencil(-1,-1, 0, center=1),
190 Stencil( 1, 1,-1, center=2),
191 ),
192 (1,-2,1) => (
193 Stencil( 1, 1, 0, 0, center=1),
194 Stencil(-2,-2, 1, 0, center=2),
195 Stencil( 1, 1,-2, 1, center=3),
196 ),
197 (-1,3,-3,1) => (
198 Stencil(-1,-1,-1, 0, 0, 0, center=1),
199 Stencil( 3, 3, 3,-1, 0, 0, center=2),
200 Stencil(-3,-3,-3, 3,-1, 0, center=3),
201 Stencil( 1, 1, 1,-3, 3,-1, center=4),
202 ),
203 )
204 @testset "interior_weights = $w" for (w, closure_stencils) ∈ cases
205 @test dissipation_transpose_lower_closure_stencils(w) == closure_stencils
206 end
207 end
208
209 @testset "dissipation_transpose_upper_closure_stencils" begin
210 cases = (
211 (-1,1) => (
212 Stencil( 1,-1, center = 1),
213 Stencil( 0, 1, center = 2),
214 ),
215 (1,-2,1) => (
216 Stencil( 1,-2, 1, 1, center=2),
217 Stencil( 0, 1,-2,-2, center=3),
218 Stencil( 0, 0, 1, 1, center=4),
219 ),
220 (-1,3,-3,1) => (
221 Stencil( 1,-3, 3,-1,-1, center=2),
222 Stencil( 0, 1,-3, 3, 3, center=3),
223 Stencil( 0, 0, 1,-3,-3, center=4),
224 Stencil( 0, 0, 0, 1, 1, center=5),
225 ),
226 )
227 @testset "interior_weights = $w" for (w, closure_stencils) ∈ cases
228 @test dissipation_transpose_upper_closure_stencils(w) == closure_stencils
229 end
230 end