Mercurial > repos > public > sbplib_julia

comparison test/BoundaryConditions/sat_test.jl @ 1479:b96858a50e35 feature/boundary_conditions

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add tests for SAT
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Sat, 23 Dec 2023 23:01:28 +0100
parents bdcdbd4ea9cd
children 329720b9ba0d
comparison
equal deleted inserted replaced
1478:fefc81a9c155 1479:b96858a50e35
1 using Test 1 using Test
2 2
3 3
4 using Sbplib.BoundaryConditions 4 using Sbplib.BoundaryConditions
5 using Sbplib.Grids 5 using Sbplib.Grids
6 using Sbplib.RegionIndices
7 using Sbplib.LazyTensors 6 using Sbplib.LazyTensors
7 using Sbplib.SbpOperators
8 8
9 grid = equidistant_grid(11, 0.0, 1.0) 9 stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order = 4)
10 (id_l,id_r) = boundary_identifiers(grid) 10
11 struct MockOp 11 struct MockOp
12 end 12 end
13 13
14 function BoundaryConditions.sat_tensors(op::MockOp, grid, bc::DirichletCondition) 14 function BoundaryConditions.sat_tensors(op::MockOp, g::TensorGrid, bc::DirichletCondition)
15 sz = size(grid) 15 e = boundary_restriction(g, stencil_set, id(bc))
16 m = sz[1] 16 L = e
17 ind = (region(bc.id) == Lower()) ? 1 : m 17 sat_op = e'
18 e = zeros(m); 18 return sat_op, L
19 e[ind] = 1.
20 eᵀ = ones(Float64,m,0);
21 e[ind] = 1.
22 c_tensor = LazyTensors.DiagonalTensor(e)
23 p_tensor = DenseTensor(eᵀ, (1,), (2,))
24 closure(u) = c_tensor*u
25 function penalty(g)
26 @show g
27 return p_tensor*g
28 end
29 return closure, penalty
30 end 19 end
31 20
21 function BoundaryConditions.sat_tensors(op::MockOp, g::TensorGrid, bc::DirichletCondition, a)
22 e = boundary_restriction(g, stencil_set, id(bc))
23 L = a*e
24 sat_op = e'
25 return sat_op, L
26 end
32 27
33 # @testset "sat" begin 28 function BoundaryConditions.sat_tensors(op::MockOp, g::TensorGrid, bc::NeumannCondition)
34 # g = ConstantBoundaryData(2.0) 29 e = boundary_restriction(g, stencil_set, id(bc))
35 # dc = DirichletCondition(g,id_l) 30 d = normal_derivative(g, stencil_set, id(bc))
36 # op = MockOp() 31 L = d
37 # f = sat(op, grid, dc) 32 sat_op = e'
38 # u = eval_on(grid, x-> -1/2 + x^2) 33 return sat_op, L
39 # @show f(0.,u) 34 end
40 # end 35
36 @testset "sat" begin
37 op = MockOp()
38 grid = equidistant_grid((11,13), (0.,0.), (1.,1.))
39 W, E, S, N = boundary_identifiers(grid)
40 u = eval_on(grid, (x,y) -> x+y^2)
41
42
43 dc_W = DirichletCondition(1.0, W)
44 SAT_W = sat(op, grid, dc_W)
45 g_W = discretize_data(grid, dc_W)
46 r_W = zeros(size(grid))
47 r_W[1,:] .= map(y -> (y^2-1.), range(0., 1., length=13))
48 @test SAT_W(u, g_W) ≈ r_W atol = 1e-13
49
50 dc_E = DirichletCondition(2, E)
51 SAT_E = sat(op, grid, dc_E, 2.)
52 g_E = discretize_data(grid, dc_E)
53 r_E = zeros(size(grid))
54 r_E[end,:] .= map(y -> (2*(1. + y^2)-2.), range(0., 1., length=13))
55 @test SAT_E(u, g_E) ≈ r_E atol = 1e-13
56
57 nc_S = NeumannCondition(.0, S)
58 SAT_S = sat(op, grid, nc_S)
59 g_S = discretize_data(grid, nc_S)
60 @test SAT_S(u, g_S) ≈ zeros(size(grid)) atol = 1e-13
61
62 nc_N = NeumannCondition(2.0, N)
63 SAT_S = sat(op, grid, nc_N)
64 g_N = discretize_data(grid, nc_N)
65 @test SAT_S(u, g_N) ≈ zeros(size(grid)) atol = 1e-13
66 end