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

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

Broke up test_swb_domain so that we can concentrate on new functionality
File size: 31.3 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        domain.default_order = 2
145
146
147        # IC
148        def x_slope(x, y):
149            z = 0*x
150            for i in range(len(x)):
151                if x[i] < 0.3:
152                    z[i] = x[i]/3
153                if 0.3 <= x[i] < 0.5:
154                    z[i] = -0.5
155                if 0.5 <= x[i] < 0.7:
156                    z[i] = 0.39
157                if 0.7 <= x[i]:
158                    z[i] = x[i]/3
159            return z
160
161        domain.set_quantity('elevation', x_slope)
162        domain.set_quantity('friction', 0)
163        domain.set_quantity('stage', 0.4) #Steady
164
165        # Boundary conditions (reflective everywhere)
166        Br = Reflective_boundary(domain)
167        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
168
169        domain.check_integrity()
170
171        initial_volume = domain.quantities['stage'].get_integral()
172        initial_xmom = domain.quantities['xmomentum'].get_integral()
173
174        import copy
175
176        ref_centroid_values = copy.copy(domain.quantities['stage'].\
177                                            centroid_values)
178
179        domain.distribute_to_vertices_and_edges()
180
181        assert num.allclose(domain.quantities['stage'].centroid_values,
182                            ref_centroid_values)
183
184        # Check that initial limiter doesn't violate cons quan
185        assert num.allclose(domain.quantities['stage'].get_integral(),
186                            initial_volume)
187
188        # Evolution
189        for t in domain.evolve(yieldstep=0.05, finaltime=10):
190            volume =  domain.quantities['stage'].get_integral()
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        initial_volume = domain.quantities['stage'].get_integral()
240        initial_xmom = domain.quantities['xmomentum'].get_integral()
241
242        import copy
243
244        ref_centroid_values = copy.copy(domain.quantities['stage'].\
245                                            centroid_values)
246
247        # Test limiter by itself
248        domain.distribute_to_vertices_and_edges()
249
250        # Check that initial limiter doesn't violate cons quan
251        assert num.allclose(domain.quantities['stage'].get_integral(),
252                            initial_volume)
253        # NOTE: This would fail if any initial stage was less than the
254        # corresponding bed elevation - but that is reasonable.
255
256        #Evolution
257        for t in domain.evolve(yieldstep=0.05, finaltime=10.0):
258            volume =  domain.quantities['stage'].get_integral()
259            assert num.allclose (volume, initial_volume)
260
261        os.remove(domain.get_name() + '.sww')
262
263    def test_conservation_5(self):
264        """Test that momentum is conserved globally in steady state scenario
265
266        This one uses a slopy bed, dirichlet and reflective bdries
267        """
268
269        from mesh_factory import rectangular
270
271        # Create basic mesh
272        points, vertices, boundary = rectangular(6, 6)
273
274        # Create shallow water domain
275        domain = Domain(points, vertices, boundary)
276        domain.smooth = False
277        domain.default_order = 2
278
279        # IC
280        def x_slope(x, y):
281            return x/3
282
283        domain.set_quantity('elevation', x_slope)
284        domain.set_quantity('friction', 0)
285        domain.set_quantity('stage', 0.4) # Steady
286
287        # Boundary conditions (reflective everywhere)
288        Br = Reflective_boundary(domain)
289        Bleft = Dirichlet_boundary([0.5, 0, 0])
290        Bright = Dirichlet_boundary([0.1, 0, 0])
291        domain.set_boundary({'left': Bleft, 'right': Bright,
292                             'top': Br, 'bottom': Br})
293
294        domain.check_integrity()
295
296        initial_volume = domain.quantities['stage'].get_integral()
297        initial_xmom = domain.quantities['xmomentum'].get_integral()
298
299        # Evolution
300        for t in domain.evolve(yieldstep=0.05, finaltime=15.0):
301            stage =  domain.quantities['stage'].get_integral()
302            xmom = domain.quantities['xmomentum'].get_integral()
303            ymom = domain.quantities['ymomentum'].get_integral()
304
305            if num.allclose(t, 10):    # Steady state reached
306                steady_xmom = domain.quantities['xmomentum'].get_integral()
307                steady_ymom = domain.quantities['ymomentum'].get_integral()
308                steady_stage = domain.quantities['stage'].get_integral()
309
310            if t > 10:
311                msg = 'time=%.2f, xmom=%.10f, steady_xmom=%.10f' % (t, xmom, steady_xmom)
312                assert num.allclose(xmom, steady_xmom), msg
313                msg = 'time=%.2f, ymom=%.10f, steady_ymom=%.10f' % (t, ymom, steady_ymom)               
314                assert num.allclose(ymom, steady_ymom), msg
315                assert num.allclose(stage, steady_stage)
316
317        os.remove(domain.get_name() + '.sww')
318
319    def test_conservation_real(self):
320        """Test that momentum is conserved globally
321
322        Stephen finally made a test that revealed the problem.
323        This test failed with code prior to 25 July 2005
324        """
325
326        import sys
327        import os.path
328        sys.path.append(os.path.join('..', 'abstract_2d_finite_volumes'))
329        from mesh_factory import rectangular
330
331        yieldstep = 0.01
332        finaltime = 0.05
333        min_depth = 1.0e-2
334
335        #Create shallow water domain
336        points, vertices, boundary = rectangular(10, 10, len1=500, len2=500)
337        domain = Domain(points, vertices, boundary)
338        domain.smooth = False
339        domain.default_order = 1
340        domain.minimum_allowed_height = min_depth
341
342        # Set initial condition
343        class Set_IC:
344            """Set an initial condition with a constant value, for x0<x<x1"""
345
346            def __init__(self, x0=0.25, x1=0.5, h=1.0):
347                self.x0 = x0
348                self.x1 = x1
349                self.h  = h
350
351            def __call__(self, x, y):
352                return self.h*((x > self.x0) & (x < self.x1))
353
354        domain.set_quantity('stage', Set_IC(200.0, 300.0, 5.0))
355
356        # Boundaries
357        R = Reflective_boundary(domain)
358        domain.set_boundary({'left': R, 'right': R, 'top':R, 'bottom': R})
359
360        ref = domain.quantities['stage'].get_integral()
361
362        # Evolution
363        for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
364            pass
365
366        now = domain.quantities['stage'].get_integral()
367
368        msg = 'Stage not conserved: was %f, now %f' % (ref, now)
369        assert num.allclose(ref, now), msg
370
371        os.remove(domain.get_name() + '.sww')
372
373
374    def test_total_volume(self):       
375        """test_total_volume
376       
377        Test that total volume can be computed correctly
378        """           
379
380        #----------------------------------------------------------------------
381        # Import necessary modules
382        #----------------------------------------------------------------------
383        from anuga.abstract_2d_finite_volumes.mesh_factory \
384                import rectangular_cross
385        from anuga.shallow_water import Domain
386
387        #----------------------------------------------------------------------
388        # Setup computational domain
389        #----------------------------------------------------------------------
390
391        length = 100.
392        width  = 20.
393        dx = dy = 5       # Resolution: of grid on both axes
394       
395        points, vertices, boundary = rectangular_cross(int(length/dx),
396                                                       int(width/dy),
397                                                       len1=length,
398                                                       len2=width)
399        domain = Domain(points, vertices, boundary)   
400
401        #----------------------------------------------------------------------
402        # Simple flat bottom bathtub
403        #----------------------------------------------------------------------
404
405        d = 1.0
406        domain.set_quantity('elevation', 0.0)
407        domain.set_quantity('stage', d)
408       
409        assert num.allclose(domain.compute_total_volume(), length*width*d)
410
411        #----------------------------------------------------------------------
412        # Slope
413        #----------------------------------------------------------------------
414               
415        slope = 1.0/10          # RHS drops to -10m
416        def topography(x, y):
417            return -x * slope
418
419        domain.set_quantity('elevation', topography)
420        domain.set_quantity('stage', 0.0)       # Domain full
421       
422        V = domain.compute_total_volume()
423        assert num.allclose(V, length*width*10/2)
424
425        domain.set_quantity('stage', -5.0)      # Domain 'half' full
426       
427        # IMPORTANT: Adjust stage to match elevation
428        domain.distribute_to_vertices_and_edges()
429       
430        V = domain.compute_total_volume()
431        assert num.allclose(V, width*(length/2)*5.0/2)
432
433
434    def test_volumetric_balance_computation(self):
435        """test_volumetric_balance_computation
436       
437        Test that total in and out flows are computed correctly
438        in a steady state situation
439        """
440
441        # Set to True if volumetric output is sought
442        verbose = False
443
444        #----------------------------------------------------------------------
445        # Import necessary modules
446        #----------------------------------------------------------------------
447
448        from anuga.abstract_2d_finite_volumes.mesh_factory \
449                import rectangular_cross
450        from anuga.shallow_water import Domain
451        from anuga.shallow_water.shallow_water_domain import Reflective_boundary
452        from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
453        from anuga.shallow_water.shallow_water_domain import Inflow
454        from anuga.shallow_water.data_manager \
455                import get_flow_through_cross_section
456
457        #----------------------------------------------------------------------
458        # Setup computational domain
459        #----------------------------------------------------------------------
460
461        finaltime = 500.0
462        length = 300.
463        width  = 20.
464        dx = dy = 5       # Resolution: of grid on both axes
465       
466        # Input parameters
467        uh = 1.0
468        vh = 0.0
469        d = 1.0
470       
471        # 20 m^3/s in the x direction across entire domain
472        ref_flow = uh*d*width
473
474        points, vertices, boundary = rectangular_cross(int(length/dx),
475                                                       int(width/dy),
476                                                       len1=length,
477                                                       len2=width)
478
479        domain = Domain(points, vertices, boundary)   
480        domain.set_name('Inflow_flowline_test')              # Output name
481
482        #----------------------------------------------------------------------
483        # Setup initial conditions
484        #----------------------------------------------------------------------
485
486        domain.set_quantity('elevation', 0.0)  # Flat bed
487        domain.set_quantity('friction', 0.0)   # Constant zero friction
488               
489        domain.set_quantity('stage', expression='elevation + %d' % d) 
490
491        #----------------------------------------------------------------------
492        # Setup boundary conditions
493        #----------------------------------------------------------------------
494
495        Br = Reflective_boundary(domain)      # Solid reflective wall
496               
497        # Constant flow in and out of domain
498        # Depth = 1m, uh=1 m/s, i.e. a flow of 20 m^3/s
499        Bi = Dirichlet_boundary([d, uh, vh]) 
500        Bo = Dirichlet_boundary([d, uh, vh])
501
502        domain.set_boundary({'left': Bi, 'right': Bo, 'top': Br, 'bottom': Br})
503
504        #----------------------------------------------------------------------
505        # Evolve system through time
506        #----------------------------------------------------------------------
507
508        for t in domain.evolve(yieldstep=50.0, finaltime=finaltime):
509            S = domain.volumetric_balance_statistics()
510            if verbose :
511                print domain.timestepping_statistics()
512                print S
513               
514            if t > 300:
515                # Steady state reached
516               
517                # Square on flowline at 200m
518                q = domain.get_flow_through_cross_section([[200.0,  0.0],
519                                                           [200.0, 20.0]])
520               
521                assert num.allclose(q, ref_flow)
522
523        os.remove('Inflow_flowline_test.sww') 
524
525    def test_volume_conservation_inflow(self):
526        """test_volume_conservation
527       
528        Test that total volume in domain is as expected, based on questions
529        raised by Petar Milevski in May 2009.
530       
531        This test adds inflow at a known rate and verifies that the total
532        terminal volume is as expected.
533       
534        """
535
536        verbose = False
537       
538
539        #---------------------------------------------------------------------
540        # Import necessary modules
541        #---------------------------------------------------------------------
542        from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
543        from anuga.shallow_water import Domain
544        from anuga.shallow_water.shallow_water_domain import Reflective_boundary
545        from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
546        from anuga.shallow_water.shallow_water_domain import Inflow
547        from anuga.shallow_water.data_manager import get_flow_through_cross_section
548
549        #----------------------------------------------------------------------
550        # Setup computational domain
551        #----------------------------------------------------------------------
552        finaltime = 200.0
553
554        length = 300.
555        width  = 20.
556        dx = dy = 5       # Resolution: of grid on both axes
557       
558
559        points, vertices, boundary = rectangular_cross(int(length/dx), 
560                                                       int(width/dy),
561                                                       len1=length, len2=width)
562
563
564        domain = Domain(points, vertices, boundary)   
565        domain.set_name('Inflow_volume_test')              # Output name
566               
567
568        #----------------------------------------------------------------------
569        # Setup initial conditions
570        #----------------------------------------------------------------------
571        slope = 0.0
572        def topography(x, y):
573            z=-x * slope
574            return z
575
576        domain.set_quantity('elevation', topography) # Use function for elevation
577        domain.set_quantity('friction', 0.0)         # Constant friction
578               
579        domain.set_quantity('stage',
580                            expression='elevation')  # Dry initially
581                           
582
583        #--------------------------------------------------------------
584        # Setup Inflow
585        #--------------------------------------------------------------
586
587        # Fixed Flowrate onto Area
588        fixed_inflow = Inflow(domain,
589                              center=(10.0, 10.0),
590                              radius=5.00,
591                              rate=10.00)                               
592                           
593        domain.forcing_terms.append(fixed_inflow)                           
594       
595        #----------------------------------------------------------------------
596        # Setup boundary conditions
597        #----------------------------------------------------------------------
598
599        Br = Reflective_boundary(domain) # Solid reflective wall
600        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
601
602       
603        #----------------------------------------------------------------------
604        # Evolve system through time
605        #----------------------------------------------------------------------
606        ref_volume = 0.0
607        ys = 10.0  # Yieldstep
608        for t in domain.evolve(yieldstep=ys, finaltime=finaltime):
609       
610            # Check volume
611            assert num.allclose(domain.compute_total_volume(), ref_volume)
612       
613            if verbose :
614                print domain.timestepping_statistics()
615                print domain.volumetric_balance_statistics()
616                print 'reference volume', ref_volume
617           
618           
619            # Update reference volume
620            ref_volume += ys * fixed_inflow.rate
621
622
623        os.remove('Inflow_volume_test.sww')
624
625
626       
627    def test_volume_conservation_rain(self):
628        """test_volume_conservation
629       
630        Test that total volume in domain is as expected, based on questions
631        raised by Petar Milevski in May 2009.
632       
633        This test adds rain at a known rate and verifies that the total
634        terminal volume is as expected.
635       
636        """
637
638        verbose = False
639       
640
641        #---------------------------------------------------------------------
642        # Import necessary modules
643        #---------------------------------------------------------------------
644        from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
645        from anuga.shallow_water import Domain
646        from anuga.shallow_water.shallow_water_domain import Reflective_boundary
647        from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
648        from anuga.shallow_water.shallow_water_domain import Rainfall
649        from anuga.shallow_water.data_manager import get_flow_through_cross_section
650
651        #----------------------------------------------------------------------
652        # Setup computational domain
653        #----------------------------------------------------------------------
654        finaltime = 200.0
655
656        length = 300.
657        width  = 20.
658        dx = dy = 5       # Resolution: of grid on both axes
659       
660
661        points, vertices, boundary = rectangular_cross(int(length/dx), 
662                                                       int(width/dy),
663                                                       len1=length, len2=width)
664
665
666        domain = Domain(points, vertices, boundary)   
667        domain.set_name('Rain_volume_test')              # Output name
668               
669
670        #----------------------------------------------------------------------
671        # Setup initial conditions
672        #----------------------------------------------------------------------
673        slope = 0.0
674        def topography(x, y):
675            z=-x * slope
676            return z
677
678        domain.set_quantity('elevation', topography) # Use function for elevation
679        domain.set_quantity('friction', 0.0)         # Constant friction
680               
681        domain.set_quantity('stage',
682                            expression='elevation')  # Dry initially
683                           
684
685        #--------------------------------------------------------------
686        # Setup rain
687        #--------------------------------------------------------------
688
689        # Fixed rain onto small circular area
690        fixed_rain = Rainfall(domain,
691                              center=(10.0, 10.0),
692                              radius=5.00,
693                              rate=10.00)   # 10 mm/s                           
694                           
695        domain.forcing_terms.append(fixed_rain)                           
696       
697        #----------------------------------------------------------------------
698        # Setup boundary conditions
699        #----------------------------------------------------------------------
700
701        Br = Reflective_boundary(domain) # Solid reflective wall
702        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
703
704       
705        #----------------------------------------------------------------------
706        # Evolve system through time
707        #----------------------------------------------------------------------
708        ref_volume = 0.0
709        ys = 10.0  # Yieldstep
710        for t in domain.evolve(yieldstep=ys, finaltime=finaltime):
711       
712            # Check volume
713            V = domain.compute_total_volume()
714            msg = 'V = %e, Ref = %e' % (V, ref_volume)
715            assert num.allclose(V, ref_volume), msg
716       
717            if verbose :
718                print domain.timestepping_statistics()
719                print domain.volumetric_balance_statistics()
720                print 'reference volume', ref_volume
721                print V
722           
723           
724            # Update reference volume.
725            # FIXME: Note that rate has now been redefined
726            # as m/s internally. This is a little confusing
727            # when it was specfied as mm/s.
728           
729            delta_V = fixed_rain.rate*fixed_rain.exchange_area
730            ref_volume += ys * delta_V
731
732        os.remove('Rain_volume_test.sww')
733
734    def Xtest_rain_conservation_and_runoff(self):
735        """test_rain_conservation_and_runoff
736       
737        Test that total volume in domain is as expected, based on questions
738        raised by Petar Milevski in May 2009.
739       
740        This test adds rain at a known rate and verifies that the total
741        volume and outflows are as expected.
742       
743        """
744
745        # FIXME (Ole): Does not work yet. Investigate boundary flows
746       
747        verbose = True #False
748       
749
750        #---------------------------------------------------------------------
751        # Import necessary modules
752        #---------------------------------------------------------------------
753        from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
754        from anuga.shallow_water import Domain
755        from anuga.shallow_water.shallow_water_domain import Reflective_boundary
756        from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
757        from anuga.shallow_water.shallow_water_domain import Rainfall
758        from anuga.shallow_water.data_manager import get_flow_through_cross_section
759
760        #----------------------------------------------------------------------
761        # Setup computational domain
762        #----------------------------------------------------------------------
763        finaltime = 500.0
764
765        length = 300.
766        width  = 20.
767        dx = dy = 5       # Resolution: of grid on both axes
768       
769
770        points, vertices, boundary = rectangular_cross(int(length/dx), 
771                                                       int(width/dy),
772                                                       len1=length, len2=width)
773
774
775        domain = Domain(points, vertices, boundary)   
776        domain.set_name('Rain_volume_runoff_test')         # Output name
777               
778
779        #----------------------------------------------------------------------
780        # Setup initial conditions
781        #----------------------------------------------------------------------
782        slope = 0.0
783        def topography(x, y):
784            z=-x * slope
785            return z
786
787        domain.set_quantity('elevation', topography) # Use function for elevation
788        domain.set_quantity('friction', 0.0)         # Constant friction
789               
790        domain.set_quantity('stage',
791                            expression='elevation')  # Dry initially
792                           
793
794        #--------------------------------------------------------------
795        # Setup rain
796        #--------------------------------------------------------------
797
798        # Fixed rain onto small circular area
799        fixed_rain = Rainfall(domain,
800                              center=(10.0, 10.0),
801                              radius=5.00,
802                              rate=10.00)   # 10 mm/s                           
803                           
804        domain.forcing_terms.append(fixed_rain)                           
805       
806        #----------------------------------------------------------------------
807        # Setup boundary conditions
808        #----------------------------------------------------------------------
809
810        Br = Reflective_boundary(domain) # Solid reflective wall
811        Bt = Transmissive_stage_zero_momentum_boundary(domain)
812        Bd = Dirichlet_boundary([-10, 0, 0])
813        domain.set_boundary({'left': Bt, 'right': Bd, 'top': Bt, 'bottom': Bt})
814
815       
816        #----------------------------------------------------------------------
817        # Evolve system through time
818        #----------------------------------------------------------------------
819        ref_volume = 0.0
820        ys = 10.0  # Yieldstep
821        for t in domain.evolve(yieldstep=ys, finaltime=finaltime):
822       
823            # Check volume
824            V = domain.compute_total_volume()
825            msg = 'V = %e, Ref = %e' % (V, ref_volume)
826            #assert num.allclose(V, ref_volume) or V < ref_volume, msg
827       
828            if verbose:
829                print domain.timestepping_statistics()
830                print domain.volumetric_balance_statistics()
831                print 'reference volume', ref_volume
832                print V
833           
834           
835            # Update reference volume.
836            # FIXME: Note that rate has now been redefined
837            # as m/s internally. This is a little confusing
838            # when it was specfied as mm/s.
839           
840            delta_V = fixed_rain.rate*fixed_rain.exchange_area
841            ref_volume += ys * delta_V
842       
843            # Compute outflow at right hand downstream boundary
844            boundary_flows, inflow , outflow = domain.compute_boundary_flows()
845            net_outflow = outflow - inflow
846       
847            outflow = boundary_flows['right']
848            if verbose:
849                print 'Outflow', outflow
850                print 'Net outflow', net_outflow
851       
852            # Update reference volume
853            ref_volume += ys * outflow           
854
855
856
857#################################################################################
858
859if __name__ == "__main__":
860    suite = unittest.makeSuite(Test_swb_conservation, 'test')
861    runner = unittest.TextTestRunner(verbosity=1)
862    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.