source: branches/numpy/anuga/culvert_flows/test_culvert_class.py @ 6982

Last change on this file since 6982 was 6304, checked in by rwilson, 16 years ago

Initial commit of numpy changes. Still a long way to go.

File size: 21.3 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
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                               use_velocity_head=True,
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*num.sin(2*pi*(t-4)/10), 0.0, 0.0])
121        Btds = Time_boundary(domain, lambda t: [-5*(num.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 num.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.1') # Shallow initial condition
311                           
312
313        filename = os.path.join(path, 'example_rating_curve.csv')
314        culvert = Culvert_flow(domain,
315                               culvert_description_filename=filename,
316                               end_point0=[9.0, 2.5], 
317                               end_point1=[13.0, 2.5],
318                               trigger_depth=0.01,
319                               verbose=False)
320
321        domain.forcing_terms.append(culvert)
322       
323
324        #-----------------------------------------------------------------------
325        # Setup boundary conditions
326        #-----------------------------------------------------------------------
327
328        # Inflow based on Flow Depth and Approaching Momentum
329
330        Br = Reflective_boundary(domain)              # Solid reflective wall
331        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
332
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 = 0.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            if not num.allclose(new_volume, ref_volume):
347                print msg
348            assert num.allclose(new_volume, ref_volume), msg       
349       
350       
351       
352        # Now try this for a range of other depths
353        for depth in [1.0, 0.5, 0.2, 0.1, 0.05]:
354            domain.set_time(0.0)
355           
356            domain.set_quantity('stage',
357                                expression='elevation + %f' % depth)
358           
359       
360            ref_volume = domain.get_quantity('stage').get_integral()
361            for t in domain.evolve(yieldstep = 0.1, finaltime = 25):
362                #print domain.timestepping_statistics()
363                new_volume = domain.get_quantity('stage').get_integral()
364           
365                #print new_volume, ref_volume, new_volume-ref_volume
366                msg = 'Total volume has changed: Is %.2f m^3 should have been %.2f m^3'\
367                    % (new_volume, ref_volume)
368
369                assert num.allclose(new_volume, ref_volume), msg
370   
371   
372   
373    def test_that_culvert_dry_bed_boyd_does_not_produce_flow(self):
374        """test_that_culvert_in_dry_bed_boyd_does_not_produce_flow(self):
375       
376        Test that culvert on a sloping dry bed doesn't produce flows
377        although there will be a 'pressure' head due to delta_w > 0
378
379        This one is using the 'Boyd' variant       
380        """
381
382        path = get_pathname_from_package('anuga.culvert_flows')   
383       
384        length = 40.
385        width = 5.
386
387        dx = dy = 1           # Resolution: Length of subdivisions on both axes
388
389        points, vertices, boundary = rectangular_cross(int(length/dx),
390                                                       int(width/dy),
391                                                       len1=length, 
392                                                       len2=width)
393        domain = Domain(points, vertices, boundary)   
394        domain.set_name('Test_culvert_dry') # Output name
395        domain.set_default_order(2)
396
397
398        #----------------------------------------------------------------------
399        # Setup initial conditions
400        #----------------------------------------------------------------------
401
402        def topography(x, y):
403            """Set up a weir
404           
405            A culvert will connect either side
406            """
407            # General Slope of Topography
408            z=-x/1000
409           
410            N = len(x)
411            for i in range(N):
412
413               # Sloping Embankment Across Channel
414                if 5.0 < x[i] < 10.1:
415                    # Cut Out Segment for Culvert face               
416                    if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0: 
417                       z[i]=z[i]
418                    else:
419                       z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
420                if 10.0 < x[i] < 12.1:
421                   z[i] +=  2.5                    # Flat Crest of Embankment
422                if 12.0 < x[i] < 14.5:
423                    # Cut Out Segment for Culvert face               
424                    if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5:
425                       z[i]=z[i]
426                    else:
427                       z[i] +=  2.5-1.0*(x[i] -12.0) # Sloping D/S Face
428                                   
429                       
430            return z
431
432
433        domain.set_quantity('elevation', topography) 
434        domain.set_quantity('friction', 0.01)         # Constant friction
435        domain.set_quantity('stage',
436                            expression='elevation')   # Dry initial condition
437
438
439        filename = os.path.join(path, 'example_rating_curve.csv')
440
441
442        culvert = Culvert_flow(domain,
443                               label='Culvert No. 1',
444                               description='This culvert is a test unit 1.2m Wide by 0.75m High',   
445                               end_point0=[9.0, 2.5], 
446                               end_point1=[13.0, 2.5],
447                               width=1.20,height=0.75,
448                               culvert_routine=boyd_generalised_culvert_model,
449                               number_of_barrels=1,
450                               update_interval=2,
451                               verbose=True)
452       
453        domain.forcing_terms.append(culvert)
454       
455
456        #-----------------------------------------------------------------------
457        # Setup boundary conditions
458        #-----------------------------------------------------------------------
459
460        # Inflow based on Flow Depth and Approaching Momentum
461
462        Br = Reflective_boundary(domain)              # Solid reflective wall
463        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
464
465
466        #-----------------------------------------------------------------------
467        # Evolve system through time
468        #-----------------------------------------------------------------------
469
470        ref_volume = domain.get_quantity('stage').get_integral()
471        for t in domain.evolve(yieldstep = 1, finaltime = 25):
472           
473            #print domain.timestepping_statistics()
474            new_volume = domain.get_quantity('stage').get_integral()
475
476            msg = 'Total volume has changed'
477            assert num.allclose(new_volume, ref_volume), msg
478            pass
479   
480
481   
482   
483
484    def test_predicted_boyd_flow(self):
485        """test_predicted_boyd_flow
486       
487        Test that flows predicted by the boyd method are consistent with what what
488        is calculated in engineering codes.
489        The data was supplied by Petar Milevski
490        """
491
492        # FIXME(Ole) this is nowhere near finished
493        path = get_pathname_from_package('anuga.culvert_flows')   
494       
495        length = 12.
496        width = 5.
497
498        dx = dy = 0.5           # Resolution: Length of subdivisions on both axes
499
500        points, vertices, boundary = rectangular_cross(int(length/dx),
501                                                       int(width/dy),
502                                                       len1=length, 
503                                                       len2=width)
504        domain = Domain(points, vertices, boundary)   
505       
506        domain.set_name('test_culvert')                 # Output name
507        domain.set_default_order(2)
508
509
510        #----------------------------------------------------------------------
511        # Setup initial conditions
512        #----------------------------------------------------------------------
513
514        def topography(x, y):
515            # General Slope of Topography
516            z=-x/10
517           
518            return z
519
520
521        domain.set_quantity('elevation', topography) 
522        domain.set_quantity('friction', 0.01)         # Constant friction
523        domain.set_quantity('stage', expression='elevation')
524
525
526        Q0 = domain.get_quantity('stage')
527        Q1 = Quantity(domain)
528       
529        # Add depths to stage       
530        head_water_depth = 0.169
531        tail_water_depth = 0.089
532       
533        inlet_poly = [[0,0], [6,0], [6,5], [0,5]]
534        outlet_poly = [[6,0], [12,0], [12,5], [6,5]]       
535       
536        Q1.set_values(Polygon_function([(inlet_poly, head_water_depth),
537                                        (outlet_poly, tail_water_depth)]))
538       
539        domain.set_quantity('stage', Q0 + Q1)
540
541
542
543        culvert = Culvert_flow(domain,
544                               label='Test culvert',
545                               description='4 m test culvert',   
546                               end_point0=[4.0, 2.5], 
547                               end_point1=[8.0, 2.5],
548                               width=1.20, 
549                               height=0.75,
550                               culvert_routine=boyd_generalised_culvert_model,       
551                               number_of_barrels=1,
552                               verbose=True)
553                               
554       
555        domain.forcing_terms.append(culvert)
556       
557        # Call
558        culvert(domain)
559   
560        #print 'Inlet flow', culvert.inlet.rate
561        #print 'Outlet flow', culvert.outlet.rate       
562       
563   
564               
565#-------------------------------------------------------------
566if __name__ == "__main__":
567    suite = unittest.makeSuite(Test_Culvert, 'test')
568    runner = unittest.TextTestRunner() #verbosity=2)
569    runner.run(suite)
570       
Note: See TracBrowser for help on using the repository browser.