source: anuga_work/development/anuga_1d/test_quantity.py @ 7825

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

anuga_1d works with numeric on 64bit machine
File size: 18.3 KB
Line 
1import quantity
2#!/usr/bin/env python
3
4import unittest
5from math import sqrt, pi
6
7
8from quantity import *
9from domain import *
10
11
12from Numeric import allclose, array, ones, Float, maximum, zeros
13
14
15class Test_Quantity(unittest.TestCase):
16    def setUp(self):
17        self.points3        = [0.0, 1.0, 2.0, 3.0]
18        self.vertex_values3 = [[1.0,2.0],[4.0,5.0],[-1.0,2.0]]
19        self.domain3        = Domain(self.points3)
20
21
22       
23        self.points4          = [0.0, 1.0, 2.5, 3.0, 5.0]
24        self.vertex_values4   = [[1.0,2.0],[4.0,5.0],[-1.0,2.0],[3.0,6.0]]
25        self.centroid_values4 = [1.5, 4.5, 0.5, 4.5]
26        self.boundary4        = {(0, 0): 'left', (3, 1): 'right'}
27        self.domain4          = Domain(self.points4,self.boundary4)
28
29        self.points10 = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]
30        self.domain10 = Domain(self.points10)
31
32        self.points6 = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0]
33        self.domain6 = Domain(self.points6)
34
35
36    def tearDown(self):
37        pass
38        #print "  Tearing down"
39
40
41    def test_creat_with_boundary(self):
42
43        assert self.domain4.boundary  == {(0, 0): 'left', (3, 1): 'right'}
44       
45    def test_creation(self):
46
47        quantity = Quantity(self.domain3)
48        assert allclose(quantity.vertex_values, [[0.0,0.0],[0.0,0.0],[0.0,0.0]])
49
50       
51        try:
52            quantity = Quantity()
53        except:
54            pass
55        else:
56            raise 'Should have raised empty quantity exception'
57
58
59        try:
60            quantity = Quantity([1,2,3])
61        except AssertionError:
62            pass
63        except:
64            raise 'Should have raised "mising domain object" error'
65
66
67    def test_creation_zeros(self):
68
69        quantity = Quantity(self.domain3)
70        assert allclose(quantity.centroid_values, [[0.,0.,0.]])
71
72
73        quantity = Quantity(self.domain4)
74        assert allclose(quantity.vertex_values, [[0.,0.], [0.,0.],
75                                                 [0.,0.], [0.,0.]])
76
77
78    def test_interpolation(self):
79        quantity = Quantity(self.domain4, self.vertex_values4)
80        assert allclose(quantity.centroid_values, self.centroid_values4) #Centroid
81
82
83
84    def test_interpolation2(self):
85        quantity = Quantity(self.domain4, self.vertex_values4)
86        assert allclose(quantity.centroid_values, self.centroid_values4) #Centroid
87
88        quantity.extrapolate_second_order()
89
90        #print quantity.vertex_values
91        assert allclose(quantity.vertex_values,[[ 0.3,  2.7],
92                                                [ 4.5,  4.5],
93                                                [ 0.5,  0.5],
94                                                [ 1.3,  7.7]])
95
96 
97
98 
99    def test_boundary_allocation(self):
100        quantity = Quantity(self.domain4,
101                            [[1,2], [5,5], [0,9], [-6, 3]])
102
103
104        assert quantity.boundary_values.shape[0] == len(self.domain4.boundary)
105
106
107    def test_set_values(self):
108        quantity = Quantity(self.domain4)
109
110        # Test vertices
111        quantity.set_values([[1,2], [5,5], [0,9], [-6, 3]], location = 'vertices')
112        assert allclose(quantity.vertex_values,
113                        [[1,2], [5,5], [0,9], [-6, 3]])
114        assert allclose(quantity.centroid_values, [1.5, 5., 4.5, -1.5]) #Centroid
115
116
117
118        # Test unique_vertices
119        quantity.set_values([1,2,4,-5,6], location='unique_vertices')
120        assert allclose(quantity.vertex_values,
121                        [[1,2], [2,4], [4,-5], [-5,6]])
122        assert allclose(quantity.centroid_values, [1.5, 3., -.5, .5]) #Centroid
123
124
125        # Test centroids
126        quantity.set_values([1,2,3,4], location = 'centroids')
127        assert allclose(quantity.centroid_values, [1., 2., 3., 4.]) #Centroid
128
129        # Test exceptions
130        try:
131            quantity.set_values([[1,3], [5,5], [0,9], [-6, 3]],
132                                location = 'bas kamel tuba')
133        except:
134            pass
135
136
137        try:
138            quantity.set_values([[1,3], [0,9]])
139        except AssertionError:
140            pass
141        except:
142            raise 'should have raised Assertionerror'
143
144
145
146    def test_set_values_const(self):
147        quantity = Quantity(self.domain4)
148
149        quantity.set_values(1.0, location = 'vertices')
150
151        assert allclose(quantity.vertex_values, [[1,1], [1,1], [1,1], [1, 1]])
152        assert allclose(quantity.centroid_values, [1, 1, 1, 1]) #Centroid
153
154
155        quantity.set_values(2.0, location = 'centroids')
156
157        assert allclose(quantity.centroid_values, [2, 2, 2, 2])
158
159
160    def test_set_values_func(self):
161        quantity = Quantity(self.domain4)
162
163        def f(x):
164            return x*x
165
166
167
168        quantity.set_values(f, location = 'vertices')
169
170
171        assert allclose(quantity.vertex_values,
172                        [[0,1], [1,6.25], [6.25,9], [9,25]])
173
174        assert allclose(quantity.centroid_values,
175                        [0.5, 3.625, 7.625, 34*0.5])
176
177
178        quantity.set_values(f, location = 'centroids')
179
180       
181        assert allclose(quantity.centroid_values,
182                        [0.25, 3.0625, 7.5625, 16.0])
183
184
185    def test_integral(self):
186        quantity = Quantity(self.domain4)
187
188        #Try constants first
189        const = 5
190        quantity.set_values(const, location = 'vertices')
191        #print 'Q', quantity.get_integral()
192
193        assert allclose(quantity.get_integral(), self.domain4.get_area() * const)
194
195        #Try with a linear function
196        def f(x):
197            return x
198
199        quantity.set_values(f, location = 'vertices')
200
201 
202        assert allclose (quantity.centroid_values,
203        [ 0.5,   1.75,  2.75,  4.  ])
204
205        assert allclose (quantity.vertex_values, [[ 0.,   1. ],
206                                                  [ 1.,   2.5],
207                                                  [ 2.5,  3. ],
208                                                  [ 3.,   5. ]])
209
210
211        ref_integral = 0.5 + 1.5*1.75 + 0.5*2.75 + 2.0*4.0
212
213        assert allclose (quantity.get_integral(), ref_integral)
214
215
216
217
218
219    def test_set_values_from_quantity(self):
220
221        quantity1 = Quantity(self.domain4)
222        quantity1.set_values([0,1,2,3,4], location='unique_vertices')
223
224        assert allclose(quantity1.vertex_values,
225                        [[0,1], [1,2], [2,3], [3,4]])
226
227
228        quantity2 = Quantity(self.domain4)
229        quantity2.set_values(quantity1)
230        assert allclose(quantity2.vertex_values,
231                        [[0,1], [1,2], [2,3], [3,4]])
232
233        quantity2.set_values(2*quantity1)
234
235        assert allclose(quantity2.vertex_values,
236                        [[0,2], [2,4], [4,6], [6,8]])
237
238        quantity2.set_values(2*quantity1 + 3)
239        assert allclose(quantity2.vertex_values,
240                        [[3,5], [5,7], [7,9], [9,11]])
241
242
243
244    def test_overloading(self):
245
246        quantity1 = Quantity(self.domain4)
247        quantity1.set_values( [[0,1],[1,2],[2,3],[3,4]],
248                            location = 'vertices')
249
250        assert allclose(quantity1.vertex_values,
251                        [[0,1], [1,2], [2,3], [3,4]])
252
253
254        quantity2 = Quantity(self.domain4)
255        quantity2.set_values([[1,2], [5,5], [0,9], [-6, 3]],
256                             location = 'vertices')
257
258
259
260        quantity3 = Quantity(self.domain4)
261        quantity3.set_values([[2,2], [7,8], [7,6], [3, -8]],
262                             location = 'vertices')
263
264
265        # Negation
266        Q = -quantity1
267        assert allclose(Q.vertex_values, -quantity1.vertex_values)
268        assert allclose(Q.centroid_values, -quantity1.centroid_values)
269
270
271        # Addition
272        Q = quantity1 + 7
273        assert allclose(Q.vertex_values, quantity1.vertex_values + 7)
274        assert allclose(Q.centroid_values, quantity1.centroid_values + 7)
275
276        Q = 7 + quantity1
277        assert allclose(Q.vertex_values, quantity1.vertex_values + 7)
278        assert allclose(Q.centroid_values, quantity1.centroid_values + 7)
279
280        Q = quantity1 + quantity2
281        assert allclose(Q.vertex_values,
282                        quantity1.vertex_values + quantity2.vertex_values)
283        assert allclose(Q.centroid_values,
284                        quantity1.centroid_values + quantity2.centroid_values)
285
286        Q = quantity1 + quantity2 - 3
287        assert allclose(Q.vertex_values,
288                        quantity1.vertex_values + quantity2.vertex_values - 3)
289
290        Q = quantity1 - quantity2
291        assert allclose(Q.vertex_values,
292                        quantity1.vertex_values - quantity2.vertex_values)
293
294        #Scaling
295        Q = quantity1*3
296        assert allclose(Q.vertex_values, quantity1.vertex_values*3)
297        assert allclose(Q.centroid_values, quantity1.centroid_values*3)
298
299        Q = 3*quantity1
300        assert allclose(Q.vertex_values, quantity1.vertex_values*3)
301
302        #Multiplication
303        Q = quantity1 * quantity2
304        assert allclose(Q.vertex_values,
305                        quantity1.vertex_values * quantity2.vertex_values)
306
307        #Linear combinations
308        Q = 4*quantity1 + 2
309        assert allclose(Q.vertex_values,
310                        4*quantity1.vertex_values + 2)
311
312        Q = quantity1*quantity2 + 2
313        assert allclose(Q.vertex_values,
314                        quantity1.vertex_values * quantity2.vertex_values + 2)
315
316        Q = quantity1*quantity2 + quantity3
317        assert allclose(Q.vertex_values,
318                        quantity1.vertex_values * quantity2.vertex_values +
319                        quantity3.vertex_values)
320        Q = quantity1*quantity2 + 3*quantity3
321        assert allclose(Q.vertex_values,
322                        quantity1.vertex_values * quantity2.vertex_values +
323                        3*quantity3.vertex_values)
324        Q = quantity1*quantity2 + 3*quantity3 + 5.0
325        assert allclose(Q.vertex_values,
326                        quantity1.vertex_values * quantity2.vertex_values +
327                        3*quantity3.vertex_values + 5)
328
329        Q = quantity1*quantity2 - quantity3
330        assert allclose(Q.vertex_values,
331                        quantity1.vertex_values * quantity2.vertex_values -
332                        quantity3.vertex_values)
333        Q = 1.5*quantity1*quantity2 - 3*quantity3 + 5.0
334        assert allclose(Q.vertex_values,
335                        1.5*quantity1.vertex_values * quantity2.vertex_values -
336                        3*quantity3.vertex_values + 5)
337
338        #Try combining quantities and arrays and scalars
339        Q = 1.5*quantity1*quantity2.vertex_values -\
340            3*quantity3.vertex_values + 5.0
341        assert allclose(Q.vertex_values,
342                        1.5*quantity1.vertex_values * quantity2.vertex_values -
343                        3*quantity3.vertex_values + 5)
344
345
346        #Powers
347        Q = quantity1**2
348        assert allclose(Q.vertex_values, quantity1.vertex_values**2)
349
350        Q = quantity1**2 +quantity2**2
351        assert allclose(Q.vertex_values,
352                        quantity1.vertex_values**2 + \
353                        quantity2.vertex_values**2)
354
355        Q = (quantity1**2 +quantity2**2)**0.5
356        assert allclose(Q.vertex_values,
357                        (quantity1.vertex_values**2 + \
358                        quantity2.vertex_values**2)**0.5)
359
360    def test_compute_gradient(self):
361        quantity = Quantity(self.domain6)
362
363        #Set up for a gradient of (2,0) at mid triangle
364        quantity.set_values([2.0, 4.0, 4.0, 5.0, 10.0, 12.0],
365                            location = 'centroids')
366
367
368        #Gradients
369        quantity.compute_gradients()
370
371        a = quantity.gradients
372
373        assert allclose(a, [ 3., 1., 0.5, 3., 3.5, 0.5])
374       
375        quantity.extrapolate_second_order()
376
377
378        assert allclose(quantity.vertex_values, [[ 1., 3. ],
379                                                 [ 4., 4. ],
380                                                 [ 4., 4. ],
381                                                 [ 4., 6.],
382                                                 [ 8.25, 11.75],
383                                                 [ 11.,  13. ]])
384
385                                                 
386
387    def test_second_order_extrapolation2(self):
388        quantity = Quantity(self.domain4)
389
390        #Set up for a gradient of (3,1), f(x) = 3x+y
391        quantity.set_values([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3],
392                            location = 'centroids')
393
394        #Gradients
395        quantity.compute_gradients()
396
397        a = quantity.gradients
398       
399
400        assert allclose(a[1], 2.8)
401
402        #Work out the others
403
404        quantity.extrapolate_second_order()
405       
406
407        assert allclose(quantity.vertex_values[1,0], 3.33333333)
408        assert allclose(quantity.vertex_values[1,1], 7.33333333)
409
410        assert allclose(quantity.centroid_values[1], 0.5*(7.33333333+3.33333333) )
411
412
413
414
415    def test_backup_saxpy_centroid_values(self):
416        quantity = Quantity(self.domain4)
417
418        #Set up for a gradient of (3,1), f(x) = 3x+y
419        c_values = array([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3])
420        d_values = array([1.0, 2.0, 3.0, 4.0])
421        quantity.set_values(c_values, location = 'centroids')
422
423        #Backup
424        quantity.backup_centroid_values()
425
426        #print quantity.vertex_values
427        assert allclose(quantity.centroid_values, quantity.centroid_backup_values)
428
429
430        quantity.set_values(d_values, location = 'centroids')
431
432        quantity.saxpy_centroid_values(2.0, 3.0)
433
434        assert allclose(quantity.centroid_values, 2.0*d_values + 3.0*c_values)
435
436
437
438    def test_first_order_extrapolator(self):
439        quantity = Quantity(self.domain4)
440
441        centroid_values = array([1.,2.,3.,4.])
442        quantity.set_values(centroid_values, location = 'centroids')
443        assert allclose(quantity.centroid_values, centroid_values) #Centroid
444
445        #Extrapolate
446        quantity.extrapolate_first_order()
447
448        #Check that gradient is zero
449        a = quantity.gradients
450        assert allclose(a, [0,0,0,0])
451       
452
453        #Check vertices but not edge values
454        assert allclose(quantity.vertex_values,
455                        [[1,1], [2,2], [3,3], [4,4]])
456
457
458    def test_second_order_extrapolator(self):
459        quantity = Quantity(self.domain4)
460
461        #Set up for a gradient of (3,0) at mid triangle
462        quantity.set_values([2.0, 4.0, 8.0, 2.0],
463                            location = 'centroids')
464
465
466
467        quantity.extrapolate_second_order()
468
469
470        #Assert that quantities are conserved
471        from Numeric import sum
472        for k in range(quantity.centroid_values.shape[0]):
473            assert allclose (quantity.centroid_values[k],
474                             sum(quantity.vertex_values[k,:])/2.0)
475
476
477    def test_limit(self):
478        quantity = Quantity(self.domain4)
479
480        #Create a deliberate overshoot (e.g. from gradient computation)
481        quantity.set_values([[0,0], [2,20], [-20,3], [8,3]])
482
483        #Limit
484        quantity.limit_minmod()
485
486       
487        #cells 1 and 2 are local max and min
488        assert quantity.vertex_values[1][0] == quantity.centroid_values[1]
489        assert quantity.vertex_values[1][1] == quantity.centroid_values[1]
490
491        assert quantity.vertex_values[2][0] == quantity.centroid_values[2]
492        assert quantity.vertex_values[2][1] == quantity.centroid_values[2]
493
494
495
496    def test_distribute_first_order(self):
497        quantity = Quantity(self.domain4)
498
499        #Test centroids
500        centroid_values = array([1.,2.,3.,4.])
501        quantity.set_values(centroid_values, location = 'centroids')
502        assert allclose(quantity.centroid_values, centroid_values) #Centroid
503
504
505        #Extrapolate from centroid to vertices and edges
506        quantity.extrapolate_first_order()
507       
508        assert allclose(quantity.vertex_values,[[ 1.,  1.],
509                                                [ 2.,  2.],
510                                                [ 3.,  3.],
511                                                [ 4.,  4.]])
512 
513
514
515    def test_distribute_second_order(self):
516        quantity = Quantity(self.domain4)
517
518        #Test centroids
519        centroid_values = array([2.,4.,8.,2.])
520        quantity.set_values(centroid_values, location = 'centroids')
521        assert allclose(quantity.centroid_values, centroid_values) #Centroid
522
523
524        #Extrapolate
525        quantity.extrapolate_second_order()
526
527        assert allclose(quantity.vertex_values, [[ 1.2,  2.8],
528                                                 [ 2.,  6. ],
529                                                 [ 8.,   8. ],
530                                                 [ 6.8, -2.8]])
531
532
533    def test_update_explicit(self):
534        quantity = Quantity(self.domain4)
535
536        #Test centroids
537        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
538        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
539
540        #Set explicit_update
541        quantity.explicit_update = array( [1.,1.,1.,1.] )
542
543        #Update with given timestep
544        quantity.update(0.1)
545
546        x = array([1, 2, 3, 4]) + array( [.1,.1,.1,.1] )
547        assert allclose( quantity.centroid_values, x)
548
549    def test_update_semi_implicit(self):
550        quantity = Quantity(self.domain4)
551
552        #Test centroids
553        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
554        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
555
556        #Set semi implicit update
557        quantity.semi_implicit_update = array([1.,1.,1.,1.])
558
559        #Update with given timestep
560        timestep = 0.1
561        quantity.update(timestep)
562
563        sem = array([1.,1.,1.,1.])
564        denom = ones(4, Float)-timestep*sem
565
566        x = array([1, 2, 3, 4])/denom
567        assert allclose( quantity.centroid_values, x)
568
569
570    def test_both_updates(self):
571        quantity = Quantity(self.domain4)
572
573        #Test centroids
574        centroid_values = array( [1, 2, 3, 4] )
575        quantity.set_values(centroid_values, location = 'centroids')
576        assert allclose(quantity.centroid_values, centroid_values) #Centroid
577
578        #Set explicit_update
579        explicit_update = array( [4.,3.,2.,1.] )
580        quantity.explicit_update[:,] = explicit_update
581
582        #Set semi implicit update
583        semi_implicit_update = array( [1.,1.,1.,1.] )
584        quantity.semi_implicit_update[:,] = semi_implicit_update
585
586        #Update with given timestep
587        timestep = 0.1
588        quantity.update(0.1)
589
590        denom = 1.0-timestep*semi_implicit_update
591        x = (centroid_values + timestep*explicit_update)/denom
592 
593        assert allclose( quantity.centroid_values, x)
594
595
596 
597#-------------------------------------------------------------
598if __name__ == "__main__":
599    suite = unittest.makeSuite(Test_Quantity, 'test')
600    runner = unittest.TextTestRunner()
601    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.