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

Last change on this file since 321 was 297, checked in by ole, 21 years ago

Adapted original pmesh_to_domain
File size: 17.9 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    def test_gradient(self):
200        quantity = Conserved_quantity(self.mesh4)
201
202        #Set up for a gradient of (3,0) at mid triangle
203        quantity.set_values([2.0, 4.0, 8.0, 2.0],
204                            location = 'centroids')
205
206
207
208        #Gradients
209        a, b = quantity.compute_gradients()
210
211        #gradient bewteen t0 and t1 is undefined as det==0
212        assert a[0] == 0.0
213        assert b[0] == 0.0
214        #The others are OK
215        for i in range(1,4):
216            assert a[i] == 3.0
217            assert b[i] == 0.0
218
219
220        quantity.extrapolate_second_order()
221
222        #print quantity.vertex_values
223        assert allclose(quantity.vertex_values, [[2., 2.,  2.],
224                                                 [0., 6.,  6.],
225                                                 [6., 6., 12.],
226                                                 [0., 0.,  6.]])
227
228               
229    def test_second_order_extrapolation2(self):
230        quantity = Conserved_quantity(self.mesh4)       
231
232        #Set up for a gradient of (3,1), f(x) = 3x+y
233        quantity.set_values([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3],
234                            location = 'centroids')
235       
236        #Gradients
237        a, b = quantity.compute_gradients()
238
239        #gradient bewteen t0 and t1 is undefined as det==0
240        assert a[0] == 0.0
241        assert b[0] == 0.0
242        #The others are OK
243        for i in range(1,4):
244            assert allclose(a[i], 3.0)
245            assert allclose(b[i], 1.0)
246
247        quantity.extrapolate_second_order()
248       
249        #print quantity.vertex_values
250        assert allclose(quantity.vertex_values[1,0], 2.0)
251        assert allclose(quantity.vertex_values[1,1], 6.0)       
252        assert allclose(quantity.vertex_values[1,2], 8.0)
253
254
255
256    def test_first_order_extrapolator(self):
257        quantity = Conserved_quantity(self.mesh4)               
258
259        #Test centroids
260        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
261        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
262
263        #Extrapolate
264        quantity.extrapolate_first_order()
265
266        #Check vertices but not edge values
267        assert allclose(quantity.vertex_values,
268                        [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
269
270
271    def test_second_order_extrapolator(self):
272        quantity = Conserved_quantity(self.mesh4)                       
273
274        #Set up for a gradient of (3,0) at mid triangle
275        quantity.set_values([2.0, 4.0, 8.0, 2.0],
276                            location = 'centroids')
277       
278
279
280        quantity.extrapolate_second_order()
281        quantity.limit()
282
283
284        #Assert that central triangle is limited by neighbours
285        assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0]
286        assert quantity.vertex_values[1,0] >= quantity.vertex_values[3,1]
287       
288        assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1]
289        assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2]
290       
291        assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0]
292        assert quantity.vertex_values[1,2] >= quantity.vertex_values[3,1]
293
294       
295        #Assert that quantities are conserved
296        from Numeric import sum
297        for k in range(quantity.centroid_values.shape[0]):
298            assert allclose (quantity.centroid_values[k],
299                             sum(quantity.vertex_values[k,:])/3)
300       
301       
302
303       
304
305    def test_limiter(self):
306        quantity = Conserved_quantity(self.mesh4)
307
308        #Create a deliberate overshoot (e.g. from gradient computation)
309        quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
310
311
312        #Limit
313        quantity.limit()
314
315        #Assert that central triangle is limited by neighbours
316        assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0]
317        assert quantity.vertex_values[1,0] <= quantity.vertex_values[3,1]
318       
319        assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1]
320        assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2]
321       
322        assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0]
323        assert quantity.vertex_values[1,2] <= quantity.vertex_values[3,1]
324
325
326       
327        #Assert that quantities are conserved
328        from Numeric import sum
329        for k in range(quantity.centroid_values.shape[0]):
330            assert allclose (quantity.centroid_values[k],
331                             sum(quantity.vertex_values[k,:])/3)
332       
333
334
335    def test_distribute_first_order(self):
336        quantity = Conserved_quantity(self.mesh4)       
337
338        #Test centroids
339        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
340        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
341
342
343        #Extrapolate
344        quantity.extrapolate_first_order()
345
346        #Interpolate
347        quantity.interpolate_from_vertices_to_edges()       
348
349        assert allclose(quantity.vertex_values,
350                        [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
351        assert allclose(quantity.edge_values, [[1,1,1], [2,2,2],
352                                               [3,3,3], [4, 4, 4]])
353
354
355
356    def test_update_explicit(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        #Set explicit_update
364        quantity.explicit_update = array( [1.,1.,1.,1.] )
365
366        #Update with given timestep
367        quantity.update(0.1)
368
369        x = array([1, 2, 3, 4]) + array( [.1,.1,.1,.1] )
370        assert allclose( quantity.centroid_values, x)
371
372    def test_update_semi_implicit(self):
373        quantity = Conserved_quantity(self.mesh4)
374
375        #Test centroids
376        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
377        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
378
379        #Set semi implicit update
380        quantity.semi_implicit_update = array( [1.,1.,1.,1.] )
381
382        #Update with given timestep
383        quantity.update(0.1)
384
385        x = array([1, 2, 3, 4])/array( [.9,.9,.9,.9] )
386        assert allclose( quantity.centroid_values, x)
387
388    def test_both_updates(self):
389        quantity = Conserved_quantity(self.mesh4)
390
391        #Test centroids
392        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
393        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
394
395        #Set explicit_update
396        quantity.explicit_update = array( [4.,3.,2.,1.] )
397       
398        #Set semi implicit update
399        quantity.semi_implicit_update = array( [1.,1.,1.,1.] )
400
401        #Update with given timestep
402        quantity.update(0.1)
403
404        x = array([1, 2, 3, 4]) + array( [.4,.3,.2,.1] )
405        x /= array( [.9,.9,.9,.9] )
406        assert allclose( quantity.centroid_values, x)       
407
408       
409
410
411    #Test smoothing   
412    def test_smoothing(self):
413
414        from mesh_factory import rectangular
415        from shallow_water import Domain, Transmissive_boundary
416        from Numeric import zeros, Float
417        from util import mean
418
419        #Create basic mesh
420        points, vertices, boundary = rectangular(2, 2)
421
422        #Create shallow water domain
423        domain = Domain(points, vertices, boundary)
424        domain.default_order=2
425        domain.reduction = mean
426
427       
428        #Set some field values
429        domain.set_quantity('elevation', lambda x,y: x)       
430        domain.set_quantity('friction', 0.03)
431
432
433        ######################
434        # Boundary conditions
435        B = Transmissive_boundary(domain)
436        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
437       
438
439        ######################
440        #Initial condition - with jumps
441
442        bed = domain.quantities['elevation'].vertex_values
443        level = zeros(bed.shape, Float)
444
445        h = 0.03
446        for i in range(level.shape[0]):
447            if i % 2 == 0:           
448                level[i,:] = bed[i,:] + h
449            else:
450                level[i,:] = bed[i,:]
451               
452        domain.set_quantity('level', level)
453
454        level = domain.quantities['level']
455
456        #Get smoothed level
457        A, V = level.get_vertex_values(xy=False, smooth=True)
458        Q = level.vertex_values
459
460       
461        assert A.shape[0] == 9 
462        assert V.shape[0] == 8 
463        assert V.shape[1] == 3                 
464       
465        #First four points
466        assert allclose(A[0], (Q[0,2] + Q[1,1])/2) 
467        assert allclose(A[1], (Q[1,0] + Q[3,1] + Q[2,2])/3)
468        assert allclose(A[2], Q[3,0])
469        assert allclose(A[3], (Q[0,0] + Q[5,1] + Q[4,2])/3)
470
471        #Center point
472        assert allclose(A[4], (Q[0,1] + Q[1,2] + Q[2,0] +\
473                               Q[5,0] + Q[6,2] + Q[7,1])/6)
474                                 
475
476        #Check V
477        assert allclose(V[0,:], [3,4,0])
478        assert allclose(V[1,:], [1,0,4])
479        assert allclose(V[2,:], [4,5,1])                       
480        assert allclose(V[3,:], [2,1,5])
481        assert allclose(V[4,:], [6,7,3])       
482        assert allclose(V[5,:], [4,3,7])
483        assert allclose(V[6,:], [7,8,4])
484        assert allclose(V[7,:], [5,4,8])                       
485
486        #Get smoothed level with XY
487        X, Y, A1, V1 = level.get_vertex_values(xy=True, smooth=True)
488
489        assert allclose(A, A1)
490        assert allclose(V, V1)       
491
492        #Check XY
493        assert allclose(X[4], 0.5)
494        assert allclose(Y[4], 0.5)
495
496        assert allclose(X[7], 1.0)
497        assert allclose(Y[7], 0.5)                 
498
499
500
501
502    def test_vertex_values_no_smoothing(self):
503
504        from mesh_factory import rectangular
505        from shallow_water import Domain, Transmissive_boundary
506        from Numeric import zeros, Float
507        from util import mean
508
509       
510        #Create basic mesh
511        points, vertices, boundary = rectangular(2, 2)
512
513        #Create shallow water domain
514        domain = Domain(points, vertices, boundary)
515        domain.default_order=2
516        domain.reduction = mean
517
518       
519        #Set some field values
520        domain.set_quantity('elevation', lambda x,y: x)       
521        domain.set_quantity('friction', 0.03)
522
523
524        ######################
525        #Initial condition - with jumps
526
527        bed = domain.quantities['elevation'].vertex_values
528        level = zeros(bed.shape, Float)
529
530        h = 0.03
531        for i in range(level.shape[0]):
532            if i % 2 == 0:           
533                level[i,:] = bed[i,:] + h
534            else:
535                level[i,:] = bed[i,:]
536               
537        domain.set_quantity('level', level)
538
539        #Get level
540        level = domain.quantities['level']       
541        A, V = level.get_vertex_values(xy=False, smooth=False)
542        Q = level.vertex_values.flat
543
544        for k in range(8):
545            assert allclose(A[k], Q[k])
546
547           
548        for k in range(8):
549            assert V[k, 0] == 3*k
550            assert V[k, 1] == 3*k+1
551            assert V[k, 2] == 3*k+2           
552           
553
554
555        X, Y, A1, V1 = level.get_vertex_values(xy=True, smooth=False)
556
557       
558        assert allclose(A, A1)
559        assert allclose(V, V1)       
560
561        #Check XY
562        assert allclose(X[1], 0.5)
563        assert allclose(Y[1], 0.5)
564        assert allclose(X[4], 0.0)
565        assert allclose(Y[4], 0.0)
566        assert allclose(X[12], 1.0)
567        assert allclose(Y[12], 0.0)               
568       
569       
570                                 
571
572       
573#-------------------------------------------------------------
574if __name__ == "__main__":
575    suite = unittest.makeSuite(TestCase,'test')
576    runner = unittest.TextTestRunner()
577    runner.run(suite)
578
579
Note: See TracBrowser for help on using the repository browser.