Mercurial > repos > public > sbplib_julia

comparison test/SbpOperators/volumeops/laplace/laplace_test.jl @ 1623:07e8080c9f1a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge feature/boundary_conditions
author Vidar Stiernström <vidar.stiernstrom@gmail.com>
date Mon, 10 Jun 2024 22:39:22 -0700
parents b74e1a21265f
children 707fc9761c2b 471a948cd2b2
comparison
equal deleted inserted replaced
1610:bb4c119e91fa 1623:07e8080c9f1a
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