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

Last change on this file since 5570 was 5570, checked in by ole, 11 years ago

Added input tests regarding polygons and points

File size: 16.9 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            for point in P[key]:
137
138                print 'Passing in:', point
139                msg = 'Point %s in polygon %s for culvert %s did not'\
140                      %(str(point), key, self.label)
141                msg += 'fall within the domain boundary.'
142                assert is_inside_polygon(point, bounding_polygon), msg
143       
144
145        # Create inflow object at each end of the culvert.
146        self.openings = []
147        self.openings.append(Inflow(domain,
148                                    polygon=P['exchange_polygon0']))
149
150        self.openings.append(Inflow(domain,
151                                    polygon=P['exchange_polygon1']))                                   
152
153
154        # Assume two openings for now: Referred to as 0 and 1
155        assert len(self.openings) == 2
156       
157        # Store basic geometry
158        self.end_points = [end_point0, end_point1]
159        self.invert_levels = [invert_level0, invert_level1]               
160        self.enquiry_polygons = [P['enquiry_polygon0'], P['enquiry_polygon1']]
161        self.vector = P['vector']
162        self.length = P['length']; assert self.length > 0.0
163        self.verbose = verbose
164        self.last_time = 0.0       
165       
166
167        # Store hydraulic parameters
168        self.manning = manning
169        self.loss_exit = loss_exit
170        self.loss_entry = loss_entry
171        self.loss_bend = loss_bend
172        self.loss_special = loss_special
173        self.sum_loss = loss_exit + loss_entry + loss_bend + loss_special
174        self.blockage_topdwn = blockage_topdwn
175        self.blockage_bottup = blockage_bottup
176
177        # Store culvert routine
178        self.culvert_routine = culvert_routine
179
180       
181        # Create objects to update momentum (a bit crude at this stage)
182
183       
184        xmom0 = General_forcing(domain, 'xmomentum',
185                                polygon=P['exchange_polygon0'])
186
187        xmom1 = General_forcing(domain, 'xmomentum',
188                                polygon=P['exchange_polygon1'])
189
190        ymom0 = General_forcing(domain, 'ymomentum',
191                                polygon=P['exchange_polygon0'])
192
193        ymom1 = General_forcing(domain, 'ymomentum',
194                                polygon=P['exchange_polygon1'])
195
196        self.opening_momentum = [ [xmom0, ymom0], [xmom1, ymom1] ]
197       
198
199        # Print something
200        s = 'Culvert Effective Length = %.2f m' %(self.length)
201        log_to_file(self.log_filename, s)
202   
203        s = 'Culvert Direction is %s\n' %str(self.vector)
204        log_to_file(self.log_filename, s)       
205       
206    def __call__(self, domain):
207        from anuga.utilities.numerical_tools import mean
208       
209        from anuga.config import g, epsilon
210        from Numeric import take, sqrt
211        from anuga.config import velocity_protection       
212
213
214        log_filename = self.log_filename
215         
216        # Time stuff
217        time = domain.get_time()
218        delta_t = time-self.last_time
219        s = '\nTime = %.2f, delta_t = %f' %(time, delta_t)
220        log_to_file(log_filename, s)
221       
222        msg = 'Time did not advance'
223        if time > 0.0: assert delta_t > 0.0, msg
224
225
226        # Get average water depths at each opening       
227        openings = self.openings   # There are two Opening [0] and [1]
228        for i, opening in enumerate(openings):
229            dq = domain.quantities
230           
231            stage = dq['stage'].get_values(location='centroids',
232                                           indices=opening.exchange_indices)
233            elevation = dq['elevation'].get_values(location='centroids',
234                                                   indices=opening.exchange_indices)
235
236            # Indices corresponding to energy enquiry field for this opening
237            coordinates = domain.get_centroid_coordinates(absolute=True) # Get all centroid points (x,y)
238            enquiry_indices = inside_polygon(coordinates,
239                                             self.enquiry_polygons[i]) 
240
241            if len(enquiry_indices) == 0:
242                msg = 'No triangles have been identified in specified region: %s' %str(self.enquiry_polygons[i])
243                raise Exception, msg               
244           
245            # Get model values for points in enquiry polygon for this opening
246            stage = dq['stage'].get_values(location='centroids',
247                                           indices=enquiry_indices)
248            xmomentum = dq['xmomentum'].get_values(location='centroids',
249                                                   indices=enquiry_indices)
250            ymomentum = dq['ymomentum'].get_values(location='centroids',
251                                                   indices=enquiry_indices)
252            elevation = dq['elevation'].get_values(location='centroids',
253                                                   indices=enquiry_indices)
254            depth = stage - elevation
255
256            # Compute mean values of selected quantitites in the enquiry area in front of the culvert
257            # Epsilon handles a dry cell case
258            ux = xmomentum/(depth+velocity_protection/depth)   # Velocity (x-direction)
259            uy = ymomentum/(depth+velocity_protection/depth)   # Velocity (y-direction)
260            print 'Velocity in culvert:', ux, uy, depth, xmomentum, ymomentum
261            v = mean(sqrt(ux**2+uy**2))      # Average velocity
262            w = mean(stage)                  # Average stage
263
264            # Store values at enquiry field
265            opening.velocity = v
266
267
268            # Compute mean values of selected quantitites in the exchange area in front of the culvert
269            # Stage and velocity comes from enquiry area and elevation from exchange area
270           
271            # Use invert level instead of elevation if specified
272            invert_level = self.invert_levels[i]
273            if invert_level is not None:
274                z = invert_level
275            else:
276                elevation = dq['elevation'].get_values(location='centroids', indices=opening.exchange_indices)
277                z = mean(elevation)                   # Average elevation
278
279            # Estimated depth above the culvert inlet
280            d = w - z  # Used for calculations involving depth
281            if d < 0.0:
282                # This is possible since w and z are taken at different locations
283                #msg = 'D < 0.0: %f' %d
284                #raise msg
285                d = 0.0
286           
287
288
289            # Depth at exchange area used to trigger calculations
290            stage = dq['stage'].get_values(location='centroids', indices=enquiry_indices)
291            elevation = dq['elevation'].get_values(location='centroids', indices=enquiry_indices)
292            depth = stage - elevation
293            d_trigger = mean(depth)
294
295
296
297            # Ratio of depth to culvert height.
298            # If ratio > 1 then culvert is running full
299            if self.culvert_type == 'circle':
300                ratio = d/self.diameter
301            else:   
302                ratio = d/self.height 
303            opening.ratio = ratio
304               
305            # Average measures of energy in front of this opening
306            Es = d + 0.5*v**2/#  Specific energy in exchange area
307            Et = w + 0.5*v**2/#  Total energy in the enquiry area
308            opening.total_energy = Et
309            opening.specific_energy = Es           
310           
311            # Store current average stage and depth with each opening object
312            opening.depth = d
313            opening.depth_trigger = d_trigger           
314            opening.stage = w
315            opening.elevation = z
316           
317
318        #################  End of the FOR loop ################################################
319
320           
321        # We now need to deal with each opening individually
322
323        # Determine flow direction based on total energy difference
324        delta_Et = openings[0].total_energy - openings[1].total_energy
325
326        if delta_Et > 0:
327            #print 'Flow U/S ---> D/S'
328            inlet=openings[0]
329            outlet=openings[1]
330
331            inlet.momentum = self.opening_momentum[0]
332            outlet.momentum = self.opening_momentum[1]
333
334        else:
335            #print 'Flow D/S ---> U/S'
336            inlet=openings[1]
337            outlet=openings[0]
338
339            inlet.momentum = self.opening_momentum[1]
340            outlet.momentum = self.opening_momentum[0]
341           
342            delta_Et = -delta_Et
343
344        msg = 'Total energy difference is negative'
345        assert delta_Et > 0.0, msg
346
347        delta_h = inlet.stage - outlet.stage
348        delta_z = inlet.elevation - outlet.elevation
349        culvert_slope = (delta_z/self.length)
350
351        if culvert_slope < 0.0:
352            # Adverse gradient - flow is running uphill
353            # Flow will be purely controlled by uphill outlet face
354            print 'WARNING: Flow is running uphill. Watch Out!', inlet.elevation, outlet.elevation
355
356
357        s = 'Delta total energy = %.3f' %(delta_Et)
358        log_to_file(log_filename, s)
359
360       
361        Q, barrel_velocity, culvert_outlet_depth = self.culvert_routine(self, inlet, outlet, delta_Et, g)
362        #####################################################
363        barrel_momentum = barrel_velocity*culvert_outlet_depth
364               
365        s = 'Barrel velocity = %f' %barrel_velocity
366        log_to_file(log_filename, s)
367
368        # Compute momentum vector at outlet
369        outlet_mom_x, outlet_mom_y = self.vector * barrel_momentum
370           
371        s = 'Directional momentum = (%f, %f)' %(outlet_mom_x, outlet_mom_y)
372        log_to_file(log_filename, s)
373
374        delta_t = time - self.last_time
375        if delta_t > 0.0:
376            xmomentum_rate = outlet_mom_x - outlet.momentum[0].value
377            xmomentum_rate /= delta_t
378               
379            ymomentum_rate = outlet_mom_y - outlet.momentum[1].value
380            ymomentum_rate /= delta_t
381                   
382            s = 'X Y MOM_RATE = (%f, %f) ' %(xmomentum_rate, ymomentum_rate)
383            log_to_file(log_filename, s)                   
384        else:
385            xmomentum_rate = ymomentum_rate = 0.0
386
387
388        # Set momentum rates for outlet jet
389        outlet.momentum[0].rate = xmomentum_rate
390        outlet.momentum[1].rate = ymomentum_rate
391
392        # Remember this value for next step (IMPORTANT)
393        outlet.momentum[0].value = outlet_mom_x
394        outlet.momentum[1].value = outlet_mom_y                   
395
396        if int(domain.time*100) % 100 == 0:
397            s = 'T=%.5f, Culvert Discharge = %.3f f'\
398                %(time, Q)
399            s += ' Depth= %0.3f  Momentum = (%0.3f, %0.3f)'\
400                 %(culvert_outlet_depth, outlet_mom_x,outlet_mom_y)
401            s += ' Momentum rate: (%.4f, %.4f)'\
402                 %(xmomentum_rate, ymomentum_rate)                   
403            s+='Outlet Vel= %.3f'\
404                %(barrel_velocity)
405            log_to_file(log_filename, s)
406           
407               
408
409
410       
411        # Execute flow term for each opening
412        # This is where Inflow objects are evaluated and update the domain
413        for opening in self.openings:
414            opening(domain)
415
416        # Execute momentum terms
417        # This is where Inflow objects are evaluated and update the domain
418        outlet.momentum[0](domain)
419        outlet.momentum[1](domain)       
420           
421        # Store value of time   
422        self.last_time = time
423
424
Note: See TracBrowser for help on using the repository browser.