Mercurial > repos > public > sbplib_julia

comparison test/SbpOperators/volumeops/laplace/laplace_test.jl @ 1858:4a9be96f2569 feature/documenter_logo

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 Sun, 12 Jan 2025 21:18:44 +0100
parents 471a948cd2b2
children f3d7e2d7a43f
comparison
equal deleted inserted replaced
1857:ffde7dad9da5 1858:4a9be96f2569
1 using Test 1 using Test
2 2
3 using Sbplib.SbpOperators 3 using Diffinitive.SbpOperators
4 using Sbplib.Grids 4 using Diffinitive.Grids
5 using Sbplib.LazyTensors 5 using Diffinitive.LazyTensors
6 6
7 @testset "Laplace" begin 7 @testset "Laplace" begin
8 g_1D = EquidistantGrid(101, 0.0, 1.) 8 # Default stencils (4th order)
9 g_3D = EquidistantGrid((51,101,52), (0.0, -1.0, 0.0), (1., 1., 1.)) 9 operator_path = sbp_operators_path()*"standard_diagonal.toml"
10 stencil_set = read_stencil_set(operator_path; order=4)
11 g_1D = equidistant_grid(0.0, 1., 101)
12 g_3D = equidistant_grid((0.0, -1.0, 0.0), (1., 1., 1.), 51, 101, 52)
13
10 @testset "Constructors" begin 14 @testset "Constructors" begin
11 stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4)
12 inner_stencil = parse_stencil(stencil_set["D2"]["inner_stencil"])
13 closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"])
14 @testset "1D" begin 15 @testset "1D" begin
15 L = laplace(g_1D, inner_stencil, closure_stencils) 16 @test Laplace(g_1D, stencil_set) == Laplace(laplace(g_1D, stencil_set), stencil_set)
16 @test L == second_derivative(g_1D, inner_stencil, closure_stencils) 17 @test Laplace(g_1D, stencil_set) isa LazyTensor{Float64,1,1}
17 @test L isa TensorMapping{T,1,1} where T
18 end 18 end
19 @testset "3D" begin 19 @testset "3D" begin
20 L = laplace(g_3D, inner_stencil, closure_stencils) 20 @test Laplace(g_3D, stencil_set) == Laplace(laplace(g_3D, stencil_set),stencil_set)
21 @test L isa TensorMapping{T,3,3} where T 21 @test Laplace(g_3D, stencil_set) isa LazyTensor{Float64,3,3}
22 Dxx = second_derivative(g_3D, inner_stencil, closure_stencils, 1)
23 Dyy = second_derivative(g_3D, inner_stencil, closure_stencils, 2)
24 Dzz = second_derivative(g_3D, inner_stencil, closure_stencils, 3)
25 @test L == Dxx + Dyy + Dzz
26 end 22 end
27 end 23 end
28 24
29 # Exact differentiation is measured point-wise. In other cases 25 # Exact differentiation is measured point-wise. In other cases
30 # the error is measured in the l2-norm. 26 # the error is measured in the l2-norm.
31 @testset "Accuracy" begin 27 @testset "Accuracy" begin
32 l2(v) = sqrt(prod(spacing(g_3D))*sum(v.^2)); 28 l2(v) = sqrt(prod(spacing.(g_3D.grids))*sum(v.^2));
33 polynomials = () 29 polynomials = ()
34 maxOrder = 4; 30 maxOrder = 4;
35 for i = 0:maxOrder-1 31 for i = 0:maxOrder-1
36 f_i(x,y,z) = 1/factorial(i)*(y^i + x^i + z^i) 32 f_i(x,y,z) = 1/factorial(i)*(y^i + x^i + z^i)
37 polynomials = (polynomials...,evalOn(g_3D,f_i)) 33 polynomials = (polynomials...,eval_on(g_3D,f_i))
38 end 34 end
39 v = evalOn(g_3D, (x,y,z) -> sin(x) + cos(y) + exp(z)) 35 # v = eval_on(g_3D, (x,y,z) -> sin(x) + cos(y) + exp(z))
40 Δv = evalOn(g_3D,(x,y,z) -> -sin(x) - cos(y) + exp(z)) 36 # Δv = eval_on(g_3D,(x,y,z) -> -sin(x) - cos(y) + exp(z))
37
38 v = eval_on(g_3D, x̄ -> sin(x̄[1]) + cos(x̄[2]) + exp(x̄[3]))
39 Δv = eval_on(g_3D, x̄ -> -sin(x̄[1]) - cos(x̄[2]) + exp(x̄[3]))
40 @inferred v[1,2,3]
41 41
42 # 2nd order interior stencil, 1st order boundary stencil, 42 # 2nd order interior stencil, 1st order boundary stencil,
43 # implies that L*v should be exact for binomials up to order 2. 43 # implies that L*v should be exact for binomials up to order 2.
44 @testset "2nd order" begin 44 @testset "2nd order" begin
45 stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=2) 45 stencil_set = read_stencil_set(operator_path; order=2)
46 inner_stencil = parse_stencil(stencil_set["D2"]["inner_stencil"]) 46 Δ = Laplace(g_3D, stencil_set)
47 closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"]) 47 @test Δ*polynomials[1] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
48 L = laplace(g_3D, inner_stencil, closure_stencils) 48 @test Δ*polynomials[2] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
49 @test L*polynomials[1] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9 49 @test Δ*polynomials[3] ≈ polynomials[1] atol = 5e-9
50 @test L*polynomials[2] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9 50 @test Δ*v ≈ Δv rtol = 5e-2 norm = l2
51 @test L*polynomials[3] ≈ polynomials[1] atol = 5e-9
52 @test L*v ≈ Δv rtol = 5e-2 norm = l2
53 end 51 end
54 52
55 # 4th order interior stencil, 2nd order boundary stencil, 53 # 4th order interior stencil, 2nd order boundary stencil,
56 # implies that L*v should be exact for binomials up to order 3. 54 # implies that L*v should be exact for binomials up to order 3.
57 @testset "4th order" begin 55 @testset "4th order" begin
58 stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4) 56 stencil_set = read_stencil_set(operator_path; order=4)
59 inner_stencil = parse_stencil(stencil_set["D2"]["inner_stencil"]) 57 Δ = Laplace(g_3D, stencil_set)
60 closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"])
61 L = laplace(g_3D, inner_stencil, closure_stencils)
62 # NOTE: high tolerances for checking the "exact" differentiation 58 # NOTE: high tolerances for checking the "exact" differentiation
63 # due to accumulation of round-off errors/cancellation errors? 59 # due to accumulation of round-off errors/cancellation errors?
64 @test L*polynomials[1] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9 60 @test Δ*polynomials[1] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
65 @test L*polynomials[2] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9 61 @test Δ*polynomials[2] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
66 @test L*polynomials[3] ≈ polynomials[1] atol = 5e-9 62 @test Δ*polynomials[3] ≈ polynomials[1] atol = 5e-9
67 @test L*polynomials[4] ≈ polynomials[2] atol = 5e-9 63 @test Δ*polynomials[4] ≈ polynomials[2] atol = 5e-9
68 @test L*v ≈ Δv rtol = 5e-4 norm = l2 64 @test Δ*v ≈ Δv rtol = 5e-4 norm = l2
69 end 65 end
70 end 66 end
71 end 67 end
68
69 @testset "laplace" begin
70 operator_path = sbp_operators_path()*"standard_diagonal.toml"
71 stencil_set = read_stencil_set(operator_path; order=4)
72 g_1D = equidistant_grid(0.0, 1., 101)
73 g_3D = equidistant_grid((0.0, -1.0, 0.0), (1., 1., 1.), 51, 101, 52)
74
75 @testset "1D" begin
76 Δ = laplace(g_1D, stencil_set)
77 @test Δ == second_derivative(g_1D, stencil_set)
78 @test Δ isa LazyTensor{Float64,1,1}
79 end
80 @testset "3D" begin
81 Δ = laplace(g_3D, stencil_set)
82 @test Δ isa LazyTensor{Float64,3,3}
83 Dxx = second_derivative(g_3D, stencil_set, 1)
84 Dyy = second_derivative(g_3D, stencil_set, 2)
85 Dzz = second_derivative(g_3D, stencil_set, 3)
86 @test Δ == Dxx + Dyy + Dzz
87 @test Δ isa LazyTensor{Float64,3,3}
88 end
89 end
90
91 @testset "sat_tensors" begin
92 # TODO: The following tests should be implemented
93 # 1. Symmetry D'H == H'D (test_broken below)
94 # 2. Test eigenvalues of and/or solution to Poisson
95 # 3. Test tuning of Dirichlet conditions
96 #
97 # These tests are likely easiest to implement once
98 # we have support for generating matrices from tensors.
99
100 operator_path = sbp_operators_path()*"standard_diagonal.toml"
101 orders = (2,4)
102 tols = (5e-2,5e-4)
103 sz = (201,401)
104 g = equidistant_grid((0.,0.), (1.,1.), sz...)
105
106 # Verify implementation of sat_tensors by testing accuracy and symmetry (TODO)
107 # of the operator D = Δ + SAT, where SAT is the tensor composition of the
108 # operators from sat_tensor. Note that SAT*u should approximate 0 for the
109 # conditions chosen.
110
111 @testset "Dirichlet" begin
112 for (o, tol) ∈ zip(orders,tols)
113 stencil_set = read_stencil_set(operator_path; order=o)
114 Δ = Laplace(g, stencil_set)
115 H = inner_product(g, stencil_set)
116 u = collect(eval_on(g, (x,y) -> sin(π*x)sin(2*π*y)))
117 Δu = collect(eval_on(g, (x,y) -> -5*π^2*sin(π*x)sin(2*π*y)))
118 D = Δ
119 for id ∈ boundary_identifiers(g)
120 D = D + foldl(∘, sat_tensors(Δ, g, DirichletCondition(0., id)))
121 end
122 e = D*u .- Δu
123 # Accuracy
124 @test sqrt(sum(H*e.^2)) ≈ 0 atol = tol
125 # Symmetry
126 r = randn(size(u))
127 @test_broken (D'∘H - H∘D)*r .≈ 0 atol = 1e-13 # TODO: Need to implement apply_transpose for D.
128 end
129 end
130
131 @testset "Neumann" begin
132 @testset "Dirichlet" begin
133 for (o, tol) ∈ zip(orders,tols)
134 stencil_set = read_stencil_set(operator_path; order=o)
135 Δ = Laplace(g, stencil_set)
136 H = inner_product(g, stencil_set)
137 u = collect(eval_on(g, (x,y) -> cos(π*x)cos(2*π*y)))
138 Δu = collect(eval_on(g, (x,y) -> -5*π^2*cos(π*x)cos(2*π*y)))
139 D = Δ
140 for id ∈ boundary_identifiers(g)
141 D = D + foldl(∘, sat_tensors(Δ, g, NeumannCondition(0., id)))
142 end
143 e = D*u .- Δu
144 # Accuracy
145 @test sqrt(sum(H*e.^2)) ≈ 0 atol = tol
146 # Symmetry
147 r = randn(size(u))
148 @test_broken (D'∘H - H∘D)*r .≈ 0 atol = 1e-13 # TODO: Need to implement apply_transpose for D.
149 end
150 end
151 end
152 end
153