1 | from anuga.abstract_2d_finite_volumes.generic_domain import Generic_Domain |
---|
2 | from anuga.abstract_2d_finite_volumes.generic_boundary_conditions import Dirichlet_boundary |
---|
3 | from kv_solver_new import Kinematic_Viscosity_Operator |
---|
4 | import numpy as num |
---|
5 | from math import sqrt |
---|
6 | import unittest |
---|
7 | |
---|
8 | class Test_kv_solver_new(unittest.TestCase): |
---|
9 | def setUp(self): |
---|
10 | pass |
---|
11 | |
---|
12 | def tearDown(self): |
---|
13 | pass |
---|
14 | |
---|
15 | #First test operator class (1 triangle) |
---|
16 | def operator1(self): |
---|
17 | points = num.array([[0.0,0.0],[1.0,0.0],[0.0,1.0]]) |
---|
18 | elements = num.array([[0,1,2]]) |
---|
19 | boundary_map = {} |
---|
20 | boundary_map[(0,0)] = Dirichlet_boundary([1,2,1]) |
---|
21 | boundary_map[(0,1)] = Dirichlet_boundary([1,1,2]) |
---|
22 | boundary_map[(0,2)] = Dirichlet_boundary([1,1,0]) |
---|
23 | #We get boundary_may.keys() = [(0, 1), (0, 0), (0, 2)] |
---|
24 | domain = Generic_Domain(source=points,triangles=elements,boundary=boundary_map) |
---|
25 | return Kinematic_Viscosity_Operator(domain) |
---|
26 | |
---|
27 | #Second test operator class (2 triangles) |
---|
28 | def operator2(self): |
---|
29 | points = num.array([[0.0,0.0],[1.0,0.0],[1.0,1.0],[0.0,1.0]]) |
---|
30 | elements = num.array([[0,1,3],[1,2,3]]) |
---|
31 | boundary_map = {} |
---|
32 | boundary_map[(0,1)] = Dirichlet_boundary([1,1,2]) |
---|
33 | boundary_map[(0,2)] = Dirichlet_boundary([1,2,2]) |
---|
34 | boundary_map[(1,0)] = Dirichlet_boundary([1,1,0]) |
---|
35 | boundary_map[(1,2)] = Dirichlet_boundary([1,2,1]) |
---|
36 | #We get boundary_map.keys() = [(0, 1), (1, 2), (1, 0), (0, 2)] |
---|
37 | domain = Generic_Domain(source=points,triangles=elements,boundary=boundary_map) |
---|
38 | return Kinematic_Viscosity_Operator(domain) |
---|
39 | |
---|
40 | def test_build_index_boundary_map(self): |
---|
41 | operator1 = self.operator1() |
---|
42 | boundary_map_index = operator1.boundary_enum |
---|
43 | assert boundary_map_index[(0,1)] == 0 |
---|
44 | assert boundary_map_index[(0,0)] == 1 |
---|
45 | assert boundary_map_index[(0,2)] == 2 |
---|
46 | |
---|
47 | operator2 = self.operator2() |
---|
48 | boundary_map_index = operator2.boundary_enum |
---|
49 | assert boundary_map_index[(0,1)] == 0 |
---|
50 | assert boundary_map_index[(1,2)] == 1 |
---|
51 | assert boundary_map_index[(1,0)] == 2 |
---|
52 | assert boundary_map_index[(0,2)] == 3 |
---|
53 | |
---|
54 | def test_kv_system_geometry(self): |
---|
55 | operator1 = self.operator1() |
---|
56 | operator2 = self.operator2() |
---|
57 | |
---|
58 | kv_geo_structure = operator1.geo_structure |
---|
59 | assert num.allclose(kv_geo_structure[0].todense(), num.array([[-6.0,0.0,6.0,0.0]])) |
---|
60 | assert num.allclose(kv_geo_structure[1].todense(), \ |
---|
61 | num.array([[-6.0/sqrt(5.0),6.0/sqrt(5.0),0.0,0.0]])) |
---|
62 | assert num.allclose(kv_geo_structure[2].todense(),\ |
---|
63 | num.array([[-6.0/sqrt(5.0),0.0,0.0,6.0/sqrt(5.0)]])) |
---|
64 | |
---|
65 | kv_geo_structure = operator2.geo_structure |
---|
66 | assert num.allclose(kv_geo_structure[0].todense(), num.array([[-3.0,3.0,0.0,0.0,0.0,0.0],\ |
---|
67 | [0.0,-6.0/sqrt(5.0),0.0,0.0,6.0/sqrt(5.0),0.0]])) |
---|
68 | assert num.allclose(kv_geo_structure[1].todense(), num.array([[-6.0/sqrt(5.0),0.0,6.0/sqrt(5.0),0.0,0.0,0.0],\ |
---|
69 | [3.0,-3.0,0.0,0.0,0.0,0.0]])) |
---|
70 | assert num.allclose(kv_geo_structure[2].todense(), num.array([[-6.0/sqrt(5.0),0.0,0.0,0.0,0.0,6.0/sqrt(5.0)],\ |
---|
71 | [0.0, -6.0/sqrt(5.0), 0.0, 6.0/sqrt(5.0), 0.0, 0.0]])) |
---|
72 | |
---|
73 | def test_apply_heights(self): |
---|
74 | operator1 = self.operator1() |
---|
75 | operator2 = self.operator2() |
---|
76 | |
---|
77 | h = num.array([1.0]) |
---|
78 | A = operator1.apply_stage_heights(h) |
---|
79 | assert num.allclose(A.todense(), num.array([-6.0-12.0/sqrt(5), 6.0/sqrt(5), 6.0, 6.0/sqrt(5)])) |
---|
80 | |
---|
81 | h = num.array([2.0]) |
---|
82 | A = operator1.apply_stage_heights(h) |
---|
83 | assert num.allclose(A.todense(), 1.5*num.array([-6.0-12.0/sqrt(5), 6.0/sqrt(5), 6.0, 6.0/sqrt(5)])) |
---|
84 | |
---|
85 | h = num.array([1.0, 1.0]) |
---|
86 | A = operator2.apply_stage_heights(h) |
---|
87 | #From test_kv_system_geometry |
---|
88 | A0 = num.array([[-3.0,3.0,0.0,0.0,0.0,0.0], |
---|
89 | [0.0,-6.0/sqrt(5.0),0.0,0.0,6.0/sqrt(5.0),0.0]]) |
---|
90 | A1 = num.array([[-6.0/sqrt(5.0),0.0,6.0/sqrt(5.0),0.0,0.0,0.0],\ |
---|
91 | [3.0,-3.0,0.0,0.0,0.0,0.0]]) |
---|
92 | A2 = num.array([[-6.0/sqrt(5.0),0.0,0.0,0.0,0.0,6.0/sqrt(5.0)],\ |
---|
93 | [0.0, -6.0/sqrt(5.0), 0.0, 6.0/sqrt(5.0), 0.0, 0.0]]) |
---|
94 | assert num.allclose(A.todense(), A0+A1+A2) |
---|
95 | |
---|
96 | h = num.array([2.0, 1.0]) |
---|
97 | A = operator2.apply_stage_heights(h) |
---|
98 | assert num.allclose(A.todense()[0,:], 1.5*A0[0,:]+1.5*A1[0,:]+1.5*A2[0,:]) |
---|
99 | assert num.allclose(A.todense()[1,:], A0[1,:]+1.5*A1[1,:]+A2[1,:]) |
---|
100 | |
---|
101 | h = num.array([-2.0, -2.0]) |
---|
102 | A = operator2.apply_stage_heights(h) |
---|
103 | assert num.allclose(A.todense()[0,:], -2*A0[0,:]-0.5*A1[0,:]-0.5*A2[0,:]) |
---|
104 | assert num.allclose(A.todense()[1,:], -0.5*A0[1,:]-2*A1[1,:]-0.5*A2[1,:]) |
---|
105 | |
---|
106 | def test_apply_kv_operator(self): |
---|
107 | operator1 = self.operator1() |
---|
108 | operator1.apply_stage_heights(num.array([[1.0]])) #h=1 |
---|
109 | operator1.build_boundary_vector() |
---|
110 | V1 = num.array([2.0]) #(uh)=2 |
---|
111 | V2 = num.array([2.0]) #(vh)=2 |
---|
112 | A = num.array([-6.0-12.0/sqrt(5), 6.0/sqrt(5), 6.0, 6.0/sqrt(5)]) |
---|
113 | U1 = num.array([[2.0],[1.0],[2.0],[1.0]]) |
---|
114 | U2 = num.array([[2.0],[2.0],[1.0],[0.0]]) |
---|
115 | |
---|
116 | operator1.set_qty_considered('u') |
---|
117 | operator1.build_boundary_vector() |
---|
118 | D1 = operator1.apply_kv_operator(V1) |
---|
119 | assert num.allclose(D1, 2*num.array(num.mat(A)*num.mat(U1)).reshape(1,)) |
---|
120 | |
---|
121 | operator1.set_qty_considered(2) |
---|
122 | operator1.build_boundary_vector() |
---|
123 | D2 = operator1.apply_kv_operator(V2) |
---|
124 | assert num.allclose(D2, 2*num.array(num.mat(A)*num.mat(U2)).reshape(1,)) |
---|
125 | |
---|
126 | def test_mul(self): |
---|
127 | operator1 = self.operator1() |
---|
128 | operator1.apply_stage_heights(num.array([[1.0]])) #h=1 |
---|
129 | operator1.set_qty_considered(1) |
---|
130 | operator1.build_boundary_vector() |
---|
131 | A = num.array([-6.0-12.0/sqrt(5), 6.0/sqrt(5), 6.0, 6.0/sqrt(5)]) |
---|
132 | V1 = num.array([2.0]) #(uh)=2 |
---|
133 | U1 = num.array([[2.0],[0.0],[0.0],[0.0]]) |
---|
134 | assert num.allclose(operator1 * V1, 2*num.array(num.mat(A)*num.mat(U1)).reshape(1,)) |
---|
135 | |
---|
136 | def test_cg_solve(self): |
---|
137 | #cf self.test_mul() |
---|
138 | operator1 = self.operator1() |
---|
139 | operator1.apply_stage_heights(num.array([[1.0]])) #h=1 |
---|
140 | operator1.set_qty_considered('u') |
---|
141 | V = num.array([2.0]) #h=1, (uh)=2 |
---|
142 | A = num.array([-6.0-12.0/sqrt(5), 6.0/sqrt(5), 6.0, 6.0/sqrt(5)]) |
---|
143 | U = num.array([[2.0,2.0],[1.0,2.0],[2.0,1.0],[1.0,0.0]]) |
---|
144 | test = 2*num.mat(A)*num.mat(U[:,0].reshape(4,1)) |
---|
145 | X = operator1.cg_solve(num.array(test).reshape(1,)) |
---|
146 | assert num.allclose(V, X) |
---|
147 | |
---|
148 | ################################################################################ |
---|
149 | |
---|
150 | if __name__ == "__main__": |
---|
151 | suite = unittest.makeSuite(Test_kv_solver_new, 'test') |
---|
152 | runner = unittest.TextTestRunner() |
---|
153 | runner.run(suite) |
---|