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

Last change on this file since 7989 was 7989, checked in by steve, 14 years ago

Hopefully merged Narimans mannings and Steve velocity and momentum switches

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