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

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

Clean up of code. Removed references to enquiry point.

File size: 8.8 KB
Line 
1from anuga.geometry.polygon import inside_polygon, polygon_area
2from anuga.config import g
3import anuga.utilities.log as log
4import inlet
5import numpy as num
6import math
7
8class Below_interval(Exception): pass 
9class Above_interval(Exception): pass
10
11
12class Generic_box_culvert:
13    """Culvert flow - transfer water from one rectangular box to another.
14    Sets up the geometry of problem
15   
16    This is the base class for culverts. Inherit from this class (and overwrite
17    compute_discharge method for specific subclasses)
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       
33        self.domain = domain
34
35        self.domain.set_fractional_step_operator(self)
36
37        self.end_points = [end_point0, end_point1]
38       
39        if height is None:
40            height = width
41
42        self.width = width
43        self.height = height
44       
45        self.verbose=verbose
46       
47        # Create the fundamental culvert polygons and create inlet objects
48        self.create_culvert_polygons()
49
50        #FIXME (SR) Put this into a foe loop to deal with more inlets
51        self.inlets = []
52        polygon0 = self.inlet_polygons[0]
53        inlet0_vector = self.culvert_vector
54        self.inlets.append(inlet.Inlet(self.domain, polygon0))
55
56        polygon1 = self.inlet_polygons[1]
57        inlet1_vector = - self.culvert_vector
58        self.inlets.append(inlet.Inlet(self.domain, polygon1))
59   
60        self.print_stats()
61
62
63    def __call__(self):
64
65        # Time stuff
66        time     = self.domain.get_time()
67        timestep = self.domain.get_timestep()
68
69
70
71        inflow  = self.inlets[0]
72        outflow = self.inlets[1]
73
74        # Determine flow direction based on total energy difference
75        delta_total_energy = inflow.get_average_total_energy() - outflow.get_average_total_energy()
76
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
94        # Transfer
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)
109
110        transfer_water = Q*timestep
111
112
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
124    def print_stats(self):
125
126        print '====================================='
127        print 'Generic Culvert Operator'
128        print '====================================='
129       
130        for i, inlet in enumerate(self.inlets):
131            print '-------------------------------------'
132            print 'Inlet %i' % i
133            print '-------------------------------------'
134
135            print 'inlet triangle indices and centres'
136            print inlet.triangle_indices[i]
137            print self.domain.get_centroid_coordinates()[inlet.triangle_indices[i]]
138       
139            print 'polygon'
140            print inlet.polygon
141
142        print '====================================='
143
144
145
146
147
148    def create_culvert_polygons(self):
149
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
156        x0, y0 = self.end_points[0]
157        x1, y1 = self.end_points[1]
158
159        dx = x1 - x0
160        dy = y1 - y0
161
162        self.culvert_vector = num.array([dx, dy])
163        self.culvert_length = math.sqrt(num.sum(self.culvert_vector**2))
164        assert self.culvert_length > 0.0, 'The length of culvert is less than 0'
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
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
173                             # direction of the culvert
174
175        self.inlet_polygons = []
176
177        # Build exchange polygon and enquiry points 0 and 1
178        for i in [0, 1]:
179            i0 = (2*i-1)
180            p0 = self.end_points[i] + w
181            p1 = self.end_points[i] - w
182            p2 = p1 + i0*h
183            p3 = p0 + i0*h
184            self.inlet_polygons.append(num.array([p0, p1, p2, p3]))
185
186        # Check that enquiry points are outside inlet polygons
187        for i in [0,1]:
188            polygon = self.inlet_polygons[i]
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
192            area = polygon_area(polygon)
193           
194            msg = 'Polygon %s ' %(polygon)
195            msg += ' has area = %f' % area
196            assert area > 0.0, msg
197   
198
199                       
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.