Mercurial > repos > public > sbplib_julia

comparison test/SbpOperators/volumeops/laplace/laplace_test.jl @ 1598:19cdec9c21cb feature/boundary_conditions

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Implement and test sat_tensors for Dirichlet and Neumann conditions
author Vidar Stiernström <vidar.stiernstrom@gmail.com>
date Sun, 26 May 2024 18:19:02 -0700
parents d68d02dd882f
children fca4a01d60c9
comparison
equal deleted inserted replaced
1597:330c39505a94 1598:19cdec9c21cb
89 end 89 end
90 end 90 end
91 91
92 @testset "sat_tensors" begin 92 @testset "sat_tensors" begin
93 operator_path = sbp_operators_path()*"standard_diagonal.toml" 93 operator_path = sbp_operators_path()*"standard_diagonal.toml"
94 stencil_set = read_stencil_set(operator_path; order=4) 94 orders = (2,4)
95 g = equidistant_grid((101,102), (-1.,-1.), (1.,1.)) 95 tols = (5e-2,5e-4)
96 W,E,S,N = boundary_identifiers(g) 96 sz = (201,401)
97 g = equidistant_grid((0.,0.), (1.,1.), sz...)
97 98
98 u = eval_on(g, (x,y) -> sin(x+y)) 99 # Verify implementation of sat_tesnors by testing accuracy and symmetry (TODO)
99 uWx = eval_on(boundary_grid(g,W), (x,y) -> -cos(x+y)) 100 # of the operator D = Δ + SAT, where SAT is the tensor composition of the
100 uEx = eval_on(boundary_grid(g,E), (x,y) -> cos(x+y)) 101 # operators from sat_tensor. Note that SAT*u should approximate 0 for the
101 uSy = eval_on(boundary_grid(g,S), (x,y) -> -cos(x+y)) 102 # conditions chosen.
102 uNy = eval_on(boundary_grid(g,N), (x,y) -> cos(x+y))
103
104
105 v = eval_on(g, (x,y) -> cos(x+y))
106 vW = eval_on(boundary_grid(g,W), (x,y) -> cos(x+y))
107 vE = eval_on(boundary_grid(g,E), (x,y) -> cos(x+y))
108 vS = eval_on(boundary_grid(g,S), (x,y) -> cos(x+y))
109 vN = eval_on(boundary_grid(g,N), (x,y) -> cos(x+y))
110 103
104 @testset "Dirichlet" begin
105 for (o, tol) ∈ zip(orders,tols)
106 stencil_set = read_stencil_set(operator_path; order=o)
107 Δ = Laplace(g, stencil_set)
108 H = inner_product(g, stencil_set)
109 u = collect(eval_on(g, (x,y) -> sin(π*x)sin(2*π*y)))
110 Δu = collect(eval_on(g, (x,y) -> -5*π^2*sin(π*x)sin(2*π*y)))
111 D = Δ
112 for id ∈ boundary_identifiers(g)
113 D = D + foldl(∘, sat_tensors(Δ, g, DirichletCondition(0., id)))
114 end
115 e = D*u .- Δu
116 # Accuracy
117 @test sqrt(sum(H*e.^2)) ≈ 0 atol = tol
118 # Symmetry
119 # TODO: # Consider generating the matrices to H and D and test D'H == H'D
120 r = randn(size(u))
121 @test_broken (D'∘H - H∘D)*r .≈ 0 atol = 1e-13 # TODO: Need to implement apply_transpose for D.
122 end
123 end
111 124
112 @testset "Neumann" begin 125 @testset "Neumann" begin
113 Δ = Laplace(g, stencil_set) 126 @testset "Dirichlet" begin
114 H = inner_product(g, stencil_set) 127 for (o, tol) ∈ zip(orders,tols)
115 HW = inner_product(boundary_grid(g,W), stencil_set) 128 stencil_set = read_stencil_set(operator_path; order=o)
116 HE = inner_product(boundary_grid(g,E), stencil_set) 129 Δ = Laplace(g, stencil_set)
117 HS = inner_product(boundary_grid(g,S), stencil_set) 130 H = inner_product(g, stencil_set)
118 HN = inner_product(boundary_grid(g,N), stencil_set) 131 u = collect(eval_on(g, (x,y) -> cos(π*x)cos(2*π*y)))
119 132 Δu = collect(eval_on(g, (x,y) -> -5*π^2*cos(π*x)cos(2*π*y)))
120 ncW = NeumannCondition(0., W) 133 op = Δ
121 ncE = NeumannCondition(0., E) 134 for id ∈ boundary_identifiers(g)
122 ncS = NeumannCondition(0., S) 135 op = op + foldl(∘, sat_tensors(Δ, g, NeumannCondition(0., id)))
123 ncN = NeumannCondition(0., N) 136 end
124 137 e = op*u .- Δu
125 SATW = foldl(∘,sat_tensors(Δ, g, ncW)) 138 # Accuracy
126 SATE = foldl(∘,sat_tensors(Δ, g, ncE)) 139 @test sqrt(sum(H*e.^2)) ≈ 0 atol = tol
127 SATS = foldl(∘,sat_tensors(Δ, g, ncS)) 140 # Symmetry
128 SATN = foldl(∘,sat_tensors(Δ, g, ncN)) 141 # TODO: # Consider generating the matrices to H and D and test D'H == H'D
129 142 r = randn(size(u))
130 143 @test_broken (D'∘H - H∘D)*r .≈ 0 atol = 1e-13 # TODO: Need to implement apply_transpose for D.
131 @test sum((H*SATW*u).*v) ≈ sum((HW*uWx).*vW) rtol = 1e-6 144 end
132 @test sum((H*SATE*u).*v) ≈ sum((HE*uEx).*vE) rtol = 1e-6 145 end
133 @test sum((H*SATS*u).*v) ≈ sum((HS*uSy).*vS) rtol = 1e-6
134 @test sum((H*SATN*u).*v) ≈ sum((HN*uNy).*vN) rtol = 1e-6
135 end 146 end
136 end 147 end
137 148