source: trunk/anuga_core/source/anuga/structures/test_culvert_class.py @ 7939

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

Adding in culvert routines into structures directory

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