source: trunk/anuga_work/development/kv/test_kinematic_viscosity.py @ 8932

Last change on this file since 8932 was 8052, checked in by steve, 15 years ago

Modified C code to take long arrays instead of int (caused error on 64 bit machine)
File size: 6.1 KB
Line 
1from anuga.abstract_2d_finite_volumes.generic_domain import Generic_Domain
2from anuga import Dirichlet_boundary
3from kinematic_viscosity import Kinematic_Viscosity_Operator
4import numpy as num
5from math import sqrt
6import unittest
7
8class Test_Kinematic_Viscosity(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                domain = Generic_Domain(source=points,triangles=elements,boundary=boundary_map)
24                return Kinematic_Viscosity_Operator(domain)
25
26        #Second test operator class (2 triangles)
27        def operator2(self):
28                points = num.array([[0.0,0.0],[1.0,0.0],[1.0,1.0],[0.0,1.0]])
29                elements = num.array([[0,1,3],[1,2,3]])
30                boundary_map = {}
31                boundary_map[(0,1)] = Dirichlet_boundary([1,1,2])
32                boundary_map[(0,2)] = Dirichlet_boundary([1,2,2])
33                boundary_map[(1,0)] = Dirichlet_boundary([1,1,0])
34                boundary_map[(1,2)] = Dirichlet_boundary([1,2,1])
35                domain = Generic_Domain(source=points,triangles=elements,boundary=boundary_map)
36                return Kinematic_Viscosity_Operator(domain)
37
38        def test_enumerate_boundary(self):
39                operator1 = self.operator1()
40                boundary_enum = operator1.boundary_enum
41                assert boundary_enum[(0,0)] == 0
42                assert boundary_enum[(0,1)] == 1
43                assert boundary_enum[(0,2)] == 2
44               
45                operator2 = self.operator2()
46                boundary_enum = operator2.boundary_enum
47                assert boundary_enum[(0,1)] == 0
48                assert boundary_enum[(0,2)] == 1
49                assert boundary_enum[(1,0)] == 2
50                assert boundary_enum[(1,2)] == 3
51
52        def test_geo_structure(self):
53                operator1 = self.operator1()
54                indices = operator1.geo_structure_indices
55                values = operator1.geo_structure_values
56                assert num.allclose(indices, num.array([[1, 2, 3]]))
57                assert num.allclose(values, num.array([[-6.0, -6.0/sqrt(5), -6.0/sqrt(5)]]))
58               
59                operator2 = self.operator2()
60                indices = operator2.geo_structure_indices
61                values = operator2.geo_structure_values
62                assert num.allclose(indices, num.array([[1,2,3],[4,0,5]]))
63                assert num.allclose(values, num.array([[-3.0,-6.0/sqrt(5),-6.0/sqrt(5)],[-6.0/sqrt(5),-3.0,-6.0/sqrt(5)]]))
64               
65        def test_apply_heights(self):
66                operator1 = self.operator1()
67                operator2 = self.operator2()
68               
69                h = num.array([1.0])
70                A = operator1.apply_stage_heights(h)
71                assert num.allclose(A.todense(), num.array([-6.0-12.0/sqrt(5), 6.0,  6.0/sqrt(5), 6.0/sqrt(5)]))
72               
73                h = num.array([2.0])
74                A = operator1.apply_stage_heights(h)
75                assert num.allclose(A.todense(), 1.5*num.array([-6.0-12.0/sqrt(5), 6.0, 6.0/sqrt(5), 6.0/sqrt(5)]))
76               
77                h = num.array([1.0, 1.0])
78                A = operator2.apply_stage_heights(h)
79                #From test_kv_system_geometry
80                A0 = num.array([[-3.0,3.0,0.0,0.0,0.0,0.0],
81                                                [0.0,-6.0/sqrt(5.0),0.0,0.0,6.0/sqrt(5.0),0.0]])
82                A1 = num.array([[-6.0/sqrt(5.0),0.0,6.0/sqrt(5.0),0.0,0.0,0.0],\
83                                                [3.0,-3.0,0.0,0.0,0.0,0.0]])
84                A2 = num.array([[-6.0/sqrt(5.0),0.0,0.0,6.0/sqrt(5.0),0.0,0.0],\
85                                                [0.0, -6.0/sqrt(5.0), 0.0, 0.0, 0.0, 6.0/sqrt(5.0)]])
86                assert num.allclose(A.todense(), A0+A1+A2)
87               
88                h = num.array([2.0, 1.0])
89                A = operator2.apply_stage_heights(h)
90                assert num.allclose(A.todense()[0,:], 1.5*A0[0,:]+1.5*A1[0,:]+1.5*A2[0,:])
91                assert num.allclose(A.todense()[1,:], A0[1,:]+1.5*A1[1,:]+A2[1,:])
92               
93                h = num.array([-2.0, -2.0])
94                A = operator2.apply_stage_heights(h)
95                assert num.allclose(A.todense()[0,:], -2*A0[0,:]-0.5*A1[0,:]-0.5*A2[0,:])
96                assert num.allclose(A.todense()[1,:], -0.5*A0[1,:]-2*A1[1,:]-0.5*A2[1,:])
97
98        def test_elliptic_multiply(self):
99                operator1 = self.operator1()
100                operator1.apply_stage_heights(num.array([[1.0]])) #h=1
101                operator1.build_boundary_vector()
102                V1 = num.array([2.0]) #(uh)=2 <- Centriod value
103                V2 = num.array([2.0]) #(vh)=2 <- Centroid value
104                A = num.array([-6.0-12.0/sqrt(5), 6.0,  6.0/sqrt(5), 6.0/sqrt(5)])
105                U1 = num.array([[2.0],[2.0],[1.0],[1.0]]) #(uh) for centroid, 3 edges
106                U2 = num.array([[2.0],[1.0],[2.0],[0.0]]) #(vh) for centroid, 3 edge
107               
108                operator1.set_qty_considered('u')
109                D1 = operator1.elliptic_multiply(V1)
110                assert num.allclose(D1, 2*num.array(num.mat(A)*num.mat(U1)).reshape(1,)) #2* for triangle_areas
111               
112                operator1.set_qty_considered(2)
113                D2 = operator1.elliptic_multiply(V2)
114                assert num.allclose(D2, 2*num.array(num.mat(A)*num.mat(U2)).reshape(1,)) #2* for triangle_areas
115       
116        def test_mul(self):
117                operator1 = self.operator1()
118                operator1.apply_stage_heights(num.array([[1.0]])) #h=1
119                operator1.set_qty_considered(1)
120                operator1.build_boundary_vector()
121                A = num.array([-6.0-12.0/sqrt(5), 6.0,  6.0/sqrt(5), 6.0/sqrt(5)])
122                V1 = num.array([2.0]) #(uh)=2
123                U1 = num.array([[2.0],[0.0],[0.0],[0.0]])
124                assert num.allclose(operator1 * V1, 2*num.array(num.mat(A)*num.mat(U1)).reshape(1,))
125       
126        def test_cg_solve(self):
127                #cf self.test_mul()
128                operator1 = self.operator1()
129                operator1.apply_stage_heights(num.array([[1.0]])) #h=1
130                operator1.set_qty_considered('u')
131                V = num.array([2.0]) #h=1, (uh)=2
132                A = num.array([-6.0-12.0/sqrt(5), 6.0,  6.0/sqrt(5), 6.0/sqrt(5)])
133                U = num.array([[2.0,2.0],[2.0,1.0],[1.0,2.0],[1.0,0.0]])
134                test = 2*num.mat(A)*num.mat(U[:,0].reshape(4,1))
135                X = operator1.cg_solve(num.array(test).reshape(1,))
136                assert num.allclose(V, X)
137       
138        def test_parabolic_solve(self):
139                operator1 = self.operator1()
140                operator1.apply_stage_heights(num.array([[1.0]])) #h=1
141                A = num.array([-6.0-12.0/sqrt(5), 6.0,  6.0/sqrt(5), 6.0/sqrt(5)])
142                U = num.array([[2.0,1.0],[2.0,1.0],[1.0,2.0],[1.0,0.0]])
143                u = num.array([[2.0,1.0]])
144                U_new = operator1.parabolic_solver(u)
145                U_mod = num.array([[0.0,0.0],[2.0,1.0],[1.0,2.0],[1.0,0.0]])
146                U_mod[0,:] = U_new
147                assert num.allclose(U_new - operator1.dt * 2 * num.mat(A)*num.mat(U_mod), U[0,:])
148
149################################################################################
150
151if __name__ == "__main__":
152    suite = unittest.makeSuite(Test_Kinematic_Viscosity, 'test', verbose=True)
153    runner = unittest.TextTestRunner()
154    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.