source: trunk/anuga_core/source/anuga/culvert_flows/test_new_culvert_class.py @ 7858

Last change on this file since 7858 was 7858, checked in by hudson, 13 years ago

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