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

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

Fixed georeferencing problem in culvert class

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