source: inundation/pyvolution/test_domain.py @ 3447

Last change on this file since 3447 was 3447, checked in by duncan, 18 years ago

Making set_region easier to use.

File size: 18.6 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').get_values(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 = domain.create_quantity_from_expression(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 = domain.create_quantity_from_expression(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    def test_set_quantity_from_expression(self):
197        """Quantity set using arbitrary expression
198
199        """
200
201
202        a = [0.0, 0.0]
203        b = [0.0, 2.0]
204        c = [2.0,0.0]
205        d = [0.0, 4.0]
206        e = [2.0, 2.0]
207        f = [4.0,0.0]
208
209        points = [a, b, c, d, e, f]
210        #bac, bce, ecf, dbe, daf, dae
211        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
212
213        domain = Domain(points, vertices, boundary=None,
214                        conserved_quantities =\
215                        ['stage', 'xmomentum', 'ymomentum'],
216                        other_quantities = ['elevation', 'friction', 'depth'])
217
218
219        domain.set_quantity('elevation', -1)
220
221
222        domain.set_quantity('stage', [[1,2,3], [5,5,5],
223                                      [0,0,9], [-6, 3, 3]])
224
225        domain.set_quantity('xmomentum', [[1,2,3], [5,5,5],
226                                          [0,0,9], [-6, 3, 3]])
227
228        domain.set_quantity('ymomentum', [[3,3,3], [4,2,1],
229                                          [2,4,-1], [1, 0, 1]])
230
231
232
233
234        domain.set_quantity('depth', expression = 'stage - elevation')
235
236        domain.check_integrity()
237
238
239
240
241        Q = domain.quantities['depth']
242
243        assert allclose(Q.vertex_values, [[2,3,4], [6,6,6],
244                                      [1,1,10], [-5, 4, 4]])
245
246
247
248
249
250
251
252
253
254    def test_boundary_indices(self):
255
256        from config import default_boundary_tag
257
258
259        a = [0.0, 0.5]
260        b = [0.0, 0.0]
261        c = [0.5, 0.5]
262
263        points = [a, b, c]
264        vertices = [ [0,1,2] ]
265        domain = Domain(points, vertices)
266
267        domain.set_boundary( {default_boundary_tag: Dirichlet_boundary([5,2,1])} )
268
269
270        domain.check_integrity()
271
272        assert allclose(domain.neighbours, [[-1,-2,-3]])
273
274
275
276    def test_boundary_conditions(self):
277
278        a = [0.0, 0.0]
279        b = [0.0, 2.0]
280        c = [2.0,0.0]
281        d = [0.0, 4.0]
282        e = [2.0, 2.0]
283        f = [4.0,0.0]
284
285        points = [a, b, c, d, e, f]
286        #bac, bce, ecf, dbe
287        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
288        boundary = { (0, 0): 'First',
289                     (0, 2): 'First',
290                     (2, 0): 'Second',
291                     (2, 1): 'Second',
292                     (3, 1): 'Second',
293                     (3, 2): 'Second'}
294
295
296        domain = Domain(points, vertices, boundary,
297                        conserved_quantities =\
298                        ['stage', 'xmomentum', 'ymomentum'])
299        domain.check_integrity()
300
301
302
303        domain.set_quantity('stage', [[1,2,3], [5,5,5],
304                                      [0,0,9], [-6, 3, 3]])
305
306
307        domain.set_boundary( {'First': Dirichlet_boundary([5,2,1]),
308                              'Second': Transmissive_boundary(domain)} )
309
310        domain.update_boundary()
311
312        assert domain.quantities['stage'].boundary_values[0] == 5. #Dirichlet
313        assert domain.quantities['stage'].boundary_values[1] == 5. #Dirichlet
314        assert domain.quantities['stage'].boundary_values[2] ==\
315               domain.get_conserved_quantities(2, edge=0)[0] #Transmissive (4.5)
316        assert domain.quantities['stage'].boundary_values[3] ==\
317               domain.get_conserved_quantities(2, edge=1)[0] #Transmissive (4.5)
318        assert domain.quantities['stage'].boundary_values[4] ==\
319               domain.get_conserved_quantities(3, edge=1)[0] #Transmissive (-1.5)
320        assert domain.quantities['stage'].boundary_values[5] ==\
321               domain.get_conserved_quantities(3, edge=2)[0] #Transmissive (-1.5)
322
323        #Check enumeration
324        for k, ((vol_id, edge_id), _) in enumerate(domain.boundary_objects):
325            assert domain.neighbours[vol_id, edge_id] == -k-1
326
327
328
329
330    def test_distribute_first_order(self):
331        """Domain implements a default first order gradient limiter
332        """
333
334        a = [0.0, 0.0]
335        b = [0.0, 2.0]
336        c = [2.0,0.0]
337        d = [0.0, 4.0]
338        e = [2.0, 2.0]
339        f = [4.0,0.0]
340
341        points = [a, b, c, d, e, f]
342        #bac, bce, ecf, dbe
343        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
344        boundary = { (0, 0): 'Third',
345                     (0, 2): 'First',
346                     (2, 0): 'Second',
347                     (2, 1): 'Second',
348                     (3, 1): 'Second',
349                     (3, 2): 'Third'}
350
351
352        domain = Domain(points, vertices, boundary,
353                        conserved_quantities =\
354                        ['stage', 'xmomentum', 'ymomentum'])
355        domain.check_integrity()
356
357
358        domain.set_quantity('stage', [[1,2,3], [5,5,5],
359                                      [0,0,9], [-6, 3, 3]])
360
361        assert allclose( domain.quantities['stage'].centroid_values,
362                         [2,5,3,0] )
363
364        domain.set_quantity('xmomentum', [[1,1,1], [2,2,2],
365                                          [3,3,3], [4, 4, 4]])
366
367        domain.set_quantity('ymomentum', [[10,10,10], [20,20,20],
368                                          [30,30,30], [40, 40, 40]])
369
370
371        domain.distribute_to_vertices_and_edges()
372
373        #First order extrapolation
374        assert allclose( domain.quantities['stage'].vertex_values,
375                         [[ 2.,  2.,  2.],
376                          [ 5.,  5.,  5.],
377                          [ 3.,  3.,  3.],
378                          [ 0.,  0.,  0.]])
379
380
381
382
383    def test_update_conserved_quantities(self):
384        a = [0.0, 0.0]
385        b = [0.0, 2.0]
386        c = [2.0,0.0]
387        d = [0.0, 4.0]
388        e = [2.0, 2.0]
389        f = [4.0,0.0]
390
391        points = [a, b, c, d, e, f]
392        #bac, bce, ecf, dbe
393        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
394        boundary = { (0, 0): 'Third',
395                     (0, 2): 'First',
396                     (2, 0): 'Second',
397                     (2, 1): 'Second',
398                     (3, 1): 'Second',
399                     (3, 2): 'Third'}
400
401
402        domain = Domain(points, vertices, boundary,
403                        conserved_quantities =\
404                        ['stage', 'xmomentum', 'ymomentum'])
405        domain.check_integrity()
406
407
408        domain.set_quantity('stage', [1,2,3,4], location='centroids')
409        domain.set_quantity('xmomentum', [1,2,3,4], location='centroids')
410        domain.set_quantity('ymomentum', [1,2,3,4], location='centroids')
411
412
413        #Assign some values to update vectors
414        #Set explicit_update
415
416        for name in domain.conserved_quantities:
417            domain.quantities[name].explicit_update = array([4.,3.,2.,1.])
418            domain.quantities[name].semi_implicit_update = array([1.,1.,1.,1.])
419
420
421        #Update with given timestep (assuming no other forcing terms)
422        domain.timestep = 0.1
423        domain.update_conserved_quantities()
424
425        sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4])
426        denom = ones(4, Float)-domain.timestep*sem
427
428#        x = array([1, 2, 3, 4]) + array( [.4,.3,.2,.1] )
429#        x /= denom
430
431        x = array([1., 2., 3., 4.])
432        x /= denom
433        x += domain.timestep*array( [4,3,2,1] )
434
435        for name in domain.conserved_quantities:
436            assert allclose(domain.quantities[name].centroid_values, x)
437
438
439    def test_set_region(self):
440        """Set quantities for sub region
441        """
442
443        a = [0.0, 0.0]
444        b = [0.0, 2.0]
445        c = [2.0,0.0]
446        d = [0.0, 4.0]
447        e = [2.0, 2.0]
448        f = [4.0,0.0]
449
450        points = [a, b, c, d, e, f]
451        #bac, bce, ecf, dbe
452        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
453        boundary = { (0, 0): 'Third',
454                     (0, 2): 'First',
455                     (2, 0): 'Second',
456                     (2, 1): 'Second',
457                     (3, 1): 'Second',
458                     (3, 2): 'Third'}
459
460        domain = Domain(points, vertices, boundary,
461                        conserved_quantities =\
462                        ['stage', 'xmomentum', 'ymomentum'])
463        domain.check_integrity()
464
465        domain.set_quantity('stage', [[1,2,3], [5,5,5],
466                                      [0,0,9], [-6, 3, 3]])
467
468        assert allclose( domain.quantities['stage'].centroid_values,
469                         [2,5,3,0] )
470
471        domain.set_quantity('xmomentum', [[1,1,1], [2,2,2],
472                                          [3,3,3], [4, 4, 4]])
473
474        domain.set_quantity('ymomentum', [[10,10,10], [20,20,20],
475                                          [30,30,30], [40, 40, 40]])
476
477
478        domain.distribute_to_vertices_and_edges()
479
480        #First order extrapolation
481        assert allclose( domain.quantities['stage'].vertex_values,
482                         [[ 2.,  2.,  2.],
483                          [ 5.,  5.,  5.],
484                          [ 3.,  3.,  3.],
485                          [ 0.,  0.,  0.]])
486
487        domain.build_tagged_elements_dictionary({'mound':[0,1]})
488        domain.set_region([add_to_verts])
489
490        self.failUnless(domain.test == "Mound",
491                        'set region failed')
492
493
494
495    def test_region_tags(self):
496        """
497        get values based on triangle lists.
498        """
499        from mesh_factory import rectangular
500        from shallow_water import Domain
501        from Numeric import zeros, Float
502
503        #Create basic mesh
504        points, vertices, boundary = rectangular(1, 3)
505
506        #Create shallow water domain
507        domain = Domain(points, vertices, boundary)
508        domain.build_tagged_elements_dictionary({'bottom':[0,1],
509                                                 'top':[4,5],
510                                                 'all':[0,1,2,3,4,5]})
511
512
513        #Set friction
514        manning = 0.07
515        domain.set_quantity('friction', manning)
516
517        domain.set_region([set_bottom_friction, set_top_friction])
518        #print domain.quantities['friction'].get_values()
519        assert allclose(domain.quantities['friction'].get_values(),\
520                        [[ 0.09,  0.09,  0.09],
521                         [ 0.09,  0.09,  0.09],
522                         [ 0.07,  0.07,  0.07],
523                         [ 0.07,  0.07,  0.07],
524                         [ 1.0,  1.0,  1.0],
525                         [ 1.0,  1.0,  1.0]])
526
527        domain.set_region([set_all_friction])
528        #print domain.quantities['friction'].get_values()
529        assert allclose(domain.quantities['friction'].get_values(),
530                        [[ 10.09, 10.09, 10.09],
531                         [ 10.09, 10.09, 10.09],
532                         [ 10.07, 10.07, 10.07],
533                         [ 10.07, 10.07, 10.07],
534                         [ 11.0,  11.0,  11.0],
535                         [ 11.0,  11.0,  11.0]])
536
537
538    def test_region_tags2(self):
539        """
540        get values based on triangle lists.
541        """
542        from mesh_factory import rectangular
543        from shallow_water import Domain
544        from Numeric import zeros, Float
545
546        #Create basic mesh
547        points, vertices, boundary = rectangular(1, 3)
548
549        #Create shallow water domain
550        domain = Domain(points, vertices, boundary)
551        domain.build_tagged_elements_dictionary({'bottom':[0,1],
552                                                 'top':[4,5],
553                                                 'all':[0,1,2,3,4,5]})
554
555
556        #Set friction
557        manning = 0.07
558        domain.set_quantity('friction', manning)
559
560        domain.set_region('top', 'friction', 1.0)
561        domain.set_region('bottom', 'friction', 0.09)
562       
563        #print domain.quantities['friction'].get_values()
564        assert allclose(domain.quantities['friction'].get_values(),\
565                        [[ 0.09,  0.09,  0.09],
566                         [ 0.09,  0.09,  0.09],
567                         [ 0.07,  0.07,  0.07],
568                         [ 0.07,  0.07,  0.07],
569                         [ 1.0,  1.0,  1.0],
570                         [ 1.0,  1.0,  1.0]])
571       
572        domain.set_region([set_bottom_friction, set_top_friction])
573        #print domain.quantities['friction'].get_values()
574        assert allclose(domain.quantities['friction'].get_values(),\
575                        [[ 0.09,  0.09,  0.09],
576                         [ 0.09,  0.09,  0.09],
577                         [ 0.07,  0.07,  0.07],
578                         [ 0.07,  0.07,  0.07],
579                         [ 1.0,  1.0,  1.0],
580                         [ 1.0,  1.0,  1.0]])
581
582        domain.set_region([set_all_friction])
583        #print domain.quantities['friction'].get_values()
584        assert allclose(domain.quantities['friction'].get_values(),
585                        [[ 10.09, 10.09, 10.09],
586                         [ 10.09, 10.09, 10.09],
587                         [ 10.07, 10.07, 10.07],
588                         [ 10.07, 10.07, 10.07],
589                         [ 11.0,  11.0,  11.0],
590                         [ 11.0,  11.0,  11.0]])
591
592#-------------------------------------------------------------
593if __name__ == "__main__":
594    suite = unittest.makeSuite(Test_Domain,'test')
595    runner = unittest.TextTestRunner()
596    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.