source: trunk/anuga_work/development/2010-projects/kv/test_kv_solver_new.py @ 7884

Last change on this file since 7884 was 7884, checked in by steve, 14 years ago

Moving 2010 project
  • Property svn:executable set to *
File size: 6.9 KB
Line 
1from anuga.abstract_2d_finite_volumes.generic_domain import Generic_Domain
2from anuga.abstract_2d_finite_volumes.generic_boundary_conditions import Dirichlet_boundary
3from kv_solver_new import Kinematic_Viscosity_Operator
4import numpy as num
5from math import sqrt
6import unittest
7
8class 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
150if __name__ == "__main__":
151    suite = unittest.makeSuite(Test_kv_solver_new, 'test')
152    runner = unittest.TextTestRunner()
153    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.