source: trunk/anuga_work/shallow_water_balanced_steve/test_swb_conservation.py @ 8335

Last change on this file since 8335 was 7866, checked in by James Hudson, 15 years ago

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