Mercurial > repos > public > sbplib_julia

comparison test/SbpOperators/volumeops/laplace/laplace_test.jl @ 1485:e96ee7d7ac9c feature/boundary_conditions

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Test sat_tensors for Laplace w. Neumann conditions
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Mon, 25 Dec 2023 23:39:56 +0100
parents 356ec6a72974
children d68d02dd882f
comparison
equal deleted inserted replaced
1484:8d60d045c2a2 1485:e96ee7d7ac9c
1 using Test 1 using Test
2 2
3 using Sbplib.SbpOperators 3 using Sbplib.SbpOperators
4 using Sbplib.Grids 4 using Sbplib.Grids
5 using Sbplib.LazyTensors 5 using Sbplib.LazyTensors
6 using Sbplib.BoundaryConditions
6 7
7 @testset "Laplace" begin 8 @testset "Laplace" begin
8 # Default stencils (4th order) 9 # Default stencils (4th order)
9 operator_path = sbp_operators_path()*"standard_diagonal.toml" 10 operator_path = sbp_operators_path()*"standard_diagonal.toml"
10 stencil_set = read_stencil_set(operator_path; order=4) 11 stencil_set = read_stencil_set(operator_path; order=4)
86 @test Δ == Dxx + Dyy + Dzz 87 @test Δ == Dxx + Dyy + Dzz
87 @test Δ isa LazyTensor{Float64,3,3} 88 @test Δ isa LazyTensor{Float64,3,3}
88 end 89 end
89 end 90 end
90 91
92 @testset "sat_tensors" begin
93 operator_path = sbp_operators_path()*"standard_diagonal.toml"
94 stencil_set = read_stencil_set(operator_path; order=4)
95 g = equidistant_grid((101,102), (-1.,-1.), (1.,1.))
96 W,E,S,N = boundary_identifiers(g)
97
98 u = eval_on(g, (x,y) -> sin(x+y))
99 uWx = eval_on(boundary_grid(g,W), (x,y) -> -cos(x+y))
100 uEx = eval_on(boundary_grid(g,E), (x,y) -> cos(x+y))
101 uSy = eval_on(boundary_grid(g,S), (x,y) -> -cos(x+y))
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
111
112 @testset "Neumann" begin
113 Δ = Laplace(g, stencil_set)
114 H = inner_product(g, stencil_set)
115 HW = inner_product(boundary_grid(g,W), stencil_set)
116 HE = inner_product(boundary_grid(g,E), stencil_set)
117 HS = inner_product(boundary_grid(g,S), stencil_set)
118 HN = inner_product(boundary_grid(g,N), stencil_set)
119
120 ncW = NeumannCondition(0., W)
121 ncE = NeumannCondition(0., E)
122 ncS = NeumannCondition(0., S)
123 ncN = NeumannCondition(0., N)
124
125 SATW = foldl(∘,sat_tensors(Δ, g, ncW))
126 SATE = foldl(∘,sat_tensors(Δ, g, ncE))
127 SATS = foldl(∘,sat_tensors(Δ, g, ncS))
128 SATN = foldl(∘,sat_tensors(Δ, g, ncN))
129
130
131 @test sum((H*SATW*u).*v) ≈ sum((HW*uWx).*vW) rtol = 1e-6
132 @test sum((H*SATE*u).*v) ≈ sum((HE*uEx).*vE) rtol = 1e-6
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
136 end
137