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

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

Implented general inlet/outlet control.
Works with rating curve but not Boyd.
File size: 20.5 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
22from Numeric import choose, greater, ones, sin, cos, exp, cosh, allclose
23
24
25
26class Test_Culvert(unittest.TestCase):
27    def setUp(self):
28        pass
29
30    def tearDown(self):
31        pass
32
33
34    def test_that_culvert_runs_rating(self):
35        """test_that_culvert_runs_rating
36       
37        This test exercises the culvert and checks values outside rating curve
38        are dealt with       
39        """
40
41        path = get_pathname_from_package('anuga.culvert_flows')   
42       
43        length = 40.
44        width = 5.
45
46        dx = dy = 1           # Resolution: Length of subdivisions on both axes
47
48        points, vertices, boundary = rectangular_cross(int(length/dx),
49                                                       int(width/dy),
50                                                       len1=length, 
51                                                       len2=width)
52        domain = Domain(points, vertices, boundary)   
53        domain.set_name('Test_culvert')                 # Output name
54        domain.set_default_order(2)
55
56
57        #----------------------------------------------------------------------
58        # Setup initial conditions
59        #----------------------------------------------------------------------
60
61        def topography(x, y):
62            """Set up a weir
63           
64            A culvert will connect either side
65            """
66            # General Slope of Topography
67            z=-x/1000
68           
69            N = len(x)
70            for i in range(N):
71
72               # Sloping Embankment Across Channel
73                if 5.0 < x[i] < 10.1:
74                    # Cut Out Segment for Culvert FACE               
75                    if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0: 
76                       z[i]=z[i]
77                    else:
78                       z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
79                if 10.0 < x[i] < 12.1:
80                   z[i] +=  2.5                    # Flat Crest of Embankment
81                if 12.0 < x[i] < 14.5:
82                    # Cut Out Segment for Culvert FACE               
83                    if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5:
84                       z[i]=z[i]
85                    else:
86                       z[i] +=  2.5-1.0*(x[i] -12.0) # Sloping D/S Face
87                                   
88                       
89            return z
90
91
92        domain.set_quantity('elevation', topography) 
93        domain.set_quantity('friction', 0.01)         # Constant friction
94        domain.set_quantity('stage',
95                            expression='elevation')   # Dry initial condition
96
97        filename=os.path.join(path, 'example_rating_curve.csv')
98        culvert = Culvert_flow(domain,
99                               culvert_description_filename=filename,       
100                               end_point0=[9.0, 2.5], 
101                               end_point1=[13.0, 2.5],
102                               verbose=False)
103
104        domain.forcing_terms.append(culvert)
105       
106
107        #-----------------------------------------------------------------------
108        # Setup boundary conditions
109        #-----------------------------------------------------------------------
110
111        # Inflow based on Flow Depth and Approaching Momentum
112        Bi = Dirichlet_boundary([0.0, 0.0, 0.0])
113        Br = Reflective_boundary(domain)              # Solid reflective wall
114        Bo = Dirichlet_boundary([-5, 0, 0])           # Outflow
115       
116        # Upstream and downstream conditions that will exceed the rating curve
117        # I.e produce delta_h outside the range [0, 10] specified in the the
118        # file example_rating_curve.csv
119        Btus = Time_boundary(domain, lambda t: [100*sin(2*pi*(t-4)/10), 0.0, 0.0])
120        Btds = Time_boundary(domain, lambda t: [-5*(cos(2*pi*(t-4)/20)), 0.0, 0.0])
121        domain.set_boundary({'left': Btus, 'right': Btds, 'top': Br, 'bottom': Br})
122
123
124        #-----------------------------------------------------------------------
125        # Evolve system through time
126        #-----------------------------------------------------------------------
127
128        min_delta_w = sys.maxint
129        max_delta_w = -min_delta_w
130        for t in domain.evolve(yieldstep = 1, finaltime = 25):
131            delta_w = culvert.inlet.stage - culvert.outlet.stage
132           
133            if delta_w > max_delta_w: max_delta_w = delta_w
134            if delta_w < min_delta_w: min_delta_w = delta_w           
135           
136            #print domain.timestepping_statistics()
137            pass
138
139        # Check that extreme values in rating curve have been exceeded
140        # so that we know that condition has been exercised
141        assert min_delta_w < 0
142        assert max_delta_w > 10       
143       
144
145    def test_that_culvert_dry_bed_rating_does_not_produce_flow(self):
146        """test_that_culvert_in_dry_bed_does_not_produce_flow(self):
147       
148        Test that culvert on a sloping dry bed doesn't produce flows
149        although there will be a 'pressure' head due to delta_w > 0
150
151        This one is using the rating curve variant
152        """
153
154        path = get_pathname_from_package('anuga.culvert_flows')   
155       
156        length = 40.
157        width = 5.
158
159        dx = dy = 1           # Resolution: Length of subdivisions on both axes
160
161        points, vertices, boundary = rectangular_cross(int(length/dx),
162                                                       int(width/dy),
163                                                       len1=length, 
164                                                       len2=width)
165        domain = Domain(points, vertices, boundary)   
166        domain.set_name('Test_culvert_dry') # Output name
167        domain.set_default_order(2)
168
169
170        #----------------------------------------------------------------------
171        # Setup initial conditions
172        #----------------------------------------------------------------------
173
174        def topography(x, y):
175            """Set up a weir
176           
177            A culvert will connect either side
178            """
179            # General Slope of Topography
180            z=-x/1000
181           
182            N = len(x)
183            for i in range(N):
184
185               # Sloping Embankment Across Channel
186                if 5.0 < x[i] < 10.1:
187                    # Cut Out Segment for Culvert FACE               
188                    if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0: 
189                       z[i]=z[i]
190                    else:
191                       z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
192                if 10.0 < x[i] < 12.1:
193                   z[i] +=  2.5                    # Flat Crest of Embankment
194                if 12.0 < x[i] < 14.5:
195                    # Cut Out Segment for Culvert FACE               
196                    if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5:
197                       z[i]=z[i]
198                    else:
199                       z[i] +=  2.5-1.0*(x[i] -12.0) # Sloping D/S Face
200                                   
201                       
202            return z
203
204
205        domain.set_quantity('elevation', topography) 
206        domain.set_quantity('friction', 0.01)         # Constant friction
207        domain.set_quantity('stage',
208                            expression='elevation')   # Dry initial condition
209
210
211        filename = os.path.join(path, 'example_rating_curve.csv')
212        culvert = Culvert_flow(domain,
213                               culvert_description_filename=filename,
214                               end_point0=[9.0, 2.5], 
215                               end_point1=[13.0, 2.5],
216                               verbose=False)
217
218        domain.forcing_terms.append(culvert)
219       
220
221        #-----------------------------------------------------------------------
222        # Setup boundary conditions
223        #-----------------------------------------------------------------------
224
225        # Inflow based on Flow Depth and Approaching Momentum
226
227        Br = Reflective_boundary(domain)              # Solid reflective wall
228        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
229
230
231        #-----------------------------------------------------------------------
232        # Evolve system through time
233        #-----------------------------------------------------------------------
234
235        ref_volume = domain.get_quantity('stage').get_integral()
236        for t in domain.evolve(yieldstep = 1, finaltime = 25):
237            #print domain.timestepping_statistics()
238            new_volume = domain.get_quantity('stage').get_integral()
239           
240            msg = 'Total volume has changed'
241            assert allclose(new_volume, ref_volume), msg
242            pass
243   
244
245   
246    def Xtest_that_culvert_rating_limits_flow_in_shallow_inlet_condition(self):
247        """test_that_culvert_rating_limits_flow_in_shallow_inlet_condition
248       
249        Test that culvert on a sloping dry bed limits flows when very little water
250        is present at inlet
251
252        This one is using the rating curve variant
253        """
254
255        path = get_pathname_from_package('anuga.culvert_flows')   
256       
257        length = 40.
258        width = 5.
259
260        dx = dy = 1           # Resolution: Length of subdivisions on both axes
261
262        points, vertices, boundary = rectangular_cross(int(length/dx),
263                                                       int(width/dy),
264                                                       len1=length, 
265                                                       len2=width)
266        domain = Domain(points, vertices, boundary)   
267        domain.set_name('Test_culvert_shallow') # Output name
268        domain.set_default_order(2)
269
270
271        #----------------------------------------------------------------------
272        # Setup initial conditions
273        #----------------------------------------------------------------------
274
275        def topography(x, y):
276            """Set up a weir
277           
278            A culvert will connect either side
279            """
280            # General Slope of Topography
281            z=-x/1000
282           
283            N = len(x)
284            for i in range(N):
285
286               # Sloping Embankment Across Channel
287                if 5.0 < x[i] < 10.1:
288                    # Cut Out Segment for Culvert FACE               
289                    if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0: 
290                       z[i]=z[i]
291                    else:
292                       z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
293                if 10.0 < x[i] < 12.1:
294                   z[i] +=  2.5                    # Flat Crest of Embankment
295                if 12.0 < x[i] < 14.5:
296                    # Cut Out Segment for Culvert FACE               
297                    if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5:
298                       z[i]=z[i]
299                    else:
300                       z[i] +=  2.5-1.0*(x[i] -12.0) # Sloping D/S Face
301                                   
302                       
303            return z
304
305
306        domain.set_quantity('elevation', topography) 
307        domain.set_quantity('friction', 0.01)         # Constant friction
308        domain.set_quantity('stage',
309                            expression='elevation + 0.2') # Shallow initial condition
310                           
311        # NOTE: Shallow values may cause this test to fail regardless of the
312        # culvert due to initial adjustments. A good value is 0.2
313
314
315        filename = os.path.join(path, 'example_rating_curve.csv')
316        culvert = Culvert_flow(domain,
317                               culvert_description_filename=filename,
318                               end_point0=[9.0, 2.5], 
319                               end_point1=[13.0, 2.5],
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 = 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            print domain.timestepping_statistics()
342            new_volume = domain.get_quantity('stage').get_integral()
343           
344            msg = 'Total volume has changed: Is %.2f m^3 should have been %.2f m^3'\
345                % (new_volume, ref_volume)
346            assert allclose(new_volume, ref_volume), msg
347   
348   
349   
350    def test_that_culvert_dry_bed_boyd_does_not_produce_flow(self):
351        """test_that_culvert_in_dry_bed_boyd_does_not_produce_flow(self):
352       
353        Test that culvert on a sloping dry bed doesn't produce flows
354        although there will be a 'pressure' head due to delta_w > 0
355
356        This one is using the 'Boyd' variant       
357        """
358
359        path = get_pathname_from_package('anuga.culvert_flows')   
360       
361        length = 40.
362        width = 5.
363
364        dx = dy = 1           # Resolution: Length of subdivisions on both axes
365
366        points, vertices, boundary = rectangular_cross(int(length/dx),
367                                                       int(width/dy),
368                                                       len1=length, 
369                                                       len2=width)
370        domain = Domain(points, vertices, boundary)   
371        domain.set_name('Test_culvert_dry') # Output name
372        domain.set_default_order(2)
373
374
375        #----------------------------------------------------------------------
376        # Setup initial conditions
377        #----------------------------------------------------------------------
378
379        def topography(x, y):
380            """Set up a weir
381           
382            A culvert will connect either side
383            """
384            # General Slope of Topography
385            z=-x/1000
386           
387            N = len(x)
388            for i in range(N):
389
390               # Sloping Embankment Across Channel
391                if 5.0 < x[i] < 10.1:
392                    # Cut Out Segment for Culvert FACE               
393                    if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0: 
394                       z[i]=z[i]
395                    else:
396                       z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
397                if 10.0 < x[i] < 12.1:
398                   z[i] +=  2.5                    # Flat Crest of Embankment
399                if 12.0 < x[i] < 14.5:
400                    # Cut Out Segment for Culvert FACE               
401                    if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5:
402                       z[i]=z[i]
403                    else:
404                       z[i] +=  2.5-1.0*(x[i] -12.0) # Sloping D/S Face
405                                   
406                       
407            return z
408
409
410        domain.set_quantity('elevation', topography) 
411        domain.set_quantity('friction', 0.01)         # Constant friction
412        domain.set_quantity('stage',
413                            expression='elevation')   # Dry initial condition
414
415
416        filename = os.path.join(path, 'example_rating_curve.csv')
417
418
419        culvert = Culvert_flow(domain,
420                               label='Culvert No. 1',
421                               description='This culvert is a test unit 1.2m Wide by 0.75m High',   
422                               end_point0=[9.0, 2.5], 
423                               end_point1=[13.0, 2.5],
424                               width=1.20,height=0.75,
425                               culvert_routine=boyd_generalised_culvert_model,
426                               number_of_barrels=1,
427                               update_interval=2,
428                               verbose=True)
429       
430        domain.forcing_terms.append(culvert)
431       
432
433        #-----------------------------------------------------------------------
434        # Setup boundary conditions
435        #-----------------------------------------------------------------------
436
437        # Inflow based on Flow Depth and Approaching Momentum
438
439        Br = Reflective_boundary(domain)              # Solid reflective wall
440        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
441
442
443        #-----------------------------------------------------------------------
444        # Evolve system through time
445        #-----------------------------------------------------------------------
446
447        ref_volume = domain.get_quantity('stage').get_integral()
448        for t in domain.evolve(yieldstep = 1, finaltime = 25):
449            #print domain.timestepping_statistics()
450            new_volume = domain.get_quantity('stage').get_integral()
451           
452            msg = 'Total volume has changed'
453            assert allclose(new_volume, ref_volume), msg
454            pass
455   
456
457   
458   
459
460    def Xtest_predicted_boyd_flow(self):
461        """test_predicted_boyd_flow
462       
463        Test that flows predicted by the boyd method are consistent with what what
464        is calculated in engineering codes.
465        The data was supplied by Petar Milevski
466        """
467
468        path = get_pathname_from_package('anuga.culvert_flows')   
469       
470        length = 12.
471        width = 5.
472
473        dx = dy = 0.5           # Resolution: Length of subdivisions on both axes
474
475        points, vertices, boundary = rectangular_cross(int(length/dx),
476                                                       int(width/dy),
477                                                       len1=length, 
478                                                       len2=width)
479        domain = Domain(points, vertices, boundary)   
480       
481        domain.set_name('test_culvert')                 # Output name
482        domain.set_default_order(2)
483
484
485        #----------------------------------------------------------------------
486        # Setup initial conditions
487        #----------------------------------------------------------------------
488
489        def topography(x, y):
490            # General Slope of Topography
491            z=-x/10
492           
493            return z
494
495
496        domain.set_quantity('elevation', topography) 
497        domain.set_quantity('friction', 0.01)         # Constant friction
498        domain.set_quantity('stage', expression='elevation')
499
500
501        Q0 = domain.get_quantity('stage')
502        Q1 = Quantity(domain)
503       
504        # Add depths to stage       
505        head_water_depth = 0.169
506        tail_water_depth = 0.089
507       
508        inlet_poly = [[0,0], [6,0], [6,5], [0,5]]
509        outlet_poly = [[6,0], [12,0], [12,5], [6,5]]       
510       
511        Q1.set_values(Polygon_function([(inlet_poly, head_water_depth),
512                                        (outlet_poly, tail_water_depth)]))
513       
514        domain.set_quantity('stage', Q0 + Q1)
515
516
517
518        culvert = Culvert_flow(domain,
519                               label='Test culvert',
520                               description='4 m test culvert',   
521                               end_point0=[4.0, 2.5], 
522                               end_point1=[8.0, 2.5],
523                               width=1.20, 
524                               height=0.75,
525                               culvert_routine=boyd_generalised_culvert_model,       
526                               number_of_barrels=1,
527                               verbose=True)
528                               
529       
530        domain.forcing_terms.append(culvert)
531       
532        # Call
533        culvert(domain)
534   
535        #print 'Inlet flow', culvert.inlet.rate
536        #print 'Outlet flow', culvert.outlet.rate       
537       
538   
539               
540#-------------------------------------------------------------
541if __name__ == "__main__":
542    #suite = unittest.makeSuite(Test_Culvert, 'test_that_culvert_rating_limits_flow_in_shallow_inlet_condition')
543    suite = unittest.makeSuite(Test_Culvert, 'test')
544    runner = unittest.TextTestRunner()
545    runner.run(suite)
546       
Note: See TracBrowser for help on using the repository browser.