Mercurial > repos > public > sbplib_julia

comparison test/SbpOperators/volumeops/laplace/laplace_test.jl @ 1651:707fc9761c2b feature/sbp_operators/laplace_curvilinear

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge feature/grids/manifolds
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 26 Jun 2024 12:47:26 +0200
parents 0685d97ebcb0 b74e1a21265f
children a63278c25c40
comparison
equal deleted inserted replaced
1648:b7dcd3dd3181 1651:707fc9761c2b
108 @test collect(Δ*gf) isa Array{<:Any,2} 108 @test collect(Δ*gf) isa Array{<:Any,2}
109 @test Δ*gf ≈ map(Δf, g) 109 @test Δ*gf ≈ map(Δf, g)
110 end 110 end
111 end 111 end
112 112
113 @testset "sat_tensors" begin
114 # TODO: The following tests should be implemented
115 # 1. Symmetry D'H == H'D (test_broken below)
116 # 2. Test eigenvalues of and/or solution to Poisson
117 # 3. Test tuning of Dirichlet conditions
118 #
119 # These tests are likely easiest to implement once
120 # we have support for generating matrices from tensors.
121
122 operator_path = sbp_operators_path()*"standard_diagonal.toml"
123 orders = (2,4)
124 tols = (5e-2,5e-4)
125 sz = (201,401)
126 g = equidistant_grid((0.,0.), (1.,1.), sz...)
127
128 # Verify implementation of sat_tensors by testing accuracy and symmetry (TODO)
129 # of the operator D = Δ + SAT, where SAT is the tensor composition of the
130 # operators from sat_tensor. Note that SAT*u should approximate 0 for the
131 # conditions chosen.
132
133 @testset "Dirichlet" begin
134 for (o, tol) ∈ zip(orders,tols)
135 stencil_set = read_stencil_set(operator_path; order=o)
136 Δ = Laplace(g, stencil_set)
137 H = inner_product(g, stencil_set)
138 u = collect(eval_on(g, (x,y) -> sin(π*x)sin(2*π*y)))
139 Δu = collect(eval_on(g, (x,y) -> -5*π^2*sin(π*x)sin(2*π*y)))
140 D = Δ
141 for id ∈ boundary_identifiers(g)
142 D = D + foldl(∘, sat_tensors(Δ, g, DirichletCondition(0., id)))
143 end
144 e = D*u .- Δu
145 # Accuracy
146 @test sqrt(sum(H*e.^2)) ≈ 0 atol = tol
147 # Symmetry
148 r = randn(size(u))
149 @test_broken (D'∘H - H∘D)*r .≈ 0 atol = 1e-13 # TODO: Need to implement apply_transpose for D.
150 end
151 end
152
153 @testset "Neumann" begin
154 @testset "Dirichlet" begin
155 for (o, tol) ∈ zip(orders,tols)
156 stencil_set = read_stencil_set(operator_path; order=o)
157 Δ = Laplace(g, stencil_set)
158 H = inner_product(g, stencil_set)
159 u = collect(eval_on(g, (x,y) -> cos(π*x)cos(2*π*y)))
160 Δu = collect(eval_on(g, (x,y) -> -5*π^2*cos(π*x)cos(2*π*y)))
161 D = Δ
162 for id ∈ boundary_identifiers(g)
163 D = D + foldl(∘, sat_tensors(Δ, g, NeumannCondition(0., id)))
164 end
165 e = D*u .- Δu
166 # Accuracy
167 @test sqrt(sum(H*e.^2)) ≈ 0 atol = tol
168 # Symmetry
169 r = randn(size(u))
170 @test_broken (D'∘H - H∘D)*r .≈ 0 atol = 1e-13 # TODO: Need to implement apply_transpose for D.
171 end
172 end
173 end
174 end
175