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

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

Added new tests

File size: 5.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                 verbose=False):
25       
26        self.domain = domain
27        self.domain.set_fractional_step_operator(self)
28        end_points = [end_point0, end_point1]
29       
30        if height is None:
31            height = width
32
33        self.width = width
34        self.height = height
35       
36        self.culvert = Boyd_box_culvert(self.domain, end_points, self.width, self.height)
37        print self.culvert
38        self.routine = self.culvert.routine
39        print self.routine
40        self.inlets = self.culvert.get_inlets()
41   
42        self.print_stats()
43
44
45    def __call__(self):
46
47        timestep = self.domain.get_timestep()
48       
49        Q, barrel_speed, culvert_outlet_depth = self.routine()
50
51        inflow  = self.routine.get_inflow()
52        outflow = self.routine.get_outflow()
53
54        outflow_direction = - outflow.outward_culvert_vector
55
56        outflow_momentum_flux = barrel_speed**2*culvert_outlet_depth*outflow_direction
57
58
59        print Q, barrel_speed, culvert_outlet_depth, outflow_momentum_flux
60
61        #FIXME (SR) Check whether we need to mult/divide by inlet area
62        inflow_transfer =  Q*timestep/inflow.get_area()
63
64        outflow_transfer = Q*timestep/outflow.get_area()
65
66
67
68        inflow.set_heights(inflow.get_average_height() - inflow_transfer)
69
70        inflow.set_xmoms(0.0)
71        inflow.set_ymoms(0.0)
72
73        #u = outflow.get_xvelocities()
74        #v = outflow.get_yvelocities()
75
76        outflow.set_heights(outflow.get_average_height() + outflow_transfer)
77        #outflow.set_xmoms(outflow.get_xmoms() + timestep*outflow_momentum_flux[0] )
78        #outflow.set_ymoms(outflow.get_ymoms() + timestep*outflow_momentum_flux[1] )
79
80    def print_stats(self):
81
82        print '====================================='
83        print 'Generic Culvert Operator'
84        print '====================================='
85       
86        for i, inlet in enumerate(self.inlets):
87            print '-------------------------------------'
88            print 'Inlet %i' % i
89            print '-------------------------------------'
90
91            print 'inlet triangle indices and centres'
92            print inlet.triangle_indices
93            print self.domain.get_centroid_coordinates()[inlet.triangle_indices]
94       
95            print 'polygon'
96            print inlet.polygon
97
98        print '====================================='
99
100                       
101# FIXME(Ole): Write in C and reuse this function by similar code
102# in interpolate.py
103def interpolate_linearly(x, xvec, yvec):
104
105    msg = 'Input to function interpolate_linearly could not be converted '
106    msg += 'to numerical scalar: x = %s' % str(x)
107    try:
108        x = float(x)
109    except:
110        raise Exception, msg
111
112
113    # Check bounds
114    if x < xvec[0]:
115        msg = 'Value provided = %.2f, interpolation minimum = %.2f.'\
116            % (x, xvec[0])
117        raise Below_interval, msg
118
119    if x > xvec[-1]:
120        msg = 'Value provided = %.2f, interpolation maximum = %.2f.'\
121            %(x, xvec[-1])
122        raise Above_interval, msg
123
124
125    # Find appropriate slot within bounds
126    i = 0
127    while x > xvec[i]: i += 1
128
129
130    x0 = xvec[i-1]
131    x1 = xvec[i]
132    alpha = (x - x0)/(x1 - x0)
133
134    y0 = yvec[i-1]
135    y1 = yvec[i]
136    y = alpha*y1 + (1-alpha)*y0
137
138    return y
139
140
141
142def read_culvert_description(culvert_description_filename):
143
144    # Read description file
145    fid = open(culvert_description_filename)
146
147    read_rating_curve_data = False
148    rating_curve = []
149    for i, line in enumerate(fid.readlines()):
150
151        if read_rating_curve_data is True:
152            fields = line.split(',')
153            head_difference = float(fields[0].strip())
154            flow_rate = float(fields[1].strip())
155            barrel_velocity = float(fields[2].strip())
156
157            rating_curve.append([head_difference, flow_rate, barrel_velocity])
158
159        if i == 0:
160            # Header
161            continue
162        if i == 1:
163            # Metadata
164            fields = line.split(',')
165            label=fields[0].strip()
166            type=fields[1].strip().lower()
167            assert type in ['box', 'pipe']
168
169            width=float(fields[2].strip())
170            height=float(fields[3].strip())
171            length=float(fields[4].strip())
172            number_of_barrels=int(fields[5].strip())
173            #fields[6] refers to losses
174            description=fields[7].strip()
175
176        if line.strip() == '': continue # Skip blanks
177
178        if line.startswith('Rating'):
179            read_rating_curve_data = True
180            # Flow data follows
181
182    fid.close()
183
184    return label, type, width, height, length, number_of_barrels, description, rating_curve
185
Note: See TracBrowser for help on using the repository browser.