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

Last change on this file since 7985 was 7985, checked in by habili, 9 years ago

deals with manning

File size: 7.2 KB
Line 
1from anuga.geometry.polygon import inside_polygon, polygon_area
2from anuga.config import g
3import anuga.utilities.log as log
4
5from boyd_box_culvert import Boyd_box_culvert
6
7class Culvert_operator:
8    """Culvert flow - transfer water from one rectangular box to another.
9    Sets up the geometry of problem
10   
11    This is the base class for culverts. Inherit from this class (and overwrite
12    compute_discharge method for specific subclasses)
13   
14    Input: Two points, pipe_size (either diameter or width, height),
15    mannings_rougness,
16    """ 
17
18    def __init__(self,
19                 domain,
20                 end_point0, 
21                 end_point1,
22                 width,
23                 height=None,
24                 apron=None,
25                 manning=0.013,
26                 enquiry_gap=0.2,
27                 verbose=False):
28       
29        self.domain = domain
30        self.domain.set_fractional_step_operator(self)
31        self.end_points = [end_point0, end_point1]
32       
33        if height is None:
34            height = width
35
36        if apron is None:
37            apron = width
38
39        self.width  = width
40        self.height = height
41        self.apron  = apron
42        self.manning = manning
43        self.enquiry_gap = enquiry_gap
44        self.verbose = verbose
45       
46        self.culvert = Boyd_box_culvert(self.domain,
47                                        self.end_points,
48                                        self.width,
49                                        self.height,
50                                        self.apron,
51                                        self.manning,
52                                        self.enquiry_gap,
53                                        self.verbose)
54       
55        self.routine = self.culvert.routine
56
57        self.inlets  = self.culvert.get_inlets()
58
59        if self.verbose:
60            self.print_stats()
61
62
63    def __call__(self):
64
65        timestep = self.domain.get_timestep()
66       
67        Q, barrel_speed, outlet_depth = self.routine()
68
69
70        inflow  = self.routine.get_inflow()
71        outflow = self.routine.get_outflow()
72
73
74        old_inflow_height = inflow.get_average_height()
75        old_inflow_xmom = inflow.get_average_xmom()
76        old_inflow_ymom = inflow.get_average_ymom()
77           
78        if old_inflow_height > 0.0 :
79                Qstar = Q/old_inflow_height
80        else:
81                Qstar = 0.0
82
83        factor = 1.0/(1.0 + Qstar*timestep/inflow.get_area())
84
85           
86           
87        new_inflow_height = old_inflow_height*factor
88        new_inflow_xmom = old_inflow_xmom*factor
89        new_inflow_ymom = old_inflow_ymom*factor
90           
91
92        inflow.set_heights(new_inflow_height)
93
94        #inflow.set_xmoms(Q/inflow.get_area())
95        #inflow.set_ymoms(0.0)
96
97
98        inflow.set_xmoms(new_inflow_xmom)
99        inflow.set_ymoms(new_inflow_ymom)
100
101
102        loss = (old_inflow_height - new_inflow_height)*inflow.get_area()
103
104           
105        # set outflow
106        if old_inflow_height > 0.0 :
107                timestep_star = timestep*new_inflow_height/old_inflow_height
108        else:
109            timestep_star = 0.0
110
111           
112            outflow_extra_height = Q*timestep_star/outflow.get_area()
113            outflow_direction = - outflow.outward_culvert_vector
114            outflow_extra_momentum = outflow_extra_height*barrel_speed*outflow_direction
115           
116
117            gain = outflow_extra_height*outflow.get_area()
118           
119            #print Q, Q*timestep, barrel_speed, outlet_depth, Qstar, factor, timestep_star
120            #print '  ', loss, gain
121
122
123            new_outflow_height = outflow.get_average_height() + outflow_extra_height
124            new_outflow_xmom = outflow.get_average_xmom() + outflow_extra_momentum[0]
125            new_outflow_ymom = outflow.get_average_ymom() + outflow_extra_momentum[1]
126
127            outflow.set_heights(new_outflow_height)
128
129            outflow.set_xmoms(barrel_speed*new_outflow_height*outflow_direction[0])
130            outflow.set_ymoms(barrel_speed*new_outflow_height*outflow_direction[1])
131
132            #outflow.set_xmoms(new_outflow_xmom)
133            #outflow.set_ymoms(new_outflow_ymom)
134           
135            #print '   outflow volume ',outflow.get_total_water_volume()
136
137    def print_stats(self):
138
139        print '====================================='
140        print 'Generic Culvert Operator'
141        print '====================================='
142
143        print 'Culvert'
144        print self.culvert
145
146        print 'Culvert Routine'
147        print self.routine
148       
149        for i, inlet in enumerate(self.inlets):
150            print '-------------------------------------'
151            print 'Inlet %i' % i
152            print '-------------------------------------'
153
154            print 'inlet triangle indices and centres'
155            print inlet.triangle_indices
156            print self.domain.get_centroid_coordinates()[inlet.triangle_indices]
157       
158            print 'polygon'
159            print inlet.polygon
160
161        print '====================================='
162
163                       
164# FIXME(Ole): Write in C and reuse this function by similar code
165# in interpolate.py
166def interpolate_linearly(x, xvec, yvec):
167
168    msg = 'Input to function interpolate_linearly could not be converted '
169    msg += 'to numerical scalar: x = %s' % str(x)
170    try:
171        x = float(x)
172    except:
173        raise Exception, msg
174
175
176    # Check bounds
177    if x < xvec[0]:
178        msg = 'Value provided = %.2f, interpolation minimum = %.2f.'\
179            % (x, xvec[0])
180        raise Below_interval, msg
181
182    if x > xvec[-1]:
183        msg = 'Value provided = %.2f, interpolation maximum = %.2f.'\
184            %(x, xvec[-1])
185        raise Above_interval, msg
186
187
188    # Find appropriate slot within bounds
189    i = 0
190    while x > xvec[i]: i += 1
191
192
193    x0 = xvec[i-1]
194    x1 = xvec[i]
195    alpha = (x - x0)/(x1 - x0)
196
197    y0 = yvec[i-1]
198    y1 = yvec[i]
199    y = alpha*y1 + (1-alpha)*y0
200
201    return y
202
203
204
205def read_culvert_description(culvert_description_filename):
206
207    # Read description file
208    fid = open(culvert_description_filename)
209
210    read_rating_curve_data = False
211    rating_curve = []
212    for i, line in enumerate(fid.readlines()):
213
214        if read_rating_curve_data is True:
215            fields = line.split(',')
216            head_difference = float(fields[0].strip())
217            flow_rate = float(fields[1].strip())
218            barrel_velocity = float(fields[2].strip())
219
220            rating_curve.append([head_difference, flow_rate, barrel_velocity])
221
222        if i == 0:
223            # Header
224            continue
225        if i == 1:
226            # Metadata
227            fields = line.split(',')
228            label=fields[0].strip()
229            type=fields[1].strip().lower()
230            assert type in ['box', 'pipe']
231
232            width=float(fields[2].strip())
233            height=float(fields[3].strip())
234            length=float(fields[4].strip())
235            number_of_barrels=int(fields[5].strip())
236            #fields[6] refers to losses
237            description=fields[7].strip()
238
239        if line.strip() == '': continue # Skip blanks
240
241        if line.startswith('Rating'):
242            read_rating_curve_data = True
243            # Flow data follows
244
245    fid.close()
246
247    return label, type, width, height, length, number_of_barrels, description, rating_curve
248
Note: See TracBrowser for help on using the repository browser.