source: anuga_core/source/anuga/shallow_water_balanced/test_swb_conservation.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: 31.5 KB
RevLine 
[7559]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()
[7573]190                       
[7559]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
[7573]239
[7559]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
[7573]258        #print domain.get_time(), initial_volume
[7559]259        for t in domain.evolve(yieldstep=0.05, finaltime=10.0):
260            volume =  domain.quantities['stage'].get_integral()
[7573]261
262            #print domain.get_time(), volume
[7559]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)
[7573]316                assert num.allclose(xmom, steady_xmom,atol=1.0e-4), msg
[7559]317
[7573]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
[7559]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.shallow_water_domain 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.shallow_water_domain 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.