source: inundation/ga/storm_surge/pyvolution/test_quantity.py @ 333

Last change on this file since 333 was 324, checked in by ole, 20 years ago

Played with file boundary.

Added set_vertex_values w/ Duncan
File size: 18.5 KB
Line 
1#!/usr/bin/env python
2
3import unittest
4from math import sqrt, pi
5
6
7from quantity import *
8from config import epsilon
9from Numeric import allclose, array
10
11       
12class TestCase(unittest.TestCase):
13    def setUp(self):
14        from domain import Domain
15       
16        a = [0.0, 0.0]
17        b = [0.0, 2.0]
18        c = [2.0,0.0]
19        d = [0.0, 4.0]
20        e = [2.0, 2.0]
21        f = [4.0,0.0]
22
23        points = [a, b, c, d, e, f]
24
25        #bac, bce, ecf, dbe
26        elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
27
28        self.mesh1 = Domain(points[:3], [elements[0]])       
29        self.mesh1.check_integrity()
30       
31        self.mesh4 = Domain(points, elements)       
32        self.mesh4.check_integrity()
33       
34    def tearDown(self):
35        pass
36        #print "  Tearing down"
37
38
39    def test_creation(self):
40       
41        quantity = Quantity(self.mesh1, [[1,2,3]])
42        assert allclose(quantity.vertex_values, [[1.,2.,3.]])
43
44        try:
45            quantity = Quantity()
46        except:
47            pass
48        else:
49            raise 'Should have raised empty quantity exception'
50
51
52        try:
53            quantity = Quantity([1,2,3])
54        except AssertionError:
55            pass
56        except:
57            raise 'Should have raised "mising mesh object" error'       
58       
59
60    def test_creation_zeros(self):
61       
62        quantity = Quantity(self.mesh1)
63        assert allclose(quantity.vertex_values, [[0.,0.,0.]])
64
65
66        quantity = Quantity(self.mesh4)
67        assert allclose(quantity.vertex_values, [[0.,0.,0.], [0.,0.,0.],
68                                                 [0.,0.,0.], [0.,0.,0.]])       
69
70       
71    def test_interpolation(self):
72        quantity = Quantity(self.mesh1, [[1,2,3]])
73        assert allclose(quantity.centroid_values, [2.0]) #Centroid
74
75        assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5]]) 
76
77
78    def test_interpolation2(self):
79        quantity = Conserved_quantity(self.mesh4,
80                            [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
81        assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
82
83
84        quantity.extrapolate_second_order()
85
86        assert allclose(quantity.vertex_values, [[2., 2., 2.],
87                                               [3.+2./3, 6.+2./3, 4.+2./3],
88                                               [7.5, 0.5, 1.],
89                                               [-5, -2.5, 7.5]])
90       
91        #print quantity.edge_values       
92        assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
93                                               [5., 5., 5.],
94                                               [4.5, 4.5, 0.],
95                                               [3.0, -1.5, -1.5]])
96
97    def test_boundary_allocation(self):
98        quantity = Conserved_quantity(self.mesh4,
99                            [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
100
101        assert quantity.boundary_values.shape[0] == len(self.mesh4.boundary)
102
103
104    def test_set_values(self):
105        quantity = Quantity(self.mesh4)
106
107
108        quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
109                            location = 'vertices')
110        assert allclose(quantity.vertex_values,
111                        [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
112        assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
113        assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
114                                               [5., 5., 5.],
115                                               [4.5, 4.5, 0.],
116                                               [3.0, -1.5, -1.5]])
117
118
119        #Test default
120        quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
121        assert allclose(quantity.vertex_values,
122                        [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])       
123        assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
124        assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
125                                               [5., 5., 5.],
126                                               [4.5, 4.5, 0.],
127                                               [3.0, -1.5, -1.5]])
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 edges
134        quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
135                            location = 'edges')
136        assert allclose(quantity.edge_values,
137                        [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])       
138
139        #Test exceptions
140        try:
141            quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
142                                location = 'bas kamel tuba')
143        except:
144            pass
145
146
147        try:
148            quantity.set_values([[1,2,3], [0,0,9]])
149        except AssertionError:
150            pass
151        except:
152            raise 'should have raised Assertionerror'
153
154
155
156    def test_set_values_const(self):
157        quantity = Quantity(self.mesh4)
158
159        quantity.set_values(1.0, location = 'vertices')
160        assert allclose(quantity.vertex_values,
161                        [[1,1,1], [1,1,1], [1,1,1], [1, 1, 1]])
162        assert allclose(quantity.centroid_values, [1, 1, 1, 1]) #Centroid
163        assert allclose(quantity.edge_values, [[1, 1, 1],
164                                               [1, 1, 1],
165                                               [1, 1, 1],
166                                               [1, 1, 1]])
167
168
169        quantity.set_values(2.0, location = 'centroids')
170        assert allclose(quantity.centroid_values, [2, 2, 2, 2])
171
172        quantity.set_values(3.0, location = 'edges')
173        assert allclose(quantity.edge_values, [[3, 3, 3],
174                                               [3, 3, 3],
175                                               [3, 3, 3],
176                                               [3, 3, 3]])       
177
178
179    def test_set_values_func(self):
180        quantity = Quantity(self.mesh4)
181
182        def f(x, y):
183            return x+y
184
185        quantity.set_values(f, location = 'vertices')
186        assert allclose(quantity.vertex_values,
187                        [[2,0,2], [2,2,4], [4,2,4], [4,2,4]])       
188        assert allclose(quantity.centroid_values,
189                        [4.0/3, 8.0/3, 10.0/3, 10.0/3])
190        assert allclose(quantity.edge_values,
191                        [[1,2,1], [3,3,2], [3,4,3], [3,4,3]])               
192
193       
194        quantity.set_values(f, location = 'centroids')
195        assert allclose(quantity.centroid_values,
196                        [4.0/3, 8.0/3, 10.0/3, 10.0/3])       
197
198
199
200
201    def test_set_vertex_values(self):
202        quantity = Quantity(self.mesh4)
203
204
205        quantity.set_vertex_values([0,1,2,3,4,5])
206
207
208        assert allclose(quantity.vertex_values,
209                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
210       
211        assert allclose(quantity.centroid_values, [1., 7./3, 11./3, 8./3]) #Centroid
212       
213        assert allclose(quantity.edge_values, [[1., 1.5, 0.5],
214                                               [3., 2.5, 1.5],
215                                               [3.5, 4.5, 3.],
216                                               [2.5, 3.5, 2]])
217
218
219
220    def test_gradient(self):
221        quantity = Conserved_quantity(self.mesh4)
222
223        #Set up for a gradient of (3,0) at mid triangle
224        quantity.set_values([2.0, 4.0, 8.0, 2.0],
225                            location = 'centroids')
226
227
228
229        #Gradients
230        a, b = quantity.compute_gradients()
231
232        #gradient bewteen t0 and t1 is undefined as det==0
233        assert a[0] == 0.0
234        assert b[0] == 0.0
235        #The others are OK
236        for i in range(1,4):
237            assert a[i] == 3.0
238            assert b[i] == 0.0
239
240
241        quantity.extrapolate_second_order()
242
243        #print quantity.vertex_values
244        assert allclose(quantity.vertex_values, [[2., 2.,  2.],
245                                                 [0., 6.,  6.],
246                                                 [6., 6., 12.],
247                                                 [0., 0.,  6.]])
248
249               
250    def test_second_order_extrapolation2(self):
251        quantity = Conserved_quantity(self.mesh4)       
252
253        #Set up for a gradient of (3,1), f(x) = 3x+y
254        quantity.set_values([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3],
255                            location = 'centroids')
256       
257        #Gradients
258        a, b = quantity.compute_gradients()
259
260        #gradient bewteen t0 and t1 is undefined as det==0
261        assert a[0] == 0.0
262        assert b[0] == 0.0
263        #The others are OK
264        for i in range(1,4):
265            assert allclose(a[i], 3.0)
266            assert allclose(b[i], 1.0)
267
268        quantity.extrapolate_second_order()
269       
270        #print quantity.vertex_values
271        assert allclose(quantity.vertex_values[1,0], 2.0)
272        assert allclose(quantity.vertex_values[1,1], 6.0)       
273        assert allclose(quantity.vertex_values[1,2], 8.0)
274
275
276
277    def test_first_order_extrapolator(self):
278        quantity = Conserved_quantity(self.mesh4)               
279
280        #Test centroids
281        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
282        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
283
284        #Extrapolate
285        quantity.extrapolate_first_order()
286
287        #Check vertices but not edge values
288        assert allclose(quantity.vertex_values,
289                        [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
290
291
292    def test_second_order_extrapolator(self):
293        quantity = Conserved_quantity(self.mesh4)                       
294
295        #Set up for a gradient of (3,0) at mid triangle
296        quantity.set_values([2.0, 4.0, 8.0, 2.0],
297                            location = 'centroids')
298       
299
300
301        quantity.extrapolate_second_order()
302        quantity.limit()
303
304
305        #Assert that central triangle is limited by neighbours
306        assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0]
307        assert quantity.vertex_values[1,0] >= quantity.vertex_values[3,1]
308       
309        assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1]
310        assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2]
311       
312        assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0]
313        assert quantity.vertex_values[1,2] >= quantity.vertex_values[3,1]
314
315       
316        #Assert that quantities are conserved
317        from Numeric import sum
318        for k in range(quantity.centroid_values.shape[0]):
319            assert allclose (quantity.centroid_values[k],
320                             sum(quantity.vertex_values[k,:])/3)
321       
322       
323
324       
325
326    def test_limiter(self):
327        quantity = Conserved_quantity(self.mesh4)
328
329        #Create a deliberate overshoot (e.g. from gradient computation)
330        quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
331
332
333        #Limit
334        quantity.limit()
335
336        #Assert that central triangle is limited by neighbours
337        assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0]
338        assert quantity.vertex_values[1,0] <= quantity.vertex_values[3,1]
339       
340        assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1]
341        assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2]
342       
343        assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0]
344        assert quantity.vertex_values[1,2] <= quantity.vertex_values[3,1]
345
346
347       
348        #Assert that quantities are conserved
349        from Numeric import sum
350        for k in range(quantity.centroid_values.shape[0]):
351            assert allclose (quantity.centroid_values[k],
352                             sum(quantity.vertex_values[k,:])/3)
353       
354
355
356    def test_distribute_first_order(self):
357        quantity = Conserved_quantity(self.mesh4)       
358
359        #Test centroids
360        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
361        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
362
363
364        #Extrapolate
365        quantity.extrapolate_first_order()
366
367        #Interpolate
368        quantity.interpolate_from_vertices_to_edges()       
369
370        assert allclose(quantity.vertex_values,
371                        [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
372        assert allclose(quantity.edge_values, [[1,1,1], [2,2,2],
373                                               [3,3,3], [4, 4, 4]])
374
375
376
377    def test_update_explicit(self):
378        quantity = Conserved_quantity(self.mesh4)
379
380        #Test centroids
381        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
382        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
383
384        #Set explicit_update
385        quantity.explicit_update = array( [1.,1.,1.,1.] )
386
387        #Update with given timestep
388        quantity.update(0.1)
389
390        x = array([1, 2, 3, 4]) + array( [.1,.1,.1,.1] )
391        assert allclose( quantity.centroid_values, x)
392
393    def test_update_semi_implicit(self):
394        quantity = Conserved_quantity(self.mesh4)
395
396        #Test centroids
397        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
398        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
399
400        #Set semi implicit update
401        quantity.semi_implicit_update = array( [1.,1.,1.,1.] )
402
403        #Update with given timestep
404        quantity.update(0.1)
405
406        x = array([1, 2, 3, 4])/array( [.9,.9,.9,.9] )
407        assert allclose( quantity.centroid_values, x)
408
409    def test_both_updates(self):
410        quantity = Conserved_quantity(self.mesh4)
411
412        #Test centroids
413        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
414        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
415
416        #Set explicit_update
417        quantity.explicit_update = array( [4.,3.,2.,1.] )
418       
419        #Set semi implicit update
420        quantity.semi_implicit_update = array( [1.,1.,1.,1.] )
421
422        #Update with given timestep
423        quantity.update(0.1)
424
425        x = array([1, 2, 3, 4]) + array( [.4,.3,.2,.1] )
426        x /= array( [.9,.9,.9,.9] )
427        assert allclose( quantity.centroid_values, x)       
428
429       
430
431
432    #Test smoothing   
433    def test_smoothing(self):
434
435        from mesh_factory import rectangular
436        from shallow_water import Domain, Transmissive_boundary
437        from Numeric import zeros, Float
438        from util import mean
439
440        #Create basic mesh
441        points, vertices, boundary = rectangular(2, 2)
442
443        #Create shallow water domain
444        domain = Domain(points, vertices, boundary)
445        domain.default_order=2
446        domain.reduction = mean
447
448       
449        #Set some field values
450        domain.set_quantity('elevation', lambda x,y: x)       
451        domain.set_quantity('friction', 0.03)
452
453
454        ######################
455        # Boundary conditions
456        B = Transmissive_boundary(domain)
457        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
458       
459
460        ######################
461        #Initial condition - with jumps
462
463        bed = domain.quantities['elevation'].vertex_values
464        level = zeros(bed.shape, Float)
465
466        h = 0.03
467        for i in range(level.shape[0]):
468            if i % 2 == 0:           
469                level[i,:] = bed[i,:] + h
470            else:
471                level[i,:] = bed[i,:]
472               
473        domain.set_quantity('level', level)
474
475        level = domain.quantities['level']
476
477        #Get smoothed level
478        A, V = level.get_vertex_values(xy=False, smooth=True)
479        Q = level.vertex_values
480
481       
482        assert A.shape[0] == 9 
483        assert V.shape[0] == 8 
484        assert V.shape[1] == 3                 
485       
486        #First four points
487        assert allclose(A[0], (Q[0,2] + Q[1,1])/2) 
488        assert allclose(A[1], (Q[1,0] + Q[3,1] + Q[2,2])/3)
489        assert allclose(A[2], Q[3,0])
490        assert allclose(A[3], (Q[0,0] + Q[5,1] + Q[4,2])/3)
491
492        #Center point
493        assert allclose(A[4], (Q[0,1] + Q[1,2] + Q[2,0] +\
494                               Q[5,0] + Q[6,2] + Q[7,1])/6)
495                                 
496
497        #Check V
498        assert allclose(V[0,:], [3,4,0])
499        assert allclose(V[1,:], [1,0,4])
500        assert allclose(V[2,:], [4,5,1])                       
501        assert allclose(V[3,:], [2,1,5])
502        assert allclose(V[4,:], [6,7,3])       
503        assert allclose(V[5,:], [4,3,7])
504        assert allclose(V[6,:], [7,8,4])
505        assert allclose(V[7,:], [5,4,8])                       
506
507        #Get smoothed level with XY
508        X, Y, A1, V1 = level.get_vertex_values(xy=True, smooth=True)
509
510        assert allclose(A, A1)
511        assert allclose(V, V1)       
512
513        #Check XY
514        assert allclose(X[4], 0.5)
515        assert allclose(Y[4], 0.5)
516
517        assert allclose(X[7], 1.0)
518        assert allclose(Y[7], 0.5)                 
519
520
521
522
523    def test_vertex_values_no_smoothing(self):
524
525        from mesh_factory import rectangular
526        from shallow_water import Domain, Transmissive_boundary
527        from Numeric import zeros, Float
528        from util import mean
529
530       
531        #Create basic mesh
532        points, vertices, boundary = rectangular(2, 2)
533
534        #Create shallow water domain
535        domain = Domain(points, vertices, boundary)
536        domain.default_order=2
537        domain.reduction = mean
538
539       
540        #Set some field values
541        domain.set_quantity('elevation', lambda x,y: x)       
542        domain.set_quantity('friction', 0.03)
543
544
545        ######################
546        #Initial condition - with jumps
547
548        bed = domain.quantities['elevation'].vertex_values
549        level = zeros(bed.shape, Float)
550
551        h = 0.03
552        for i in range(level.shape[0]):
553            if i % 2 == 0:           
554                level[i,:] = bed[i,:] + h
555            else:
556                level[i,:] = bed[i,:]
557               
558        domain.set_quantity('level', level)
559
560        #Get level
561        level = domain.quantities['level']       
562        A, V = level.get_vertex_values(xy=False, smooth=False)
563        Q = level.vertex_values.flat
564
565        for k in range(8):
566            assert allclose(A[k], Q[k])
567
568           
569        for k in range(8):
570            assert V[k, 0] == 3*k
571            assert V[k, 1] == 3*k+1
572            assert V[k, 2] == 3*k+2           
573           
574
575
576        X, Y, A1, V1 = level.get_vertex_values(xy=True, smooth=False)
577
578       
579        assert allclose(A, A1)
580        assert allclose(V, V1)       
581
582        #Check XY
583        assert allclose(X[1], 0.5)
584        assert allclose(Y[1], 0.5)
585        assert allclose(X[4], 0.0)
586        assert allclose(Y[4], 0.0)
587        assert allclose(X[12], 1.0)
588        assert allclose(Y[12], 0.0)               
589       
590       
591                                 
592
593       
594#-------------------------------------------------------------
595if __name__ == "__main__":
596    suite = unittest.makeSuite(TestCase,'test')
597    runner = unittest.TextTestRunner()
598    runner.run(suite)
599
600
Note: See TracBrowser for help on using the repository browser.