Mercurial > repos > public > sbplib_julia

annotate test/SbpOperators/volumeops/laplace/laplace_test.jl @ 2055:274f4c1ce4b5 feature/sbp_operators/laplace_curvilinear tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add test for SBP property of laplace in mapped grid
author Jonatan Werpers <jonatan@werpers.com>
date Sun, 08 Feb 2026 00:05:40 +0100
parents f3d7e2d7a43f
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
728
45966c77cb20 Split tests for SbpOperators over several files
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
1 using Test
45966c77cb20 Split tests for SbpOperators over several files
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
2
1726
471a948cd2b2 Rename project from Sbplib to Diffinitive
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1615
diff changeset
3 using Diffinitive.SbpOperators
471a948cd2b2 Rename project from Sbplib to Diffinitive
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1615
diff changeset
4 using Diffinitive.Grids
471a948cd2b2 Rename project from Sbplib to Diffinitive
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1615
diff changeset
5 using Diffinitive.LazyTensors
728
45966c77cb20 Split tests for SbpOperators over several files
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
6
1588
f6774e98d223 Add basic test
Jonatan Werpers <jonatan@werpers.com>
parents: 1559
diff changeset
7 using StaticArrays
2055
274f4c1ce4b5 Add test for SBP property of laplace in mapped grid
Jonatan Werpers <jonatan@werpers.com>
parents: 1751
diff changeset
8 using SparseArrays
274f4c1ce4b5 Add test for SBP property of laplace in mapped grid
Jonatan Werpers <jonatan@werpers.com>
parents: 1751
diff changeset
9 using Tokens
274f4c1ce4b5 Add test for SBP property of laplace in mapped grid
Jonatan Werpers <jonatan@werpers.com>
parents: 1751
diff changeset
10 using LinearAlgebra
1588
f6774e98d223 Add basic test
Jonatan Werpers <jonatan@werpers.com>
parents: 1559
diff changeset
11
1291
356ec6a72974 Implement changes in SbpOperators
Jonatan Werpers <jonatan@werpers.com>
parents: 1285
diff changeset
12 @testset "Laplace" begin
1285
7d52c4835d15 Skip broken testsets
Jonatan Werpers <jonatan@werpers.com>
parents: 1040
diff changeset
13 # Default stencils (4th order)
7d52c4835d15 Skip broken testsets
Jonatan Werpers <jonatan@werpers.com>
parents: 1040
diff changeset
14 operator_path = sbp_operators_path()*"standard_diagonal.toml"
7d52c4835d15 Skip broken testsets
Jonatan Werpers <jonatan@werpers.com>
parents: 1040
diff changeset
15 stencil_set = read_stencil_set(operator_path; order=4)
1529
43aaf710463e Change to signature of equidistant_grid to same style as many array methods.
Jonatan Werpers <jonatan@werpers.com>
parents: 1291
diff changeset
16 g_1D = equidistant_grid(0.0, 1., 101)
43aaf710463e Change to signature of equidistant_grid to same style as many array methods.
Jonatan Werpers <jonatan@werpers.com>
parents: 1291
diff changeset
17 g_3D = equidistant_grid((0.0, -1.0, 0.0), (1., 1., 1.), 51, 101, 52)
728
45966c77cb20 Split tests for SbpOperators over several files
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
18
45966c77cb20 Split tests for SbpOperators over several files
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
19 @testset "Constructors" begin
750
f88b2117dc69 Merge in default
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 732
diff changeset
20 @testset "1D" begin
1291
356ec6a72974 Implement changes in SbpOperators
Jonatan Werpers <jonatan@werpers.com>
parents: 1285
diff changeset
21 @test Laplace(g_1D, stencil_set) == Laplace(laplace(g_1D, stencil_set), stencil_set)
356ec6a72974 Implement changes in SbpOperators
Jonatan Werpers <jonatan@werpers.com>
parents: 1285
diff changeset
22 @test Laplace(g_1D, stencil_set) isa LazyTensor{Float64,1,1}
750
f88b2117dc69 Merge in default
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 732
diff changeset
23 end
f88b2117dc69 Merge in default
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 732
diff changeset
24 @testset "3D" begin
1291
356ec6a72974 Implement changes in SbpOperators
Jonatan Werpers <jonatan@werpers.com>
parents: 1285
diff changeset
25 @test Laplace(g_3D, stencil_set) == Laplace(laplace(g_3D, stencil_set),stencil_set)
356ec6a72974 Implement changes in SbpOperators
Jonatan Werpers <jonatan@werpers.com>
parents: 1285
diff changeset
26 @test Laplace(g_3D, stencil_set) isa LazyTensor{Float64,3,3}
728
45966c77cb20 Split tests for SbpOperators over several files
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
27 end
45966c77cb20 Split tests for SbpOperators over several files
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
28 end
45966c77cb20 Split tests for SbpOperators over several files
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
29
45966c77cb20 Split tests for SbpOperators over several files
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
30 # Exact differentiation is measured point-wise. In other cases
45966c77cb20 Split tests for SbpOperators over several files
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
31 # the error is measured in the l2-norm.
45966c77cb20 Split tests for SbpOperators over several files
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
32 @testset "Accuracy" begin
1291
356ec6a72974 Implement changes in SbpOperators
Jonatan Werpers <jonatan@werpers.com>
parents: 1285
diff changeset
33 l2(v) = sqrt(prod(spacing.(g_3D.grids))*sum(v.^2));
728
45966c77cb20 Split tests for SbpOperators over several files
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
34 polynomials = ()
45966c77cb20 Split tests for SbpOperators over several files
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
35 maxOrder = 4;
45966c77cb20 Split tests for SbpOperators over several files
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
36 for i = 0:maxOrder-1
45966c77cb20 Split tests for SbpOperators over several files
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
37 f_i(x,y,z) = 1/factorial(i)*(y^i + x^i + z^i)
1291
356ec6a72974 Implement changes in SbpOperators
Jonatan Werpers <jonatan@werpers.com>
parents: 1285
diff changeset
38 polynomials = (polynomials...,eval_on(g_3D,f_i))
728
45966c77cb20 Split tests for SbpOperators over several files
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
39 end
1291
356ec6a72974 Implement changes in SbpOperators
Jonatan Werpers <jonatan@werpers.com>
parents: 1285
diff changeset
40 # v = eval_on(g_3D, (x,y,z) -> sin(x) + cos(y) + exp(z))
356ec6a72974 Implement changes in SbpOperators
Jonatan Werpers <jonatan@werpers.com>
parents: 1285
diff changeset
41 # Δv = eval_on(g_3D,(x,y,z) -> -sin(x) - cos(y) + exp(z))
356ec6a72974 Implement changes in SbpOperators
Jonatan Werpers <jonatan@werpers.com>
parents: 1285
diff changeset
42
356ec6a72974 Implement changes in SbpOperators
Jonatan Werpers <jonatan@werpers.com>
parents: 1285
diff changeset
43 v = eval_on(g_3D, x̄ -> sin(x̄[1]) + cos(x̄[2]) + exp(x̄[3]))
356ec6a72974 Implement changes in SbpOperators
Jonatan Werpers <jonatan@werpers.com>
parents: 1285
diff changeset
44 Δv = eval_on(g_3D, x̄ -> -sin(x̄[1]) - cos(x̄[2]) + exp(x̄[3]))
356ec6a72974 Implement changes in SbpOperators
Jonatan Werpers <jonatan@werpers.com>
parents: 1285
diff changeset
45 @inferred v[1,2,3]
728
45966c77cb20 Split tests for SbpOperators over several files
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
46
45966c77cb20 Split tests for SbpOperators over several files
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
47 # 2nd order interior stencil, 1st order boundary stencil,
45966c77cb20 Split tests for SbpOperators over several files
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
48 # implies that L*v should be exact for binomials up to order 2.
45966c77cb20 Split tests for SbpOperators over several files
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
49 @testset "2nd order" begin
1018
5ec49dd2c7c4 Reintroduce read_stencil_set
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 990
diff changeset
50 stencil_set = read_stencil_set(operator_path; order=2)
924
12e8e431b43c Start restructuring Laplace making it more minimal.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 872
diff changeset
51 Δ = Laplace(g_3D, stencil_set)
12e8e431b43c Start restructuring Laplace making it more minimal.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 872
diff changeset
52 @test Δ*polynomials[1] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
12e8e431b43c Start restructuring Laplace making it more minimal.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 872
diff changeset
53 @test Δ*polynomials[2] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
12e8e431b43c Start restructuring Laplace making it more minimal.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 872
diff changeset
54 @test Δ*polynomials[3] ≈ polynomials[1] atol = 5e-9
12e8e431b43c Start restructuring Laplace making it more minimal.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 872
diff changeset
55 @test Δ*v ≈ Δv rtol = 5e-2 norm = l2
728
45966c77cb20 Split tests for SbpOperators over several files
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
56 end
45966c77cb20 Split tests for SbpOperators over several files
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
57
45966c77cb20 Split tests for SbpOperators over several files
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
58 # 4th order interior stencil, 2nd order boundary stencil,
45966c77cb20 Split tests for SbpOperators over several files
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
59 # implies that L*v should be exact for binomials up to order 3.
45966c77cb20 Split tests for SbpOperators over several files
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
60 @testset "4th order" begin
1018
5ec49dd2c7c4 Reintroduce read_stencil_set
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 990
diff changeset
61 stencil_set = read_stencil_set(operator_path; order=4)
924
12e8e431b43c Start restructuring Laplace making it more minimal.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 872
diff changeset
62 Δ = Laplace(g_3D, stencil_set)
728
45966c77cb20 Split tests for SbpOperators over several files
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
63 # NOTE: high tolerances for checking the "exact" differentiation
45966c77cb20 Split tests for SbpOperators over several files
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
64 # due to accumulation of round-off errors/cancellation errors?
924
12e8e431b43c Start restructuring Laplace making it more minimal.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 872
diff changeset
65 @test Δ*polynomials[1] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
12e8e431b43c Start restructuring Laplace making it more minimal.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 872
diff changeset
66 @test Δ*polynomials[2] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
12e8e431b43c Start restructuring Laplace making it more minimal.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 872
diff changeset
67 @test Δ*polynomials[3] ≈ polynomials[1] atol = 5e-9
12e8e431b43c Start restructuring Laplace making it more minimal.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 872
diff changeset
68 @test Δ*polynomials[4] ≈ polynomials[2] atol = 5e-9
12e8e431b43c Start restructuring Laplace making it more minimal.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 872
diff changeset
69 @test Δ*v ≈ Δv rtol = 5e-4 norm = l2
728
45966c77cb20 Split tests for SbpOperators over several files
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
70 end
45966c77cb20 Split tests for SbpOperators over several files
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
71 end
45966c77cb20 Split tests for SbpOperators over several files
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
72 end
924
12e8e431b43c Start restructuring Laplace making it more minimal.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 872
diff changeset
73
1291
356ec6a72974 Implement changes in SbpOperators
Jonatan Werpers <jonatan@werpers.com>
parents: 1285
diff changeset
74 @testset "laplace" begin
1285
7d52c4835d15 Skip broken testsets
Jonatan Werpers <jonatan@werpers.com>
parents: 1040
diff changeset
75 operator_path = sbp_operators_path()*"standard_diagonal.toml"
7d52c4835d15 Skip broken testsets
Jonatan Werpers <jonatan@werpers.com>
parents: 1040
diff changeset
76 stencil_set = read_stencil_set(operator_path; order=4)
1529
43aaf710463e Change to signature of equidistant_grid to same style as many array methods.
Jonatan Werpers <jonatan@werpers.com>
parents: 1291
diff changeset
77 g_1D = equidistant_grid(0.0, 1., 101)
43aaf710463e Change to signature of equidistant_grid to same style as many array methods.
Jonatan Werpers <jonatan@werpers.com>
parents: 1291
diff changeset
78 g_3D = equidistant_grid((0.0, -1.0, 0.0), (1., 1., 1.), 51, 101, 52)
1285
7d52c4835d15 Skip broken testsets
Jonatan Werpers <jonatan@werpers.com>
parents: 1040
diff changeset
79
1559
43e6acbefdd1 Rename testsets and add one for mapped grids
Jonatan Werpers <jonatan@werpers.com>
parents: 1529
diff changeset
80 @testset "EquidistantGrid" begin
1291
356ec6a72974 Implement changes in SbpOperators
Jonatan Werpers <jonatan@werpers.com>
parents: 1285
diff changeset
81 Δ = laplace(g_1D, stencil_set)
356ec6a72974 Implement changes in SbpOperators
Jonatan Werpers <jonatan@werpers.com>
parents: 1285
diff changeset
82 @test Δ == second_derivative(g_1D, stencil_set)
356ec6a72974 Implement changes in SbpOperators
Jonatan Werpers <jonatan@werpers.com>
parents: 1285
diff changeset
83 @test Δ isa LazyTensor{Float64,1,1}
924
12e8e431b43c Start restructuring Laplace making it more minimal.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 872
diff changeset
84 end
1559
43e6acbefdd1 Rename testsets and add one for mapped grids
Jonatan Werpers <jonatan@werpers.com>
parents: 1529
diff changeset
85 @testset "TensorGrid" begin
1291
356ec6a72974 Implement changes in SbpOperators
Jonatan Werpers <jonatan@werpers.com>
parents: 1285
diff changeset
86 Δ = laplace(g_3D, stencil_set)
356ec6a72974 Implement changes in SbpOperators
Jonatan Werpers <jonatan@werpers.com>
parents: 1285
diff changeset
87 @test Δ isa LazyTensor{Float64,3,3}
356ec6a72974 Implement changes in SbpOperators
Jonatan Werpers <jonatan@werpers.com>
parents: 1285
diff changeset
88 Dxx = second_derivative(g_3D, stencil_set, 1)
356ec6a72974 Implement changes in SbpOperators
Jonatan Werpers <jonatan@werpers.com>
parents: 1285
diff changeset
89 Dyy = second_derivative(g_3D, stencil_set, 2)
356ec6a72974 Implement changes in SbpOperators
Jonatan Werpers <jonatan@werpers.com>
parents: 1285
diff changeset
90 Dzz = second_derivative(g_3D, stencil_set, 3)
924
12e8e431b43c Start restructuring Laplace making it more minimal.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 872
diff changeset
91 @test Δ == Dxx + Dyy + Dzz
1291
356ec6a72974 Implement changes in SbpOperators
Jonatan Werpers <jonatan@werpers.com>
parents: 1285
diff changeset
92 @test Δ isa LazyTensor{Float64,3,3}
924
12e8e431b43c Start restructuring Laplace making it more minimal.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 872
diff changeset
93 end
1559
43e6acbefdd1 Rename testsets and add one for mapped grids
Jonatan Werpers <jonatan@werpers.com>
parents: 1529
diff changeset
94
43e6acbefdd1 Rename testsets and add one for mapped grids
Jonatan Werpers <jonatan@werpers.com>
parents: 1529
diff changeset
95 @testset "MappedGrid" begin
1588
f6774e98d223 Add basic test
Jonatan Werpers <jonatan@werpers.com>
parents: 1559
diff changeset
96 c = Chart(unitsquare()) do (ξ,η)
f6774e98d223 Add basic test
Jonatan Werpers <jonatan@werpers.com>
parents: 1559
diff changeset
97 @SVector[2ξ + η*(1-η), 3η+(1+η/2)*ξ^2]
f6774e98d223 Add basic test
Jonatan Werpers <jonatan@werpers.com>
parents: 1559
diff changeset
98 end
1590
62cb622cbe6b Fix typo in test jacobian
Jonatan Werpers <jonatan@werpers.com>
parents: 1588
diff changeset
99 Grids.jacobian(c::typeof(c), (ξ,η)) = @SMatrix[2 1-2η; (2+η)*ξ 3+ξ^2/2]
1588
f6774e98d223 Add basic test
Jonatan Werpers <jonatan@werpers.com>
parents: 1559
diff changeset
100
1697
a63278c25c40 Fix test of laplace on mapped grid
Jonatan Werpers <jonatan@werpers.com>
parents: 1651
diff changeset
101 g = equidistant_grid(c, 60,60)
1588
f6774e98d223 Add basic test
Jonatan Werpers <jonatan@werpers.com>
parents: 1559
diff changeset
102
f6774e98d223 Add basic test
Jonatan Werpers <jonatan@werpers.com>
parents: 1559
diff changeset
103 @test laplace(g, stencil_set) isa LazyTensor{<:Any,2,2}
f6774e98d223 Add basic test
Jonatan Werpers <jonatan@werpers.com>
parents: 1559
diff changeset
104
1647
0685d97ebcb0 Start adding test for laplace on mapped grid
Jonatan Werpers <jonatan@werpers.com>
parents: 1590
diff changeset
105 f((x,y)) = sin(4(x + y))
1697
a63278c25c40 Fix test of laplace on mapped grid
Jonatan Werpers <jonatan@werpers.com>
parents: 1651
diff changeset
106 Δf((x,y)) = -32sin(4(x + y))
1647
0685d97ebcb0 Start adding test for laplace on mapped grid
Jonatan Werpers <jonatan@werpers.com>
parents: 1590
diff changeset
107 gf = map(f,g)
1588
f6774e98d223 Add basic test
Jonatan Werpers <jonatan@werpers.com>
parents: 1559
diff changeset
108
f6774e98d223 Add basic test
Jonatan Werpers <jonatan@werpers.com>
parents: 1559
diff changeset
109 Δ = laplace(g, stencil_set)
f6774e98d223 Add basic test
Jonatan Werpers <jonatan@werpers.com>
parents: 1559
diff changeset
110
f6774e98d223 Add basic test
Jonatan Werpers <jonatan@werpers.com>
parents: 1559
diff changeset
111 @test collect(Δ*gf) isa Array{<:Any,2}
1697
a63278c25c40 Fix test of laplace on mapped grid
Jonatan Werpers <jonatan@werpers.com>
parents: 1651
diff changeset
112 @test Δ*gf ≈ map(Δf, g) rtol=2e-2
2055
274f4c1ce4b5 Add test for SBP property of laplace in mapped grid
Jonatan Werpers <jonatan@werpers.com>
parents: 1751
diff changeset
113
274f4c1ce4b5 Add test for SBP property of laplace in mapped grid
Jonatan Werpers <jonatan@werpers.com>
parents: 1751
diff changeset
114
274f4c1ce4b5 Add test for SBP property of laplace in mapped grid
Jonatan Werpers <jonatan@werpers.com>
parents: 1751
diff changeset
115 @testset "SBP property" begin
274f4c1ce4b5 Add test for SBP property of laplace in mapped grid
Jonatan Werpers <jonatan@werpers.com>
parents: 1751
diff changeset
116 g = equidistant_grid(c, 20,20)
274f4c1ce4b5 Add test for SBP property of laplace in mapped grid
Jonatan Werpers <jonatan@werpers.com>
parents: 1751
diff changeset
117 Δ = laplace(g, stencil_set)
274f4c1ce4b5 Add test for SBP property of laplace in mapped grid
Jonatan Werpers <jonatan@werpers.com>
parents: 1751
diff changeset
118 H = inner_product(g, stencil_set)
274f4c1ce4b5 Add test for SBP property of laplace in mapped grid
Jonatan Werpers <jonatan@werpers.com>
parents: 1751
diff changeset
119 es = map(boundary_identifiers(g)) do id
274f4c1ce4b5 Add test for SBP property of laplace in mapped grid
Jonatan Werpers <jonatan@werpers.com>
parents: 1751
diff changeset
120 boundary_restriction(g, stencil_set, id)
274f4c1ce4b5 Add test for SBP property of laplace in mapped grid
Jonatan Werpers <jonatan@werpers.com>
parents: 1751
diff changeset
121 end
274f4c1ce4b5 Add test for SBP property of laplace in mapped grid
Jonatan Werpers <jonatan@werpers.com>
parents: 1751
diff changeset
122 ds = map(boundary_identifiers(g)) do id
274f4c1ce4b5 Add test for SBP property of laplace in mapped grid
Jonatan Werpers <jonatan@werpers.com>
parents: 1751
diff changeset
123 normal_derivative(g, stencil_set, id)
274f4c1ce4b5 Add test for SBP property of laplace in mapped grid
Jonatan Werpers <jonatan@werpers.com>
parents: 1751
diff changeset
124 end
274f4c1ce4b5 Add test for SBP property of laplace in mapped grid
Jonatan Werpers <jonatan@werpers.com>
parents: 1751
diff changeset
125 Hᵧs = map(boundary_identifiers(g)) do id
274f4c1ce4b5 Add test for SBP property of laplace in mapped grid
Jonatan Werpers <jonatan@werpers.com>
parents: 1751
diff changeset
126 inner_product(boundary_grid(g, id), stencil_set)
274f4c1ce4b5 Add test for SBP property of laplace in mapped grid
Jonatan Werpers <jonatan@werpers.com>
parents: 1751
diff changeset
127 end
274f4c1ce4b5 Add test for SBP property of laplace in mapped grid
Jonatan Werpers <jonatan@werpers.com>
parents: 1751
diff changeset
128
274f4c1ce4b5 Add test for SBP property of laplace in mapped grid
Jonatan Werpers <jonatan@werpers.com>
parents: 1751
diff changeset
129 BT = mapreduce(+, es, ds,Hᵧs) do e, d, Hᵧ
274f4c1ce4b5 Add test for SBP property of laplace in mapped grid
Jonatan Werpers <jonatan@werpers.com>
parents: 1751
diff changeset
130 e'∘Hᵧ∘d
274f4c1ce4b5 Add test for SBP property of laplace in mapped grid
Jonatan Werpers <jonatan@werpers.com>
parents: 1751
diff changeset
131 end
274f4c1ce4b5 Add test for SBP property of laplace in mapped grid
Jonatan Werpers <jonatan@werpers.com>
parents: 1751
diff changeset
132 M = -H∘Δ + BT
274f4c1ce4b5 Add test for SBP property of laplace in mapped grid
Jonatan Werpers <jonatan@werpers.com>
parents: 1751
diff changeset
133
274f4c1ce4b5 Add test for SBP property of laplace in mapped grid
Jonatan Werpers <jonatan@werpers.com>
parents: 1751
diff changeset
134 M = sparse(M)
274f4c1ce4b5 Add test for SBP property of laplace in mapped grid
Jonatan Werpers <jonatan@werpers.com>
parents: 1751
diff changeset
135 @test M ≈ M'
274f4c1ce4b5 Add test for SBP property of laplace in mapped grid
Jonatan Werpers <jonatan@werpers.com>
parents: 1751
diff changeset
136
274f4c1ce4b5 Add test for SBP property of laplace in mapped grid
Jonatan Werpers <jonatan@werpers.com>
parents: 1751
diff changeset
137 function issemiposdef(A, tol=1e-8)
274f4c1ce4b5 Add test for SBP property of laplace in mapped grid
Jonatan Werpers <jonatan@werpers.com>
parents: 1751
diff changeset
138 return isposdef(A+tol*I)
274f4c1ce4b5 Add test for SBP property of laplace in mapped grid
Jonatan Werpers <jonatan@werpers.com>
parents: 1751
diff changeset
139 end
274f4c1ce4b5 Add test for SBP property of laplace in mapped grid
Jonatan Werpers <jonatan@werpers.com>
parents: 1751
diff changeset
140
274f4c1ce4b5 Add test for SBP property of laplace in mapped grid
Jonatan Werpers <jonatan@werpers.com>
parents: 1751
diff changeset
141 @test issemiposdef(Symmetric(M))
274f4c1ce4b5 Add test for SBP property of laplace in mapped grid
Jonatan Werpers <jonatan@werpers.com>
parents: 1751
diff changeset
142 end
1559
43e6acbefdd1 Rename testsets and add one for mapped grids
Jonatan Werpers <jonatan@werpers.com>
parents: 1529
diff changeset
143 end
924
12e8e431b43c Start restructuring Laplace making it more minimal.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 872
diff changeset
144 end
12e8e431b43c Start restructuring Laplace making it more minimal.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 872
diff changeset
145
1485
e96ee7d7ac9c Test sat_tensors for Laplace w. Neumann conditions
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1291
diff changeset
146 @testset "sat_tensors" begin
1615
b74e1a21265f Add todos to test
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1613
diff changeset
147 # TODO: The following tests should be implemented
b74e1a21265f Add todos to test
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1613
diff changeset
148 # 1. Symmetry D'H == H'D (test_broken below)
b74e1a21265f Add todos to test
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1613
diff changeset
149 # 2. Test eigenvalues of and/or solution to Poisson
b74e1a21265f Add todos to test
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1613
diff changeset
150 # 3. Test tuning of Dirichlet conditions
b74e1a21265f Add todos to test
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1613
diff changeset
151 #
b74e1a21265f Add todos to test
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1613
diff changeset
152 # These tests are likely easiest to implement once
b74e1a21265f Add todos to test
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1613
diff changeset
153 # we have support for generating matrices from tensors.
b74e1a21265f Add todos to test
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1613
diff changeset
154
1485
e96ee7d7ac9c Test sat_tensors for Laplace w. Neumann conditions
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1291
diff changeset
155 operator_path = sbp_operators_path()*"standard_diagonal.toml"
1598
19cdec9c21cb Implement and test sat_tensors for Dirichlet and Neumann conditions
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1594
diff changeset
156 orders = (2,4)
19cdec9c21cb Implement and test sat_tensors for Dirichlet and Neumann conditions
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1594
diff changeset
157 tols = (5e-2,5e-4)
19cdec9c21cb Implement and test sat_tensors for Dirichlet and Neumann conditions
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1594
diff changeset
158 sz = (201,401)
19cdec9c21cb Implement and test sat_tensors for Dirichlet and Neumann conditions
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1594
diff changeset
159 g = equidistant_grid((0.,0.), (1.,1.), sz...)
1485
e96ee7d7ac9c Test sat_tensors for Laplace w. Neumann conditions
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1291
diff changeset
160
1615
b74e1a21265f Add todos to test
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1613
diff changeset
161 # Verify implementation of sat_tensors by testing accuracy and symmetry (TODO)
1598
19cdec9c21cb Implement and test sat_tensors for Dirichlet and Neumann conditions
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1594
diff changeset
162 # of the operator D = Δ + SAT, where SAT is the tensor composition of the
19cdec9c21cb Implement and test sat_tensors for Dirichlet and Neumann conditions
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1594
diff changeset
163 # operators from sat_tensor. Note that SAT*u should approximate 0 for the
19cdec9c21cb Implement and test sat_tensors for Dirichlet and Neumann conditions
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1594
diff changeset
164 # conditions chosen.
1485
e96ee7d7ac9c Test sat_tensors for Laplace w. Neumann conditions
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1291
diff changeset
165
1598
19cdec9c21cb Implement and test sat_tensors for Dirichlet and Neumann conditions
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1594
diff changeset
166 @testset "Dirichlet" begin
19cdec9c21cb Implement and test sat_tensors for Dirichlet and Neumann conditions
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1594
diff changeset
167 for (o, tol) ∈ zip(orders,tols)
19cdec9c21cb Implement and test sat_tensors for Dirichlet and Neumann conditions
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1594
diff changeset
168 stencil_set = read_stencil_set(operator_path; order=o)
19cdec9c21cb Implement and test sat_tensors for Dirichlet and Neumann conditions
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1594
diff changeset
169 Δ = Laplace(g, stencil_set)
19cdec9c21cb Implement and test sat_tensors for Dirichlet and Neumann conditions
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1594
diff changeset
170 H = inner_product(g, stencil_set)
19cdec9c21cb Implement and test sat_tensors for Dirichlet and Neumann conditions
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1594
diff changeset
171 u = collect(eval_on(g, (x,y) -> sin(π*x)sin(2*π*y)))
19cdec9c21cb Implement and test sat_tensors for Dirichlet and Neumann conditions
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1594
diff changeset
172 Δu = collect(eval_on(g, (x,y) -> -5*π^2*sin(π*x)sin(2*π*y)))
19cdec9c21cb Implement and test sat_tensors for Dirichlet and Neumann conditions
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1594
diff changeset
173 D = Δ
19cdec9c21cb Implement and test sat_tensors for Dirichlet and Neumann conditions
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1594
diff changeset
174 for id ∈ boundary_identifiers(g)
19cdec9c21cb Implement and test sat_tensors for Dirichlet and Neumann conditions
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1594
diff changeset
175 D = D + foldl(∘, sat_tensors(Δ, g, DirichletCondition(0., id)))
19cdec9c21cb Implement and test sat_tensors for Dirichlet and Neumann conditions
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1594
diff changeset
176 end
19cdec9c21cb Implement and test sat_tensors for Dirichlet and Neumann conditions
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1594
diff changeset
177 e = D*u .- Δu
19cdec9c21cb Implement and test sat_tensors for Dirichlet and Neumann conditions
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1594
diff changeset
178 # Accuracy
19cdec9c21cb Implement and test sat_tensors for Dirichlet and Neumann conditions
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1594
diff changeset
179 @test sqrt(sum(H*e.^2)) ≈ 0 atol = tol
19cdec9c21cb Implement and test sat_tensors for Dirichlet and Neumann conditions
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1594
diff changeset
180 # Symmetry
19cdec9c21cb Implement and test sat_tensors for Dirichlet and Neumann conditions
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1594
diff changeset
181 r = randn(size(u))
19cdec9c21cb Implement and test sat_tensors for Dirichlet and Neumann conditions
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1594
diff changeset
182 @test_broken (D'∘H - H∘D)*r .≈ 0 atol = 1e-13 # TODO: Need to implement apply_transpose for D.
19cdec9c21cb Implement and test sat_tensors for Dirichlet and Neumann conditions
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1594
diff changeset
183 end
19cdec9c21cb Implement and test sat_tensors for Dirichlet and Neumann conditions
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1594
diff changeset
184 end
1485
e96ee7d7ac9c Test sat_tensors for Laplace w. Neumann conditions
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1291
diff changeset
185
e96ee7d7ac9c Test sat_tensors for Laplace w. Neumann conditions
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1291
diff changeset
186 @testset "Neumann" begin
1598
19cdec9c21cb Implement and test sat_tensors for Dirichlet and Neumann conditions
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1594
diff changeset
187 @testset "Dirichlet" begin
19cdec9c21cb Implement and test sat_tensors for Dirichlet and Neumann conditions
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1594
diff changeset
188 for (o, tol) ∈ zip(orders,tols)
19cdec9c21cb Implement and test sat_tensors for Dirichlet and Neumann conditions
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1594
diff changeset
189 stencil_set = read_stencil_set(operator_path; order=o)
19cdec9c21cb Implement and test sat_tensors for Dirichlet and Neumann conditions
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1594
diff changeset
190 Δ = Laplace(g, stencil_set)
19cdec9c21cb Implement and test sat_tensors for Dirichlet and Neumann conditions
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1594
diff changeset
191 H = inner_product(g, stencil_set)
19cdec9c21cb Implement and test sat_tensors for Dirichlet and Neumann conditions
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1594
diff changeset
192 u = collect(eval_on(g, (x,y) -> cos(π*x)cos(2*π*y)))
19cdec9c21cb Implement and test sat_tensors for Dirichlet and Neumann conditions
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1594
diff changeset
193 Δu = collect(eval_on(g, (x,y) -> -5*π^2*cos(π*x)cos(2*π*y)))
1613
15488c889a50 Fix variable name
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1603
diff changeset
194 D = Δ
1598
19cdec9c21cb Implement and test sat_tensors for Dirichlet and Neumann conditions
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1594
diff changeset
195 for id ∈ boundary_identifiers(g)
1613
15488c889a50 Fix variable name
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1603
diff changeset
196 D = D + foldl(∘, sat_tensors(Δ, g, NeumannCondition(0., id)))
1598
19cdec9c21cb Implement and test sat_tensors for Dirichlet and Neumann conditions
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1594
diff changeset
197 end
1613
15488c889a50 Fix variable name
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1603
diff changeset
198 e = D*u .- Δu
1598
19cdec9c21cb Implement and test sat_tensors for Dirichlet and Neumann conditions
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1594
diff changeset
199 # Accuracy
19cdec9c21cb Implement and test sat_tensors for Dirichlet and Neumann conditions
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1594
diff changeset
200 @test sqrt(sum(H*e.^2)) ≈ 0 atol = tol
19cdec9c21cb Implement and test sat_tensors for Dirichlet and Neumann conditions
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1594
diff changeset
201 # Symmetry
19cdec9c21cb Implement and test sat_tensors for Dirichlet and Neumann conditions
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1594
diff changeset
202 r = randn(size(u))
19cdec9c21cb Implement and test sat_tensors for Dirichlet and Neumann conditions
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1594
diff changeset
203 @test_broken (D'∘H - H∘D)*r .≈ 0 atol = 1e-13 # TODO: Need to implement apply_transpose for D.
19cdec9c21cb Implement and test sat_tensors for Dirichlet and Neumann conditions
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1594
diff changeset
204 end
19cdec9c21cb Implement and test sat_tensors for Dirichlet and Neumann conditions
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1594
diff changeset
205 end
1485
e96ee7d7ac9c Test sat_tensors for Laplace w. Neumann conditions
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1291
diff changeset
206 end
e96ee7d7ac9c Test sat_tensors for Laplace w. Neumann conditions
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1291
diff changeset
207 end
e96ee7d7ac9c Test sat_tensors for Laplace w. Neumann conditions
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 1291
diff changeset
208