source: storm_surge/pyvolution-1d/test_quantity.py @ 865

Last change on this file since 865 was 1, checked in by duncan, 20 years ago

initial import
File size: 13.1 KB
Line 
1#!/usr/bin/env python
2
3import unittest
4from math import sqrt, pi
5
6
7from quantity import *
8#from config import epsilon
9from Numeric import allclose, array, zeros, Float
10
11       
12class TestCase(unittest.TestCase):
13    def setUp(self):
14        from domain import Domain
15       
16        a = 0.0
17        b = 1.0
18        c = 2.0
19        d = 2.5
20        e = 3.1
21        f = 4.0
22
23        self.points          = [a, b, c, d, e, f]
24        self.centroids       = [(a+b)/2,(b+c)/2,(c+d)/2,(d+e)/2,(e+f)/2]
25        self.vertex_values   = [[1.0,2.0],[2.0,3.0],[3.0,4.0],[4.5,5],[5.5,5.6]]
26        self.centroid_values = [[1.5,      2.5,      3.5,      4.75,   5.55]]
27
28        self.domain1 = Domain(self.points[0:2])
29        self.domain5 = Domain(self.points)
30       
31    def tearDown(self):
32        pass
33        #print "  Tearing down"
34
35
36    def test_creation(self):
37       
38        quantity = Quantity(self.domain5, self.vertex_values)
39        assert allclose(quantity.centroid_values, self.centroid_values)
40
41        try:
42            quantity = Quantity()
43        except:
44            pass
45        else:
46            raise 'Should have raised empty quantity exception'
47
48
49        try:
50            quantity = Quantity([1,2])
51        except AssertionError:
52            pass
53        except:
54            raise 'Should have raised "missing domain object" error'       
55       
56
57    def test_creation_zeros(self):
58       
59        quantity = Quantity(self.domain1)
60        assert allclose(quantity.vertex_values, [[0.,0.]])
61
62
63        quantity = Quantity(self.domain5)
64        assert allclose(quantity.vertex_values, [[0.,0.], [0.,0.],
65                                                 [0.,0.], [0.,0.],
66                                                 [0.,0.]])       
67
68       
69    def test_interpolation(self):
70        quantity = Quantity(self.domain1, [[1,2]])
71        assert allclose(quantity.centroid_values, 1.5) #Centroid
72
73
74    def test_interpolation2(self):
75        quantity = Quantity(self.domain5,
76                            [[1,2], [5,5], [0,9], [-6, 3], [3,4]])
77        assert allclose(quantity.centroid_values, [1.5, 5., 4.5, -1.5, 3.5 ]) #Centroid
78
79    def test_set_values(self):
80        quantity = Quantity(self.domain5)
81
82
83        quantity.set_values([[1,2], [5,5], [0,0], [-6, 3], [-2,4]],
84                            location = 'vertices')
85        assert allclose(quantity.vertex_values,
86                        [[1,2], [5,5], [0,0], [-6, 3], [-2,4]])
87        assert allclose(quantity.centroid_values, [1.5, 5., 0., -1.5, 1.0]) #Centroid
88
89        #Test default
90        quantity.set_values([[1,2], [5,5], [0,0], [-6, 3], [-2,4]])
91        assert allclose(quantity.vertex_values,
92                        [[1,2], [5,5], [0,0], [-6, 3], [-2,4]])       
93        assert allclose(quantity.centroid_values, [1.5, 5., 0., -1.5, 1.0]) #Centroid
94
95        #Test centroids
96        quantity.set_values([1,2,3,4,5], location = 'centroids')
97        assert allclose(quantity.centroid_values, [1., 2., 3., 4., 5.]) #Centroid
98
99        #Test exceptionsdatamining.anu.edu.au/svn/
100        try:
101            quantity.set_values([[1,2], [5,5], [0,0], [-6, 3], [-2,4]],
102                                location = 'bas kamel tuba')
103        except:
104            pass
105
106
107        try:
108            quantity.set_values([[1,2], [0,0]])
109        except AssertionError:
110            pass
111        except:
112            raise 'should have raised AssertionError'
113
114    def test_set_values_const(self):
115        quantity = Quantity(self.domain5)
116
117        quantity.set_values(1.0, location = 'vertices')
118        assert allclose(quantity.vertex_values,
119                        [[1,1], [1,1], [1,1], [1,1], [1,1]])
120        assert allclose(quantity.centroid_values, [1, 1, 1, 1, 1]) #Centroid
121
122
123        quantity.set_values(2.0, location = 'centroids')
124        assert allclose(quantity.centroid_values, [2, 2, 2, 2, 2])
125
126
127    def test_set_values_func(self):
128        quantity = Quantity(self.domain5)
129
130        def f(x):
131            return x**2
132
133        quantity.set_values(f, location = 'vertices')
134        assert allclose(quantity.vertex_values,
135                        [[0,1], [1,4], [4,6.25], [6.25,9.61], [9.61,16]])       
136        assert allclose(quantity.centroid_values,
137                        [0.5, 2.5, 5.125, 7.93, 12.805 ])
138       
139        quantity.set_values(f, location = 'centroids')
140        assert allclose(quantity.centroid_values,
141                        [0.25, 1.5**2, 2.25**2, 2.8**2, 3.55**2])       
142
143
144    def test_gradient(self):
145       
146        quantity = Conserved_quantity(self.domain5)
147
148        #Set up for a gradient of (3,0) at mid triangle
149        quantity.set_values([2.0, 4.0, 8.0, 2.0, 5.0],
150                            location = 'centroids')
151       
152        #Gradients
153        quantity.compute_gradients()
154        N = quantity.gradients.shape[0]
155       
156        G = quantity.gradients
157        G0 = zeros(N, Float)
158        Q = quantity.centroid_values
159        X = quantity.domain.centroids
160
161        def grad0(x0,x1,x2,q0,q1,q2):
162            return ((q0-q1)/(x0-x1)*(x2-x1) - (q2-q1)/(x2-x1)*(x0-x1))/(x2-x0)
163
164        def grad1(x0,x1,q0,q1):
165            return  (q1 - q0)/(x1 - x0)
166
167        #Check Gradients
168        G0[0] = grad1(X[0],X[1],Q[0],Q[1])
169        for i in range(1,4):
170            G0[i] = grad0(X[i-1],X[i],X[i+1],Q[i-1],Q[i],Q[i+1])
171        G0[4] = grad1(X[3],X[4],Q[3],Q[4])
172
173        assert allclose(G,G0)
174
175        quantity.extrapolate_first_order()
176       
177        assert allclose(quantity.vertex_values, [[2., 2.],
178                                                 [4., 4.],
179                                                 [8., 8.],
180                                                 [2., 2.],
181                                                 [5., 5.]])
182
183        quantity.extrapolate_second_order()
184
185        V = quantity.domain.vertices
186        for k in range(quantity.domain.number_of_elements):
187            G0[k] = G0[k]*(V[k,1]-V[k,0])*0.5
188
189        assert allclose(quantity.vertex_values, [[2.-G0[0], 2.+G0[0]],
190                                                 [4.-G0[1], 4.+G0[1]],
191                                                 [8.-G0[2], 8.+G0[2]],
192                                                 [2.-G0[3], 2.+G0[3]],
193                                                 [5.-G0[4], 5.+G0[4]]])
194
195
196
197    def test_second_order_extrapolation2(self):
198        quantity = Conserved_quantity(self.domain5)       
199
200        #Set up for a gradient of (2), f(x) = 2x
201        quantity.set_values([1., 3., 4.5, 5.6, 7.1],
202                            location = 'centroids')
203       
204        #Gradients
205        quantity.compute_gradients()
206        G = quantity.gradients
207       
208        allclose(G, 2.0)
209
210        quantity.extrapolate_second_order()
211
212        V = quantity.domain.vertices
213        Q = quantity.vertex_values
214       
215        assert allclose(Q[:,0], 2.0*V[:,0])
216        assert allclose(Q[:,1], 2.0*V[:,1])       
217
218
219
220
221
222       
223
224
225##     def test_first_order_extrapolator(self):
226##         quantity = Conserved_quantity(self.mesh4)               
227
228##         #Test centroids
229##         quantity.set_values([1.,2.,3.,4.], location = 'centroids')
230##         assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
231
232##         #Extrapolate
233##         quantity.extrapolate_first_order()
234
235##         #Check vertices but not edge values
236##         assert allclose(quantity.vertex_values,
237##                         [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
238
239
240    def test_second_order_limit_extrapolator(self):
241        quantity = Conserved_quantity(self.domain5)                       
242
243        #Set up for a gradient of (3,0) at mid triangle
244        quantity.set_values([2.0, 4.0, 8.0, 2.0, 0.0],
245                            location = 'centroids')
246       
247
248
249        quantity.extrapolate_second_order()
250        quantity.limit()
251
252
253        #Assert that central triangle is limited by neighbours
254        assert quantity.vertex_values[0,0] >= 2.0
255        assert quantity.vertex_values[0,0] <= 4.0
256        assert quantity.vertex_values[0,1] >= 2.0
257        assert quantity.vertex_values[0,1] <= 4.0
258       
259        assert quantity.vertex_values[1,0] >= 2.0
260        assert quantity.vertex_values[1,0] <= 8.0
261        assert quantity.vertex_values[1,1] >= 2.0
262        assert quantity.vertex_values[1,1] <= 8.0
263
264        assert quantity.vertex_values[2,0] >= 2.0
265        assert quantity.vertex_values[2,0] <= 8.0
266        assert quantity.vertex_values[2,1] >= 2.0
267        assert quantity.vertex_values[2,1] <= 8.0
268
269        assert quantity.vertex_values[3,0] >= 0.0
270        assert quantity.vertex_values[3,0] <= 8.0
271        assert quantity.vertex_values[3,1] >= 0.0
272        assert quantity.vertex_values[3,1] <= 8.0
273
274        assert quantity.vertex_values[4,0] >= 0.0
275        assert quantity.vertex_values[4,0] <= 2.0
276        assert quantity.vertex_values[4,1] >= 0.0
277        assert quantity.vertex_values[4,1] <= 2.0
278       
279
280       
281        #Assert that quantities are conserved
282        from Numeric import sum
283        for k in range(quantity.centroid_values.shape[0]):
284            assert allclose (quantity.centroid_values[k],
285                             sum(quantity.vertex_values[k,:])/2)
286       
287       
288       
289
290    def test_limiter(self):
291        quantity = Conserved_quantity(self.domain5)
292
293        #Create a deliberate overshoot (e.g. from gradient computation)
294        quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
295
296
297        #Limit
298        quantity.limit()
299
300        #Assert that central triangle is limited by neighbours
301        assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0]
302        assert quantity.vertex_values[1,0] <= quantity.vertex_values[3,1]
303       
304        assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1]
305        assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2]
306       
307        assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0]
308        assert quantity.vertex_values[1,2] <= quantity.vertex_values[3,1]
309
310
311       
312        #Assert that quantities are conserved
313        from Numeric import sum
314        for k in range(quantity.centroid_values.shape[0]):
315            assert allclose (quantity.centroid_values[k],
316                             sum(quantity.vertex_values[k,:])/3)
317       
318
319
320##     def test_distribute_first_order(self):
321##         quantity = Conserved_quantity(self.mesh4)       
322
323##         #Test centroids
324##         quantity.set_values([1.,2.,3.,4.], location = 'centroids')
325##         assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
326
327
328##         #Extrapolate
329##         quantity.extrapolate_first_order()
330
331##         #Interpolate
332##         quantity.interpolate_from_vertices_to_edges()       
333
334##         assert allclose(quantity.vertex_values,
335##                         [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
336##         assert allclose(quantity.edge_values, [[1,1,1], [2,2,2],
337##                                                [3,3,3], [4, 4, 4]])
338
339
340
341    def test_update_explicit(self):
342        quantity = Conserved_quantity(self.domain5)
343
344        #Test centroids
345        quantity.set_values([1.,2.,3.,4.,5.], location = 'centroids')
346        assert allclose(quantity.centroid_values, [1, 2, 3, 4, 5]) #Centroid
347
348        #Set explicit_update
349        quantity.explicit_update = array( [1.,1.,1.,1.,1.] )
350
351        #Update with given timestep
352        quantity.update(0.1)
353
354        x = array([1, 2, 3, 4, 5]) + array( [.1,.1,.1,.1,.1] )
355        assert allclose( quantity.centroid_values, x)
356
357    def test_update_semi_implicit(self):
358        quantity = Conserved_quantity(self.domain5)
359
360        #Test centroids
361        quantity.set_values([1.,2.,3.,4.,5.], location = 'centroids')
362        assert allclose(quantity.centroid_values, [1, 2, 3, 4, 5]) #Centroid
363
364        #Set semi implicit update
365        quantity.semi_implicit_update = array( [1.,1.,1.,1.,1.] )
366
367        #Update with given timestep
368        quantity.update(0.1)
369
370        x = array([1, 2, 3, 4, 5])/array( [.9,.9,.9,.9,.9] )
371        assert allclose( quantity.centroid_values, x)
372
373    def test_both_updates(self):
374        quantity = Conserved_quantity(self.domain5)
375
376        #Test centroids
377        quantity.set_values([1.,2.,3.,4.,5.], location = 'centroids')
378        assert allclose(quantity.centroid_values, [1, 2, 3, 4, 5]) #Centroid
379
380        #Set explicit_update
381        quantity.explicit_update = array( [4.,3.,2.,1.,0.0] )
382       
383        #Set semi implicit update
384        quantity.semi_implicit_update = array( [1.,1.,1.,1.,1.] )
385
386        #Update with given timestep
387        quantity.update(0.1)
388
389        x = array([1, 2, 3, 4, 5]) + array( [.4,.3,.2,.1,0.0] )
390        x /= array( [.9,.9,.9,.9,.9] )
391        assert allclose( quantity.centroid_values, x)       
392
393       
394       
395#-------------------------------------------------------------
396if __name__ == "__main__":
397    suite = unittest.makeSuite(TestCase,'test')
398    runner = unittest.TextTestRunner()
399    runner.run(suite)
400
401
Note: See TracBrowser for help on using the repository browser.