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

Last change on this file since 7976 was 7976, checked in by habili, 13 years ago

Clean up of code. Removed references to enquiry point.

File size: 8.8 KB
RevLine 
[7976]1from anuga.geometry.polygon import inside_polygon, polygon_area
2from anuga.config import g
[7939]3import anuga.utilities.log as log
[7975]4import inlet
[7939]5import numpy as num
[7971]6import math
[7939]7
8class Below_interval(Exception): pass 
9class Above_interval(Exception): pass
10
11
[7955]12class Generic_box_culvert:
[7958]13    """Culvert flow - transfer water from one rectangular box to another.
[7955]14    Sets up the geometry of problem
[7939]15   
[7955]16    This is the base class for culverts. Inherit from this class (and overwrite
17    compute_discharge method for specific subclasses)
[7939]18   
19    Input: Two points, pipe_size (either diameter or width, height),
20    mannings_rougness,
21    """ 
22
23    def __init__(self,
24                 domain,
25                 end_point0=None, 
26                 end_point1=None,
27                 width=None,
28                 height=None,
29                 verbose=False):
30       
31        # Input check
32       
[7958]33        self.domain = domain
34
[7962]35        self.domain.set_fractional_step_operator(self)
36
[7971]37        self.end_points = [end_point0, end_point1]
[7958]38       
[7955]39        if height is None:
40            height = width
[7939]41
[7958]42        self.width = width
[7939]43        self.height = height
44       
[7955]45        self.verbose=verbose
46       
[7968]47        # Create the fundamental culvert polygons and create inlet objects
[7955]48        self.create_culvert_polygons()
[7939]49
[7968]50        #FIXME (SR) Put this into a foe loop to deal with more inlets
51        self.inlets = []
[7971]52        polygon0 = self.inlet_polygons[0]
[7968]53        inlet0_vector = self.culvert_vector
[7976]54        self.inlets.append(inlet.Inlet(self.domain, polygon0))
[7962]55
[7971]56        polygon1 = self.inlet_polygons[1]
[7968]57        inlet1_vector = - self.culvert_vector
[7976]58        self.inlets.append(inlet.Inlet(self.domain, polygon1))
[7971]59   
[7968]60        self.print_stats()
61
62
[7962]63    def __call__(self):
64
65        # Time stuff
66        time     = self.domain.get_time()
67        timestep = self.domain.get_timestep()
68
69
70
[7975]71        inflow  = self.inlets[0]
72        outflow = self.inlets[1]
[7962]73
[7975]74        # Determine flow direction based on total energy difference
75        delta_total_energy = inflow.get_average_total_energy() - outflow.get_average_total_energy()
[7968]76
[7975]77        if delta_total_energy < 0:
78            inflow  = self.inlets[1]
79            outflow = self.inlets[0]
80            delta_total_energy = -delta_total_energy
81
82        delta_z = inflow.get_average_elevation() - outflow.get_average_elevation()
83        culvert_slope = delta_z/self.culvert_length
84
85        # Determine controlling energy (driving head) for culvert
86        if inflow.get_average_specific_energy() > delta_total_energy:
87            # Outlet control
88            driving_head = delta_total_energy
89        else:
90            # Inlet control
91            driving_head = inflow.get_average_specific_energy()
92           
93
[7962]94        # Transfer
[7975]95        from culvert_routines import boyd_generalised_culvert_model
96        Q, barrel_velocity, culvert_outlet_depth =\
97                    boyd_generalised_culvert_model(inflow.get_average_height(),
98                                         outflow.get_average_height(),
99                                         inflow.get_average_speed(),
100                                         outflow.get_average_speed(),
101                                         inflow.get_average_specific_energy(),
102                                         delta_total_energy,
103                                         g,
104                                         culvert_length=self.culvert_length,
105                                         culvert_width=self.width,
106                                         culvert_height=self.height,
107                                         culvert_type='box',
108                                         manning=0.01)
[7962]109
[7975]110        transfer_water = Q*timestep
[7962]111
112
[7975]113        inflow.set_heights(inflow.get_average_height() - transfer_water)
114        inflow.set_xmoms(0.0)
115        inflow.set_ymoms(0.0)
116
117
118        outflow.set_heights(outflow.get_average_height() + transfer_water)
119        outflow.set_xmoms(0.0)
120        outflow.set_ymoms(0.0)
121
122
123
[7962]124    def print_stats(self):
125
126        print '====================================='
127        print 'Generic Culvert Operator'
128        print '====================================='
129       
[7968]130        for i, inlet in enumerate(self.inlets):
[7962]131            print '-------------------------------------'
[7968]132            print 'Inlet %i' % i
[7962]133            print '-------------------------------------'
134
[7968]135            print 'inlet triangle indices and centres'
136            print inlet.triangle_indices[i]
137            print self.domain.get_centroid_coordinates()[inlet.triangle_indices[i]]
[7962]138       
[7968]139            print 'polygon'
140            print inlet.polygon
[7962]141
142        print '====================================='
143
144
[7939]145
[7955]146
147
148    def create_culvert_polygons(self):
[7971]149
[7955]150        """Create polygons at the end of a culvert inlet and outlet.
151        At either end two polygons will be created; one for the actual flow to pass through and one a little further away
152        for enquiring the total energy at both ends of the culvert and transferring flow.
153        """
154
155        # Calculate geometry
[7957]156        x0, y0 = self.end_points[0]
157        x1, y1 = self.end_points[1]
[7955]158
[7958]159        dx = x1 - x0
160        dy = y1 - y0
[7955]161
162        self.culvert_vector = num.array([dx, dy])
[7971]163        self.culvert_length = math.sqrt(num.sum(self.culvert_vector**2))
[7958]164        assert self.culvert_length > 0.0, 'The length of culvert is less than 0'
[7955]165
166        # Unit direction vector and normal
167        self.culvert_vector /= self.culvert_length                      # Unit vector in culvert direction
168        self.culvert_normal = num.array([-dy, dx])/self.culvert_length  # Normal vector
169
170        # Short hands
[7958]171        w = 0.5*self.width*self.culvert_normal # Perpendicular vector of 1/2 width
172        h = self.height*self.culvert_vector    # Vector of length=height in the
[7955]173                             # direction of the culvert
174
[7971]175        self.inlet_polygons = []
[7955]176
[7957]177        # Build exchange polygon and enquiry points 0 and 1
[7958]178        for i in [0, 1]:
[7962]179            i0 = (2*i-1)
[7957]180            p0 = self.end_points[i] + w
[7958]181            p1 = self.end_points[i] - w
[7962]182            p2 = p1 + i0*h
183            p3 = p0 + i0*h
[7971]184            self.inlet_polygons.append(num.array([p0, p1, p2, p3]))
[7955]185
[7971]186        # Check that enquiry points are outside inlet polygons
[7957]187        for i in [0,1]:
[7971]188            polygon = self.inlet_polygons[i]
[7962]189            # FIXME (SR) Probably should calculate the area of all the triangles
190            # associated with this polygon, as there is likely to be some
191            # inconsistency between triangles and ploygon
[7957]192            area = polygon_area(polygon)
[7962]193           
[7955]194            msg = 'Polygon %s ' %(polygon)
195            msg += ' has area = %f' % area
196            assert area > 0.0, msg
[7971]197   
[7939]198
199                       
[7957]200# FIXME(Ole): Write in C and reuse this function by similar code
201# in interpolate.py
202def interpolate_linearly(x, xvec, yvec):
203
204    msg = 'Input to function interpolate_linearly could not be converted '
205    msg += 'to numerical scalar: x = %s' % str(x)
206    try:
207        x = float(x)
208    except:
209        raise Exception, msg
210
211
212    # Check bounds
213    if x < xvec[0]:
214        msg = 'Value provided = %.2f, interpolation minimum = %.2f.'\
215            % (x, xvec[0])
216        raise Below_interval, msg
217
218    if x > xvec[-1]:
219        msg = 'Value provided = %.2f, interpolation maximum = %.2f.'\
220            %(x, xvec[-1])
221        raise Above_interval, msg
222
223
224    # Find appropriate slot within bounds
225    i = 0
226    while x > xvec[i]: i += 1
227
228
229    x0 = xvec[i-1]
230    x1 = xvec[i]
231    alpha = (x - x0)/(x1 - x0)
232
233    y0 = yvec[i-1]
234    y1 = yvec[i]
235    y = alpha*y1 + (1-alpha)*y0
236
237    return y
238
239
240
241def read_culvert_description(culvert_description_filename):
242
243    # Read description file
244    fid = open(culvert_description_filename)
245
246    read_rating_curve_data = False
247    rating_curve = []
248    for i, line in enumerate(fid.readlines()):
249
250        if read_rating_curve_data is True:
251            fields = line.split(',')
252            head_difference = float(fields[0].strip())
253            flow_rate = float(fields[1].strip())
254            barrel_velocity = float(fields[2].strip())
255
256            rating_curve.append([head_difference, flow_rate, barrel_velocity])
257
258        if i == 0:
259            # Header
260            continue
261        if i == 1:
262            # Metadata
263            fields = line.split(',')
264            label=fields[0].strip()
265            type=fields[1].strip().lower()
266            assert type in ['box', 'pipe']
267
268            width=float(fields[2].strip())
269            height=float(fields[3].strip())
270            length=float(fields[4].strip())
271            number_of_barrels=int(fields[5].strip())
272            #fields[6] refers to losses
273            description=fields[7].strip()
274
275        if line.strip() == '': continue # Skip blanks
276
277        if line.startswith('Rating'):
278            read_rating_curve_data = True
279            # Flow data follows
280
281    fid.close()
282
283    return label, type, width, height, length, number_of_barrels, description, rating_curve
284
Note: See TracBrowser for help on using the repository browser.