source: anuga_core/source/anuga/shallow_water_balanced/test_swb_conservation.py @ 7733

Last change on this file since 7733 was 7733, checked in by hudson, 14 years ago

Fixed unit test failures.

File size: 31.5 KB
Line 
1#!/usr/bin/env python
2
3import unittest, os
4import os.path
5from math import pi, sqrt
6import tempfile
7
8from anuga.config import g, epsilon
9from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
10from anuga.utilities.numerical_tools import mean
11from anuga.utilities.polygon import is_inside_polygon
12from anuga.coordinate_transforms.geo_reference import Geo_reference
13from anuga.abstract_2d_finite_volumes.quantity import Quantity
14from anuga.geospatial_data.geospatial_data import Geospatial_data
15from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
16
17from anuga.utilities.system_tools import get_pathname_from_package
18from swb_domain import *
19
20import numpy as num
21
22# Get gateway to C implementation of flux function for direct testing
23from shallow_water_ext import flux_function_central as flux_function
24
25
26
27
28class Test_swb_conservation(unittest.TestCase):
29    def setUp(self):
30        pass
31
32    def tearDown(self):
33        pass
34
35
36    def test_conservation_1(self):
37        """Test that stage is conserved globally
38
39        This one uses a flat bed, reflective bdries and a suitable
40        initial condition
41        """
42
43        from mesh_factory import rectangular
44
45        # Create basic mesh
46        points, vertices, boundary = rectangular(6, 6)
47
48        # Create shallow water domain
49        domain = Domain(points, vertices, boundary)
50        domain.smooth = False
51        domain.default_order = 2
52
53        # IC
54        def x_slope(x, y):
55            return x/3
56
57        domain.set_quantity('elevation', 0)
58        domain.set_quantity('friction', 0)
59        domain.set_quantity('stage', x_slope)
60
61        # Boundary conditions (reflective everywhere)
62        Br = Reflective_boundary(domain)
63        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
64
65        domain.check_integrity()
66
67        initial_volume = domain.quantities['stage'].get_integral()
68        initial_xmom = domain.quantities['xmomentum'].get_integral()
69
70        # Evolution
71        for t in domain.evolve(yieldstep=0.05, finaltime=5.0):
72            volume = domain.quantities['stage'].get_integral()
73            assert num.allclose(volume, initial_volume)
74
75            #I don't believe that the total momentum should be the same
76            #It starts with zero and ends with zero though
77            #xmom = domain.quantities['xmomentum'].get_integral()
78            #print xmom
79            #assert allclose (xmom, initial_xmom)
80
81        os.remove(domain.get_name() + '.sww')
82
83    def test_conservation_2(self):
84        """Test that stage is conserved globally
85
86        This one uses a slopy bed, reflective bdries and a suitable
87        initial condition
88        """
89
90        from mesh_factory import rectangular
91
92        # Create basic mesh
93        points, vertices, boundary = rectangular(6, 6)
94
95        # Create shallow water domain
96        domain = Domain(points, vertices, boundary)
97        domain.smooth = False
98        domain.default_order = 2
99
100        # IC
101        def x_slope(x, y):
102            return x/3
103
104        domain.set_quantity('elevation', x_slope)
105        domain.set_quantity('friction', 0)
106        domain.set_quantity('stage', 0.4)    # Steady
107
108        # Boundary conditions (reflective everywhere)
109        Br = Reflective_boundary(domain)
110        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
111
112        domain.check_integrity()
113
114        initial_volume = domain.quantities['stage'].get_integral()
115        initial_xmom = domain.quantities['xmomentum'].get_integral()
116
117        # Evolution
118        for t in domain.evolve(yieldstep=0.05, finaltime=5.0):
119            volume = domain.quantities['stage'].get_integral()
120            assert num.allclose(volume, initial_volume)
121
122            #FIXME: What would we expect from momentum
123            #xmom = domain.quantities['xmomentum'].get_integral()
124            #print xmom
125            #assert allclose (xmom, initial_xmom)
126
127        os.remove(domain.get_name() + '.sww')
128
129    def test_conservation_3(self):
130        """Test that stage is conserved globally
131
132        This one uses a larger grid, convoluted bed, reflective boundaries
133        and a suitable initial condition
134        """
135
136        from mesh_factory import rectangular
137
138        # Create basic mesh
139        points, vertices, boundary = rectangular(2, 1)
140
141        # Create shallow water domain
142        domain = Domain(points, vertices, boundary)
143        domain.smooth = False
144
145
146        # IC
147        def x_slope(x, y):
148            z = 0*x
149            for i in range(len(x)):
150                if x[i] < 0.3:
151                    z[i] = x[i]/3
152                if 0.3 <= x[i] < 0.5:
153                    z[i] = -0.5
154                if 0.5 <= x[i] < 0.7:
155                    z[i] = 0.39
156                if 0.7 <= x[i]:
157                    z[i] = x[i]/3
158            return z
159
160        domain.set_quantity('elevation', x_slope)
161        domain.set_quantity('friction', 0)
162        domain.set_quantity('stage', 0.4) #Steady
163
164        # Boundary conditions (reflective everywhere)
165        Br = Reflective_boundary(domain)
166        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
167
168        domain.check_integrity()
169
170        initial_volume = domain.quantities['stage'].get_integral()
171        initial_xmom = domain.quantities['xmomentum'].get_integral()
172
173        import copy
174
175        ref_centroid_values = copy.copy(domain.quantities['stage'].\
176                                            centroid_values)
177
178        domain.distribute_to_vertices_and_edges()
179
180        assert num.allclose(domain.quantities['stage'].centroid_values,
181                            ref_centroid_values)
182
183        # Check that initial limiter doesn't violate cons quan
184        assert num.allclose(domain.quantities['stage'].get_integral(),
185                            initial_volume)
186
187        # Evolution
188        for t in domain.evolve(yieldstep=0.05, finaltime=10):
189            volume =  domain.quantities['stage'].get_integral()
190                       
191            assert num.allclose (volume, initial_volume)
192
193        os.remove(domain.get_name() + '.sww')
194
195    def test_conservation_4(self):
196        """Test that stage is conserved globally
197
198        This one uses a larger grid, convoluted bed, reflective boundaries
199        and a suitable initial condition
200        """
201
202        from mesh_factory import rectangular
203
204        # Create basic mesh
205        points, vertices, boundary = rectangular(6, 6)
206
207        # Create shallow water domain
208        domain = Domain(points, vertices, boundary)
209        domain.smooth = False
210        domain.default_order = 2
211
212        # IC
213        def x_slope(x, y):
214            z = 0*x
215            for i in range(len(x)):
216                if x[i] < 0.3:
217                    z[i] = x[i]/3
218                if 0.3 <= x[i] < 0.5:
219                    z[i] = -0.5
220                if 0.5 <= x[i] < 0.7:
221                    #z[i] = 0.3     # OK with beta == 0.2
222                    z[i] = 0.34     # OK with beta == 0.0
223                    #z[i] = 0.35    # Fails after 80 timesteps with an error
224                                    # of the order 1.0e-5
225                if 0.7 <= x[i]:
226                    z[i] = x[i]/3
227            return z
228
229        domain.set_quantity('elevation', x_slope)
230        domain.set_quantity('friction', 0)
231        domain.set_quantity('stage', 0.4) #Steady
232
233        # Boundary conditions (reflective everywhere)
234        Br = Reflective_boundary(domain)
235        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
236
237        domain.check_integrity()
238
239
240        initial_volume = domain.quantities['stage'].get_integral()
241        initial_xmom = domain.quantities['xmomentum'].get_integral()
242
243        import copy
244
245        ref_centroid_values = copy.copy(domain.quantities['stage'].\
246                                            centroid_values)
247
248        # Test limiter by itself
249        domain.distribute_to_vertices_and_edges()
250
251        # Check that initial limiter doesn't violate cons quan
252        assert num.allclose(domain.quantities['stage'].get_integral(),
253                            initial_volume)
254        # NOTE: This would fail if any initial stage was less than the
255        # corresponding bed elevation - but that is reasonable.
256
257        #Evolution
258        #print domain.get_time(), initial_volume
259        for t in domain.evolve(yieldstep=0.05, finaltime=10.0):
260            volume =  domain.quantities['stage'].get_integral()
261
262            #print domain.get_time(), volume
263            assert num.allclose (volume, initial_volume)
264
265        os.remove(domain.get_name() + '.sww')
266
267    def test_conservation_5(self):
268        """Test that momentum is conserved globally in steady state scenario
269
270        This one uses a slopy bed, dirichlet and reflective bdries
271        """
272
273        from mesh_factory import rectangular
274
275        # Create basic mesh
276        points, vertices, boundary = rectangular(6, 6)
277
278        # Create shallow water domain
279        domain = Domain(points, vertices, boundary)
280        domain.smooth = False
281        domain.default_order = 2
282
283        # IC
284        def x_slope(x, y):
285            return x/3
286
287        domain.set_quantity('elevation', x_slope)
288        domain.set_quantity('friction', 0)
289        domain.set_quantity('stage', 0.4) # Steady
290
291        # Boundary conditions (reflective everywhere)
292        Br = Reflective_boundary(domain)
293        Bleft = Dirichlet_boundary([0.5, 0, 0])
294        Bright = Dirichlet_boundary([0.1, 0, 0])
295        domain.set_boundary({'left': Bleft, 'right': Bright,
296                             'top': Br, 'bottom': Br})
297
298        domain.check_integrity()
299
300        initial_volume = domain.quantities['stage'].get_integral()
301        initial_xmom = domain.quantities['xmomentum'].get_integral()
302
303        # Evolution
304        for t in domain.evolve(yieldstep=0.05, finaltime=15.0):
305            stage =  domain.quantities['stage'].get_integral()
306            xmom = domain.quantities['xmomentum'].get_integral()
307            ymom = domain.quantities['ymomentum'].get_integral()
308
309            if num.allclose(t, 10):    # Steady state reached
310                steady_xmom = domain.quantities['xmomentum'].get_integral()
311                steady_ymom = domain.quantities['ymomentum'].get_integral()
312                steady_stage = domain.quantities['stage'].get_integral()
313
314            if t > 10:
315                msg = 'time=%.2f, xmom=%.10f, steady_xmom=%.10f' % (t, xmom, steady_xmom)
316                assert num.allclose(xmom, steady_xmom,atol=1.0e-4), msg
317
318                msg = 'time=%.2f, ymom=%.10f, steady_ymom=%.10f' % (t, ymom, steady_ymom)
319                assert num.allclose(ymom, steady_ymom,atol=1.0e-4), msg
320
321                msg = 'time=%.2f, stage=%.10f, steady_stage=%.10f' % (t, stage, steady_stage)
322                assert num.allclose(stage, steady_stage,atol=1.0e-4)
323
324        os.remove(domain.get_name() + '.sww')
325
326    def test_conservation_real(self):
327        """Test that momentum is conserved globally
328
329        Stephen finally made a test that revealed the problem.
330        This test failed with code prior to 25 July 2005
331        """
332
333        import sys
334        import os.path
335        sys.path.append(os.path.join('..', 'abstract_2d_finite_volumes'))
336        from mesh_factory import rectangular
337
338        yieldstep = 0.01
339        finaltime = 0.05
340        min_depth = 1.0e-2
341
342        #Create shallow water domain
343        points, vertices, boundary = rectangular(10, 10, len1=500, len2=500)
344        domain = Domain(points, vertices, boundary)
345        domain.smooth = False
346        domain.default_order = 1
347        domain.minimum_allowed_height = min_depth
348
349        # Set initial condition
350        class Set_IC:
351            """Set an initial condition with a constant value, for x0<x<x1"""
352
353            def __init__(self, x0=0.25, x1=0.5, h=1.0):
354                self.x0 = x0
355                self.x1 = x1
356                self.= h
357
358            def __call__(self, x, y):
359                return self.h*((x > self.x0) & (x < self.x1))
360
361        domain.set_quantity('stage', Set_IC(200.0, 300.0, 5.0))
362
363        # Boundaries
364        R = Reflective_boundary(domain)
365        domain.set_boundary({'left': R, 'right': R, 'top':R, 'bottom': R})
366
367        ref = domain.quantities['stage'].get_integral()
368
369        # Evolution
370        for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
371            pass
372
373        now = domain.quantities['stage'].get_integral()
374
375        msg = 'Stage not conserved: was %f, now %f' % (ref, now)
376        assert num.allclose(ref, now), msg
377
378        os.remove(domain.get_name() + '.sww')
379
380
381    def test_total_volume(self):       
382        """test_total_volume
383       
384        Test that total volume can be computed correctly
385        """           
386
387        #----------------------------------------------------------------------
388        # Import necessary modules
389        #----------------------------------------------------------------------
390        from anuga.abstract_2d_finite_volumes.mesh_factory \
391                import rectangular_cross
392        from anuga.shallow_water import Domain
393
394        #----------------------------------------------------------------------
395        # Setup computational domain
396        #----------------------------------------------------------------------
397
398        length = 100.
399        width  = 20.
400        dx = dy = 5       # Resolution: of grid on both axes
401       
402        points, vertices, boundary = rectangular_cross(int(length/dx),
403                                                       int(width/dy),
404                                                       len1=length,
405                                                       len2=width)
406        domain = Domain(points, vertices, boundary)   
407
408        #----------------------------------------------------------------------
409        # Simple flat bottom bathtub
410        #----------------------------------------------------------------------
411
412        d = 1.0
413        domain.set_quantity('elevation', 0.0)
414        domain.set_quantity('stage', d)
415       
416        assert num.allclose(domain.compute_total_volume(), length*width*d)
417
418        #----------------------------------------------------------------------
419        # Slope
420        #----------------------------------------------------------------------
421               
422        slope = 1.0/10          # RHS drops to -10m
423        def topography(x, y):
424            return -x * slope
425
426        domain.set_quantity('elevation', topography)
427        domain.set_quantity('stage', 0.0)       # Domain full
428       
429        V = domain.compute_total_volume()
430        assert num.allclose(V, length*width*10/2)
431
432        domain.set_quantity('stage', -5.0)      # Domain 'half' full
433       
434        # IMPORTANT: Adjust stage to match elevation
435        domain.distribute_to_vertices_and_edges()
436       
437        V = domain.compute_total_volume()
438        assert num.allclose(V, width*(length/2)*5.0/2)
439
440
441    def test_volumetric_balance_computation(self):
442        """test_volumetric_balance_computation
443       
444        Test that total in and out flows are computed correctly
445        in a steady state situation
446        """
447
448        # Set to True if volumetric output is sought
449        verbose = False
450
451        #----------------------------------------------------------------------
452        # Import necessary modules
453        #----------------------------------------------------------------------
454
455        from anuga.abstract_2d_finite_volumes.mesh_factory \
456                import rectangular_cross
457        from anuga.shallow_water import Domain
458        from anuga.shallow_water.shallow_water_domain import Reflective_boundary
459        from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
460        from anuga.shallow_water.forcing import Inflow
461        from anuga.shallow_water.data_manager \
462                import get_flow_through_cross_section
463
464        #----------------------------------------------------------------------
465        # Setup computational domain
466        #----------------------------------------------------------------------
467
468        finaltime = 500.0
469        length = 300.
470        width  = 20.
471        dx = dy = 5       # Resolution: of grid on both axes
472       
473        # Input parameters
474        uh = 1.0
475        vh = 0.0
476        d = 1.0
477       
478        # 20 m^3/s in the x direction across entire domain
479        ref_flow = uh*d*width
480
481        points, vertices, boundary = rectangular_cross(int(length/dx),
482                                                       int(width/dy),
483                                                       len1=length,
484                                                       len2=width)
485
486        domain = Domain(points, vertices, boundary)   
487        domain.set_name('Inflow_flowline_test')              # Output name
488
489        #----------------------------------------------------------------------
490        # Setup initial conditions
491        #----------------------------------------------------------------------
492
493        domain.set_quantity('elevation', 0.0)  # Flat bed
494        domain.set_quantity('friction', 0.0)   # Constant zero friction
495               
496        domain.set_quantity('stage', expression='elevation + %d' % d) 
497
498        #----------------------------------------------------------------------
499        # Setup boundary conditions
500        #----------------------------------------------------------------------
501
502        Br = Reflective_boundary(domain)      # Solid reflective wall
503               
504        # Constant flow in and out of domain
505        # Depth = 1m, uh=1 m/s, i.e. a flow of 20 m^3/s
506        Bi = Dirichlet_boundary([d, uh, vh]) 
507        Bo = Dirichlet_boundary([d, uh, vh])
508
509        domain.set_boundary({'left': Bi, 'right': Bo, 'top': Br, 'bottom': Br})
510
511        #----------------------------------------------------------------------
512        # Evolve system through time
513        #----------------------------------------------------------------------
514
515        for t in domain.evolve(yieldstep=50.0, finaltime=finaltime):
516            S = domain.volumetric_balance_statistics()
517            if verbose :
518                print domain.timestepping_statistics()
519                print S
520               
521            if t > 300:
522                # Steady state reached
523               
524                # Square on flowline at 200m
525                q = domain.get_flow_through_cross_section([[200.0,  0.0],
526                                                           [200.0, 20.0]])
527               
528                assert num.allclose(q, ref_flow)
529
530        os.remove('Inflow_flowline_test.sww') 
531
532    def test_volume_conservation_inflow(self):
533        """test_volume_conservation
534       
535        Test that total volume in domain is as expected, based on questions
536        raised by Petar Milevski in May 2009.
537       
538        This test adds inflow at a known rate and verifies that the total
539        terminal volume is as expected.
540       
541        """
542
543        verbose = False
544       
545
546        #---------------------------------------------------------------------
547        # Import necessary modules
548        #---------------------------------------------------------------------
549        from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
550        from anuga.shallow_water import Domain
551        from anuga.shallow_water.shallow_water_domain import Reflective_boundary
552        from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
553        from anuga.shallow_water.forcing import Inflow
554        from anuga.shallow_water.data_manager import get_flow_through_cross_section
555
556        #----------------------------------------------------------------------
557        # Setup computational domain
558        #----------------------------------------------------------------------
559        finaltime = 200.0
560
561        length = 300.
562        width  = 20.
563        dx = dy = 5       # Resolution: of grid on both axes
564       
565
566        points, vertices, boundary = rectangular_cross(int(length/dx), 
567                                                       int(width/dy),
568                                                       len1=length, len2=width)
569
570
571        domain = Domain(points, vertices, boundary)   
572        domain.set_name('Inflow_volume_test')              # Output name
573               
574
575        #----------------------------------------------------------------------
576        # Setup initial conditions
577        #----------------------------------------------------------------------
578        slope = 0.0
579        def topography(x, y):
580            z=-x * slope
581            return z
582
583        domain.set_quantity('elevation', topography) # Use function for elevation
584        domain.set_quantity('friction', 0.0)         # Constant friction
585               
586        domain.set_quantity('stage',
587                            expression='elevation')  # Dry initially
588                           
589
590        #--------------------------------------------------------------
591        # Setup Inflow
592        #--------------------------------------------------------------
593
594        # Fixed Flowrate onto Area
595        fixed_inflow = Inflow(domain,
596                              center=(10.0, 10.0),
597                              radius=5.00,
598                              rate=10.00)                               
599                           
600        domain.forcing_terms.append(fixed_inflow)                           
601       
602        #----------------------------------------------------------------------
603        # Setup boundary conditions
604        #----------------------------------------------------------------------
605
606        Br = Reflective_boundary(domain) # Solid reflective wall
607        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
608
609       
610        #----------------------------------------------------------------------
611        # Evolve system through time
612        #----------------------------------------------------------------------
613        ref_volume = 0.0
614        ys = 10.0  # Yieldstep
615        for t in domain.evolve(yieldstep=ys, finaltime=finaltime):
616       
617            # Check volume
618            assert num.allclose(domain.compute_total_volume(), ref_volume)
619       
620            if verbose :
621                print domain.timestepping_statistics()
622                print domain.volumetric_balance_statistics()
623                print 'reference volume', ref_volume
624           
625           
626            # Update reference volume
627            ref_volume += ys * fixed_inflow.rate
628
629
630        os.remove('Inflow_volume_test.sww')
631
632
633       
634    def test_volume_conservation_rain(self):
635        """test_volume_conservation
636       
637        Test that total volume in domain is as expected, based on questions
638        raised by Petar Milevski in May 2009.
639       
640        This test adds rain at a known rate and verifies that the total
641        terminal volume is as expected.
642       
643        """
644
645        verbose = False
646       
647
648        #---------------------------------------------------------------------
649        # Import necessary modules
650        #---------------------------------------------------------------------
651        from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
652        from anuga.shallow_water import Domain
653        from anuga.shallow_water.shallow_water_domain import Reflective_boundary
654        from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
655        from anuga.shallow_water.shallow_water_domain import Rainfall
656        from anuga.shallow_water.data_manager import get_flow_through_cross_section
657
658        #----------------------------------------------------------------------
659        # Setup computational domain
660        #----------------------------------------------------------------------
661        finaltime = 200.0
662
663        length = 300.
664        width  = 20.
665        dx = dy = 5       # Resolution: of grid on both axes
666       
667
668        points, vertices, boundary = rectangular_cross(int(length/dx), 
669                                                       int(width/dy),
670                                                       len1=length, len2=width)
671
672
673        domain = Domain(points, vertices, boundary)   
674        domain.set_name('Rain_volume_test')              # Output name
675               
676
677        #----------------------------------------------------------------------
678        # Setup initial conditions
679        #----------------------------------------------------------------------
680        slope = 0.0
681        def topography(x, y):
682            z=-x * slope
683            return z
684
685        domain.set_quantity('elevation', topography) # Use function for elevation
686        domain.set_quantity('friction', 0.0)         # Constant friction
687               
688        domain.set_quantity('stage',
689                            expression='elevation')  # Dry initially
690                           
691
692        #--------------------------------------------------------------
693        # Setup rain
694        #--------------------------------------------------------------
695
696        # Fixed rain onto small circular area
697        fixed_rain = Rainfall(domain,
698                              center=(10.0, 10.0),
699                              radius=5.00,
700                              rate=10.00)   # 10 mm/s                           
701                           
702        domain.forcing_terms.append(fixed_rain)                           
703       
704        #----------------------------------------------------------------------
705        # Setup boundary conditions
706        #----------------------------------------------------------------------
707
708        Br = Reflective_boundary(domain) # Solid reflective wall
709        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
710
711       
712        #----------------------------------------------------------------------
713        # Evolve system through time
714        #----------------------------------------------------------------------
715        ref_volume = 0.0
716        ys = 10.0  # Yieldstep
717        for t in domain.evolve(yieldstep=ys, finaltime=finaltime):
718       
719            # Check volume
720            V = domain.compute_total_volume()
721            msg = 'V = %e, Ref = %e' % (V, ref_volume)
722            assert num.allclose(V, ref_volume), msg
723       
724            if verbose :
725                print domain.timestepping_statistics()
726                print domain.volumetric_balance_statistics()
727                print 'reference volume', ref_volume
728                print V
729           
730           
731            # Update reference volume.
732            # FIXME: Note that rate has now been redefined
733            # as m/s internally. This is a little confusing
734            # when it was specfied as mm/s.
735           
736            delta_V = fixed_rain.rate*fixed_rain.exchange_area
737            ref_volume += ys * delta_V
738
739        os.remove('Rain_volume_test.sww')
740
741    def Xtest_rain_conservation_and_runoff(self):
742        """test_rain_conservation_and_runoff
743       
744        Test that total volume in domain is as expected, based on questions
745        raised by Petar Milevski in May 2009.
746       
747        This test adds rain at a known rate and verifies that the total
748        volume and outflows are as expected.
749       
750        """
751
752        # FIXME (Ole): Does not work yet. Investigate boundary flows
753       
754        verbose = True #False
755       
756
757        #---------------------------------------------------------------------
758        # Import necessary modules
759        #---------------------------------------------------------------------
760        from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
761        from anuga.shallow_water import Domain
762        from anuga.shallow_water.shallow_water_domain import Reflective_boundary
763        from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
764        from anuga.shallow_water.shallow_water_domain import Rainfall
765        from anuga.shallow_water.data_manager import get_flow_through_cross_section
766
767        #----------------------------------------------------------------------
768        # Setup computational domain
769        #----------------------------------------------------------------------
770        finaltime = 500.0
771
772        length = 300.
773        width  = 20.
774        dx = dy = 5       # Resolution: of grid on both axes
775       
776
777        points, vertices, boundary = rectangular_cross(int(length/dx), 
778                                                       int(width/dy),
779                                                       len1=length, len2=width)
780
781
782        domain = Domain(points, vertices, boundary)   
783        domain.set_name('Rain_volume_runoff_test')         # Output name
784               
785
786        #----------------------------------------------------------------------
787        # Setup initial conditions
788        #----------------------------------------------------------------------
789        slope = 0.0
790        def topography(x, y):
791            z=-x * slope
792            return z
793
794        domain.set_quantity('elevation', topography) # Use function for elevation
795        domain.set_quantity('friction', 0.0)         # Constant friction
796               
797        domain.set_quantity('stage',
798                            expression='elevation')  # Dry initially
799                           
800
801        #--------------------------------------------------------------
802        # Setup rain
803        #--------------------------------------------------------------
804
805        # Fixed rain onto small circular area
806        fixed_rain = Rainfall(domain,
807                              center=(10.0, 10.0),
808                              radius=5.00,
809                              rate=10.00)   # 10 mm/s                           
810                           
811        domain.forcing_terms.append(fixed_rain)                           
812       
813        #----------------------------------------------------------------------
814        # Setup boundary conditions
815        #----------------------------------------------------------------------
816
817        Br = Reflective_boundary(domain) # Solid reflective wall
818        Bt = Transmissive_stage_zero_momentum_boundary(domain)
819        Bd = Dirichlet_boundary([-10, 0, 0])
820        domain.set_boundary({'left': Bt, 'right': Bd, 'top': Bt, 'bottom': Bt})
821
822       
823        #----------------------------------------------------------------------
824        # Evolve system through time
825        #----------------------------------------------------------------------
826        ref_volume = 0.0
827        ys = 10.0  # Yieldstep
828        for t in domain.evolve(yieldstep=ys, finaltime=finaltime):
829       
830            # Check volume
831            V = domain.compute_total_volume()
832            msg = 'V = %e, Ref = %e' % (V, ref_volume)
833            #assert num.allclose(V, ref_volume) or V < ref_volume, msg
834       
835            if verbose:
836                print domain.timestepping_statistics()
837                print domain.volumetric_balance_statistics()
838                print 'reference volume', ref_volume
839                print V
840           
841           
842            # Update reference volume.
843            # FIXME: Note that rate has now been redefined
844            # as m/s internally. This is a little confusing
845            # when it was specfied as mm/s.
846           
847            delta_V = fixed_rain.rate*fixed_rain.exchange_area
848            ref_volume += ys * delta_V
849       
850            # Compute outflow at right hand downstream boundary
851            boundary_flows, inflow , outflow = domain.compute_boundary_flows()
852            net_outflow = outflow - inflow
853       
854            outflow = boundary_flows['right']
855            if verbose:
856                print 'Outflow', outflow
857                print 'Net outflow', net_outflow
858       
859            # Update reference volume
860            ref_volume += ys * outflow           
861
862
863
864#################################################################################
865
866if __name__ == "__main__":
867    suite = unittest.makeSuite(Test_swb_conservation, 'test')
868    runner = unittest.TextTestRunner(verbosity=1)
869    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.