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

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

Refactoring of routines.

File size: 6.1 KB
Line 
1from anuga.geometry.polygon import inside_polygon, polygon_area
2from anuga.config import g
3import anuga.utilities.log as log
4import box_culvert
5
6class Culvert_operator:
7    """Culvert flow - transfer water from one rectangular box to another.
8    Sets up the geometry of problem
9   
10    This is the base class for culverts. Inherit from this class (and overwrite
11    compute_discharge method for specific subclasses)
12   
13    Input: Two points, pipe_size (either diameter or width, height),
14    mannings_rougness,
15    """ 
16
17    def __init__(self,
18                 domain,
19                 end_point0=None, 
20                 end_point1=None,
21                 width=None,
22                 height=None,
23                 verbose=False):
24       
25        self.domain = domain
26        self.domain.set_fractional_step_operator(self)
27        end_points = [end_point0, end_point1]
28       
29        if height is None:
30            height = width
31
32        self.width = width
33        self.height = height
34       
35        self.culvert = box_culvert.Box_culvert(self.domain, end_points, self.width, self.height)
36        self.inlets = self.culvert.get_inlets()
37   
38        self.print_stats()
39
40
41    def __call__(self):
42
43        timestep = self.domain.get_timestep()
44
45        inflow  = self.inlets[0]
46        outflow = self.inlets[1]
47
48        # Determine flow direction based on total energy difference
49        delta_total_energy = inflow.get_average_total_energy() - outflow.get_average_total_energy()
50
51        if delta_total_energy < 0:
52            inflow  = self.inlets[1]
53            outflow = self.inlets[0]
54            delta_total_energy = -delta_total_energy
55
56        delta_z = inflow.get_average_elevation() - outflow.get_average_elevation()
57        culvert_slope = delta_z/self.culvert.get_culvert_length()
58
59        # Determine controlling energy (driving head) for culvert
60        if inflow.get_average_specific_energy() > delta_total_energy:
61            # Outlet control
62            driving_head = delta_total_energy
63        else:
64            # Inlet control
65            driving_head = inflow.get_average_specific_energy()
66           
67        # Transfer
68        from culvert_routines import boyd_box, boyd_circle
69        Q, barrel_velocity, culvert_outlet_depth =\
70                              boyd_circle(inflow.get_average_height(),
71                                         outflow.get_average_height(),
72                                         inflow.get_average_speed(),
73                                         outflow.get_average_speed(),
74                                         inflow.get_average_specific_energy(),
75                                         delta_total_energy,
76                                         culvert_length=self.culvert.get_culvert_length(),
77                                         culvert_width=self.width,
78                                         culvert_height=self.height,
79                                         manning=0.01)
80
81        transfer_water = Q*timestep
82
83        inflow.set_heights(inflow.get_average_height() - transfer_water)
84        inflow.set_xmoms(0.0)
85        inflow.set_ymoms(0.0)
86
87        outflow.set_heights(outflow.get_average_height() + transfer_water)
88        outflow.set_xmoms(0.0)
89        outflow.set_ymoms(0.0)
90
91    def print_stats(self):
92
93        print '====================================='
94        print 'Generic Culvert Operator'
95        print '====================================='
96       
97        for i, inlet in enumerate(self.inlets):
98            print '-------------------------------------'
99            print 'Inlet %i' % i
100            print '-------------------------------------'
101
102            print 'inlet triangle indices and centres'
103            print inlet.triangle_indices[i]
104            print self.domain.get_centroid_coordinates()[inlet.triangle_indices[i]]
105       
106            print 'polygon'
107            print inlet.polygon
108
109        print '====================================='
110
111                       
112# FIXME(Ole): Write in C and reuse this function by similar code
113# in interpolate.py
114def interpolate_linearly(x, xvec, yvec):
115
116    msg = 'Input to function interpolate_linearly could not be converted '
117    msg += 'to numerical scalar: x = %s' % str(x)
118    try:
119        x = float(x)
120    except:
121        raise Exception, msg
122
123
124    # Check bounds
125    if x < xvec[0]:
126        msg = 'Value provided = %.2f, interpolation minimum = %.2f.'\
127            % (x, xvec[0])
128        raise Below_interval, msg
129
130    if x > xvec[-1]:
131        msg = 'Value provided = %.2f, interpolation maximum = %.2f.'\
132            %(x, xvec[-1])
133        raise Above_interval, msg
134
135
136    # Find appropriate slot within bounds
137    i = 0
138    while x > xvec[i]: i += 1
139
140
141    x0 = xvec[i-1]
142    x1 = xvec[i]
143    alpha = (x - x0)/(x1 - x0)
144
145    y0 = yvec[i-1]
146    y1 = yvec[i]
147    y = alpha*y1 + (1-alpha)*y0
148
149    return y
150
151
152
153def read_culvert_description(culvert_description_filename):
154
155    # Read description file
156    fid = open(culvert_description_filename)
157
158    read_rating_curve_data = False
159    rating_curve = []
160    for i, line in enumerate(fid.readlines()):
161
162        if read_rating_curve_data is True:
163            fields = line.split(',')
164            head_difference = float(fields[0].strip())
165            flow_rate = float(fields[1].strip())
166            barrel_velocity = float(fields[2].strip())
167
168            rating_curve.append([head_difference, flow_rate, barrel_velocity])
169
170        if i == 0:
171            # Header
172            continue
173        if i == 1:
174            # Metadata
175            fields = line.split(',')
176            label=fields[0].strip()
177            type=fields[1].strip().lower()
178            assert type in ['box', 'pipe']
179
180            width=float(fields[2].strip())
181            height=float(fields[3].strip())
182            length=float(fields[4].strip())
183            number_of_barrels=int(fields[5].strip())
184            #fields[6] refers to losses
185            description=fields[7].strip()
186
187        if line.strip() == '': continue # Skip blanks
188
189        if line.startswith('Rating'):
190            read_rating_curve_data = True
191            # Flow data follows
192
193    fid.close()
194
195    return label, type, width, height, length, number_of_barrels, description, rating_curve
196
Note: See TracBrowser for help on using the repository browser.