Mercurial > repos > public > sbplib_julia

comparison test/SbpOperators/volumeops/laplace/laplace_test.jl @ 1854:654a2b7e6824 tooling/benchmarks

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 Sat, 11 Jan 2025 10:19:47 +0100
parents 471a948cd2b2
children f3d7e2d7a43f
comparison
equal deleted inserted replaced
1378:2b5480e2d4bf 1854:654a2b7e6824
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 # Default stencils (4th order) 8 # Default stencils (4th order)
9 operator_path = sbp_operators_path()*"standard_diagonal.toml" 9 operator_path = sbp_operators_path()*"standard_diagonal.toml"
10 stencil_set = read_stencil_set(operator_path; order=4) 10 stencil_set = read_stencil_set(operator_path; order=4)
11 g_1D = equidistant_grid(101, 0.0, 1.) 11 g_1D = equidistant_grid(0.0, 1., 101)
12 g_3D = equidistant_grid((51,101,52), (0.0, -1.0, 0.0), (1., 1., 1.)) 12 g_3D = equidistant_grid((0.0, -1.0, 0.0), (1., 1., 1.), 51, 101, 52)
13 13
14 @testset "Constructors" begin 14 @testset "Constructors" begin
15 @testset "1D" begin 15 @testset "1D" begin
16 @test Laplace(g_1D, stencil_set) == Laplace(laplace(g_1D, stencil_set), stencil_set) 16 @test Laplace(g_1D, stencil_set) == Laplace(laplace(g_1D, stencil_set), stencil_set)
17 @test Laplace(g_1D, stencil_set) isa LazyTensor{Float64,1,1} 17 @test Laplace(g_1D, stencil_set) isa LazyTensor{Float64,1,1}
67 end 67 end
68 68
69 @testset "laplace" begin 69 @testset "laplace" begin
70 operator_path = sbp_operators_path()*"standard_diagonal.toml" 70 operator_path = sbp_operators_path()*"standard_diagonal.toml"
71 stencil_set = read_stencil_set(operator_path; order=4) 71 stencil_set = read_stencil_set(operator_path; order=4)
72 g_1D = equidistant_grid(101, 0.0, 1.) 72 g_1D = equidistant_grid(0.0, 1., 101)
73 g_3D = equidistant_grid((51,101,52), (0.0, -1.0, 0.0), (1., 1., 1.)) 73 g_3D = equidistant_grid((0.0, -1.0, 0.0), (1., 1., 1.), 51, 101, 52)
74 74
75 @testset "1D" begin 75 @testset "1D" begin
76 Δ = laplace(g_1D, stencil_set) 76 Δ = laplace(g_1D, stencil_set)
77 @test Δ == second_derivative(g_1D, stencil_set) 77 @test Δ == second_derivative(g_1D, stencil_set)
78 @test Δ isa LazyTensor{Float64,1,1} 78 @test Δ isa LazyTensor{Float64,1,1}
86 @test Δ == Dxx + Dyy + Dzz 86 @test Δ == Dxx + Dyy + Dzz
87 @test Δ isa LazyTensor{Float64,3,3} 87 @test Δ isa LazyTensor{Float64,3,3}
88 end 88 end
89 end 89 end
90 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