source: inundation/pyvolution/test_domain.py @ 1916

Last change on this file since 1916 was 1916, checked in by ole, 18 years ago

Implemented operator overloading for (pow) in class Quantity and tested.
Wrote create_quantity_from_expression and test.

File size: 15.0 KB
Line 
1#!/usr/bin/env python
2
3import unittest
4from math import sqrt
5
6from domain import *
7from config import epsilon
8from Numeric import allclose, array, ones, Float
9
10
11def add_to_verts(tag, elements, domain):
12    if tag == "mound":
13        domain.test = "Mound"
14
15
16def set_bottom_friction(tag, elements, domain):
17    if tag == "bottom":
18        #print 'bottom - indices',elements
19        domain.set_quantity('friction', 0.09, indices = elements)
20
21def set_top_friction(tag, elements, domain):
22    if tag == "top":
23        #print 'top - indices',elements
24        domain.set_quantity('friction', 1., indices = elements)
25
26
27def set_all_friction(tag, elements, domain):
28    if tag == "all":
29        new_values = domain.get_quantity('friction', indices = elements) + 10.0
30
31        domain.set_quantity('friction', new_values, indices = elements)
32
33
34class Test_Domain(unittest.TestCase):
35    def setUp(self):
36        pass
37
38
39    def tearDown(self):
40        pass
41
42
43    def test_simple(self):
44        a = [0.0, 0.0]
45        b = [0.0, 2.0]
46        c = [2.0,0.0]
47        d = [0.0, 4.0]
48        e = [2.0, 2.0]
49        f = [4.0,0.0]
50
51        points = [a, b, c, d, e, f]
52        #bac, bce, ecf, dbe, daf, dae
53        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
54
55        conserved_quantities = ['stage', 'xmomentum', 'ymomentum']
56        other_quantities = ['elevation', 'friction']
57
58        domain = Domain(points, vertices, None,
59                        conserved_quantities, other_quantities)
60        domain.check_integrity()
61
62        for name in conserved_quantities + other_quantities:
63            assert domain.quantities.has_key(name)
64
65
66        assert domain.get_conserved_quantities(0, edge=1) == 0.
67
68
69    def test_conserved_quantities(self):
70
71        a = [0.0, 0.0]
72        b = [0.0, 2.0]
73        c = [2.0,0.0]
74        d = [0.0, 4.0]
75        e = [2.0, 2.0]
76        f = [4.0,0.0]
77
78        points = [a, b, c, d, e, f]
79        #bac, bce, ecf, dbe, daf, dae
80        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
81
82        domain = Domain(points, vertices, boundary=None,
83                        conserved_quantities =\
84                        ['stage', 'xmomentum', 'ymomentum'])
85
86
87        domain.set_quantity('stage', [[1,2,3], [5,5,5],
88                                      [0,0,9], [-6, 3, 3]])
89
90        domain.set_quantity('xmomentum', [[1,2,3], [5,5,5],
91                                          [0,0,9], [-6, 3, 3]])
92
93        domain.check_integrity()
94
95        #Centroids
96        q = domain.get_conserved_quantities(0)
97        assert allclose(q, [2., 2., 0.])
98
99        q = domain.get_conserved_quantities(1)
100        assert allclose(q, [5., 5., 0.])
101
102        q = domain.get_conserved_quantities(2)
103        assert allclose(q, [3., 3., 0.])
104
105        q = domain.get_conserved_quantities(3)
106        assert allclose(q, [0., 0., 0.])
107
108
109        #Edges
110        q = domain.get_conserved_quantities(0, edge=0)
111        assert allclose(q, [2.5, 2.5, 0.])
112        q = domain.get_conserved_quantities(0, edge=1)
113        assert allclose(q, [2., 2., 0.])
114        q = domain.get_conserved_quantities(0, edge=2)
115        assert allclose(q, [1.5, 1.5, 0.])
116
117        for i in range(3):
118            q = domain.get_conserved_quantities(1, edge=i)
119            assert allclose(q, [5, 5, 0.])
120
121
122        q = domain.get_conserved_quantities(2, edge=0)
123        assert allclose(q, [4.5, 4.5, 0.])
124        q = domain.get_conserved_quantities(2, edge=1)
125        assert allclose(q, [4.5, 4.5, 0.])
126        q = domain.get_conserved_quantities(2, edge=2)
127        assert allclose(q, [0., 0., 0.])
128
129
130        q = domain.get_conserved_quantities(3, edge=0)
131        assert allclose(q, [3., 3., 0.])
132        q = domain.get_conserved_quantities(3, edge=1)
133        assert allclose(q, [-1.5, -1.5, 0.])
134        q = domain.get_conserved_quantities(3, edge=2)
135        assert allclose(q, [-1.5, -1.5, 0.])
136
137
138
139    def test_create_quantity_from_expression(self):
140        """Quantity created from other quantities using arbitrary expression
141       
142        """
143
144
145        a = [0.0, 0.0]
146        b = [0.0, 2.0]
147        c = [2.0,0.0]
148        d = [0.0, 4.0]
149        e = [2.0, 2.0]
150        f = [4.0,0.0]
151
152        points = [a, b, c, d, e, f]
153        #bac, bce, ecf, dbe, daf, dae
154        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
155
156        domain = Domain(points, vertices, boundary=None,
157                        conserved_quantities =\
158                        ['stage', 'xmomentum', 'ymomentum'],
159                        other_quantities = ['elevation', 'friction'])
160
161
162        domain.set_quantity('elevation', -1)
163
164       
165        domain.set_quantity('stage', [[1,2,3], [5,5,5],
166                                      [0,0,9], [-6, 3, 3]])
167
168        domain.set_quantity('xmomentum', [[1,2,3], [5,5,5],
169                                          [0,0,9], [-6, 3, 3]])
170
171        domain.set_quantity('ymomentum', [[3,3,3], [4,2,1],
172                                          [2,4,-1], [1, 0, 1]])       
173
174        domain.check_integrity()
175
176
177
178        expression = 'stage - elevation'
179        Q = create_quantity_from_expression(domain, expression)
180
181        assert allclose(Q.vertex_values, [[2,3,4], [6,6,6],
182                                      [1,1,10], [-5, 4, 4]])
183
184        expression = '(xmomentum*xmomentum + ymomentum*ymomentum)**0.5'
185        Q = create_quantity_from_expression(domain, expression)
186
187        X = domain.quantities['xmomentum'].vertex_values
188        Y = domain.quantities['ymomentum'].vertex_values       
189
190        assert allclose(Q.vertex_values, (X**2 + Y**2)**0.5)
191                                     
192
193
194
195
196
197    def test_boundary_indices(self):
198
199        from config import default_boundary_tag
200
201
202        a = [0.0, 0.5]
203        b = [0.0, 0.0]
204        c = [0.5, 0.5]
205
206        points = [a, b, c]
207        vertices = [ [0,1,2] ]
208        domain = Domain(points, vertices)
209
210        domain.set_boundary( {default_boundary_tag: Dirichlet_boundary([5,2,1])} )
211
212
213        domain.check_integrity()
214
215        assert allclose(domain.neighbours, [[-1,-2,-3]])
216
217
218
219    def test_boundary_conditions(self):
220
221        a = [0.0, 0.0]
222        b = [0.0, 2.0]
223        c = [2.0,0.0]
224        d = [0.0, 4.0]
225        e = [2.0, 2.0]
226        f = [4.0,0.0]
227
228        points = [a, b, c, d, e, f]
229        #bac, bce, ecf, dbe
230        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
231        boundary = { (0, 0): 'First',
232                     (0, 2): 'First',
233                     (2, 0): 'Second',
234                     (2, 1): 'Second',
235                     (3, 1): 'Second',
236                     (3, 2): 'Second'}
237
238
239        domain = Domain(points, vertices, boundary,
240                        conserved_quantities =\
241                        ['stage', 'xmomentum', 'ymomentum'])
242        domain.check_integrity()
243
244
245
246        domain.set_quantity('stage', [[1,2,3], [5,5,5],
247                                      [0,0,9], [-6, 3, 3]])
248
249
250        domain.set_boundary( {'First': Dirichlet_boundary([5,2,1]),
251                              'Second': Transmissive_boundary(domain)} )
252
253        domain.update_boundary()
254
255        assert domain.quantities['stage'].boundary_values[0] == 5. #Dirichlet
256        assert domain.quantities['stage'].boundary_values[1] == 5. #Dirichlet
257        assert domain.quantities['stage'].boundary_values[2] ==\
258               domain.get_conserved_quantities(2, edge=0)[0] #Transmissive (4.5)
259        assert domain.quantities['stage'].boundary_values[3] ==\
260               domain.get_conserved_quantities(2, edge=1)[0] #Transmissive (4.5)
261        assert domain.quantities['stage'].boundary_values[4] ==\
262               domain.get_conserved_quantities(3, edge=1)[0] #Transmissive (-1.5)
263        assert domain.quantities['stage'].boundary_values[5] ==\
264               domain.get_conserved_quantities(3, edge=2)[0] #Transmissive (-1.5)
265
266        #Check enumeration
267        for k, ((vol_id, edge_id), _) in enumerate(domain.boundary_objects):
268            assert domain.neighbours[vol_id, edge_id] == -k-1
269
270
271
272
273    def test_distribute_first_order(self):
274        """Domain implements a default first order gradient limiter
275        """
276
277        a = [0.0, 0.0]
278        b = [0.0, 2.0]
279        c = [2.0,0.0]
280        d = [0.0, 4.0]
281        e = [2.0, 2.0]
282        f = [4.0,0.0]
283
284        points = [a, b, c, d, e, f]
285        #bac, bce, ecf, dbe
286        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
287        boundary = { (0, 0): 'Third',
288                     (0, 2): 'First',
289                     (2, 0): 'Second',
290                     (2, 1): 'Second',
291                     (3, 1): 'Second',
292                     (3, 2): 'Third'}
293
294
295        domain = Domain(points, vertices, boundary,
296                        conserved_quantities =\
297                        ['stage', 'xmomentum', 'ymomentum'])
298        domain.check_integrity()
299
300
301        domain.set_quantity('stage', [[1,2,3], [5,5,5],
302                                      [0,0,9], [-6, 3, 3]])
303
304        assert allclose( domain.quantities['stage'].centroid_values,
305                         [2,5,3,0] )
306
307        domain.set_quantity('xmomentum', [[1,1,1], [2,2,2],
308                                          [3,3,3], [4, 4, 4]])
309
310        domain.set_quantity('ymomentum', [[10,10,10], [20,20,20],
311                                          [30,30,30], [40, 40, 40]])
312
313
314        domain.distribute_to_vertices_and_edges()
315
316        #First order extrapolation
317        assert allclose( domain.quantities['stage'].vertex_values,
318                         [[ 2.,  2.,  2.],
319                          [ 5.,  5.,  5.],
320                          [ 3.,  3.,  3.],
321                          [ 0.,  0.,  0.]])
322
323
324
325
326    def test_update_conserved_quantities(self):
327        a = [0.0, 0.0]
328        b = [0.0, 2.0]
329        c = [2.0,0.0]
330        d = [0.0, 4.0]
331        e = [2.0, 2.0]
332        f = [4.0,0.0]
333
334        points = [a, b, c, d, e, f]
335        #bac, bce, ecf, dbe
336        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
337        boundary = { (0, 0): 'Third',
338                     (0, 2): 'First',
339                     (2, 0): 'Second',
340                     (2, 1): 'Second',
341                     (3, 1): 'Second',
342                     (3, 2): 'Third'}
343
344
345        domain = Domain(points, vertices, boundary,
346                        conserved_quantities =\
347                        ['stage', 'xmomentum', 'ymomentum'])
348        domain.check_integrity()
349
350
351        domain.set_quantity('stage', [1,2,3,4], location='centroids')
352        domain.set_quantity('xmomentum', [1,2,3,4], location='centroids')
353        domain.set_quantity('ymomentum', [1,2,3,4], location='centroids')
354
355
356        #Assign some values to update vectors
357        #Set explicit_update
358
359        for name in domain.conserved_quantities:
360            domain.quantities[name].explicit_update = array([4.,3.,2.,1.])
361            domain.quantities[name].semi_implicit_update = array([1.,1.,1.,1.])
362
363
364        #Update with given timestep (assuming no other forcing terms)
365        domain.timestep = 0.1
366        domain.update_conserved_quantities()
367
368        sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4])
369        denom = ones(4, Float)-domain.timestep*sem
370
371        x = array([1, 2, 3, 4]) + array( [.4,.3,.2,.1] )
372        x /= denom
373
374        for name in domain.conserved_quantities:
375            assert allclose(domain.quantities[name].centroid_values, x)
376
377
378    def test_set_region(self):
379        """Set quantities for sub region
380        """
381
382        a = [0.0, 0.0]
383        b = [0.0, 2.0]
384        c = [2.0,0.0]
385        d = [0.0, 4.0]
386        e = [2.0, 2.0]
387        f = [4.0,0.0]
388
389        points = [a, b, c, d, e, f]
390        #bac, bce, ecf, dbe
391        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
392        boundary = { (0, 0): 'Third',
393                     (0, 2): 'First',
394                     (2, 0): 'Second',
395                     (2, 1): 'Second',
396                     (3, 1): 'Second',
397                     (3, 2): 'Third'}
398
399        domain = Domain(points, vertices, boundary,
400                        conserved_quantities =\
401                        ['stage', 'xmomentum', 'ymomentum'])
402        domain.check_integrity()
403
404        domain.set_quantity('stage', [[1,2,3], [5,5,5],
405                                      [0,0,9], [-6, 3, 3]])
406
407        assert allclose( domain.quantities['stage'].centroid_values,
408                         [2,5,3,0] )
409
410        domain.set_quantity('xmomentum', [[1,1,1], [2,2,2],
411                                          [3,3,3], [4, 4, 4]])
412
413        domain.set_quantity('ymomentum', [[10,10,10], [20,20,20],
414                                          [30,30,30], [40, 40, 40]])
415
416
417        domain.distribute_to_vertices_and_edges()
418
419        #First order extrapolation
420        assert allclose( domain.quantities['stage'].vertex_values,
421                         [[ 2.,  2.,  2.],
422                          [ 5.,  5.,  5.],
423                          [ 3.,  3.,  3.],
424                          [ 0.,  0.,  0.]])
425
426        domain.build_tagged_elements_dictionary({'mound':[0,1]})
427        domain.set_region([add_to_verts])
428
429        self.failUnless(domain.test == "Mound",
430                        'set region failed')
431
432
433
434    def test_region_tags(self):
435        """
436        get values based on triangle lists.
437        """
438        from mesh_factory import rectangular
439        from shallow_water import Domain
440        from Numeric import zeros, Float
441
442        #Create basic mesh
443        points, vertices, boundary = rectangular(1, 3)
444
445        #Create shallow water domain
446        domain = Domain(points, vertices, boundary)
447        domain.build_tagged_elements_dictionary({'bottom':[0,1],
448                                                 'top':[4,5],
449                                                 'all':[0,1,2,3,4,5]})
450
451
452        #Set friction
453        manning = 0.07
454        domain.set_quantity('friction', manning)
455
456        domain.set_region([set_bottom_friction, set_top_friction])
457        #print domain.quantities['friction'].get_values()
458        assert allclose(domain.quantities['friction'].get_values(),\
459                        [[ 0.09,  0.09,  0.09],
460                         [ 0.09,  0.09,  0.09],
461                         [ 0.07,  0.07,  0.07],
462                         [ 0.07,  0.07,  0.07],
463                         [ 1.0,  1.0,  1.0],
464                         [ 1.0,  1.0,  1.0]])
465
466        domain.set_region([set_all_friction])
467        #print domain.quantities['friction'].get_values()
468        assert allclose(domain.quantities['friction'].get_values(),
469                        [[ 10.09, 10.09, 10.09],
470                         [ 10.09, 10.09, 10.09],
471                         [ 10.07, 10.07, 10.07],
472                         [ 10.07, 10.07, 10.07],
473                         [ 11.0,  11.0,  11.0],
474                         [ 11.0,  11.0,  11.0]])
475
476
477#-------------------------------------------------------------
478if __name__ == "__main__":
479    suite = unittest.makeSuite(Test_Domain,'test')
480    runner = unittest.TextTestRunner()
481    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.