source: anuga_core/source/anuga/culvert_flows/test_culvert_class.py @ 7496

Last change on this file since 7496 was 7496, checked in by ole, 15 years ago

Fixed conflict

File size: 29.9 KB
Line 
1#!/usr/bin/env python
2
3
4import unittest
5import os.path
6import sys
7
8from anuga.utilities.system_tools import get_pathname_from_package
9from anuga.utilities.polygon import Polygon_function
10       
11from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
12from anuga.abstract_2d_finite_volumes.quantity import Quantity
13
14from anuga.shallow_water import Domain, Reflective_boundary,\
15    Dirichlet_boundary,\
16    Transmissive_boundary, Time_boundary
17
18from anuga.culvert_flows.culvert_class import Culvert_flow, Culvert_flow_rating, Culvert_flow_energy
19from anuga.culvert_flows.culvert_routines import boyd_generalised_culvert_model
20     
21from math import pi,pow,sqrt
22
23import numpy as num
24
25
26# Helper functions
27def run_culvert_flow_problem(depth):
28    """Run flow with culvert given depth
29    """
30
31    length = 40.
32    width = 5.
33
34    dx = dy = 1           # Resolution: Length of subdivisions on both axes
35
36    points, vertices, boundary = rectangular_cross(int(length/dx),
37                                                   int(width/dy),
38                                                   len1=length, 
39                                                   len2=width)
40    domain = Domain(points, vertices, boundary)   
41    domain.set_name('Test_culvert_shallow') # Output name
42    domain.set_default_order(2)
43
44
45    #----------------------------------------------------------------------
46    # Setup initial conditions
47    #----------------------------------------------------------------------
48
49    def topography(x, y):
50        """Set up a weir
51       
52        A culvert will connect either side
53        """
54        # General Slope of Topography
55        z=-x/1000
56       
57        N = len(x)
58        for i in range(N):
59
60           # Sloping Embankment Across Channel
61            if 5.0 < x[i] < 10.1:
62                # Cut Out Segment for Culvert face               
63                if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0: 
64                   z[i]=z[i]
65                else:
66                   z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
67            if 10.0 < x[i] < 12.1:
68               z[i] +=  2.5                    # Flat Crest of Embankment
69            if 12.0 < x[i] < 14.5:
70                # Cut Out Segment for Culvert face               
71                if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5:
72                   z[i]=z[i]
73                else:
74                   z[i] +=  2.5-1.0*(x[i] -12.0) # Sloping D/S Face
75                           
76               
77        return z
78
79
80    domain.set_quantity('elevation', topography) 
81    domain.set_quantity('friction', 0.01)         # Constant friction
82    domain.set_quantity('stage',
83                        expression='elevation + %f' % depth) # Shallow initial condition
84
85    # Boyd culvert
86    culvert = Culvert_flow(domain,
87                           label='Culvert No. 1',
88                           description='This culvert is a test unit 1.2m Wide by 0.75m High',   
89                           end_point0=[9.0, 2.5], 
90                           end_point1=[13.0, 2.5],
91                           width=1.20, height=0.75,
92                           culvert_routine=boyd_generalised_culvert_model,
93                           number_of_barrels=1,
94                           update_interval=2,
95                           verbose=False)
96   
97
98    domain.forcing_terms.append(culvert)
99   
100
101    #-----------------------------------------------------------------------
102    # Setup boundary conditions
103    #-----------------------------------------------------------------------
104
105    # Inflow based on Flow Depth and Approaching Momentum
106
107    Br = Reflective_boundary(domain)              # Solid reflective wall
108    domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
109
110
111   
112    #-----------------------------------------------------------------------
113    # Evolve system through time
114    #-----------------------------------------------------------------------
115
116    #print 'depth', depth
117    ref_volume = domain.get_quantity('stage').get_integral()
118    for t in domain.evolve(yieldstep = 0.1, finaltime = 25):
119        new_volume = domain.get_quantity('stage').get_integral()
120       
121        msg = ('Total volume has changed: Is %.8f m^3 should have been %.8f m^3'
122               % (new_volume, ref_volume))
123        assert num.allclose(new_volume, ref_volume), msg       
124       
125
126
127class Test_Culvert(unittest.TestCase):
128    def setUp(self):
129        pass
130
131    def tearDown(self):
132        pass
133
134
135    def test_that_culvert_runs_rating(self):
136        """test_that_culvert_runs_rating
137       
138        This test exercises the culvert and checks values outside rating curve
139        are dealt with       
140        """
141
142        path = get_pathname_from_package('anuga.culvert_flows')   
143       
144        length = 40.
145        width = 5.
146
147        dx = dy = 1           # Resolution: Length of subdivisions on both axes
148
149        points, vertices, boundary = rectangular_cross(int(length/dx),
150                                                       int(width/dy),
151                                                       len1=length, 
152                                                       len2=width)
153        domain = Domain(points, vertices, boundary)   
154        domain.set_name('Test_culvert')                 # Output name
155        domain.set_default_order(2)
156
157
158        #----------------------------------------------------------------------
159        # Setup initial conditions
160        #----------------------------------------------------------------------
161
162        def topography(x, y):
163            """Set up a weir
164           
165            A culvert will connect either side
166            """
167            # General Slope of Topography
168            z=-x/1000
169           
170            N = len(x)
171            for i in range(N):
172
173               # Sloping Embankment Across Channel
174                if 5.0 < x[i] < 10.1:
175                    # Cut Out Segment for Culvert face               
176                    if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0: 
177                       z[i]=z[i]
178                    else:
179                       z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
180                if 10.0 < x[i] < 12.1:
181                   z[i] +=  2.5                    # Flat Crest of Embankment
182                if 12.0 < x[i] < 14.5:
183                    # Cut Out Segment for Culvert face               
184                    if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5:
185                       z[i]=z[i]
186                    else:
187                       z[i] +=  2.5-1.0*(x[i] -12.0) # Sloping D/S Face
188                                   
189                       
190            return z
191
192
193        domain.set_quantity('elevation', topography) 
194        domain.set_quantity('friction', 0.01)         # Constant friction
195        domain.set_quantity('stage',
196                            expression='elevation')   # Dry initial condition
197
198        filename=os.path.join(path, 'example_rating_curve.csv')
199        culvert = Culvert_flow(domain,
200                               culvert_description_filename=filename,       
201                               end_point0=[9.0, 2.5], 
202                               end_point1=[13.0, 2.5],
203                               width=1.00,
204                               use_velocity_head=True,
205                               verbose=False)
206
207        domain.forcing_terms.append(culvert)
208       
209
210        #-----------------------------------------------------------------------
211        # Setup boundary conditions
212        #-----------------------------------------------------------------------
213
214        # Inflow based on Flow Depth and Approaching Momentum
215        Bi = Dirichlet_boundary([0.0, 0.0, 0.0])
216        Br = Reflective_boundary(domain)              # Solid reflective wall
217        Bo = Dirichlet_boundary([-5, 0, 0])           # Outflow
218       
219        # Upstream and downstream conditions that will exceed the rating curve
220        # I.e produce delta_h outside the range [0, 10] specified in the the
221        # file example_rating_curve.csv
222        Btus = Time_boundary(domain, lambda t: [100*num.sin(2*pi*(t-4)/10), 0.0, 0.0])
223        Btds = Time_boundary(domain, lambda t: [-5*(num.cos(2*pi*(t-4)/20)), 0.0, 0.0])
224        domain.set_boundary({'left': Btus, 'right': Btds, 'top': Br, 'bottom': Br})
225
226
227        #-----------------------------------------------------------------------
228        # Evolve system through time
229        #-----------------------------------------------------------------------
230
231        min_delta_w = sys.maxint
232        max_delta_w = -min_delta_w
233        for t in domain.evolve(yieldstep = 1, finaltime = 25):
234            delta_w = culvert.inlet.stage - culvert.outlet.stage
235           
236            if delta_w > max_delta_w: max_delta_w = delta_w
237            if delta_w < min_delta_w: min_delta_w = delta_w           
238           
239            pass
240
241        # Check that extreme values in rating curve have been exceeded
242        # so that we know that condition has been exercised
243        assert min_delta_w < 0
244        assert max_delta_w > 10       
245       
246
247    def test_that_culvert_dry_bed_rating_does_not_produce_flow(self):
248        """test_that_culvert_in_dry_bed_does_not_produce_flow(self):
249       
250        Test that culvert on a sloping dry bed doesn't produce flows
251        although there will be a 'pressure' head due to delta_w > 0
252
253        This one is using the rating curve variant
254        """
255
256        path = get_pathname_from_package('anuga.culvert_flows')   
257       
258        length = 40.
259        width = 5.
260
261        dx = dy = 1           # Resolution: Length of subdivisions on both axes
262
263        points, vertices, boundary = rectangular_cross(int(length/dx),
264                                                       int(width/dy),
265                                                       len1=length, 
266                                                       len2=width)
267        domain = Domain(points, vertices, boundary)   
268        domain.set_name('Test_culvert_dry') # Output name
269        domain.set_default_order(2)
270
271
272        #----------------------------------------------------------------------
273        # Setup initial conditions
274        #----------------------------------------------------------------------
275
276        def topography(x, y):
277            """Set up a weir
278           
279            A culvert will connect either side
280            """
281            # General Slope of Topography
282            z=-x/1000
283           
284            N = len(x)
285            for i in range(N):
286
287               # Sloping Embankment Across Channel
288                if 5.0 < x[i] < 10.1:
289                    # Cut Out Segment for Culvert face               
290                    if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0: 
291                       z[i]=z[i]
292                    else:
293                       z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
294                if 10.0 < x[i] < 12.1:
295                   z[i] +=  2.5                    # Flat Crest of Embankment
296                if 12.0 < x[i] < 14.5:
297                    # Cut Out Segment for Culvert face               
298                    if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5:
299                       z[i]=z[i]
300                    else:
301                       z[i] +=  2.5-1.0*(x[i] -12.0) # Sloping D/S Face
302                                   
303                       
304            return z
305
306
307        domain.set_quantity('elevation', topography) 
308        domain.set_quantity('friction', 0.01)         # Constant friction
309        domain.set_quantity('stage',
310                            expression='elevation')   # Dry initial condition
311
312
313        filename = os.path.join(path, 'example_rating_curve.csv')
314        culvert = Culvert_flow(domain,
315                               culvert_description_filename=filename,
316                               end_point0=[9.0, 2.5], 
317                               end_point1=[13.0, 2.5],
318                               height=0.75,
319                               verbose=False)
320
321        domain.forcing_terms.append(culvert)
322       
323
324        #-----------------------------------------------------------------------
325        # Setup boundary conditions
326        #-----------------------------------------------------------------------
327
328        # Inflow based on Flow Depth and Approaching Momentum
329
330        Br = Reflective_boundary(domain)              # Solid reflective wall
331        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
332
333
334        #-----------------------------------------------------------------------
335        # Evolve system through time
336        #-----------------------------------------------------------------------
337
338        ref_volume = domain.get_quantity('stage').get_integral()
339        for t in domain.evolve(yieldstep = 1, finaltime = 25):
340            new_volume = domain.get_quantity('stage').get_integral()
341           
342            msg = 'Total volume has changed'
343            assert num.allclose(new_volume, ref_volume, rtol=1.0e-10), msg
344            pass
345   
346
347   
348
349       
350    def test_that_culvert_flows_conserves_volume(self):
351        """test_that_culvert_flows_conserves_volume
352
353        Test that culvert on a sloping dry bed limits flows when very little water
354        is present at inlet.
355
356        Uses helper function: run_culvert_flow_problem(depth):
357
358        """
359
360        # Try this for a range of depths
361        for depth in [0.1, 0.2, 0.5, 1.0]:
362            run_culvert_flow_problem(depth)
363
364
365    def OBSOLETE_XXXtest_that_culvert_rating_limits_flow_in_shallow_inlet_condition(self):
366        """test_that_culvert_rating_limits_flow_in_shallow_inlet_condition
367       
368        Test that culvert on a sloping dry bed limits flows when very little water
369        is present at inlet
370
371        This one is using the rating curve variant
372        """
373
374       
375
376        path = get_pathname_from_package('anuga.culvert_flows')   
377       
378        length = 40.
379        width = 5.
380
381        dx = dy = 1           # Resolution: Length of subdivisions on both axes
382
383        points, vertices, boundary = rectangular_cross(int(length/dx),
384                                                       int(width/dy),
385                                                       len1=length, 
386                                                       len2=width)
387        domain = Domain(points, vertices, boundary)   
388        domain.set_name('Test_culvert_shallow') # Output name
389        domain.set_default_order(2)
390
391
392        #----------------------------------------------------------------------
393        # Setup initial conditions
394        #----------------------------------------------------------------------
395
396        def topography(x, y):
397            """Set up a weir
398           
399            A culvert will connect either side
400            """
401            # General Slope of Topography
402            z=-x/1000
403           
404            N = len(x)
405            for i in range(N):
406
407               # Sloping Embankment Across Channel
408                if 5.0 < x[i] < 10.1:
409                    # Cut Out Segment for Culvert face               
410                    if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0: 
411                       z[i]=z[i]
412                    else:
413                       z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
414                if 10.0 < x[i] < 12.1:
415                   z[i] +=  2.5                    # Flat Crest of Embankment
416                if 12.0 < x[i] < 14.5:
417                    # Cut Out Segment for Culvert face               
418                    if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5:
419                       z[i]=z[i]
420                    else:
421                       z[i] +=  2.5-1.0*(x[i] -12.0) # Sloping D/S Face
422                                   
423                       
424            return z
425
426
427        domain.set_quantity('elevation', topography) 
428        domain.set_quantity('friction', 0.01)         # Constant friction
429        domain.set_quantity('stage',
430                            expression='elevation + 0.1') # Shallow initial condition
431
432        # Boyd culvert
433        culvert = Culvert_flow(domain,
434                               label='Culvert No. 1',
435                               description='This culvert is a test unit 1.2m Wide by 0.75m High',   
436                               end_point0=[9.0, 2.5], 
437                               end_point1=[13.0, 2.5],
438                               width=1.20, height=0.75,
439                               culvert_routine=boyd_generalised_culvert_model,
440                               number_of_barrels=1,
441                               update_interval=2,
442                               verbose=False)
443       
444        # Rating curve
445        #filename = os.path.join(path, 'example_rating_curve.csv')
446        #culvert = Culvert_flow(domain,
447        #                       culvert_description_filename=filename,
448        #                       end_point0=[9.0, 2.5],
449        #                       end_point1=[13.0, 2.5],
450        #                       trigger_depth=0.01,
451        #                       verbose=False)
452
453        domain.forcing_terms.append(culvert)
454       
455
456        #-----------------------------------------------------------------------
457        # Setup boundary conditions
458        #-----------------------------------------------------------------------
459
460        # Inflow based on Flow Depth and Approaching Momentum
461
462        Br = Reflective_boundary(domain)              # Solid reflective wall
463        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
464
465
466       
467        #-----------------------------------------------------------------------
468        # Evolve system through time
469        #-----------------------------------------------------------------------
470
471        print 'depth', 0.1
472        ref_volume = domain.get_quantity('stage').get_integral()
473        for t in domain.evolve(yieldstep = 0.1, finaltime = 25):
474            new_volume = domain.get_quantity('stage').get_integral()
475
476            msg = ('Total volume has changed: Is %.8f m^3 should have been %.8f m^3'
477                   % (new_volume, ref_volume))
478            assert num.allclose(new_volume, ref_volume, rtol=1.0e-10), msg       
479       
480       
481        return
482        # Now try this again for a depth of 10 cm and for a range of other depths
483        for depth in [0.1, 0.2, 0.5, 1.0]:
484            print 'depth', depth
485            domain.set_time(0.0)
486           
487            domain.set_quantity('elevation', topography) 
488            domain.set_quantity('friction', 0.01)         # Constant friction
489            domain.set_quantity('stage',
490                                expression='elevation + %f' % depth)
491           
492       
493            ref_volume = domain.get_quantity('stage').get_integral()
494            for t in domain.evolve(yieldstep = 0.1, finaltime = 25):
495                new_volume = domain.get_quantity('stage').get_integral()
496           
497                msg = 'Total volume has changed: Is %.8f m^3 should have been %.8f m^3'\
498                    % (new_volume, ref_volume)
499
500                assert num.allclose(new_volume, ref_volume, rtol=1.0e-10), msg
501   
502   
503   
504    def test_that_culvert_dry_bed_boyd_does_not_produce_flow(self):
505        """test_that_culvert_in_dry_bed_boyd_does_not_produce_flow(self):
506       
507        Test that culvert on a sloping dry bed doesn't produce flows
508        although there will be a 'pressure' head due to delta_w > 0
509
510        This one is using the 'Boyd' variant       
511        """
512
513        path = get_pathname_from_package('anuga.culvert_flows')   
514       
515        length = 40.
516        width = 5.
517
518        dx = dy = 1           # Resolution: Length of subdivisions on both axes
519
520        points, vertices, boundary = rectangular_cross(int(length/dx),
521                                                       int(width/dy),
522                                                       len1=length, 
523                                                       len2=width)
524        domain = Domain(points, vertices, boundary)   
525        domain.set_name('Test_culvert_dry') # Output name
526        domain.set_default_order(2)
527
528
529        #----------------------------------------------------------------------
530        # Setup initial conditions
531        #----------------------------------------------------------------------
532
533        def topography(x, y):
534            """Set up a weir
535           
536            A culvert will connect either side
537            """
538            # General Slope of Topography
539            z=-x/1000
540           
541            N = len(x)
542            for i in range(N):
543
544               # Sloping Embankment Across Channel
545                if 5.0 < x[i] < 10.1:
546                    # Cut Out Segment for Culvert face               
547                    if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0: 
548                       z[i]=z[i]
549                    else:
550                       z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
551                if 10.0 < x[i] < 12.1:
552                   z[i] +=  2.5                    # Flat Crest of Embankment
553                if 12.0 < x[i] < 14.5:
554                    # Cut Out Segment for Culvert face               
555                    if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5:
556                       z[i]=z[i]
557                    else:
558                       z[i] +=  2.5-1.0*(x[i] -12.0) # Sloping D/S Face
559                                   
560                       
561            return z
562
563
564        domain.set_quantity('elevation', topography) 
565        domain.set_quantity('friction', 0.01)         # Constant friction
566        domain.set_quantity('stage',
567                            expression='elevation')   # Dry initial condition
568
569
570        filename = os.path.join(path, 'example_rating_curve.csv')
571
572
573        culvert = Culvert_flow(domain,
574                               label='Culvert No. 1',
575                               description='This culvert is a test unit 1.2m Wide by 0.75m High',   
576                               end_point0=[9.0, 2.5], 
577                               end_point1=[13.0, 2.5],
578                               width=1.20, height=0.75,
579                               culvert_routine=boyd_generalised_culvert_model,
580                               number_of_barrels=1,
581                               update_interval=2,
582                               verbose=False)
583       
584        domain.forcing_terms.append(culvert)
585       
586
587        #-----------------------------------------------------------------------
588        # Setup boundary conditions
589        #-----------------------------------------------------------------------
590
591        # Inflow based on Flow Depth and Approaching Momentum
592
593        Br = Reflective_boundary(domain)              # Solid reflective wall
594        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
595
596
597        #-----------------------------------------------------------------------
598        # Evolve system through time
599        #-----------------------------------------------------------------------
600
601        ref_volume = domain.get_quantity('stage').get_integral()
602        for t in domain.evolve(yieldstep = 1, finaltime = 25):
603           
604            new_volume = domain.get_quantity('stage').get_integral()
605
606            msg = 'Total volume has changed'
607            assert num.allclose(new_volume, ref_volume, rtol=1.0e-10), msg
608            pass
609   
610
611   
612   
613
614    def test_predicted_boyd_flow(self):
615        """test_predicted_boyd_flow
616       
617        Test that flows predicted by the boyd method are consistent with what what
618        is calculated in engineering codes.
619        The data was supplied by Petar Milevski
620        """
621
622        # FIXME(Ole) this is nowhere near finished
623        path = get_pathname_from_package('anuga.culvert_flows')   
624       
625        length = 12.
626        width = 5.
627
628        dx = dy = 0.5           # Resolution: Length of subdivisions on both axes
629
630        points, vertices, boundary = rectangular_cross(int(length/dx),
631                                                       int(width/dy),
632                                                       len1=length, 
633                                                       len2=width)
634        domain = Domain(points, vertices, boundary)   
635       
636        domain.set_name('test_culvert')                 # Output name
637        domain.set_default_order(2)
638
639
640        #----------------------------------------------------------------------
641        # Setup initial conditions
642        #----------------------------------------------------------------------
643
644        def topography(x, y):
645            # General Slope of Topography
646            z=-x/10
647           
648            return z
649
650
651        domain.set_quantity('elevation', topography) 
652        domain.set_quantity('friction', 0.01)         # Constant friction
653        domain.set_quantity('stage', expression='elevation')
654
655
656        Q0 = domain.get_quantity('stage')
657        Q1 = Quantity(domain)
658       
659        # Add depths to stage       
660        head_water_depth = 0.169
661        tail_water_depth = 0.089
662       
663        inlet_poly = [[0,0], [6,0], [6,5], [0,5]]
664        outlet_poly = [[6,0], [12,0], [12,5], [6,5]]       
665       
666        Q1.set_values(Polygon_function([(inlet_poly, head_water_depth),
667                                        (outlet_poly, tail_water_depth)]))
668       
669        domain.set_quantity('stage', Q0 + Q1)
670
671
672
673        culvert = Culvert_flow(domain,
674                               label='Test culvert',
675                               description='4 m test culvert',   
676                               end_point0=[4.0, 2.5], 
677                               end_point1=[8.0, 2.5],
678                               width=1.20, 
679                               height=0.75,
680                               culvert_routine=boyd_generalised_culvert_model,       
681                               number_of_barrels=1,
682                               verbose=False)
683                               
684       
685        domain.forcing_terms.append(culvert)
686       
687        # Call
688        culvert(domain)
689   
690
691   
692
693    def test_momentum_jet(self):
694        """test_momentum_jet
695       
696        Test that culvert_class can accept keyword use_momentum_jet
697        This does not yet imply that the values have been tested. FIXME
698        """
699
700
701        length = 40.
702        width = 5.
703
704        dx = dy = 1           # Resolution: Length of subdivisions on both axes
705
706        points, vertices, boundary = rectangular_cross(int(length/dx),
707                                                       int(width/dy),
708                                                       len1=length, 
709                                                       len2=width)
710        domain = Domain(points, vertices, boundary)   
711        domain.set_name('Test_culvert_shallow') # Output name
712        domain.set_default_order(2)
713
714
715        #----------------------------------------------------------------------
716        # Setup initial conditions
717        #----------------------------------------------------------------------
718
719        def topography(x, y):
720            """Set up a weir
721           
722            A culvert will connect either side
723            """
724            # General Slope of Topography
725            z=-x/1000
726           
727            N = len(x)
728            for i in range(N):
729
730               # Sloping Embankment Across Channel
731                if 5.0 < x[i] < 10.1:
732                    # Cut Out Segment for Culvert face               
733                    if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0: 
734                       z[i]=z[i]
735                    else:
736                       z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
737                if 10.0 < x[i] < 12.1:
738                   z[i] +=  2.5                    # Flat Crest of Embankment
739                if 12.0 < x[i] < 14.5:
740                    # Cut Out Segment for Culvert face               
741                    if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5:
742                       z[i]=z[i]
743                    else:
744                       z[i] +=  2.5-1.0*(x[i] -12.0) # Sloping D/S Face
745                                   
746                       
747            return z
748
749
750        domain.set_quantity('elevation', topography) 
751        domain.set_quantity('friction', 0.01)         # Constant friction
752        domain.set_quantity('stage',
753                            expression='elevation + 1.0') # Shallow initial condition
754
755        # Boyd culvert
756        culvert = Culvert_flow(domain,
757                               label='Culvert No. 1',
758                               description='This culvert is a test unit 1.2m Wide by 0.75m High',   
759                               end_point0=[9.0, 2.5], 
760                               end_point1=[13.0, 2.5],
761                               width=1.20, height=0.75,
762                               culvert_routine=boyd_generalised_culvert_model,
763                               number_of_barrels=1,
764                               use_momentum_jet=True,
765                               update_interval=2,
766                               verbose=False)
767   
768
769        domain.forcing_terms.append(culvert)
770
771       
772        # Call
773        culvert(domain)
774   
775
776        #-----------------------------------------------------------------------
777        # Setup boundary conditions
778        #-----------------------------------------------------------------------
779
780
781        Br = Reflective_boundary(domain)              # Solid reflective wall
782        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
783
784        #-----------------------------------------------------------------------
785        # Evolve system through time
786        #-----------------------------------------------------------------------
787
788        ref_volume = domain.get_quantity('stage').get_integral()
789        for t in domain.evolve(yieldstep = 0.1, finaltime = 25):
790            new_volume = domain.get_quantity('stage').get_integral()
791           
792            msg = ('Total volume has changed: Is %.8f m^3 should have been %.8f m^3'
793                   % (new_volume, ref_volume))
794            assert num.allclose(new_volume, ref_volume), msg       
795       
796
797       
798
799               
800#-------------------------------------------------------------
801
802if __name__ == "__main__":
803    suite = unittest.makeSuite(Test_Culvert, 'test')
804    runner = unittest.TextTestRunner() #verbosity=2)
805    runner.run(suite)
806       
Note: See TracBrowser for help on using the repository browser.