source: anuga_work/development/2010-projects/anuga_1d/generic/test_quantity.py @ 7840

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

Changing name of 1d projects so that it will be easy to moveto the trunk
File size: 18.4 KB
Line 
1import quantity
2#!/usr/bin/env python
3
4import unittest
5from math import sqrt, pi
6
7
8from generic_domain import Generic_domain as Domain
9#from shallow_water_domain import flux_function as domain_flux_function
10
11from 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
40    def tearDown(self):
41        pass
42        #print "  Tearing down"
43
44
45    def test_creat_with_boundary(self):
46
47        assert self.domain4.boundary  == {(0, 0): 'left', (3, 1): 'right'}
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_distribute_first_order(self):
500        quantity = Quantity(self.domain4)
501
502        #Test centroids
503        centroid_values = array([1.,2.,3.,4.])
504        quantity.set_values(centroid_values, location = 'centroids')
505        assert allclose(quantity.centroid_values, centroid_values) #Centroid
506
507
508        #Extrapolate from centroid to vertices and edges
509        quantity.extrapolate_first_order()
510       
511        assert allclose(quantity.vertex_values,[[ 1.,  1.],
512                                                [ 2.,  2.],
513                                                [ 3.,  3.],
514                                                [ 4.,  4.]])
515 
516
517
518    def test_distribute_second_order(self):
519        quantity = Quantity(self.domain4)
520
521        #Test centroids
522        centroid_values = array([2.,4.,8.,2.])
523        quantity.set_values(centroid_values, location = 'centroids')
524        assert allclose(quantity.centroid_values, centroid_values) #Centroid
525
526
527        #Extrapolate
528        quantity.extrapolate_second_order()
529
530        assert allclose(quantity.vertex_values, [[ 1.2,  2.8],
531                                                 [ 2.,  6. ],
532                                                 [ 8.,   8. ],
533                                                 [ 6.8, -2.8]])
534
535
536    def test_update_explicit(self):
537        quantity = Quantity(self.domain4)
538
539        #Test centroids
540        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
541        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
542
543        #Set explicit_update
544        quantity.explicit_update = array( [1.,1.,1.,1.] )
545
546        #Update with given timestep
547        quantity.update(0.1)
548
549        x = array([1, 2, 3, 4]) + array( [.1,.1,.1,.1] )
550        assert allclose( quantity.centroid_values, x)
551
552    def test_update_semi_implicit(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 semi implicit update
560        quantity.semi_implicit_update = array([1.,1.,1.,1.])
561
562        #Update with given timestep
563        timestep = 0.1
564        quantity.update(timestep)
565
566        sem = array([1.,1.,1.,1.])
567        denom = ones(4, numpy.float)-timestep*sem
568
569        x = array([1, 2, 3, 4])/denom
570        assert allclose( quantity.centroid_values, x)
571
572
573    def test_both_updates(self):
574        quantity = Quantity(self.domain4)
575
576        #Test centroids
577        centroid_values = array( [1, 2, 3, 4] )
578        quantity.set_values(centroid_values, location = 'centroids')
579        assert allclose(quantity.centroid_values, centroid_values) #Centroid
580
581        #Set explicit_update
582        explicit_update = array( [4.,3.,2.,1.] )
583        quantity.explicit_update[:,] = explicit_update
584
585        #Set semi implicit update
586        semi_implicit_update = array( [1.,1.,1.,1.] )
587        quantity.semi_implicit_update[:,] = semi_implicit_update
588
589        #Update with given timestep
590        timestep = 0.1
591        quantity.update(0.1)
592
593        denom = 1.0-timestep*semi_implicit_update
594        x = (centroid_values + timestep*explicit_update)/denom
595 
596        assert allclose( quantity.centroid_values, x)
597
598
599 
600#-------------------------------------------------------------
601if __name__ == "__main__":
602    suite = unittest.makeSuite(Test_Quantity, 'test')
603    runner = unittest.TextTestRunner()
604    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.