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

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

Added input tests regarding polygons and points

File size: 16.9 KB
RevLine
[5435]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
[5570]4from anuga.utilities.polygon import inside_polygon
5from anuga.utilities.polygon import is_inside_polygon
[5435]6
[5570]7
[5434]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
[5437]99        if label is None: label = 'culvert_flow'
100        label += '_' + str(id(self))
[5434]101
102        # Open log file for storing some specific results...
103        self.log_filename = label + '.log'
104        self.label = label
105
[5437]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)
[5434]110
[5437]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
[5434]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')
[5570]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
[5434]143
[5570]144
145        # Create inflow object at each end of the culvert.
[5434]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
[5570]208
[5434]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):
[5570]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)
[5434]235
236            # Indices corresponding to energy enquiry field for this opening
[5453]237            coordinates = domain.get_centroid_coordinates(absolute=True) # Get all centroid points (x,y)
[5570]238            enquiry_indices = inside_polygon(coordinates,
239                                             self.enquiry_polygons[i])
[5453]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
[5434]244
245            # Get model values for points in enquiry polygon for this opening
[5570]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)
[5434]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)
[5570]260            print 'Velocity in culvert:', ux, uy, depth, xmomentum, ymomentum
[5434]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
[5437]278
[5434]279            # Estimated depth above the culvert inlet
[5437]280            d = w - z  # Used for calculations involving depth
[5434]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
[5437]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
[5434]297            # Ratio of depth to culvert height.
298            # If ratio > 1 then culvert is running full
[5437]299            if self.culvert_type == 'circle':
300                ratio = d/self.diameter
301            else:
302                ratio = d/self.height
[5434]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
[5437]311            # Store current average stage and depth with each opening object
[5434]312            opening.depth = d
[5437]313            opening.depth_trigger = d_trigger
[5434]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:
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.