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
RevLine 
[7976]1from anuga.geometry.polygon import inside_polygon, polygon_area
2from anuga.config import g
[7939]3import anuga.utilities.log as log
4
[7980]5from boyd_box_culvert import Boyd_box_culvert
6
[7977]7class Culvert_operator:
[7958]8    """Culvert flow - transfer water from one rectangular box to another.
[7955]9    Sets up the geometry of problem
[7939]10   
[7955]11    This is the base class for culverts. Inherit from this class (and overwrite
12    compute_discharge method for specific subclasses)
[7939]13   
14    Input: Two points, pipe_size (either diameter or width, height),
15    mannings_rougness,
[7985]16    """ 
[7939]17
18    def __init__(self,
19                 domain,
[7978]20                 end_point0, 
21                 end_point1,
22                 width,
[7939]23                 height=None,
[7984]24                 apron=None,
[7985]25                 manning=0.013,
[7984]26                 enquiry_gap=0.2,
[7939]27                 verbose=False):
28       
[7958]29        self.domain = domain
[7962]30        self.domain.set_fractional_step_operator(self)
[7984]31        self.end_points = [end_point0, end_point1]
[7958]32       
[7955]33        if height is None:
34            height = width
[7939]35
[7984]36        if apron is None:
37            apron = width
38
39        self.width  = width
[7939]40        self.height = height
[7984]41        self.apron  = apron
[7985]42        self.manning = manning
[7984]43        self.enquiry_gap = enquiry_gap
44        self.verbose = verbose
[7955]45       
[7984]46        self.culvert = Boyd_box_culvert(self.domain,
47                                        self.end_points,
48                                        self.width,
49                                        self.height,
50                                        self.apron,
[7985]51                                        self.manning,
[7984]52                                        self.enquiry_gap,
53                                        self.verbose)
54       
[7980]55        self.routine = self.culvert.routine
[7968]56
[7984]57        self.inlets  = self.culvert.get_inlets()
[7968]58
[7984]59        if self.verbose:
60            self.print_stats()
61
62
[7962]63    def __call__(self):
64
65        timestep = self.domain.get_timestep()
[7978]66       
[7981]67        Q, barrel_speed, outlet_depth = self.routine()
[7962]68
[7984]69
[7980]70        inflow  = self.routine.get_inflow()
71        outflow = self.routine.get_outflow()
[7968]72
[7975]73
[7981]74        old_inflow_height = inflow.get_average_height()
[7985]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
[7980]82
[7985]83        factor = 1.0/(1.0 + Qstar*timestep/inflow.get_area())
[7980]84
[7985]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           
[7980]91
[7981]92        inflow.set_heights(new_inflow_height)
[7984]93
94        #inflow.set_xmoms(Q/inflow.get_area())
95        #inflow.set_ymoms(0.0)
96
97
[7981]98        inflow.set_xmoms(new_inflow_xmom)
99        inflow.set_ymoms(new_inflow_ymom)
[7980]100
[7984]101
102        loss = (old_inflow_height - new_inflow_height)*inflow.get_area()
103
[7985]104           
105        # set outflow
106        if old_inflow_height > 0.0 :
107                timestep_star = timestep*new_inflow_height/old_inflow_height
108        else:
[7984]109            timestep_star = 0.0
110
[7985]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           
[7980]116
[7985]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
[7980]121
[7984]122
[7985]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]
[7984]126
[7985]127            outflow.set_heights(new_outflow_height)
[7984]128
[7985]129            outflow.set_xmoms(barrel_speed*new_outflow_height*outflow_direction[0])
130            outflow.set_ymoms(barrel_speed*new_outflow_height*outflow_direction[1])
[7984]131
[7985]132            #outflow.set_xmoms(new_outflow_xmom)
133            #outflow.set_ymoms(new_outflow_ymom)
134           
135            #print '   outflow volume ',outflow.get_total_water_volume()
[7984]136
[7962]137    def print_stats(self):
138
139        print '====================================='
140        print 'Generic Culvert Operator'
141        print '====================================='
[7984]142
143        print 'Culvert'
144        print self.culvert
145
146        print 'Culvert Routine'
147        print self.routine
[7962]148       
[7968]149        for i, inlet in enumerate(self.inlets):
[7962]150            print '-------------------------------------'
[7968]151            print 'Inlet %i' % i
[7962]152            print '-------------------------------------'
153
[7968]154            print 'inlet triangle indices and centres'
[7980]155            print inlet.triangle_indices
156            print self.domain.get_centroid_coordinates()[inlet.triangle_indices]
[7962]157       
[7968]158            print 'polygon'
159            print inlet.polygon
[7962]160
161        print '====================================='
162
[7939]163                       
[7957]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.