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

Last change on this file since 6079 was 6079, checked in by ole, 16 years ago

Fixed dry bed culvert behaviour
File size: 12.8 KB
Line 
1#!/usr/bin/env python
2
3
4import unittest
5import os.path
6
7from anuga.utilities.system_tools import get_pathname_from_package       
8from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
9
10from anuga.shallow_water import Domain, Reflective_boundary,\
11    Dirichlet_boundary,\
12    Transmissive_boundary, Time_boundary
13
14from anuga.culvert_flows.culvert_class import Culvert_flow_rating, Culvert_flow_energy
15from anuga.culvert_flows.culvert_routines import boyd_generalised_culvert_model
16     
17from math import pi,pow,sqrt
18from Numeric import choose, greater, ones, sin, exp, cosh, allclose
19
20
21
22class Test_Culvert(unittest.TestCase):
23    def setUp(self):
24        pass
25
26    def tearDown(self):
27        pass
28
29
30    def test_that_culvert_runs(self):
31        """test_that_culvert_runs
32       
33        This test only exercises the culvert - there is no quantitative
34        diagnostics yet
35        """
36
37        path = get_pathname_from_package('anuga.culvert_flows')   
38       
39        length = 40.
40        width = 5.
41
42        dx = dy = 1           # Resolution: Length of subdivisions on both axes
43
44        points, vertices, boundary = rectangular_cross(int(length/dx),
45                                                       int(width/dy),
46                                                       len1=length, 
47                                                       len2=width)
48        domain = Domain(points, vertices, boundary)   
49        domain.set_name('Test_culvert')                 # Output name
50        domain.set_default_order(2)
51
52
53        #----------------------------------------------------------------------
54        # Setup initial conditions
55        #----------------------------------------------------------------------
56
57        def topography(x, y):
58            """Set up a weir
59           
60            A culvert will connect either side
61            """
62            # General Slope of Topography
63            z=-x/1000
64           
65            N = len(x)
66            for i in range(N):
67
68               # Sloping Embankment Across Channel
69                if 5.0 < x[i] < 10.1:
70                    # Cut Out Segment for Culvert FACE               
71                    if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0: 
72                       z[i]=z[i]
73                    else:
74                       z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
75                if 10.0 < x[i] < 12.1:
76                   z[i] +=  2.5                    # Flat Crest of Embankment
77                if 12.0 < x[i] < 14.5:
78                    # Cut Out Segment for Culvert FACE               
79                    if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5:
80                       z[i]=z[i]
81                    else:
82                       z[i] +=  2.5-1.0*(x[i] -12.0) # Sloping D/S Face
83                                   
84                       
85            return z
86
87
88        domain.set_quantity('elevation', topography) 
89        domain.set_quantity('friction', 0.01)         # Constant friction
90        domain.set_quantity('stage',
91                            expression='elevation')   # Dry initial condition
92
93        filename=os.path.join(path, 'example_rating_curve.csv')
94        culvert = Culvert_flow_rating(domain,
95                               culvert_description_filename=filename,       
96                               end_point0=[9.0, 2.5], 
97                               end_point1=[13.0, 2.5],
98                               verbose=False)
99
100        domain.forcing_terms.append(culvert)
101       
102
103        #-----------------------------------------------------------------------
104        # Setup boundary conditions
105        #-----------------------------------------------------------------------
106
107        # Inflow based on Flow Depth and Approaching Momentum
108        Bi = Dirichlet_boundary([0.0, 0.0, 0.0])
109        Br = Reflective_boundary(domain)              # Solid reflective wall
110        Bo = Dirichlet_boundary([-5, 0, 0])           # Outflow
111        Btus = Time_boundary(domain, lambda t: [0.0+ 1.25*(1+sin(2*pi*(t-4)/10)), 0.0, 0.0])
112        Btds = Time_boundary(domain, lambda t: [0.0+ 0.75*(1+sin(2*pi*(t-4)/20)), 0.0, 0.0])
113        domain.set_boundary({'left': Btus, 'right': Btds, 'top': Br, 'bottom': Br})
114
115
116        #-----------------------------------------------------------------------
117        # Evolve system through time
118        #-----------------------------------------------------------------------
119
120        for t in domain.evolve(yieldstep = 1, finaltime = 25):
121            #print domain.timestepping_statistics()
122            pass
123   
124
125
126    def test_that_culvert_dry_bed_rating(self):
127        """test_that_culvert_in_dry_bed_does_not_produce_flow(self):
128       
129        Test that culvert on a sloping dry bed doesn't produce flows
130        although there will be a 'pressure' head due to delta_w > 0
131
132        This one is using the rating curve variant
133        """
134
135        path = get_pathname_from_package('anuga.culvert_flows')   
136       
137        length = 40.
138        width = 5.
139
140        dx = dy = 1           # Resolution: Length of subdivisions on both axes
141
142        points, vertices, boundary = rectangular_cross(int(length/dx),
143                                                       int(width/dy),
144                                                       len1=length, 
145                                                       len2=width)
146        domain = Domain(points, vertices, boundary)   
147        domain.set_name('Test_culvert_dry') # Output name
148        domain.set_default_order(2)
149
150
151        #----------------------------------------------------------------------
152        # Setup initial conditions
153        #----------------------------------------------------------------------
154
155        def topography(x, y):
156            """Set up a weir
157           
158            A culvert will connect either side
159            """
160            # General Slope of Topography
161            z=-x/1000
162           
163            N = len(x)
164            for i in range(N):
165
166               # Sloping Embankment Across Channel
167                if 5.0 < x[i] < 10.1:
168                    # Cut Out Segment for Culvert FACE               
169                    if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0: 
170                       z[i]=z[i]
171                    else:
172                       z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
173                if 10.0 < x[i] < 12.1:
174                   z[i] +=  2.5                    # Flat Crest of Embankment
175                if 12.0 < x[i] < 14.5:
176                    # Cut Out Segment for Culvert FACE               
177                    if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5:
178                       z[i]=z[i]
179                    else:
180                       z[i] +=  2.5-1.0*(x[i] -12.0) # Sloping D/S Face
181                                   
182                       
183            return z
184
185
186        domain.set_quantity('elevation', topography) 
187        domain.set_quantity('friction', 0.01)         # Constant friction
188        domain.set_quantity('stage',
189                            expression='elevation')   # Dry initial condition
190
191
192        filename = os.path.join(path, 'example_rating_curve.csv')
193        culvert = Culvert_flow_rating(domain,
194                               culvert_description_filename=filename,
195                               end_point0=[9.0, 2.5], 
196                               end_point1=[13.0, 2.5],
197                               verbose=False)
198
199        domain.forcing_terms.append(culvert)
200       
201
202        #-----------------------------------------------------------------------
203        # Setup boundary conditions
204        #-----------------------------------------------------------------------
205
206        # Inflow based on Flow Depth and Approaching Momentum
207
208        Br = Reflective_boundary(domain)              # Solid reflective wall
209        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
210
211
212        #-----------------------------------------------------------------------
213        # Evolve system through time
214        #-----------------------------------------------------------------------
215
216        ref_volume = domain.get_quantity('stage').get_integral()
217        for t in domain.evolve(yieldstep = 1, finaltime = 25):
218            #print domain.timestepping_statistics()
219            new_volume = domain.get_quantity('stage').get_integral()
220           
221            msg = 'Total volume has changed'
222            assert allclose(new_volume, ref_volume), msg
223            pass
224   
225
226
227    def test_that_culvert_dry_bed_boyd(self):
228        """test_that_culvert_in_dry_bed_does_not_produce_flow(self):
229       
230        Test that culvert on a sloping dry bed doesn't produce flows
231        although there will be a 'pressure' head due to delta_w > 0
232
233        This one is using the 'Boyd' variant       
234        """
235
236        path = get_pathname_from_package('anuga.culvert_flows')   
237       
238        length = 40.
239        width = 5.
240
241        dx = dy = 1           # Resolution: Length of subdivisions on both axes
242
243        points, vertices, boundary = rectangular_cross(int(length/dx),
244                                                       int(width/dy),
245                                                       len1=length, 
246                                                       len2=width)
247        domain = Domain(points, vertices, boundary)   
248        domain.set_name('Test_culvert_dry') # Output name
249        domain.set_default_order(2)
250
251
252        #----------------------------------------------------------------------
253        # Setup initial conditions
254        #----------------------------------------------------------------------
255
256        def topography(x, y):
257            """Set up a weir
258           
259            A culvert will connect either side
260            """
261            # General Slope of Topography
262            z=-x/1000
263           
264            N = len(x)
265            for i in range(N):
266
267               # Sloping Embankment Across Channel
268                if 5.0 < x[i] < 10.1:
269                    # Cut Out Segment for Culvert FACE               
270                    if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0: 
271                       z[i]=z[i]
272                    else:
273                       z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
274                if 10.0 < x[i] < 12.1:
275                   z[i] +=  2.5                    # Flat Crest of Embankment
276                if 12.0 < x[i] < 14.5:
277                    # Cut Out Segment for Culvert FACE               
278                    if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5:
279                       z[i]=z[i]
280                    else:
281                       z[i] +=  2.5-1.0*(x[i] -12.0) # Sloping D/S Face
282                                   
283                       
284            return z
285
286
287        domain.set_quantity('elevation', topography) 
288        domain.set_quantity('friction', 0.01)         # Constant friction
289        domain.set_quantity('stage',
290                            expression='elevation')   # Dry initial condition
291
292
293        filename = os.path.join(path, 'example_rating_curve.csv')
294
295
296        culvert = Culvert_flow_energy(domain,
297                                      label='Culvert No. 1',
298                                      description='This culvert is a test unit 1.2m Wide by 0.75m High',   
299                                      end_point0=[9.0, 2.5], 
300                                      end_point1=[13.0, 2.5],
301                                      width=1.20,height=0.75,
302                                      culvert_routine=boyd_generalised_culvert_model,       
303                                      number_of_barrels=1,
304                                      update_interval=2,
305                                      verbose=True)
306       
307        domain.forcing_terms.append(culvert)
308       
309
310        #-----------------------------------------------------------------------
311        # Setup boundary conditions
312        #-----------------------------------------------------------------------
313
314        # Inflow based on Flow Depth and Approaching Momentum
315
316        Br = Reflective_boundary(domain)              # Solid reflective wall
317        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
318
319
320        #-----------------------------------------------------------------------
321        # Evolve system through time
322        #-----------------------------------------------------------------------
323
324        ref_volume = domain.get_quantity('stage').get_integral()
325        for t in domain.evolve(yieldstep = 1, finaltime = 25):
326            #print domain.timestepping_statistics()
327            new_volume = domain.get_quantity('stage').get_integral()
328           
329            msg = 'Total volume has changed'
330            assert allclose(new_volume, ref_volume), msg
331            pass
332   
333
334   
335               
336#-------------------------------------------------------------
337if __name__ == "__main__":
338    suite = unittest.makeSuite(Test_Culvert, 'test')
339    runner = unittest.TextTestRunner()
340    runner.run(suite)
341       
Note: See TracBrowser for help on using the repository browser.