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

Last change on this file since 7981 was 7981, checked in by steve, 9 years ago

Implemented semi-implicit update of culvert operator

File size: 5.6 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, outlet_depth = self.routine()
50
51        inflow  = self.routine.get_inflow()
52        outflow = self.routine.get_outflow()
53
54
55        old_inflow_height = inflow.get_average_height()
56                old_inflow_xmom = inflow.get_average_xmom()
57                old_inflow_ymom = inflow.get_average_ymom()
58               
59                if old_inflow_height > 0.0 :
60                        Qstar = Q/old_inflow_height/inflow.get_area()
61                else:
62                        Qstar = 0.0
63
64                factor = 1.0/(1.0 + Qstar*timestep)
65
66               
67               
68                new_inflow_height = old_inflow_height*factor
69                new_inflow_xmom = old_inflow_xmom*factor
70                new_inflow_ymom = old_inflow_ymom*factor
71               
72
73        inflow.set_heights(new_inflow_height)
74        inflow.set_xmoms(new_inflow_xmom)
75        inflow.set_ymoms(new_inflow_ymom)
76
77               
78                # set outflow
79                if old_inflow_height > 0.0 :
80                        timestep_star = timestep*new_inflow_height/old_inflow_height
81                else:
82                        timestep_star = 0.0
83               
84                print Q, barrel_speed, outlet_depth, Qstar, factor, timestep_star
85               
86               
87        outflow_extra_height = Q*timestep_star/outflow.get_area()
88        outflow_direction = - outflow.outward_culvert_vector
89        outflow_extra_momentum = outflow_extra_height*barrel_speed*outflow_direction
90               
91
92        outflow.set_heights(outflow.get_average_height() + outflow_extra_height)
93        outflow.set_xmoms(outflow.get_average_xmom() + outflow_extra_momentum[0] )
94        outflow.set_ymoms(outflow.get_average_ymom() + outflow_extra_momentum[1] )
95
96    def print_stats(self):
97
98        print '====================================='
99        print 'Generic Culvert Operator'
100        print '====================================='
101       
102        for i, inlet in enumerate(self.inlets):
103            print '-------------------------------------'
104            print 'Inlet %i' % i
105            print '-------------------------------------'
106
107            print 'inlet triangle indices and centres'
108            print inlet.triangle_indices
109            print self.domain.get_centroid_coordinates()[inlet.triangle_indices]
110       
111            print 'polygon'
112            print inlet.polygon
113
114        print '====================================='
115
116                       
117# FIXME(Ole): Write in C and reuse this function by similar code
118# in interpolate.py
119def interpolate_linearly(x, xvec, yvec):
120
121    msg = 'Input to function interpolate_linearly could not be converted '
122    msg += 'to numerical scalar: x = %s' % str(x)
123    try:
124        x = float(x)
125    except:
126        raise Exception, msg
127
128
129    # Check bounds
130    if x < xvec[0]:
131        msg = 'Value provided = %.2f, interpolation minimum = %.2f.'\
132            % (x, xvec[0])
133        raise Below_interval, msg
134
135    if x > xvec[-1]:
136        msg = 'Value provided = %.2f, interpolation maximum = %.2f.'\
137            %(x, xvec[-1])
138        raise Above_interval, msg
139
140
141    # Find appropriate slot within bounds
142    i = 0
143    while x > xvec[i]: i += 1
144
145
146    x0 = xvec[i-1]
147    x1 = xvec[i]
148    alpha = (x - x0)/(x1 - x0)
149
150    y0 = yvec[i-1]
151    y1 = yvec[i]
152    y = alpha*y1 + (1-alpha)*y0
153
154    return y
155
156
157
158def read_culvert_description(culvert_description_filename):
159
160    # Read description file
161    fid = open(culvert_description_filename)
162
163    read_rating_curve_data = False
164    rating_curve = []
165    for i, line in enumerate(fid.readlines()):
166
167        if read_rating_curve_data is True:
168            fields = line.split(',')
169            head_difference = float(fields[0].strip())
170            flow_rate = float(fields[1].strip())
171            barrel_velocity = float(fields[2].strip())
172
173            rating_curve.append([head_difference, flow_rate, barrel_velocity])
174
175        if i == 0:
176            # Header
177            continue
178        if i == 1:
179            # Metadata
180            fields = line.split(',')
181            label=fields[0].strip()
182            type=fields[1].strip().lower()
183            assert type in ['box', 'pipe']
184
185            width=float(fields[2].strip())
186            height=float(fields[3].strip())
187            length=float(fields[4].strip())
188            number_of_barrels=int(fields[5].strip())
189            #fields[6] refers to losses
190            description=fields[7].strip()
191
192        if line.strip() == '': continue # Skip blanks
193
194        if line.startswith('Rating'):
195            read_rating_curve_data = True
196            # Flow data follows
197
198    fid.close()
199
200    return label, type, width, height, length, number_of_barrels, description, rating_curve
201
Note: See TracBrowser for help on using the repository browser.