# source:anuga_core/source/anuga/culvert_flows/test_culvert_class.py@6128

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

Consolidated culverts - tests pass.
Outstanding issues are

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