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

Last change on this file since 7573 was 7573, checked in by steve, 14 years ago

Committing a version of shallow_water_balanced which passes it unit test
using a version of edge limiting which doesn't limit boundary edges. THis
is useful to allow linear functions to be reconstructed.

Had to play with the transmissive BC and use centroid values which is
set via domain set_centroid_transmissive_bc

File size: 34.8 KB
Line 
1#!/usr/bin/env python
2
3import unittest
4from math import sqrt
5
6from anuga.abstract_2d_finite_volumes.domain import *
7from anuga.config import epsilon
8
9import numpy as num
10
11
12def add_to_verts(tag, elements, domain):
13    if tag == "mound":
14        domain.test = "Mound"
15
16
17def set_bottom_friction(tag, elements, domain):
18    if tag == "bottom":
19        domain.set_quantity('friction', 0.09, indices = elements)
20
21def set_top_friction(tag, elements, domain):
22    if tag == "top":
23        domain.set_quantity('friction', 1., indices = elements)
24
25
26def set_all_friction(tag, elements, domain):
27    if tag == 'all':
28        new_values = domain.get_quantity('friction').get_values(indices = elements) + 10.0
29
30        domain.set_quantity('friction', new_values, indices = elements)
31
32
33class Test_Domain(unittest.TestCase):
34    def setUp(self):
35        pass
36
37
38    def tearDown(self):
39        pass
40
41
42    def test_simple(self):
43        a = [0.0, 0.0]
44        b = [0.0, 2.0]
45        c = [2.0,0.0]
46        d = [0.0, 4.0]
47        e = [2.0, 2.0]
48        f = [4.0,0.0]
49
50        points = [a, b, c, d, e, f]
51        #bac, bce, ecf, dbe, daf, dae
52        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
53
54        conserved_quantities = ['stage', 'xmomentum', 'ymomentum']
55        evolved_quantities = ['stage', 'xmomentum', 'ymomentum', 'xvelocity']
56       
57        other_quantities = ['elevation', 'friction']
58
59        domain = Domain(points, vertices, None,
60                        conserved_quantities, evolved_quantities, other_quantities)
61        domain.check_integrity()
62
63        for name in conserved_quantities + other_quantities:
64            assert domain.quantities.has_key(name)
65
66
67        assert num.alltrue(domain.get_conserved_quantities(0, edge=1) == 0.)
68
69
70
71    def test_CFL(self):
72        a = [0.0, 0.0]
73        b = [0.0, 2.0]
74        c = [2.0,0.0]
75        d = [0.0, 4.0]
76        e = [2.0, 2.0]
77        f = [4.0,0.0]
78
79        points = [a, b, c, d, e, f]
80        #bac, bce, ecf, dbe, daf, dae
81        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
82
83        conserved_quantities = ['stage', 'xmomentum', 'ymomentum']
84        evolved_quantities = ['stage', 'xmomentum', 'ymomentum', 'xvelocity']
85       
86        other_quantities = ['elevation', 'friction']
87
88        domain = Domain(points, vertices, None,
89                        conserved_quantities, evolved_quantities, other_quantities)
90
91        try:
92            domain.set_CFL(-0.1)
93        except:
94            pass
95        else:
96            msg = 'Should have caught a negative cfl'
97            raise Exception, msg
98
99
100       
101        try:
102            domain.set_CFL(2.0)
103        except:
104            pass
105        else:
106            msg = 'Should have warned of cfl>1.0'
107            raise Exception, msg
108
109        assert domain.CFL == 2.0
110       
111
112    def test_conserved_quantities(self):
113
114        a = [0.0, 0.0]
115        b = [0.0, 2.0]
116        c = [2.0,0.0]
117        d = [0.0, 4.0]
118        e = [2.0, 2.0]
119        f = [4.0,0.0]
120
121        points = [a, b, c, d, e, f]
122        #bac, bce, ecf, dbe, daf, dae
123        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
124
125        domain = Domain(points, vertices, boundary=None,
126                        conserved_quantities =\
127                        ['stage', 'xmomentum', 'ymomentum'])
128
129
130        domain.set_quantity('stage', [[1,2,3], [5,5,5],
131                                      [0,0,9], [-6, 3, 3]])
132
133        domain.set_quantity('xmomentum', [[1,2,3], [5,5,5],
134                                          [0,0,9], [-6, 3, 3]])
135
136        domain.check_integrity()
137
138        #Centroids
139        q = domain.get_conserved_quantities(0)
140        assert num.allclose(q, [2., 2., 0.])
141
142        q = domain.get_conserved_quantities(1)
143        assert num.allclose(q, [5., 5., 0.])
144
145        q = domain.get_conserved_quantities(2)
146        assert num.allclose(q, [3., 3., 0.])
147
148        q = domain.get_conserved_quantities(3)
149        assert num.allclose(q, [0., 0., 0.])
150
151
152        #Edges
153        q = domain.get_conserved_quantities(0, edge=0)
154        assert num.allclose(q, [2.5, 2.5, 0.])
155        q = domain.get_conserved_quantities(0, edge=1)
156        assert num.allclose(q, [2., 2., 0.])
157        q = domain.get_conserved_quantities(0, edge=2)
158        assert num.allclose(q, [1.5, 1.5, 0.])
159
160        for i in range(3):
161            q = domain.get_conserved_quantities(1, edge=i)
162            assert num.allclose(q, [5, 5, 0.])
163
164
165        q = domain.get_conserved_quantities(2, edge=0)
166        assert num.allclose(q, [4.5, 4.5, 0.])
167        q = domain.get_conserved_quantities(2, edge=1)
168        assert num.allclose(q, [4.5, 4.5, 0.])
169        q = domain.get_conserved_quantities(2, edge=2)
170        assert num.allclose(q, [0., 0., 0.])
171
172
173        q = domain.get_conserved_quantities(3, edge=0)
174        assert num.allclose(q, [3., 3., 0.])
175        q = domain.get_conserved_quantities(3, edge=1)
176        assert num.allclose(q, [-1.5, -1.5, 0.])
177        q = domain.get_conserved_quantities(3, edge=2)
178        assert num.allclose(q, [-1.5, -1.5, 0.])
179
180
181
182    def test_create_quantity_from_expression(self):
183        """Quantity created from other quantities using arbitrary expression
184
185        """
186
187
188        a = [0.0, 0.0]
189        b = [0.0, 2.0]
190        c = [2.0,0.0]
191        d = [0.0, 4.0]
192        e = [2.0, 2.0]
193        f = [4.0,0.0]
194
195        points = [a, b, c, d, e, f]
196        #bac, bce, ecf, dbe, daf, dae
197        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
198
199        domain = Domain(points, vertices, boundary=None,
200                        conserved_quantities =\
201                        ['stage', 'xmomentum', 'ymomentum'],
202                        other_quantities = ['elevation', 'friction'])
203
204
205        domain.set_quantity('elevation', -1)
206
207
208        domain.set_quantity('stage', [[1,2,3], [5,5,5],
209                                      [0,0,9], [-6, 3, 3]])
210
211        domain.set_quantity('xmomentum', [[1,2,3], [5,5,5],
212                                          [0,0,9], [-6, 3, 3]])
213
214        domain.set_quantity('ymomentum', [[3,3,3], [4,2,1],
215                                          [2,4,-1], [1, 0, 1]])
216
217        domain.check_integrity()
218
219
220
221        expression = 'stage - elevation'
222        Q = domain.create_quantity_from_expression(expression)
223
224        assert num.allclose(Q.vertex_values, [[2,3,4], [6,6,6],
225                                              [1,1,10], [-5, 4, 4]])
226
227        expression = '(xmomentum*xmomentum + ymomentum*ymomentum)**0.5'
228        Q = domain.create_quantity_from_expression(expression)
229
230        X = domain.quantities['xmomentum'].vertex_values
231        Y = domain.quantities['ymomentum'].vertex_values
232
233        assert num.allclose(Q.vertex_values, (X**2 + Y**2)**0.5)
234
235
236
237    def test_set_quanitities_to_be_monitored(self):
238        """test_set_quanitities_to_be_monitored
239        """
240
241        a = [0.0, 0.0]
242        b = [0.0, 2.0]
243        c = [2.0,0.0]
244        d = [0.0, 4.0]
245        e = [2.0, 2.0]
246        f = [4.0,0.0]
247
248        points = [a, b, c, d, e, f]
249        #bac, bce, ecf, dbe, daf, dae
250        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
251
252
253        domain = Domain(points, vertices, boundary=None,
254                        conserved_quantities =\
255                        ['stage', 'xmomentum', 'ymomentum'],
256                        other_quantities = ['elevation', 'friction', 'depth'])
257
258
259        assert domain.quantities_to_be_monitored is None
260        domain.set_quantities_to_be_monitored(['stage', 'stage-elevation'])
261        assert len(domain.quantities_to_be_monitored) == 2
262        assert domain.quantities_to_be_monitored.has_key('stage')
263        assert domain.quantities_to_be_monitored.has_key('stage-elevation')
264        for key in domain.quantities_to_be_monitored['stage'].keys():
265            assert domain.quantities_to_be_monitored['stage'][key] is None
266
267
268        # Check that invalid requests are dealt with
269        try:
270            domain.set_quantities_to_be_monitored(['yyyyy'])       
271        except:
272            pass
273        else:
274            msg = 'Should have caught illegal quantity'
275            raise Exception, msg
276
277        try:
278            domain.set_quantities_to_be_monitored(['stage-xx'])       
279        except NameError:
280            pass
281        else:
282            msg = 'Should have caught illegal quantity'
283            raise Exception, msg
284
285        try:
286            domain.set_quantities_to_be_monitored('stage', 'stage-elevation')
287        except:
288            pass
289        else:
290            msg = 'Should have caught too many arguments'
291            raise Exception, msg
292
293        try:
294            domain.set_quantities_to_be_monitored('stage', 'blablabla')
295        except:
296            pass
297        else:
298            msg = 'Should have caught polygon as a string'
299            raise Exception, msg       
300
301
302
303        # Now try with a polygon restriction
304        domain.set_quantities_to_be_monitored('xmomentum',
305                                              polygon=[[1,1], [1,3], [3,3], [3,1]],
306                                              time_interval = [0,3])
307        assert domain.monitor_indices[0] == 1
308        assert domain.monitor_time_interval[0] == 0
309        assert domain.monitor_time_interval[1] == 3       
310       
311
312    def test_set_quantity_from_expression(self):
313        """Quantity set using arbitrary expression
314
315        """
316
317
318        a = [0.0, 0.0]
319        b = [0.0, 2.0]
320        c = [2.0,0.0]
321        d = [0.0, 4.0]
322        e = [2.0, 2.0]
323        f = [4.0,0.0]
324
325        points = [a, b, c, d, e, f]
326        #bac, bce, ecf, dbe, daf, dae
327        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
328
329        domain = Domain(points, vertices, boundary=None,
330                        conserved_quantities =\
331                        ['stage', 'xmomentum', 'ymomentum'],
332                        other_quantities = ['elevation', 'friction', 'depth'])
333
334
335        domain.set_quantity('elevation', -1)
336
337
338        domain.set_quantity('stage', [[1,2,3], [5,5,5],
339                                      [0,0,9], [-6, 3, 3]])
340
341        domain.set_quantity('xmomentum', [[1,2,3], [5,5,5],
342                                          [0,0,9], [-6, 3, 3]])
343
344        domain.set_quantity('ymomentum', [[3,3,3], [4,2,1],
345                                          [2,4,-1], [1, 0, 1]])
346
347
348
349
350        domain.set_quantity('depth', expression = 'stage - elevation')
351
352        domain.check_integrity()
353
354
355
356
357        Q = domain.quantities['depth']
358
359        assert num.allclose(Q.vertex_values, [[2,3,4], [6,6,6],
360                                              [1,1,10], [-5, 4, 4]])
361
362
363
364                                     
365    def test_add_quantity(self):
366        """Test that quantities already set can be added to using
367        add_quantity
368
369        """
370
371
372        a = [0.0, 0.0]
373        b = [0.0, 2.0]
374        c = [2.0,0.0]
375        d = [0.0, 4.0]
376        e = [2.0, 2.0]
377        f = [4.0,0.0]
378
379        points = [a, b, c, d, e, f]
380        #bac, bce, ecf, dbe, daf, dae
381        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
382
383        domain = Domain(points, vertices, boundary=None,
384                        conserved_quantities =\
385                        ['stage', 'xmomentum', 'ymomentum'],
386                        other_quantities = ['elevation', 'friction', 'depth'])
387
388
389        A = num.array([[1,2,3], [5,5,-5], [0,0,9], [-6,3,3]], num.float)
390        B = num.array([[2,4,4], [3,2,1], [6,-3,4], [4,5,-1]], num.float)
391       
392        # Shorthands
393        stage = domain.quantities['stage']
394        elevation = domain.quantities['elevation']
395        depth = domain.quantities['depth']
396       
397        # Go testing
398        domain.set_quantity('elevation', A)
399        domain.add_quantity('elevation', B)
400        assert num.allclose(elevation.vertex_values, A+B)
401       
402        domain.add_quantity('elevation', 4)
403        assert num.allclose(elevation.vertex_values, A+B+4)       
404       
405       
406        # Test using expression
407        domain.set_quantity('stage', [[1,2,3], [5,5,5],
408                                      [0,0,9], [-6, 3, 3]])       
409        domain.set_quantity('depth', 1.0)                                     
410        domain.add_quantity('depth', expression = 'stage - elevation')       
411        assert num.allclose(depth.vertex_values, stage.vertex_values-elevation.vertex_values+1)
412               
413       
414        # Check self referential expression
415        reference = 2*stage.vertex_values - depth.vertex_values
416        domain.add_quantity('stage', expression = 'stage - depth')               
417        assert num.allclose(stage.vertex_values, reference)       
418                                     
419
420        # Test using a function
421        def f(x, y):
422            return x+y
423           
424        domain.set_quantity('elevation', f)           
425        domain.set_quantity('stage', 5.0)
426        domain.set_quantity('depth', expression = 'stage - elevation')
427       
428        domain.add_quantity('depth', f)
429        assert num.allclose(stage.vertex_values, depth.vertex_values)               
430         
431           
432       
433       
434                                     
435                                     
436    def test_setting_timestepping_method(self):
437        """test_setting_timestepping_method
438        """
439
440        a = [0.0, 0.0]
441        b = [0.0, 2.0]
442        c = [2.0,0.0]
443        d = [0.0, 4.0]
444        e = [2.0, 2.0]
445        f = [4.0,0.0]
446
447        points = [a, b, c, d, e, f]
448        #bac, bce, ecf, dbe, daf, dae
449        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
450
451
452        domain = Domain(points, vertices, boundary=None,
453                        conserved_quantities =\
454                        ['stage', 'xmomentum', 'ymomentum'],
455                        other_quantities = ['elevation', 'friction', 'depth'])
456
457
458        domain.timestepping_method = None
459
460
461        # Check that invalid requests are dealt with
462        try:
463            domain.set_timestepping_method('eee')       
464        except:
465            pass
466        else:
467            msg = 'Should have caught illegal method'
468            raise Exception, msg
469
470
471        #Should have no trouble with euler, rk2 or rk3
472        domain.set_timestepping_method('euler')
473        domain.set_timestepping_method('rk2')
474        domain.set_timestepping_method('rk3')
475
476        domain.set_timestepping_method(1)
477        domain.set_timestepping_method(2)
478        domain.set_timestepping_method(3)
479
480        #test get timestepping method
481        assert domain.get_timestepping_method() == 'rk3'
482
483
484
485    def test_boundary_indices(self):
486
487        from anuga.config import default_boundary_tag
488
489
490        a = [0.0, 0.5]
491        b = [0.0, 0.0]
492        c = [0.5, 0.5]
493
494        points = [a, b, c]
495        vertices = [ [0,1,2] ]
496        domain = Domain(points, vertices)
497
498        domain.set_boundary( {default_boundary_tag: Dirichlet_boundary([5,2,1])} )
499
500
501        domain.check_integrity()
502
503        assert num.allclose(domain.neighbours, [[-1,-2,-3]])
504
505
506
507    def test_boundary_conditions(self):
508
509        a = [0.0, 0.0]
510        b = [0.0, 2.0]
511        c = [2.0,0.0]
512        d = [0.0, 4.0]
513        e = [2.0, 2.0]
514        f = [4.0,0.0]
515
516        points = [a, b, c, d, e, f]
517        #bac, bce, ecf, dbe
518        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
519        boundary = { (0, 0): 'First',
520                     (0, 2): 'First',
521                     (2, 0): 'Second',
522                     (2, 1): 'Second',
523                     (3, 1): 'Second',
524                     (3, 2): 'Second'}
525
526
527        domain = Domain(points, vertices, boundary,
528                        conserved_quantities =\
529                        ['stage', 'xmomentum', 'ymomentum'])
530        domain.check_integrity()
531
532
533
534        domain.set_quantity('stage', [[1,2,3], [5,5,5],
535                                      [0,0,9], [-6, 3, 3]])
536
537
538        domain.set_boundary( {'First': Dirichlet_boundary([5,2,1]),
539                              'Second': Transmissive_boundary(domain)} )
540
541        domain.update_boundary()
542
543        assert domain.quantities['stage'].boundary_values[0] == 5. #Dirichlet
544        assert domain.quantities['stage'].boundary_values[1] == 5. #Dirichlet
545        assert domain.quantities['stage'].boundary_values[2] ==\
546               domain.get_conserved_quantities(2, edge=0)[0] #Transmissive (4.5)
547        assert domain.quantities['stage'].boundary_values[3] ==\
548               domain.get_conserved_quantities(2, edge=1)[0] #Transmissive (4.5)
549        assert domain.quantities['stage'].boundary_values[4] ==\
550               domain.get_conserved_quantities(3, edge=1)[0] #Transmissive (-1.5)
551        assert domain.quantities['stage'].boundary_values[5] ==\
552               domain.get_conserved_quantities(3, edge=2)[0] #Transmissive (-1.5)
553
554        #Check enumeration
555        for k, ((vol_id, edge_id), _) in enumerate(domain.boundary_objects):
556            assert domain.neighbours[vol_id, edge_id] == -k-1
557
558
559
560
561    def test_conserved_evolved_boundary_conditions(self):
562
563        a = [0.0, 0.0]
564        b = [0.0, 2.0]
565        c = [2.0,0.0]
566        d = [0.0, 4.0]
567        e = [2.0, 2.0]
568        f = [4.0,0.0]
569
570        points = [a, b, c, d, e, f]
571        #bac, bce, ecf, dbe
572        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
573        boundary = { (0, 0): 'First',
574                     (0, 2): 'First',
575                     (2, 0): 'Second',
576                     (2, 1): 'Second',
577                     (3, 1): 'Second',
578                     (3, 2): 'Second'}
579
580
581 
582        try:
583            domain = Domain(points, vertices, boundary,
584                            conserved_quantities = ['stage', 'xmomentum', 'ymomentum'],
585                            evolved_quantities =\
586                                   ['stage', 'xmomentum', 'xvelocity', 'ymomentum', 'yvelocity'])
587        except:
588            pass
589        else:
590            msg = 'Should have caught the evolved quantities not being in order'
591            raise Exception, msg           
592
593
594        domain = Domain(points, vertices, boundary,
595                        conserved_quantities = ['stage', 'xmomentum', 'ymomentum'],
596                        evolved_quantities =\
597                        ['stage', 'xmomentum', 'ymomentum', 'xvelocity', 'yvelocity'])
598
599
600        domain.set_quantity('stage', [[1,2,3], [5,5,5],
601                                      [0,0,9], [6, -3, 3]])
602
603
604        domain.set_boundary( {'First': Dirichlet_boundary([5,2,1,4,6]),
605                              'Second': Transmissive_boundary(domain)} )
606
607        try:
608            domain.update_boundary()
609        except:
610            pass
611        else:
612            msg = 'Should have caught the lack of conserved_values_to_evolved_values member function'
613            raise Exception, msg
614
615
616        def  conserved_values_to_evolved_values(q_cons, q_evol):
617
618            q_evol[0:3] = q_cons
619            q_evol[3] = q_cons[1]/q_cons[0]
620            q_evol[4] = q_cons[2]/q_cons[0]
621
622            return q_evol
623
624        domain.conserved_values_to_evolved_values = conserved_values_to_evolved_values
625
626        domain.update_boundary()
627
628
629        assert domain.quantities['stage'].boundary_values[0] == 5. #Dirichlet
630        assert domain.quantities['stage'].boundary_values[1] == 5. #Dirichlet
631        assert domain.quantities['xvelocity'].boundary_values[0] == 4. #Dirichlet
632        assert domain.quantities['yvelocity'].boundary_values[1] == 6. #Dirichlet
633
634        q_cons = domain.get_conserved_quantities(2, edge=0) #Transmissive
635        assert domain.quantities['stage'    ].boundary_values[2] == q_cons[0]
636        assert domain.quantities['xmomentum'].boundary_values[2] == q_cons[1]
637        assert domain.quantities['ymomentum'].boundary_values[2] == q_cons[2]
638        assert domain.quantities['xvelocity'].boundary_values[2] == q_cons[1]/q_cons[0]
639        assert domain.quantities['yvelocity'].boundary_values[2] == q_cons[2]/q_cons[0]
640
641        q_cons = domain.get_conserved_quantities(2, edge=1) #Transmissive
642        assert domain.quantities['stage'    ].boundary_values[3] == q_cons[0]
643        assert domain.quantities['xmomentum'].boundary_values[3] == q_cons[1]
644        assert domain.quantities['ymomentum'].boundary_values[3] == q_cons[2]
645        assert domain.quantities['xvelocity'].boundary_values[3] == q_cons[1]/q_cons[0]
646        assert domain.quantities['yvelocity'].boundary_values[3] == q_cons[2]/q_cons[0]       
647
648
649        q_cons = domain.get_conserved_quantities(3, edge=1) #Transmissive
650        assert domain.quantities['stage'    ].boundary_values[4] == q_cons[0]
651        assert domain.quantities['xmomentum'].boundary_values[4] == q_cons[1]
652        assert domain.quantities['ymomentum'].boundary_values[4] == q_cons[2]
653        assert domain.quantities['xvelocity'].boundary_values[4] == q_cons[1]/q_cons[0]
654        assert domain.quantities['yvelocity'].boundary_values[4] == q_cons[2]/q_cons[0]               
655
656
657        q_cons = domain.get_conserved_quantities(3, edge=2) #Transmissive
658        assert domain.quantities['stage'    ].boundary_values[5] == q_cons[0]
659        assert domain.quantities['xmomentum'].boundary_values[5] == q_cons[1]
660        assert domain.quantities['ymomentum'].boundary_values[5] == q_cons[2]
661        assert domain.quantities['xvelocity'].boundary_values[5] == q_cons[1]/q_cons[0]
662        assert domain.quantities['yvelocity'].boundary_values[5] == q_cons[2]/q_cons[0]
663 
664
665    def test_distribute_first_order(self):
666        """Domain implements a default first order gradient limiter
667        """
668
669        a = [0.0, 0.0]
670        b = [0.0, 2.0]
671        c = [2.0,0.0]
672        d = [0.0, 4.0]
673        e = [2.0, 2.0]
674        f = [4.0,0.0]
675
676        points = [a, b, c, d, e, f]
677        #bac, bce, ecf, dbe
678        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
679        boundary = { (0, 0): 'Third',
680                     (0, 2): 'First',
681                     (2, 0): 'Second',
682                     (2, 1): 'Second',
683                     (3, 1): 'Second',
684                     (3, 2): 'Third'}
685
686
687        domain = Domain(points, vertices, boundary,
688                        conserved_quantities =\
689                        ['stage', 'xmomentum', 'ymomentum'])
690        domain.set_default_order(1)
691        domain.check_integrity()
692
693
694        domain.set_quantity('stage', [[1,2,3], [5,5,5],
695                                      [0,0,9], [-6, 3, 3]])
696
697        assert num.allclose( domain.quantities['stage'].centroid_values,
698                             [2,5,3,0] )
699
700        domain.set_quantity('xmomentum', [[1,1,1], [2,2,2],
701                                          [3,3,3], [4, 4, 4]])
702
703        domain.set_quantity('ymomentum', [[10,10,10], [20,20,20],
704                                          [30,30,30], [40, 40, 40]])
705
706
707        domain.distribute_to_vertices_and_edges()
708
709        #First order extrapolation
710        assert num.allclose( domain.quantities['stage'].vertex_values,
711                             [[ 2.,  2.,  2.],
712                              [ 5.,  5.,  5.],
713                              [ 3.,  3.,  3.],
714                              [ 0.,  0.,  0.]])
715
716
717
718
719    def test_update_conserved_quantities(self):
720        a = [0.0, 0.0]
721        b = [0.0, 2.0]
722        c = [2.0,0.0]
723        d = [0.0, 4.0]
724        e = [2.0, 2.0]
725        f = [4.0,0.0]
726
727        points = [a, b, c, d, e, f]
728        #bac, bce, ecf, dbe
729        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
730        boundary = { (0, 0): 'Third',
731                     (0, 2): 'First',
732                     (2, 0): 'Second',
733                     (2, 1): 'Second',
734                     (3, 1): 'Second',
735                     (3, 2): 'Third'}
736
737
738        domain = Domain(points, vertices, boundary,
739                        conserved_quantities =\
740                        ['stage', 'xmomentum', 'ymomentum'])
741        domain.check_integrity()
742
743
744        domain.set_quantity('stage', [1,2,3,4], location='centroids')
745        domain.set_quantity('xmomentum', [1,2,3,4], location='centroids')
746        domain.set_quantity('ymomentum', [1,2,3,4], location='centroids')
747
748
749        #Assign some values to update vectors
750        #Set explicit_update
751
752
753        for name in domain.conserved_quantities:
754            domain.quantities[name].explicit_update = num.array([4.,3.,2.,1.])
755            domain.quantities[name].semi_implicit_update = num.array([1.,1.,1.,1.])
756
757
758        #Update with given timestep (assuming no other forcing terms)
759        domain.timestep = 0.1
760        domain.update_conserved_quantities()
761
762        sem = num.array([1.,1.,1.,1.])/num.array([1, 2, 3, 4])
763        denom = num.ones(4, num.float) - domain.timestep*sem
764
765#        x = array([1, 2, 3, 4]) + array( [.4,.3,.2,.1] )
766#        x /= denom
767
768        x = num.array([1., 2., 3., 4.])
769        x += domain.timestep*num.array( [4,3,2,1] )
770        x /= denom
771
772
773        for name in domain.conserved_quantities:
774            assert num.allclose(domain.quantities[name].centroid_values, x)
775
776
777    def test_set_region(self):
778        """Set quantities for sub region
779        """
780
781        a = [0.0, 0.0]
782        b = [0.0, 2.0]
783        c = [2.0,0.0]
784        d = [0.0, 4.0]
785        e = [2.0, 2.0]
786        f = [4.0,0.0]
787
788        points = [a, b, c, d, e, f]
789        #bac, bce, ecf, dbe
790        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
791        boundary = { (0, 0): 'Third',
792                     (0, 2): 'First',
793                     (2, 0): 'Second',
794                     (2, 1): 'Second',
795                     (3, 1): 'Second',
796                     (3, 2): 'Third'}
797
798        domain = Domain(points, vertices, boundary,
799                        conserved_quantities =\
800                        ['stage', 'xmomentum', 'ymomentum'])
801        domain.set_default_order(1)                       
802        domain.check_integrity()
803
804        domain.set_quantity('stage', [[1,2,3], [5,5,5],
805                                      [0,0,9], [-6, 3, 3]])
806
807        assert num.allclose( domain.quantities['stage'].centroid_values,
808                             [2,5,3,0] )
809
810        domain.set_quantity('xmomentum', [[1,1,1], [2,2,2],
811                                          [3,3,3], [4, 4, 4]])
812
813        domain.set_quantity('ymomentum', [[10,10,10], [20,20,20],
814                                          [30,30,30], [40, 40, 40]])
815
816
817        domain.distribute_to_vertices_and_edges()
818
819        #First order extrapolation
820        assert num.allclose( domain.quantities['stage'].vertex_values,
821                             [[ 2.,  2.,  2.],
822                              [ 5.,  5.,  5.],
823                              [ 3.,  3.,  3.],
824                              [ 0.,  0.,  0.]])
825
826        domain.build_tagged_elements_dictionary({'mound':[0,1]})
827        domain.set_region([add_to_verts])
828
829        self.failUnless(domain.test == "Mound",
830                        'set region failed')
831
832
833    def test_track_speeds(self):
834        """
835        get values based on triangle lists.
836        """
837        from mesh_factory import rectangular
838        from shallow_water import Domain
839
840        #Create basic mesh
841        points, vertices, boundary = rectangular(1, 3)
842
843        #Create shallow water domain
844        domain = Domain(points, vertices, boundary)
845        domain.timestepping_statistics(track_speeds=True)
846
847
848
849    def test_region_tags(self):
850        """
851        get values based on triangle lists.
852        """
853        from mesh_factory import rectangular
854        from shallow_water import Domain
855
856        #Create basic mesh
857        points, vertices, boundary = rectangular(1, 3)
858
859        #Create shallow water domain
860        domain = Domain(points, vertices, boundary)
861        domain.build_tagged_elements_dictionary({'bottom':[0,1],
862                                                 'top':[4,5],
863                                                 'all':[0,1,2,3,4,5]})
864
865
866        #Set friction
867        manning = 0.07
868        domain.set_quantity('friction', manning)
869
870        domain.set_region([set_bottom_friction, set_top_friction])
871        assert num.allclose(domain.quantities['friction'].get_values(),\
872                            [[ 0.09,  0.09,  0.09],
873                             [ 0.09,  0.09,  0.09],
874                             [ 0.07,  0.07,  0.07],
875                             [ 0.07,  0.07,  0.07],
876                             [ 1.0,  1.0,  1.0],
877                             [ 1.0,  1.0,  1.0]])
878
879        domain.set_region([set_all_friction])
880        assert num.allclose(domain.quantities['friction'].get_values(),
881                            [[ 10.09, 10.09, 10.09],
882                             [ 10.09, 10.09, 10.09],
883                             [ 10.07, 10.07, 10.07],
884                             [ 10.07, 10.07, 10.07],
885                             [ 11.0,  11.0,  11.0],
886                             [ 11.0,  11.0,  11.0]])
887
888
889    def test_region_tags2(self):
890        """
891        get values based on triangle lists.
892        """
893        from mesh_factory import rectangular
894        from shallow_water import Domain
895
896        #Create basic mesh
897        points, vertices, boundary = rectangular(1, 3)
898
899        #Create shallow water domain
900        domain = Domain(points, vertices, boundary)
901        domain.build_tagged_elements_dictionary({'bottom':[0,1],
902                                                 'top':[4,5],
903                                                 'all':[0,1,2,3,4,5]})
904
905
906        #Set friction
907        manning = 0.07
908        domain.set_quantity('friction', manning)
909
910        domain.set_region('top', 'friction', 1.0)
911        domain.set_region('bottom', 'friction', 0.09)
912       
913        msg = ("domain.quantities['friction'].get_values()=\n%s\n"
914               'should equal\n'
915               '[[ 0.09,  0.09,  0.09],\n'
916               ' [ 0.09,  0.09,  0.09],\n'
917               ' [ 0.07,  0.07,  0.07],\n'
918               ' [ 0.07,  0.07,  0.07],\n'
919               ' [ 1.0,  1.0,  1.0],\n'
920               ' [ 1.0,  1.0,  1.0]]'
921               % str(domain.quantities['friction'].get_values()))
922        assert num.allclose(domain.quantities['friction'].get_values(),
923                            [[ 0.09,  0.09,  0.09],
924                             [ 0.09,  0.09,  0.09],
925                             [ 0.07,  0.07,  0.07],
926                             [ 0.07,  0.07,  0.07],
927                             [ 1.0,  1.0,  1.0],
928                             [ 1.0,  1.0,  1.0]]), msg
929       
930        domain.set_region([set_bottom_friction, set_top_friction])
931        assert num.allclose(domain.quantities['friction'].get_values(),
932                            [[ 0.09,  0.09,  0.09],
933                             [ 0.09,  0.09,  0.09],
934                             [ 0.07,  0.07,  0.07],
935                             [ 0.07,  0.07,  0.07],
936                             [ 1.0,  1.0,  1.0],
937                             [ 1.0,  1.0,  1.0]])
938
939        domain.set_region([set_all_friction])
940        assert num.allclose(domain.quantities['friction'].get_values(),
941                            [[ 10.09, 10.09, 10.09],
942                             [ 10.09, 10.09, 10.09],
943                             [ 10.07, 10.07, 10.07],
944                             [ 10.07, 10.07, 10.07],
945                             [ 11.0,  11.0,  11.0],
946                             [ 11.0,  11.0,  11.0]])
947                             
948    def test_rectangular_periodic_and_ghosts(self):
949
950        from mesh_factory import rectangular_periodic
951       
952
953        M=5
954        N=2
955        points, vertices, boundary, full_send_dict, ghost_recv_dict = rectangular_periodic(M, N)
956
957        assert num.allclose(ghost_recv_dict[0][0], [24, 25, 26, 27,  0,  1,  2,  3])
958        assert num.allclose(full_send_dict[0][0] , [ 4,  5,  6,  7, 20, 21, 22, 23])
959
960        conserved_quantities = ['quant1', 'quant2']
961        domain = Domain(points, vertices, boundary, conserved_quantities,
962                        full_send_dict=full_send_dict,
963                        ghost_recv_dict=ghost_recv_dict)
964
965
966
967
968        assert num.allclose(ghost_recv_dict[0][0], [24, 25, 26, 27,  0,  1,  2,  3])
969        assert num.allclose(full_send_dict[0][0] , [ 4,  5,  6,  7, 20, 21, 22, 23])
970
971        def xylocation(x,y):
972            return 15*x + 9*y
973
974       
975        domain.set_quantity('quant1',xylocation,location='centroids')
976        domain.set_quantity('quant2',xylocation,location='centroids')
977
978
979        assert num.allclose(domain.quantities['quant1'].centroid_values,
980                            [  0.5,   1.,   5.,    5.5,   3.5,   4.,    8.,    8.5,   6.5,  7.,   11.,   11.5,   9.5,
981                               10.,   14.,   14.5,  12.5,  13.,   17.,   17.5,  15.5,  16.,   20.,   20.5,
982                               18.5,  19.,   23.,   23.5])
983
984
985
986        assert num.allclose(domain.quantities['quant2'].centroid_values,
987                            [  0.5,   1.,   5.,    5.5,   3.5,   4.,    8.,    8.5,   6.5,  7.,   11.,   11.5,   9.5,
988                               10.,   14.,   14.5,  12.5,  13.,   17.,   17.5,  15.5,  16.,   20.,   20.5,
989                               18.5,  19.,   23.,   23.5])
990
991        domain.update_ghosts()
992
993
994        assert num.allclose(domain.quantities['quant1'].centroid_values,
995                            [  15.5,  16.,   20.,   20.5,   3.5,   4.,    8.,    8.5,   6.5,  7.,   11.,   11.5,   9.5,
996                               10.,   14.,   14.5,  12.5,  13.,   17.,   17.5,  15.5,  16.,   20.,   20.5,
997                                3.5,   4.,    8.,    8.5])
998
999
1000
1001        assert num.allclose(domain.quantities['quant2'].centroid_values,
1002                            [  15.5,  16.,   20.,   20.5,   3.5,   4.,    8.,    8.5,   6.5,  7.,   11.,   11.5,   9.5,
1003                               10.,   14.,   14.5,  12.5,  13.,   17.,   17.5,  15.5,  16.,   20.,   20.5,
1004                                3.5,   4.,    8.,    8.5])
1005
1006       
1007        assert num.allclose(domain.tri_full_flag, [0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1008                                                   1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0])
1009
1010        #Test that points are arranged in a counter clock wise order
1011        domain.check_integrity()
1012
1013
1014    def test_that_mesh_methods_exist(self):
1015        """test_that_mesh_methods_exist
1016       
1017        Test that relavent mesh methods are made available in
1018        domain through composition
1019        """
1020        from mesh_factory import rectangular
1021        from shallow_water import Domain
1022
1023        # Create basic mesh
1024        points, vertices, boundary = rectangular(1, 3)
1025
1026        # Create shallow water domain
1027        domain = Domain(points, vertices, boundary)                             
1028       
1029       
1030        domain.get_centroid_coordinates()
1031        domain.get_radii()
1032        domain.get_areas()
1033        domain.get_area()
1034        domain.get_vertex_coordinates()
1035        domain.get_triangles()
1036        domain.get_nodes()
1037        domain.get_number_of_nodes()
1038        domain.get_normal(0,0)
1039        domain.get_triangle_containing_point([0.4,0.5])
1040        domain.get_intersecting_segments([[0.0, 0.0], [0.0, 1.0]])
1041        domain.get_disconnected_triangles()
1042        domain.get_boundary_tags()
1043        domain.get_boundary_polygon()
1044        #domain.get_number_of_triangles_per_node()
1045        domain.get_triangles_and_vertices_per_node()
1046        domain.get_interpolation_object()
1047        domain.get_tagged_elements()
1048        domain.get_lone_vertices()
1049        domain.get_unique_vertices()
1050        g = domain.get_georeference()
1051        domain.set_georeference(g)
1052        domain.build_tagged_elements_dictionary()
1053        domain.statistics()
1054        domain.get_extent()
1055
1056       
1057       
1058
1059#-------------------------------------------------------------
1060
1061if __name__ == "__main__":
1062    suite = unittest.makeSuite(Test_Domain,'test')
1063    runner = unittest.TextTestRunner()
1064    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.