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

Last change on this file since 7984 was 7984, checked in by steve, 13 years ago

Added enquiry points to culverts

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