source: anuga_core/source/anuga/culvert_flows/culvert_class.py @ 5584

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

Testing that culvert polygons are inside mesh

File size: 17.1 KB
Line 
1from anuga.shallow_water.shallow_water_domain import Inflow, General_forcing
2from anuga.culvert_flows.culvert_polygons import create_culvert_polygons
3from anuga.utilities.system_tools import log_to_file
4from anuga.utilities.polygon import inside_polygon
5from anuga.utilities.polygon import is_inside_polygon
6
7
8class Culvert_flow:
9    """Culvert flow - transfer water from one hole to another
10   
11    Using Momentum as Calculated by Culvert Flow !!
12    Could be Several Methods Investigated to do This !!!
13
14    2008_May_08
15    To Ole:
16    OK so here we need to get the Polygon Creating code to create polygons for the culvert Based on
17    the two input Points (X0,Y0) and (X1,Y1) - need to be passed to create polygon
18
19    The two centers are now passed on to create_polygon.
20   
21
22    Input: Two points, pipe_size (either diameter or width, height), mannings_rougness,
23    inlet/outlet energy_loss_coefficients, internal_bend_coefficent,
24    top-down_blockage_factor and bottom_up_blockage_factor
25   
26   
27    And the Delta H enquiry should be change from Openings in line 412 to the enquiry Polygons infront
28    of the culvert
29    At the moment this script uses only Depth, later we can change it to Energy...
30
31    Once we have Delta H can calculate a Flow Rate and from Flow Rate an Outlet Velocity
32    The Outlet Velocity x Outlet Depth = Momentum to be applied at the Outlet...
33       
34    """ 
35
36    def __init__(self,
37                 domain,
38                 label=None, 
39                 description=None,
40                 end_point0=None, 
41                 end_point1=None,
42                 width=None,
43                 height=None,
44                 diameter=None,
45                 manning=None,          # Mannings Roughness for Culvert
46                 invert_level0=None,    # Invert level if not the same as the Elevation on the Domain
47                 invert_level1=None,    # Invert level if not the same as the Elevation on the Domain
48                 loss_exit=None,
49                 loss_entry=None,
50                 loss_bend=None,
51                 loss_special=None,
52                 blockage_topdwn=None,
53                 blockage_bottup=None,
54                 culvert_routine=None,
55                 verbose=False):
56       
57        from Numeric import sqrt, sum
58
59        # Input check
60        if diameter is not None:
61            self.culvert_type = 'circle'
62            self.diameter = diameter
63            if height is not None or width is not None:
64                msg = 'Either diameter or width&height must be specified, but not both.'
65                raise Exception, msg
66        else:
67            if height is not None:
68                if width is None:
69                    self.culvert_type = 'square'                               
70                    width = height
71                else:
72                    self.culvert_type = 'rectangle'
73            elif width is not None:
74                if height is None:
75                    self.culvert_type = 'square'                               
76                    height = width
77            else:
78                msg = 'Either diameter or width&height must be specified.'
79                raise Exception, msg               
80               
81            if height == width:
82                self.culvert_type = 'square'                                               
83               
84            self.height = height
85            self.width = width
86
87           
88        assert self.culvert_type in ['circle', 'square', 'rectangle']
89       
90        # Set defaults
91        if manning is None: manning = 0.012   # Set a Default Mannings Roughness for Pipe
92        if loss_exit is None: loss_exit = 1.00
93        if loss_entry is None: loss_entry = 0.50
94        if loss_bend is None: loss_bend=0.00
95        if loss_special is None: loss_special=0.00
96        if blockage_topdwn is None: blockage_topdwn=0.00
97        if blockage_bottup is None: blockage_bottup=0.00
98        if culvert_routine is None: culvert_routine=boyd_generalised_culvert_model
99        if label is None: label = 'culvert_flow'
100        label += '_' + str(id(self)) 
101       
102        # Open log file for storing some specific results...
103        self.log_filename = label + '.log' 
104        self.label = label
105
106        # Print something
107        log_to_file(self.log_filename, self.label)       
108        log_to_file(self.log_filename, description)
109        log_to_file(self.log_filename, self.culvert_type)       
110
111
112        # Create the fundamental culvert polygons from POLYGON
113        if self.culvert_type == 'circle':
114            # Redefine width and height for use with create_culvert_polygons
115            width = height = diameter
116       
117        P = create_culvert_polygons(end_point0,
118                                    end_point1,
119                                    width=width,   
120                                    height=height)
121       
122        if verbose is True:
123            pass
124            #plot_polygons([[end_point0, end_point1],
125            #               P['exchange_polygon0'],
126            #               P['exchange_polygon1'],
127            #               P['enquiry_polygon0'],
128            #               P['enquiry_polygon1']],
129            #              figname='culvert_polygon_output')
130
131
132        # Check that all polygons lie within the mesh.
133        bounding_polygon = domain.get_boundary_polygon()
134        for key in P.keys():
135            print 'Key', key
136            if key in ['exchange_polygon0', 
137                       'exchange_polygon1',
138                       'enquiry_polygon0',
139                       'enquiry_polygon1']:
140                for point in P[key]:
141               
142                    print 'Passing in:', point
143                    msg = 'Point %s in polygon %s for culvert %s did not'\
144                        %(str(point), key, self.label)
145                    msg += 'fall within the domain boundary.'
146                    assert is_inside_polygon(point, bounding_polygon), msg
147       
148
149        # Create inflow object at each end of the culvert.
150        self.openings = []
151        self.openings.append(Inflow(domain,
152                                    polygon=P['exchange_polygon0']))
153
154        self.openings.append(Inflow(domain,
155                                    polygon=P['exchange_polygon1']))                                   
156
157
158        # Assume two openings for now: Referred to as 0 and 1
159        assert len(self.openings) == 2
160       
161        # Store basic geometry
162        self.end_points = [end_point0, end_point1]
163        self.invert_levels = [invert_level0, invert_level1]               
164        self.enquiry_polygons = [P['enquiry_polygon0'], P['enquiry_polygon1']]
165        self.vector = P['vector']
166        self.length = P['length']; assert self.length > 0.0
167        self.verbose = verbose
168        self.last_time = 0.0       
169       
170
171        # Store hydraulic parameters
172        self.manning = manning
173        self.loss_exit = loss_exit
174        self.loss_entry = loss_entry
175        self.loss_bend = loss_bend
176        self.loss_special = loss_special
177        self.sum_loss = loss_exit + loss_entry + loss_bend + loss_special
178        self.blockage_topdwn = blockage_topdwn
179        self.blockage_bottup = blockage_bottup
180
181        # Store culvert routine
182        self.culvert_routine = culvert_routine
183
184       
185        # Create objects to update momentum (a bit crude at this stage)
186
187       
188        xmom0 = General_forcing(domain, 'xmomentum',
189                                polygon=P['exchange_polygon0'])
190
191        xmom1 = General_forcing(domain, 'xmomentum',
192                                polygon=P['exchange_polygon1'])
193
194        ymom0 = General_forcing(domain, 'ymomentum',
195                                polygon=P['exchange_polygon0'])
196
197        ymom1 = General_forcing(domain, 'ymomentum',
198                                polygon=P['exchange_polygon1'])
199
200        self.opening_momentum = [ [xmom0, ymom0], [xmom1, ymom1] ]
201       
202
203        # Print something
204        s = 'Culvert Effective Length = %.2f m' %(self.length)
205        log_to_file(self.log_filename, s)
206   
207        s = 'Culvert Direction is %s\n' %str(self.vector)
208        log_to_file(self.log_filename, s)       
209       
210    def __call__(self, domain):
211        from anuga.utilities.numerical_tools import mean
212       
213        from anuga.config import g, epsilon
214        from Numeric import take, sqrt
215        from anuga.config import velocity_protection       
216
217
218        log_filename = self.log_filename
219         
220        # Time stuff
221        time = domain.get_time()
222        delta_t = time-self.last_time
223        s = '\nTime = %.2f, delta_t = %f' %(time, delta_t)
224        log_to_file(log_filename, s)
225       
226        msg = 'Time did not advance'
227        if time > 0.0: assert delta_t > 0.0, msg
228
229
230        # Get average water depths at each opening       
231        openings = self.openings   # There are two Opening [0] and [1]
232        for i, opening in enumerate(openings):
233            dq = domain.quantities
234           
235            stage = dq['stage'].get_values(location='centroids',
236                                           indices=opening.exchange_indices)
237            elevation = dq['elevation'].get_values(location='centroids',
238                                                   indices=opening.exchange_indices)
239
240            # Indices corresponding to energy enquiry field for this opening
241            coordinates = domain.get_centroid_coordinates(absolute=True) # Get all centroid points (x,y)
242            enquiry_indices = inside_polygon(coordinates,
243                                             self.enquiry_polygons[i]) 
244
245            if len(enquiry_indices) == 0:
246                msg = 'No triangles have been identified in specified region: %s' %str(self.enquiry_polygons[i])
247                raise Exception, msg               
248           
249            # Get model values for points in enquiry polygon for this opening
250            stage = dq['stage'].get_values(location='centroids',
251                                           indices=enquiry_indices)
252            xmomentum = dq['xmomentum'].get_values(location='centroids',
253                                                   indices=enquiry_indices)
254            ymomentum = dq['ymomentum'].get_values(location='centroids',
255                                                   indices=enquiry_indices)
256            elevation = dq['elevation'].get_values(location='centroids',
257                                                   indices=enquiry_indices)
258            depth = stage - elevation
259
260            # Compute mean values of selected quantitites in the enquiry area in front of the culvert
261            # Epsilon handles a dry cell case
262            ux = xmomentum/(depth+velocity_protection/depth)   # Velocity (x-direction)
263            uy = ymomentum/(depth+velocity_protection/depth)   # Velocity (y-direction)
264            print 'Velocity in culvert:', ux, uy, depth, xmomentum, ymomentum
265            v = mean(sqrt(ux**2+uy**2))      # Average velocity
266            w = mean(stage)                  # Average stage
267
268            # Store values at enquiry field
269            opening.velocity = v
270
271
272            # Compute mean values of selected quantitites in the exchange area in front of the culvert
273            # Stage and velocity comes from enquiry area and elevation from exchange area
274           
275            # Use invert level instead of elevation if specified
276            invert_level = self.invert_levels[i]
277            if invert_level is not None:
278                z = invert_level
279            else:
280                elevation = dq['elevation'].get_values(location='centroids', indices=opening.exchange_indices)
281                z = mean(elevation)                   # Average elevation
282
283            # Estimated depth above the culvert inlet
284            d = w - z  # Used for calculations involving depth
285            if d < 0.0:
286                # This is possible since w and z are taken at different locations
287                #msg = 'D < 0.0: %f' %d
288                #raise msg
289                d = 0.0
290           
291
292
293            # Depth at exchange area used to trigger calculations
294            stage = dq['stage'].get_values(location='centroids', indices=enquiry_indices)
295            elevation = dq['elevation'].get_values(location='centroids', indices=enquiry_indices)
296            depth = stage - elevation
297            d_trigger = mean(depth)
298
299
300
301            # Ratio of depth to culvert height.
302            # If ratio > 1 then culvert is running full
303            if self.culvert_type == 'circle':
304                ratio = d/self.diameter
305            else:   
306                ratio = d/self.height 
307            opening.ratio = ratio
308               
309            # Average measures of energy in front of this opening
310            Es = d + 0.5*v**2/#  Specific energy in exchange area
311            Et = w + 0.5*v**2/#  Total energy in the enquiry area
312            opening.total_energy = Et
313            opening.specific_energy = Es           
314           
315            # Store current average stage and depth with each opening object
316            opening.depth = d
317            opening.depth_trigger = d_trigger           
318            opening.stage = w
319            opening.elevation = z
320           
321
322        #################  End of the FOR loop ################################################
323
324           
325        # We now need to deal with each opening individually
326
327        # Determine flow direction based on total energy difference
328        delta_Et = openings[0].total_energy - openings[1].total_energy
329
330        if delta_Et > 0:
331            #print 'Flow U/S ---> D/S'
332            inlet=openings[0]
333            outlet=openings[1]
334
335            inlet.momentum = self.opening_momentum[0]
336            outlet.momentum = self.opening_momentum[1]
337
338        else:
339            #print 'Flow D/S ---> U/S'
340            inlet=openings[1]
341            outlet=openings[0]
342
343            inlet.momentum = self.opening_momentum[1]
344            outlet.momentum = self.opening_momentum[0]
345           
346            delta_Et = -delta_Et
347
348        msg = 'Total energy difference is negative'
349        assert delta_Et > 0.0, msg
350
351        delta_h = inlet.stage - outlet.stage
352        delta_z = inlet.elevation - outlet.elevation
353        culvert_slope = (delta_z/self.length)
354
355        if culvert_slope < 0.0:
356            # Adverse gradient - flow is running uphill
357            # Flow will be purely controlled by uphill outlet face
358            print 'WARNING: Flow is running uphill. Watch Out!', inlet.elevation, outlet.elevation
359
360
361        s = 'Delta total energy = %.3f' %(delta_Et)
362        log_to_file(log_filename, s)
363
364       
365        Q, barrel_velocity, culvert_outlet_depth = self.culvert_routine(self, inlet, outlet, delta_Et, g)
366        #####################################################
367        barrel_momentum = barrel_velocity*culvert_outlet_depth
368               
369        s = 'Barrel velocity = %f' %barrel_velocity
370        log_to_file(log_filename, s)
371
372        # Compute momentum vector at outlet
373        outlet_mom_x, outlet_mom_y = self.vector * barrel_momentum
374           
375        s = 'Directional momentum = (%f, %f)' %(outlet_mom_x, outlet_mom_y)
376        log_to_file(log_filename, s)
377
378        delta_t = time - self.last_time
379        if delta_t > 0.0:
380            xmomentum_rate = outlet_mom_x - outlet.momentum[0].value
381            xmomentum_rate /= delta_t
382               
383            ymomentum_rate = outlet_mom_y - outlet.momentum[1].value
384            ymomentum_rate /= delta_t
385                   
386            s = 'X Y MOM_RATE = (%f, %f) ' %(xmomentum_rate, ymomentum_rate)
387            log_to_file(log_filename, s)                   
388        else:
389            xmomentum_rate = ymomentum_rate = 0.0
390
391
392        # Set momentum rates for outlet jet
393        outlet.momentum[0].rate = xmomentum_rate
394        outlet.momentum[1].rate = ymomentum_rate
395
396        # Remember this value for next step (IMPORTANT)
397        outlet.momentum[0].value = outlet_mom_x
398        outlet.momentum[1].value = outlet_mom_y                   
399
400        if int(domain.time*100) % 100 == 0:
401            s = 'T=%.5f, Culvert Discharge = %.3f f'\
402                %(time, Q)
403            s += ' Depth= %0.3f  Momentum = (%0.3f, %0.3f)'\
404                 %(culvert_outlet_depth, outlet_mom_x,outlet_mom_y)
405            s += ' Momentum rate: (%.4f, %.4f)'\
406                 %(xmomentum_rate, ymomentum_rate)                   
407            s+='Outlet Vel= %.3f'\
408                %(barrel_velocity)
409            log_to_file(log_filename, s)
410           
411               
412
413
414       
415        # Execute flow term for each opening
416        # This is where Inflow objects are evaluated and update the domain
417        for opening in self.openings:
418            opening(domain)
419
420        # Execute momentum terms
421        # This is where Inflow objects are evaluated and update the domain
422        outlet.momentum[0](domain)
423        outlet.momentum[1](domain)       
424           
425        # Store value of time   
426        self.last_time = time
427
428
Note: See TracBrowser for help on using the repository browser.