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

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

Moving 2010 project

File size: 21.9 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        else:
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.beta = 2.0
380        quantity.limiter = "minmod_kurganov"
381        quantity.extrapolate_second_order()
382
383
384        assert allclose(quantity.vertex_values, [[ 1., 3. ],
385                                                 [ 4., 4. ],
386                                                 [ 4., 4. ],
387                                                 [ 4., 6.],
388                                                 [ 8.25, 11.75],
389                                                 [ 11.,  13. ]])
390
391                                                 
392
393    def test_second_order_extrapolation2(self):
394        quantity = Quantity(self.domain4)
395
396        #Set up for a gradient of (3,1), f(x) = 3x+y
397        quantity.set_values([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3],
398                            location = 'centroids')
399
400        #Gradients
401        quantity.compute_gradients()
402
403        a = quantity.gradients
404       
405
406        assert allclose(a[1], 2.8)
407
408        #Work out the others
409        quantity.beta = 2.0
410        quantity.limiter = "minmod_kurganov"
411        quantity.extrapolate_second_order()
412       
413
414        assert allclose(quantity.vertex_values[1,0], 3.33333333)
415        assert allclose(quantity.vertex_values[1,1], 7.33333333)
416
417        assert allclose(quantity.centroid_values[1], 0.5*(7.33333333+3.33333333) )
418
419
420
421
422    def test_backup_saxpy_centroid_values(self):
423        quantity = Quantity(self.domain4)
424
425        #Set up for a gradient of (3,1), f(x) = 3x+y
426        c_values = array([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3])
427        d_values = array([1.0, 2.0, 3.0, 4.0])
428        quantity.set_values(c_values, location = 'centroids')
429
430        #Backup
431        quantity.backup_centroid_values()
432
433        #print quantity.vertex_values
434        assert allclose(quantity.centroid_values, quantity.centroid_backup_values)
435
436
437        quantity.set_values(d_values, location = 'centroids')
438
439        quantity.saxpy_centroid_values(2.0, 3.0)
440
441        assert allclose(quantity.centroid_values, 2.0*d_values + 3.0*c_values)
442
443
444
445    def test_first_order_extrapolator(self):
446        quantity = Quantity(self.domain4)
447
448        centroid_values = array([1.,2.,3.,4.])
449        quantity.set_values(centroid_values, location = 'centroids')
450        assert allclose(quantity.centroid_values, centroid_values) #Centroid
451
452        #Extrapolate
453        quantity.extrapolate_first_order()
454
455        #Check that gradient is zero
456        a = quantity.gradients
457        assert allclose(a, [0,0,0,0])
458       
459
460        #Check vertices but not edge values
461        assert allclose(quantity.vertex_values,
462                        [[1,1], [2,2], [3,3], [4,4]])
463
464
465    def test_second_order_extrapolator(self):
466        quantity = Quantity(self.domain4)
467
468        #Set up for a gradient of (3,0) at mid triangle
469        quantity.set_values([2.0, 4.0, 8.0, 2.0],
470                            location = 'centroids')
471
472
473
474        quantity.extrapolate_second_order()
475
476
477        #Assert that quantities are conserved
478        for k in range(quantity.centroid_values.shape[0]):
479            assert allclose (quantity.centroid_values[k],
480                             numpy.sum(quantity.vertex_values[k,:])/2.0)
481
482
483    def test_limit(self):
484        quantity = Quantity(self.domain4)
485
486        #Create a deliberate overshoot (e.g. from gradient computation)
487        quantity.set_values([[0,0], [2,20], [-20,3], [8,3]])
488
489        #Limit
490        quantity.limit_minmod()
491
492       
493        #cells 1 and 2 are local max and min
494        assert quantity.vertex_values[1][0] == quantity.centroid_values[1]
495        assert quantity.vertex_values[1][1] == quantity.centroid_values[1]
496
497        assert quantity.vertex_values[2][0] == quantity.centroid_values[2]
498        assert quantity.vertex_values[2][1] == quantity.centroid_values[2]
499
500
501
502
503    def test_limit_minmod(self):
504        quantity = Quantity(self.domain4)
505
506        #Create a deliberate overshoot (e.g. from gradient computation)
507        quantity.set_values([[0,0], [2,10], [9,13], [12,14]])
508
509        #Limit
510        quantity.limiter = 'minmod'
511        quantity.limit_range()
512
513
514
515
516        assert allclose(quantity.vertex_values, [[ -2.4,   2.4],
517                                                 [  2.4,   9.6],
518                                                 [ 10.6,  11.4],
519                                                 [ 11.4,  14.6]] )
520
521        from anuga_1d.base.quantity_ext import limit_minmod_ext
522
523        quantity.set_values([[0,0], [2,10], [9,13], [12,14]])
524
525        limit_minmod_ext(quantity)
526
527
528
529
530        assert allclose(quantity.vertex_values, [[ -2.4,   2.4],
531                                                 [  2.4,   9.6],
532                                                 [ 10.6,  11.4],
533                                                 [ 11.4,  14.6]] )
534
535
536
537    def test_limit_minmod_kurganov(self):
538        quantity = Quantity(self.domain4)
539
540        #Create a deliberate overshoot (e.g. from gradient computation)
541        quantity.set_values([[0,0], [2,10], [9,13], [12,14]])
542
543        #Limit
544        quantity.limiter = 'minmod_kurganov'
545        quantity.beta = 0.5
546        quantity.limit_range()
547
548
549        #print quantity.vertex_values
550
551        assert allclose(quantity.vertex_values, [[ -2.4,   2.4],
552                                                 [  4.2,   7.8],
553                                                 [ 10.8,  11.2],
554                                                 [ 11.4,  14.6]])
555
556        from anuga_1d.base.quantity_ext import limit_minmod_kurganov_ext
557
558        quantity.set_values([[0,0], [2,10], [9,13], [12,14]])
559
560        limit_minmod_kurganov_ext(quantity)
561
562
563        #print quantity.vertex_values
564
565
566    def test_limit_vanleer(self):
567        quantity = Quantity(self.domain4)
568
569        #Create a deliberate overshoot (e.g. from gradient computation)
570        quantity.set_values([[0,0], [2,10], [9,13], [12,14]])
571
572        #Limit
573        quantity.set_limiter('vanleer')
574        quantity.limit_range()
575
576
577        #print quantity.vertex_values
578
579        assert allclose(quantity.vertex_values, [[ -2.4,          2.4       ],
580                                                 [  2.32653061,   9.67346939],
581                                                 [ 10.39393939,  11.60606061],
582                                                 [ 11.4,         14.6       ]] )
583
584        from anuga_1d.base.quantity_ext import limit_vanleer_ext
585
586        quantity.set_values([[0,0], [2,10], [9,13], [12,14]])
587
588        limit_vanleer_ext(quantity)
589
590
591        #print quantity.vertex_values
592
593        assert allclose(quantity.vertex_values, [[ -2.4,          2.4       ],
594                                                 [  2.32653061,   9.67346939],
595                                                 [ 10.39393939,  11.60606061],
596                                                 [ 11.4,         14.6       ]] )
597#
598
599    def test_find_qmax_qmin(self):
600        quantity = Quantity(self.domain4)
601
602
603        quantity.set_values([1.0, 4.0, 8.0, 2.0],
604                            location = 'centroids')
605
606
607
608        quantity.find_qmax_qmin()
609
610
611        assert allclose(quantity.qmax, [4.0, 8.0, 8.0, 8.0])
612        assert allclose(quantity.qmin, [1.0, 1.0, 2.0, 2.0])
613
614
615    def test_distribute_first_order(self):
616        quantity = Quantity(self.domain4)
617
618        #Test centroids
619        centroid_values = array([1.,2.,3.,4.])
620        quantity.set_values(centroid_values, location = 'centroids')
621        assert allclose(quantity.centroid_values, centroid_values) #Centroid
622
623
624        #Extrapolate from centroid to vertices and edges
625        quantity.extrapolate_first_order()
626       
627        assert allclose(quantity.vertex_values,[[ 1.,  1.],
628                                                [ 2.,  2.],
629                                                [ 3.,  3.],
630                                                [ 4.,  4.]])
631 
632
633
634    def test_distribute_second_order(self):
635        quantity = Quantity(self.domain4)
636
637        #Test centroids
638        centroid_values = array([2.,4.,8.,2.])
639        quantity.set_values(centroid_values, location = 'centroids')
640        assert allclose(quantity.centroid_values, centroid_values) #Centroid
641
642
643        #Extrapolate
644        quantity.beta = 2.0
645        quantity.limiter = "minmod_kurganov"
646        quantity.extrapolate_second_order()
647
648        assert allclose(quantity.vertex_values, [[ 1.2,  2.8],
649                                                 [ 2.,  6. ],
650                                                 [ 8.,   8. ],
651                                                 [ 6.8, -2.8]])
652
653
654    def test_update_explicit(self):
655        quantity = Quantity(self.domain4)
656
657        #Test centroids
658        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
659        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
660
661        #Set explicit_update
662        quantity.explicit_update = array( [1.,1.,1.,1.] )
663
664        #Update with given timestep
665        quantity.update(0.1)
666
667        x = array([1, 2, 3, 4]) + array( [.1,.1,.1,.1] )
668        assert allclose( quantity.centroid_values, x)
669
670    def test_update_semi_implicit(self):
671        quantity = Quantity(self.domain4)
672
673        #Test centroids
674        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
675        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
676
677        #Set semi implicit update
678        quantity.semi_implicit_update = array([1.,1.,1.,1.])
679
680        #Update with given timestep
681        timestep = 0.1
682        quantity.update(timestep)
683
684        sem = array([1.,1.,1.,1.])
685        denom = ones(4, numpy.float)-timestep*sem
686
687        x = array([1, 2, 3, 4])/denom
688        assert allclose( quantity.centroid_values, x)
689
690
691    def test_both_updates(self):
692        quantity = Quantity(self.domain4)
693
694        #Test centroids
695        centroid_values = array( [1, 2, 3, 4] )
696        quantity.set_values(centroid_values, location = 'centroids')
697        assert allclose(quantity.centroid_values, centroid_values) #Centroid
698
699        #Set explicit_update
700        explicit_update = array( [4.,3.,2.,1.] )
701        quantity.explicit_update[:,] = explicit_update
702
703        #Set semi implicit update
704        semi_implicit_update = array( [1.,1.,1.,1.] )
705        quantity.semi_implicit_update[:,] = semi_implicit_update
706
707        #Update with given timestep
708        timestep = 0.1
709        quantity.update(0.1)
710
711        denom = 1.0-timestep*semi_implicit_update
712        x = (centroid_values + timestep*explicit_update)/denom
713 
714        assert allclose( quantity.centroid_values, x)
715
716
717 
718#-------------------------------------------------------------
719if __name__ == "__main__":
720    suite = unittest.makeSuite(Test_Quantity, 'test')
721    runner = unittest.TextTestRunner()
722    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.