source: anuga_core/source/anuga/abstract_2d_finite_volumes/test_domain.py @ 3698

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

Hi all,
I'm doing a change in the anuga structure, moving the code to

\anuga_core\source\anuga

After you have done an svn update, the PYTHONPATH has to be changed to;
PYTHONPATH = anuga_core/source/

This is part of changes required to make installation of anuga quicker and reducing the size of our sandpits.

If any imports are broken, try fixing them. With adding anuga. to them for example. If this seems to have really broken things, email/phone me.

Cheers
Duncan

File size: 18.6 KB
Line 
1#!/usr/bin/env python
2
3import unittest
4from math import sqrt
5
6from domain import *
7from anuga.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 anuga.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.