source: trunk/anuga_core/source/anuga/structures/culvert_operator.py @ 7958

Last change on this file since 7958 was 7958, checked in by habili, 14 years ago

Corrections made to the code

File size: 26.2 KB
Line 
1import sys
2
3from anuga.shallow_water.forcing import Inflow, General_forcing
4from anuga.utilities.system_tools import log_to_file
5from anuga.geometry.polygon import inside_polygon, is_inside_polygon
6from anuga.geometry.polygon import plot_polygons, polygon_area
7
8
9from anuga.utilities.numerical_tools import mean
10from anuga.utilities.numerical_tools import ensure_numeric, sign
11       
12from anuga.config import g, epsilon
13from anuga.config import minimum_allowed_height, velocity_protection       
14import anuga.utilities.log as log
15
16import numpy as num
17from math import sqrt
18
19class Below_interval(Exception): pass 
20class Above_interval(Exception): pass
21
22   
23
24class Generic_box_culvert:
25    """Culvert flow - transfer water from one rectangular box to another.
26    Sets up the geometry of problem
27   
28    This is the base class for culverts. Inherit from this class (and overwrite
29    compute_discharge method for specific subclasses)
30   
31    Input: Two points, pipe_size (either diameter or width, height),
32    mannings_rougness,
33    """ 
34
35    def __init__(self,
36                 domain,
37                 end_point0=None, 
38                 end_point1=None,
39                 enquiry_gap_factor=0.2,
40                 width=None,
41                 height=None,
42                 verbose=False):
43       
44        # Input check
45       
46        self.domain = domain
47
48        self.end_points= [end_point0, end_point1]
49        self.enquiry_gap_factor = enquiry_gap_factor
50       
51        if height is None:
52            height = width
53
54        self.width = width
55        self.height = height
56       
57        self.verbose=verbose
58        self.filename = None
59       
60        # Create the fundamental culvert polygons from polygon
61        self.create_culvert_polygons()
62        self.compute_enquiry_indices()
63        self.check_culvert_inside_domain()
64
65        # Establish initial values at each enquiry point
66        self.enquiry_quantity_values = []
67        dq = domain.quantities
68        for i in [0, 1]:
69            idx = self.enquiry_indices[i]
70            elevation = dq['elevation'].get_values(location='centroids',
71                                                   indices=[idx])[0]
72            stage = dq['stage'].get_values(location='centroids',
73                                           indices=[idx])[0]
74            depth = stage - elevation
75           
76            quantity_values = {'stage' : stage, 'elevation' : elevation, 'depth' : depth }
77            self.enquiry_quantity_values.append(quantity_values)
78
79    def set_store_hydrograph_discharge(self, filename=None):
80
81        if filename is None:
82            self.filename = 'culvert_discharge_hydrograph'
83        else:
84            self.filename = filename
85
86        self.discharge_hydrograph = True
87       
88        self.timeseries_filename = self.filename + '_timeseries.csv'
89        fid = open(self.timeseries_filename, 'w')
90        fid.write('time, discharge\n')
91        fid.close()
92
93    def create_culvert_polygons(self):
94        """Create polygons at the end of a culvert inlet and outlet.
95        At either end two polygons will be created; one for the actual flow to pass through and one a little further away
96        for enquiring the total energy at both ends of the culvert and transferring flow.
97        """
98
99        # Calculate geometry
100        x0, y0 = self.end_points[0]
101        x1, y1 = self.end_points[1]
102
103        dx = x1 - x0
104        dy = y1 - y0
105
106        self.culvert_vector = num.array([dx, dy])
107        self.culvert_length = sqrt(num.sum(self.culvert_vector**2))
108        assert self.culvert_length > 0.0, 'The length of culvert is less than 0'
109
110        # Unit direction vector and normal
111        self.culvert_vector /= self.culvert_length                      # Unit vector in culvert direction
112        self.culvert_normal = num.array([-dy, dx])/self.culvert_length  # Normal vector
113
114        # Short hands
115        w = 0.5*self.width*self.culvert_normal # Perpendicular vector of 1/2 width
116        h = self.height*self.culvert_vector    # Vector of length=height in the
117                             # direction of the culvert
118        gap = (1 + self.enquiry_gap_factor)*h
119
120        self.exchange_polygons = []
121        self.enquiry_points = []
122
123        # Build exchange polygon and enquiry points 0 and 1
124        for i in [0, 1]:
125            p0 = self.end_points[i] + w
126            p1 = self.end_points[i] - w
127            p2 = p1 - h
128            p3 = p0 - h
129            self.exchange_polygons.append(num.array([p0, p1, p2, p3]))
130            self.enquiry_points.append(self.end_points[i] + (2*i-1)*gap)
131
132        self.polygon_areas = []
133
134        # Check that enquiry points are outside exchange polygons
135        for i in [0,1]:
136            polygon = self.exchange_polygons[i]
137            # FIXME(SR) should use area of triangles associated with polygon
138            area = polygon_area(polygon)
139            self.polygon_areas.append(area)
140
141            msg = 'Polygon %s ' %(polygon)
142            msg += ' has area = %f' % area
143            assert area > 0.0, msg
144
145            for j in [0,1]:
146                point = self.enquiry_points[j]
147                msg = 'Enquiry point falls inside a culvert polygon.'
148
149                assert not inside_polygon(point, polygon), msg
150
151       
152    def __call__(self, domain):
153
154        # Time stuff
155        time = domain.get_time()
156       
157       
158        update = False
159        if self.update_interval is None:
160            # Use next timestep as has been computed in domain.py       
161            delta_t = domain.timestep           
162            update = True
163        else:   
164            # Use update interval
165            delta_t = self.update_interval           
166            if time - self.last_update > self.update_interval or time == 0.0:
167                update = True
168           
169        if self.log_filename is not None:       
170            s = '\nTime = %.2f, delta_t = %f' %(time, delta_t)
171            log_to_file(self.log_filename, s)
172               
173                               
174        if update is True:
175            self.compute_rates(delta_t)
176           
177   
178        # Execute flow term for each opening
179        # This is where Inflow objects are evaluated using the last rate
180        # that has been calculated
181        #
182        # This will take place at every internal timestep and update the domain
183        for opening in self.openings:
184            opening(domain)
185           
186
187
188    def compute_enquiry_indices(self):
189        """Get indices for nearest centroids to self.enquiry_points
190        """
191       
192        domain = self.domain
193       
194        enquiry_indices = []                 
195        for point in self.enquiry_points:
196            # Find nearest centroid
197            N = len(domain)   
198            points = domain.get_centroid_coordinates(absolute=True)
199
200            # Calculate indices in exchange area for this forcing term
201           
202            triangle_id = min_dist = sys.maxint
203            for k in range(N):
204                x, y = points[k,:] # Centroid
205
206                c = point                               
207                distance = (x-c[0])**2+(y-c[1])**2
208                if distance < min_dist:
209                    min_dist = distance
210                    triangle_id = k
211
212                   
213            if triangle_id < sys.maxint:
214                msg = 'found triangle with centroid (%f, %f)'\
215                    %tuple(points[triangle_id, :])
216                msg += ' for point (%f, %f)' %tuple(point)
217               
218                enquiry_indices.append(triangle_id)
219            else:
220                msg = 'Triangle not found for point (%f, %f)' %point
221                raise Exception, msg
222       
223        self.enquiry_indices = enquiry_indices
224
225       
226    def check_culvert_inside_domain(self):
227        """Check that all polygons and enquiry points lie within the mesh.
228        """
229        bounding_polygon = self.domain.get_boundary_polygon()
230        for i in [0, 1]:
231            for point in list(self.exchange_polygons[i]) + self.enquiry_points:
232                msg = 'Point %s did not '\
233                    %(str(point))
234                msg += 'fall within the domain boundary.'
235                assert is_inside_polygon(point, bounding_polygon), msg
236           
237
238    def adjust_flow_for_available_water_at_inlet(self, Q, delta_t):
239        """Adjust Q downwards depending on available water at inlet
240
241           This is a critical step in modelling bridges and Culverts
242           the predicted flow through a structure based on an abstract
243           algorithm can at times request for water that is simply not
244           available due to any number of constrictions that limit the
245           flow approaching the structure In order to ensure that
246           there is adequate flow available certain checks are
247           required There needs to be a check using the Static Water
248           Volume sitting infront of the structure, In addition if the
249           water is moving the available water will be larger than the
250           static volume
251           
252           NOTE To temporarily switch this off for Debugging purposes
253           rem out line in function def compute_rates below
254        """
255   
256        if delta_t < epsilon:
257            # No need to adjust if time step is very small or zero
258            # In this case the possible flow will be very large
259            # anyway.
260            return Q
261       
262        # Short hands
263        domain = self.domain       
264        dq = domain.quantities               
265        time = domain.get_time()
266        I = self.inlet
267        idx = I.exchange_indices   
268
269        # Find triangle with the smallest depth
270        stage = dq['stage'].get_values(location='centroids', 
271                                               indices=[idx])
272        elevation = dq['elevation'].get_values(location='centroids', 
273                                               indices=[idx])       
274        depth = stage-elevation
275        min_depth = min(depth.flat)  # This may lead to errors if edge of area is at a higher level !!!!
276        avg_depth = mean(depth.flat) # Yes, but this one violates the conservation unit tests
277
278
279
280        # FIXME (Ole): If you want these, use log.critical() and
281        # make the statements depend on verbose
282        #print I.depth
283        #print I.velocity
284        #print self.width
285
286        # max_Q Based on Volume Calcs
287
288
289        depth_term = min_depth*I.exchange_area/delta_t
290        if min_depth < 0.2:
291            # Only add velocity term in shallow waters (< 20 cm)
292            # This is a little ad hoc, but maybe it is reasonable
293            velocity_term = self.width*min_depth*I.velocity
294        else:
295            velocity_term = 0.0
296
297        # This one takes approaching water into account   
298        max_Q = max(velocity_term, depth_term)
299
300        # This one preserves Volume
301        #max_Q = depth_term
302
303
304        if self.verbose is True:
305            log.critical('Max_Q = %f' % max_Q)           
306            msg = 'Width = %.2fm, Depth at inlet = %.2f m, Velocity = %.2f m/s.      ' % (self.width, I.depth, I.velocity)
307            msg += 'Max Q = %.2f m^3/s' %(max_Q)
308            log.critical(msg)
309
310        if self.log_filename is not None:               
311            log_to_file(self.log_filename, msg)
312        # New Procedure for assessing the flow available to the Culvert
313        # This routine uses the GET FLOW THROUGH CROSS SECTION
314        #   Need to check Several Polyline however
315        # Firstly 3 sides of the exchange Poly
316        # then only the Line Directly infront of the Polygon
317        # Access polygon Points from   self.inlet.polygon
318     
319        #  The Following computes the flow crossing over 3 sides of the exchange polygon for the structure
320        # Clearly the flow in the culvert can not be more than that flowing toward it through the exhange polygon
321       
322        #q1 = domain.get_flow_through_cross_section(self.culvert_polygons['exchange_polygon0'][1:3]) # First Side Segment
323        #q2 = domain.get_flow_through_cross_section(self.culvert_polygons['exchange_polygon0'][2:])   # Second Face Segment
324        #q3 =domain.get_flow_through_cross_section(self.culvert_polygons['exchange_polygon0'].take([3,0], axis=0)) # Third Side Segment
325        # q4 = domain.get_flow_through_cross_section([self.culvert_polygons['exchange_polygon0'][1:4]][0])
326        #q4=max(q1,0.0)+max(q2,0.0)+max(q3,0.0)
327        # To use only the Flow crossing the 3 sides of the Exchange Polygon use the following Line Only
328        #max_Q=max(q1,q2,q3,q4)
329        # Try Simple Smoothing using Average of 2 approaches
330        #max_Q=(max(q1,q2,q3,q4)+max_Q)/2.0
331        # Calculate the minimum in absolute terms of
332        # the requsted flow and the possible flow
333        Q_reduced = sign(Q)*min(abs(Q), abs(max_Q))
334        if self.verbose is True:
335            msg = 'Initial Q Reduced = %.2f m3/s.      ' % (Q_reduced)
336            log.critical(msg)
337
338        if self.log_filename is not None:               
339            log_to_file(self.log_filename, msg)
340        # Now Keep Rolling Average of Computed Discharge to Reduce / Remove Oscillations
341        #  can use delta_t if we want to averageover a time frame for example
342        # N = 5.0/delta_t  Will provide the average over 5 seconds
343
344        self.i=(self.i+1)%self.N
345        self.Q_list[self.i]=Q_reduced
346        Q_reduced = sum(self.Q_list)/len(self.Q_list)
347
348        if self.verbose is True:
349            msg = 'Final Q Reduced = %.2f m3/s.      ' % (Q_reduced)
350            log.critical(msg)
351
352        if self.log_filename is not None:               
353            log_to_file(self.log_filename, msg)
354
355
356        if abs(Q_reduced) < abs(Q): 
357            msg = '%.2fs: Requested flow is ' % time
358            msg += 'greater than what is supported by the smallest '
359            msg += 'depth at inlet exchange area:\n        '
360            msg += 'inlet exchange area: %.2f '% (I.exchange_area) 
361            msg += 'velocity at inlet :%.2f '% (I.velocity)
362            msg += 'Vel* Exch Area = : %.2f '% (I.velocity*avg_depth*self.width)
363            msg += 'h_min*inlet_area/delta_t = %.2f*%.2f/%.2f '\
364                % (avg_depth, I.exchange_area, delta_t)
365            msg += ' = %.2f m^3/s\n        ' % Q_reduced
366            msg += 'Q will be reduced from %.2f m^3/s to %.2f m^3/s.' % (Q, Q_reduced)
367            msg += 'Note calculate max_Q from V %.2f m^3/s ' % (max_Q)
368            if self.verbose is True:
369                log.critical(msg)
370               
371            if self.log_filename is not None:               
372                log_to_file(self.log_filename, msg)
373       
374        return Q_reduced   
375                       
376           
377    def compute_rates(self, delta_t):
378        """Compute new rates for inlet and outlet
379        """
380
381        # Short hands
382        domain = self.domain       
383        dq = domain.quantities               
384       
385        # Time stuff
386        time = domain.get_time()
387        self.last_update = time
388
389           
390        if hasattr(self, 'log_filename'):
391            log_filename = self.log_filename
392           
393        # Compute stage, energy and velocity at the
394        # enquiry points at each end of the culvert
395        openings = self.openings
396        for i, opening in enumerate(openings):
397            idx = self.enquiry_indices[i]               
398           
399            stage = dq['stage'].get_values(location='centroids',
400                                           indices=[idx])[0]
401            depth = h = stage-opening.elevation
402                                                           
403           
404            # Get velocity                                 
405            xmomentum = dq['xmomentum'].get_values(location='centroids',
406                                                   indices=[idx])[0]
407            ymomentum = dq['xmomentum'].get_values(location='centroids',
408                                                   indices=[idx])[0]
409
410            if h > minimum_allowed_height:
411                u = xmomentum/(h + velocity_protection/h)
412                v = ymomentum/(h + velocity_protection/h)
413            else:
414                u = v = 0.0
415               
416            v_squared = u*u + v*v
417           
418            if self.use_velocity_head is True:
419                velocity_head = 0.5*v_squared/g   
420            else:
421                velocity_head = 0.0
422           
423            opening.total_energy = velocity_head + stage
424            opening.specific_energy = velocity_head + depth
425            opening.stage = stage
426            opening.depth = depth
427            opening.velocity = sqrt(v_squared)
428           
429
430        # We now need to deal with each opening individually
431        # Determine flow direction based on total energy difference
432        delta_total_energy = openings[0].total_energy - openings[1].total_energy
433        if delta_total_energy > 0:
434            inlet = openings[0]
435            outlet = openings[1]
436
437            # FIXME: I think this whole momentum jet thing could be a bit more elegant
438            inlet.momentum = self.opening_momentum[0]
439            outlet.momentum = self.opening_momentum[1]
440        else:
441            inlet = openings[1]
442            outlet = openings[0]
443           
444            inlet.momentum = self.opening_momentum[1]
445            outlet.momentum = self.opening_momentum[0]
446
447            delta_total_energy = -delta_total_energy
448
449        self.inlet = inlet
450        self.outlet = outlet
451           
452        msg = 'Total energy difference is negative'
453        assert delta_total_energy >= 0.0, msg
454
455        # Recompute slope and issue warning if flow is uphill
456        # These values do not enter the computation
457        delta_z = inlet.elevation - outlet.elevation
458        culvert_slope = (delta_z/self.length)
459        if culvert_slope < 0.0:
460            # Adverse gradient - flow is running uphill
461            # Flow will be purely controlled by uphill outlet face
462            if self.verbose is True:
463                log.critical('%.2fs - WARNING: Flow is running uphill.' % time)
464           
465        if self.log_filename is not None:
466            s = 'Time=%.2f, inlet stage = %.2f, outlet stage = %.2f'\
467                %(time, self.inlet.stage, self.outlet.stage)
468            log_to_file(self.log_filename, s)
469            s = 'Delta total energy = %.3f' %(delta_total_energy)
470            log_to_file(log_filename, s)
471
472           
473        # Determine controlling energy (driving head) for culvert
474        if inlet.specific_energy > delta_total_energy:
475            # Outlet control
476            driving_head = delta_total_energy
477        else:
478            # Inlet control
479            driving_head = inlet.specific_energy
480           
481
482
483        if self.inlet.depth <= self.trigger_depth:
484            Q = 0.0
485        else:
486            # Calculate discharge for one barrel and
487            # set inlet.rate and outlet.rate
488           
489            if self.culvert_description_filename is not None:
490                try:
491                    Q = interpolate_linearly(driving_head, 
492                                             self.rating_curve[:,0], 
493                                             self.rating_curve[:,1]) 
494                except Below_interval, e:
495                    Q = self.rating_curve[0,1]             
496                    msg = '%.2fs: ' % time
497                    msg += 'Delta head smaller than rating curve minimum: '
498                    msg += str(e)
499                    msg += '\n        '
500                    msg += 'I will use minimum discharge %.2f m^3/s ' % Q
501                    msg += 'for culvert "%s"' % self.label
502                   
503                    if hasattr(self, 'log_filename'):                   
504                        log_to_file(self.log_filename, msg)
505                except Above_interval, e:
506                    Q = self.rating_curve[-1,1]         
507                    msg = '%.2fs: ' % time                 
508                    msg += 'Delta head greater than rating curve maximum: '
509                    msg += str(e)
510                    msg += '\n        '
511                    msg += 'I will use maximum discharge %.2f m^3/s ' % Q
512                    msg += 'for culvert "%s"' % self.label
513                   
514                    if self.log_filename is not None:                   
515                        log_to_file(self.log_filename, msg)
516            else:
517                # User culvert routine
518                Q, barrel_velocity, culvert_outlet_depth =\
519                    self.culvert_routine(inlet.depth,
520                                         outlet.depth,
521                                         inlet.velocity,
522                                         outlet.velocity,
523                                         inlet.specific_energy, 
524                                         delta_total_energy, 
525                                         g,
526                                         culvert_length=self.length,
527                                         culvert_width=self.width,
528                                         culvert_height=self.height,
529                                         culvert_type=self.culvert_type,
530                                         manning=self.manning,
531                                         sum_loss=self.sum_loss,
532                                         log_filename=self.log_filename)
533           
534           
535       
536        # Adjust discharge for multiple barrels
537        Q *= self.number_of_barrels
538
539        # Adjust discharge for available water at the inlet
540        Q = self.adjust_flow_for_available_water_at_inlet(Q, delta_t)
541       
542        self.inlet.rate = -Q
543        self.outlet.rate = Q
544
545
546        # Momentum jet stuff
547        if self.use_momentum_jet is True:
548
549
550            # Compute barrel momentum
551            barrel_momentum = barrel_velocity*culvert_outlet_depth
552
553            if self.log_filename is not None:                                   
554                s = 'Barrel velocity = %f' %barrel_velocity
555                log_to_file(self.log_filename, s)
556
557            # Compute momentum vector at outlet
558            outlet_mom_x, outlet_mom_y = self.vector * barrel_momentum
559               
560            if self.log_filename is not None:               
561                s = 'Directional momentum = (%f, %f)' %(outlet_mom_x, outlet_mom_y)
562                log_to_file(self.log_filename, s)
563
564
565            # Update momentum       
566            if delta_t > 0.0:
567                xmomentum_rate = outlet_mom_x - outlet.momentum[0].value
568                xmomentum_rate /= delta_t
569                   
570                ymomentum_rate = outlet_mom_y - outlet.momentum[1].value
571                ymomentum_rate /= delta_t
572                       
573                if self.log_filename is not None:               
574                    s = 'X Y MOM_RATE = (%f, %f) ' %(xmomentum_rate, ymomentum_rate)
575                    log_to_file(self.log_filename, s)                   
576            else:
577                xmomentum_rate = ymomentum_rate = 0.0
578
579
580            # Set momentum rates for outlet jet
581            outlet.momentum[0].rate = xmomentum_rate
582            outlet.momentum[1].rate = ymomentum_rate
583
584            # Remember this value for next step (IMPORTANT)
585            outlet.momentum[0].value = outlet_mom_x
586            outlet.momentum[1].value = outlet_mom_y                   
587
588            if int(domain.time*100) % 100 == 0:
589
590                if self.log_filename is not None:               
591                    s = 'T=%.5f, Culvert Discharge = %.3f f'\
592                        %(time, Q)
593                    s += ' Depth= %0.3f  Momentum = (%0.3f, %0.3f)'\
594                        %(culvert_outlet_depth, outlet_mom_x,outlet_mom_y)
595                    s += ' Momentum rate: (%.4f, %.4f)'\
596                        %(xmomentum_rate, ymomentum_rate)                   
597                    s+='Outlet Vel= %.3f'\
598                        %(barrel_velocity)
599                    log_to_file(self.log_filename, s)
600
601
602            # Execute momentum terms
603            # This is where Inflow objects are evaluated and update the domain
604                self.outlet.momentum[0](domain)
605                self.outlet.momentum[1](domain)       
606           
607
608           
609        # Log timeseries to file
610        try:
611            fid = open(self.timeseries_filename, 'a')       
612        except:
613            pass
614        else:   
615            fid.write('%.2f, %.2f\n' %(time, Q))
616            fid.close()
617
618        # Store value of time
619        self.last_time = time
620
621
622# FIXME(Ole): Write in C and reuse this function by similar code
623# in interpolate.py
624def interpolate_linearly(x, xvec, yvec):
625
626    msg = 'Input to function interpolate_linearly could not be converted '
627    msg += 'to numerical scalar: x = %s' % str(x)
628    try:
629        x = float(x)
630    except:
631        raise Exception, msg
632
633
634    # Check bounds
635    if x < xvec[0]:
636        msg = 'Value provided = %.2f, interpolation minimum = %.2f.'\
637            % (x, xvec[0])
638        raise Below_interval, msg
639
640    if x > xvec[-1]:
641        msg = 'Value provided = %.2f, interpolation maximum = %.2f.'\
642            %(x, xvec[-1])
643        raise Above_interval, msg
644
645
646    # Find appropriate slot within bounds
647    i = 0
648    while x > xvec[i]: i += 1
649
650
651    x0 = xvec[i-1]
652    x1 = xvec[i]
653    alpha = (x - x0)/(x1 - x0)
654
655    y0 = yvec[i-1]
656    y1 = yvec[i]
657    y = alpha*y1 + (1-alpha)*y0
658
659    return y
660
661
662
663def read_culvert_description(culvert_description_filename):
664
665    # Read description file
666    fid = open(culvert_description_filename)
667
668    read_rating_curve_data = False
669    rating_curve = []
670    for i, line in enumerate(fid.readlines()):
671
672        if read_rating_curve_data is True:
673            fields = line.split(',')
674            head_difference = float(fields[0].strip())
675            flow_rate = float(fields[1].strip())
676            barrel_velocity = float(fields[2].strip())
677
678            rating_curve.append([head_difference, flow_rate, barrel_velocity])
679
680        if i == 0:
681            # Header
682            continue
683        if i == 1:
684            # Metadata
685            fields = line.split(',')
686            label=fields[0].strip()
687            type=fields[1].strip().lower()
688            assert type in ['box', 'pipe']
689
690            width=float(fields[2].strip())
691            height=float(fields[3].strip())
692            length=float(fields[4].strip())
693            number_of_barrels=int(fields[5].strip())
694            #fields[6] refers to losses
695            description=fields[7].strip()
696
697        if line.strip() == '': continue # Skip blanks
698
699        if line.startswith('Rating'):
700            read_rating_curve_data = True
701            # Flow data follows
702
703    fid.close()
704
705    return label, type, width, height, length, number_of_barrels, description, rating_curve
706
Note: See TracBrowser for help on using the repository browser.