Mercurial > repos > public > sbplib_julia

comparison test/SbpOperators/volumeops/derivatives/second_derivative_test.jl @ 1207:f1c2a4fa0ee1 performance/get_region_type_inference

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, 03 Feb 2023 22:14:47 +0100
parents c94a12327737
children 7d52c4835d15
comparison
equal deleted inserted replaced
919:b41180efb6c2 1207:f1c2a4fa0ee1
4 using Sbplib.Grids 4 using Sbplib.Grids
5 using Sbplib.LazyTensors 5 using Sbplib.LazyTensors
6 6
7 import Sbplib.SbpOperators.VolumeOperator 7 import Sbplib.SbpOperators.VolumeOperator
8 8
9 # TODO: Refactor these test to look more like the tests in first_derivative_test.jl.
10
9 @testset "SecondDerivative" begin 11 @testset "SecondDerivative" begin
10 stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4) 12 operator_path = sbp_operators_path()*"standard_diagonal.toml"
13 stencil_set = read_stencil_set(operator_path; order=4)
11 inner_stencil = parse_stencil(stencil_set["D2"]["inner_stencil"]) 14 inner_stencil = parse_stencil(stencil_set["D2"]["inner_stencil"])
12 closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"]) 15 closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"])
13 Lx = 3.5 16 Lx = 3.5
14 Ly = 3. 17 Ly = 3.
15 g_1D = EquidistantGrid(121, 0.0, Lx) 18 g_1D = EquidistantGrid(121, 0.0, Lx)
16 g_2D = EquidistantGrid((121,123), (0.0, 0.0), (Lx, Ly)) 19 g_2D = EquidistantGrid((121,123), (0.0, 0.0), (Lx, Ly))
17 20
18 @testset "Constructors" begin 21 @testset "Constructors" begin
19 @testset "1D" begin 22 @testset "1D" begin
20 Dₓₓ = second_derivative(g_1D,inner_stencil,closure_stencils) 23 Dₓₓ = second_derivative(g_1D,inner_stencil,closure_stencils,1)
21 @test Dₓₓ == second_derivative(g_1D,inner_stencil,closure_stencils,1) 24 @test Dₓₓ == second_derivative(g_1D,inner_stencil,closure_stencils)
25 @test Dₓₓ == second_derivative(g_1D,stencil_set,1)
26 @test Dₓₓ == second_derivative(g_1D,stencil_set)
22 @test Dₓₓ isa VolumeOperator 27 @test Dₓₓ isa VolumeOperator
23 end 28 end
24 @testset "2D" begin 29 @testset "2D" begin
25 Dₓₓ = second_derivative(g_2D,inner_stencil,closure_stencils,1) 30 Dₓₓ = second_derivative(g_2D,inner_stencil,closure_stencils,1)
26 D2 = second_derivative(g_1D,inner_stencil,closure_stencils) 31 D2 = second_derivative(g_1D,inner_stencil,closure_stencils,1)
27 I = IdentityMapping{Float64}(size(g_2D)[2]) 32 I = IdentityTensor{Float64}(size(g_2D)[2])
28 @test Dₓₓ == D2⊗I 33 @test Dₓₓ == D2⊗I
29 @test Dₓₓ isa TensorMapping{T,2,2} where T 34 @test Dₓₓ == second_derivative(g_2D,stencil_set,1)
35 @test Dₓₓ isa LazyTensor{T,2,2} where T
30 end 36 end
31 end 37 end
32 38
33 # Exact differentiation is measured point-wise. In other cases 39 # Exact differentiation is measured point-wise. In other cases
34 # the error is measured in the l2-norm. 40 # the error is measured in the l2-norm.
45 vₓₓ = evalOn(g_1D,x -> -sin(x)) 51 vₓₓ = evalOn(g_1D,x -> -sin(x))
46 52
47 # 2nd order interior stencil, 1nd order boundary stencil, 53 # 2nd order interior stencil, 1nd order boundary stencil,
48 # implies that L*v should be exact for monomials up to order 2. 54 # implies that L*v should be exact for monomials up to order 2.
49 @testset "2nd order" begin 55 @testset "2nd order" begin
50 stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=2) 56 stencil_set = read_stencil_set(operator_path; order=2)
51 inner_stencil = parse_stencil(stencil_set["D2"]["inner_stencil"]) 57 Dₓₓ = second_derivative(g_1D,stencil_set)
52 closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"])
53 Dₓₓ = second_derivative(g_1D,inner_stencil,closure_stencils)
54 @test Dₓₓ*monomials[1] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10 58 @test Dₓₓ*monomials[1] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10
55 @test Dₓₓ*monomials[2] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10 59 @test Dₓₓ*monomials[2] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10
56 @test Dₓₓ*monomials[3] ≈ monomials[1] atol = 5e-10 60 @test Dₓₓ*monomials[3] ≈ monomials[1] atol = 5e-10
57 @test Dₓₓ*v ≈ vₓₓ rtol = 5e-2 norm = l2 61 @test Dₓₓ*v ≈ vₓₓ rtol = 5e-2 norm = l2
58 end 62 end
59 63
60 # 4th order interior stencil, 2nd order boundary stencil, 64 # 4th order interior stencil, 2nd order boundary stencil,
61 # implies that L*v should be exact for monomials up to order 3. 65 # implies that L*v should be exact for monomials up to order 3.
62 @testset "4th order" begin 66 @testset "4th order" begin
63 stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4) 67 stencil_set = read_stencil_set(operator_path; order=4)
64 inner_stencil = parse_stencil(stencil_set["D2"]["inner_stencil"]) 68 Dₓₓ = second_derivative(g_1D,stencil_set)
65 closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"])
66 Dₓₓ = second_derivative(g_1D,inner_stencil,closure_stencils)
67 # NOTE: high tolerances for checking the "exact" differentiation 69 # NOTE: high tolerances for checking the "exact" differentiation
68 # due to accumulation of round-off errors/cancellation errors? 70 # due to accumulation of round-off errors/cancellation errors?
69 @test Dₓₓ*monomials[1] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10 71 @test Dₓₓ*monomials[1] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10
70 @test Dₓₓ*monomials[2] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10 72 @test Dₓₓ*monomials[2] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10
71 @test Dₓₓ*monomials[3] ≈ monomials[1] atol = 5e-10 73 @test Dₓₓ*monomials[3] ≈ monomials[1] atol = 5e-10
86 v_yy = evalOn(g_2D,(x,y) -> -cos(y)) 88 v_yy = evalOn(g_2D,(x,y) -> -cos(y))
87 89
88 # 2nd order interior stencil, 1st order boundary stencil, 90 # 2nd order interior stencil, 1st order boundary stencil,
89 # implies that L*v should be exact for binomials up to order 2. 91 # implies that L*v should be exact for binomials up to order 2.
90 @testset "2nd order" begin 92 @testset "2nd order" begin
91 stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=2) 93 stencil_set = read_stencil_set(operator_path; order=2)
92 inner_stencil = parse_stencil(stencil_set["D2"]["inner_stencil"]) 94 Dyy = second_derivative(g_2D,stencil_set,2)
93 closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"])
94 Dyy = second_derivative(g_2D,inner_stencil,closure_stencils,2)
95 @test Dyy*binomials[1] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9 95 @test Dyy*binomials[1] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9
96 @test Dyy*binomials[2] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9 96 @test Dyy*binomials[2] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9
97 @test Dyy*binomials[3] ≈ evalOn(g_2D,(x,y)->1.) atol = 5e-9 97 @test Dyy*binomials[3] ≈ evalOn(g_2D,(x,y)->1.) atol = 5e-9
98 @test Dyy*v ≈ v_yy rtol = 5e-2 norm = l2 98 @test Dyy*v ≈ v_yy rtol = 5e-2 norm = l2
99 end 99 end
100 100
101 # 4th order interior stencil, 2nd order boundary stencil, 101 # 4th order interior stencil, 2nd order boundary stencil,
102 # implies that L*v should be exact for binomials up to order 3. 102 # implies that L*v should be exact for binomials up to order 3.
103 @testset "4th order" begin 103 @testset "4th order" begin
104 stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4) 104 stencil_set = read_stencil_set(operator_path; order=4)
105 inner_stencil = parse_stencil(stencil_set["D2"]["inner_stencil"]) 105 Dyy = second_derivative(g_2D,stencil_set,2)
106 closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"])
107 Dyy = second_derivative(g_2D,inner_stencil,closure_stencils,2)
108 # NOTE: high tolerances for checking the "exact" differentiation 106 # NOTE: high tolerances for checking the "exact" differentiation
109 # due to accumulation of round-off errors/cancellation errors? 107 # due to accumulation of round-off errors/cancellation errors?
110 @test Dyy*binomials[1] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9 108 @test Dyy*binomials[1] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9
111 @test Dyy*binomials[2] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9 109 @test Dyy*binomials[2] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9
112 @test Dyy*binomials[3] ≈ evalOn(g_2D,(x,y)->1.) atol = 5e-9 110 @test Dyy*binomials[3] ≈ evalOn(g_2D,(x,y)->1.) atol = 5e-9