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

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

Work on culvert testing. New test disabled due to negative total energy issue in Boyd routine.

File size: 16.4 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_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_rating(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(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_rating(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 test_that_culvert_dry_bed_boyd(self):
247        """test_that_culvert_in_dry_bed_does_not_produce_flow(self):
248
249        Test that culvert on a sloping dry bed doesn't produce flows
250        although there will be a 'pressure' head due to delta_w > 0
251
252        This one is using the 'Boyd' 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_dry') # 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')   # Dry initial condition
310
311
312        filename = os.path.join(path, 'example_rating_curve.csv')
313
314
315        culvert = Culvert_flow_energy(domain,
316                                      label='Culvert No. 1',
317                                      description='This culvert is a test unit 1.2m Wide by 0.75m High',
318                                      end_point0=[9.0, 2.5],
319                                      end_point1=[13.0, 2.5],
320                                      width=1.20,height=0.75,
321                                      culvert_routine=boyd_generalised_culvert_model,
322                                      number_of_barrels=1,
323                                      update_interval=2,
324                                      verbose=True)
325
326        domain.forcing_terms.append(culvert)
327
328
329        #-----------------------------------------------------------------------
330        # Setup boundary conditions
331        #-----------------------------------------------------------------------
332
333        # Inflow based on Flow Depth and Approaching Momentum
334
335        Br = Reflective_boundary(domain)              # Solid reflective wall
336        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
337
338
339        #-----------------------------------------------------------------------
340        # Evolve system through time
341        #-----------------------------------------------------------------------
342
343        ref_volume = domain.get_quantity('stage').get_integral()
344        for t in domain.evolve(yieldstep = 1, finaltime = 25):
345            #print domain.timestepping_statistics()
346            new_volume = domain.get_quantity('stage').get_integral()
347
348            msg = 'Total volume has changed'
349            assert allclose(new_volume, ref_volume), msg
350            pass
351
352
353
354
355
356    def Xtest_predicted_boyd_flow(self):
357        """test_predicted_boyd_flow
358
359        Test that flows predicted by the boyd method are consistent with what what
360        is calculated in engineering codes.
361        The data was supplied by Petar Milevski
362        """
363
364        path = get_pathname_from_package('anuga.culvert_flows')
365
366        length = 12.
367        width = 5.
368
369        dx = dy = 0.5           # Resolution: Length of subdivisions on both axes
370
371        points, vertices, boundary = rectangular_cross(int(length/dx),
372                                                       int(width/dy),
373                                                       len1=length,
374                                                       len2=width)
375        domain = Domain(points, vertices, boundary)
376
377        domain.set_name('test_culvert')                 # Output name
378        domain.set_default_order(2)
379
380
381        #----------------------------------------------------------------------
382        # Setup initial conditions
383        #----------------------------------------------------------------------
384
385        def topography(x, y):
386            # General Slope of Topography
387            z=-x/10
388
389            return z
390
391
392        domain.set_quantity('elevation', topography)
393        domain.set_quantity('friction', 0.01)         # Constant friction
394        domain.set_quantity('stage', expression='elevation')
395
396
397        Q0 = domain.get_quantity('stage')
398        Q1 = Quantity(domain)
399
400        # Add depths to stage
402        tail_water_depth = 0.089
403
404        inlet_poly = [[0,0], [6,0], [6,5], [0,5]]
405        outlet_poly = [[6,0], [12,0], [12,5], [6,5]]
406
408                                        (outlet_poly, tail_water_depth)]))
409
410        domain.set_quantity('stage', Q0 + Q1)
411
412
413
414        culvert = Culvert_flow_energy(domain,
415                                      label='Test culvert',
416                                      description='4 m test culvert',
417                                      end_point0=[4.0, 2.5],
418                                      end_point1=[8.0, 2.5],
419                                      width=1.20,
420                                      height=0.75,
421                                      culvert_routine=boyd_generalised_culvert_model,
422                                      number_of_barrels=1,
423                                      verbose=True)
424
425
426        domain.forcing_terms.append(culvert)
427
428        # Call
429        culvert(domain)
430
431        #print 'Inlet flow', culvert.inlet.rate
432        #print 'Outlet flow', culvert.outlet.rate
433
434
435
436#-------------------------------------------------------------
437if __name__ == "__main__":
438    #suite = unittest.makeSuite(Test_Culvert, 'test_predicted_boyd_flow')
439    suite = unittest.makeSuite(Test_Culvert, 'test')
440    runner = unittest.TextTestRunner()
441    runner.run(suite)
442
Note: See TracBrowser for help on using the repository browser.