source: anuga_work/development/2010-projects/anuga_1d/base/test_quantity_1d.py @ 7855

Last change on this file since 7855 was 7855, checked in by steve, 13 years ago

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