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

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

Updating the balanced and parallel code

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