Mercurial > repos > public > sbplib_julia

comparison test/SbpOperators/boundaryops/boundary_restriction_test.jl @ 1652:65b2d2c72fbc feature/sbp_operators/laplace_curvilinear

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add boundary restriction operator for mapped grid
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 26 Jun 2024 12:54:29 +0200
parents 43aaf710463e
children f3d7e2d7a43f
comparison
equal deleted inserted replaced
1651:707fc9761c2b 1652:65b2d2c72fbc
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.RegionIndices 6 using Sbplib.RegionIndices
7 using Sbplib.SbpOperators: BoundaryOperator, Stencil 7 using Sbplib.SbpOperators: BoundaryOperator, Stencil
8
9 using StaticArrays
8 10
9 @testset "boundary_restriction" begin 11 @testset "boundary_restriction" begin
10 stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order = 4) 12 stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order = 4)
11 e_closure = parse_stencil(stencil_set["e"]["closure"]) 13 e_closure = parse_stencil(stencil_set["e"]["closure"])
12 g_1D = equidistant_grid(0.0, 1.0, 11) 14 g_1D = equidistant_grid(0.0, 1.0, 11)
31 @test e_w isa LazyTensor{T,1,2} where T 33 @test e_w isa LazyTensor{T,1,2} where T
32 end 34 end
33 end 35 end
34 36
35 @testset "Application" begin 37 @testset "Application" begin
36 @testset "1D" begin 38 @testset "EquidistantGrid" begin
37 e_l, e_r = boundary_restriction.(Ref(g_1D), Ref(stencil_set), boundary_identifiers(g_1D)) 39 e_l, e_r = boundary_restriction.(Ref(g_1D), Ref(stencil_set), boundary_identifiers(g_1D))
38 v = eval_on(g_1D,x->1+x^2) 40 v = eval_on(g_1D,x->1+x^2)
39 u = fill(3.124) 41 u = fill(3.124)
40 42
41 @test (e_l*v)[] == v[1] 43 @test (e_l*v)[] == v[1]
42 @test (e_r*v)[] == v[end] 44 @test (e_r*v)[] == v[end]
43 @test (e_r*v)[1] == v[end] 45 @test (e_r*v)[1] == v[end]
44 end 46 end
45 47
46 @testset "2D" begin 48 @testset "TensorGrid" begin
47 e_w, e_e, e_s, e_n = boundary_restriction.(Ref(g_2D), Ref(stencil_set), boundary_identifiers(g_2D)) 49 e_w, e_e, e_s, e_n = boundary_restriction.(Ref(g_2D), Ref(stencil_set), boundary_identifiers(g_2D))
48 v = rand(11, 15) 50 v = rand(11, 15)
49 u = fill(3.124) 51 u = fill(3.124)
50 52
51 @test e_w*v == v[1,:] 53 @test e_w*v == v[1,:]
52 @test e_e*v == v[end,:] 54 @test e_e*v == v[end,:]
53 @test e_s*v == v[:,1] 55 @test e_s*v == v[:,1]
54 @test e_n*v == v[:,end] 56 @test e_n*v == v[:,end]
55 end 57 end
58
59 @testset "MappedGrid" begin
60 c = Chart(unitsquare()) do (ξ,η)
61 @SVector[2ξ + η*(1-η), 3η+(1+η/2)*ξ^2]
62 end
63 Grids.jacobian(c::typeof(c), (ξ,η)) = @SMatrix[2 1-2η; (2+η)*ξ 3+ξ^2/2]
64
65 mg = equidistant_grid(c, 10,13)
66
67 e_w, e_e, e_s, e_n = boundary_restriction.(Ref(mg), Ref(stencil_set), boundary_identifiers(mg))
68 v = rand(10, 13)
69
70 @test e_w*v == v[1,:]
71 @test e_e*v == v[end,:]
72 @test e_s*v == v[:,1]
73 @test e_n*v == v[:,end]
74 end
56 end 75 end
57 end 76 end