source: trunk/anuga_core/source/anuga/structures/structure_operator.py @ 8008

Last change on this file since 8008 was 8008, checked in by habili, 12 years ago

Deleting unnecessary structure test scripts

File size: 9.2 KB
Line 
1import anuga
2import numpy as num
3import math
4import inlet
5
6class Structure_operator:
7    """Structure Operator - transfer water from one rectangular box to another.
8    Sets up the geometry of problem
9   
10    This is the base class for structures (culverts, pipes, bridges etc). Inherit from this class (and overwrite
11    discharge_routine method for specific subclasses)
12   
13    Input: Two points, pipe_size (either diameter or width, height),
14    mannings_rougness,
15    """ 
16
17    def __init__(self,
18                 domain,
19                 end_point0, 
20                 end_point1,
21                 width,
22                 height,
23                 apron,
24                 manning,
25                 enquiry_gap,
26                 description,
27                 verbose):
28       
29        self.domain = domain
30        self.domain.set_fractional_step_operator(self)
31        self.end_points = [end_point0, end_point1]
32       
33        if height is None:
34            height = width
35
36        if apron is None:
37            apron = width
38
39        self.width  = width
40        self.height = height
41        self.apron  = apron
42        self.manning = manning
43        self.enquiry_gap = enquiry_gap
44        self.description = description
45        self.verbose = verbose
46
47        self.discharge = 0.0
48        self.velocity = 0.0
49        self.delta_total_energy = 0.0
50        self.driving_energy = 0.0
51       
52        self.__create_exchange_polygons()
53
54        self.inlets = []
55        polygon0 = self.inlet_polygons[0]
56        enquiry_point0 = self.inlet_equiry_points[0]
57        outward_vector0 = self.culvert_vector
58        self.inlets.append(inlet.Inlet(self.domain, polygon0, enquiry_point0, outward_vector0))
59
60        polygon1 = self.inlet_polygons[1]
61        exchange_polygon1 = self.inlet_equiry_points[1]
62        outward_vector1  = - self.culvert_vector
63        self.inlets.append(inlet.Inlet(self.domain, polygon1, exchange_polygon1, outward_vector1))
64
65
66    def __call__(self):
67
68        timestep = self.domain.get_timestep()
69       
70        self.__determine_inflow_outflow()
71       
72        Q, barrel_speed, outlet_depth = self.discharge_routine()
73
74        old_inflow_height = self.inflow.get_average_height()
75        old_inflow_xmom = self.inflow.get_average_xmom()
76        old_inflow_ymom = self.inflow.get_average_ymom()
77           
78        if old_inflow_height > 0.0 :
79                Q_star = Q/old_inflow_height
80        else:
81                Q_star = 0.0
82
83        factor = 1.0/(1.0 + Q_star*timestep/self.inflow.get_area())
84
85        new_inflow_height = old_inflow_height*factor
86        new_inflow_xmom = old_inflow_xmom*factor
87        new_inflow_ymom = old_inflow_ymom*factor
88           
89        self.inflow.set_heights(new_inflow_height)
90
91        #inflow.set_xmoms(Q/inflow.get_area())
92        #inflow.set_ymoms(0.0)
93
94        self.inflow.set_xmoms(new_inflow_xmom)
95        self.inflow.set_ymoms(new_inflow_ymom)
96
97        loss = (old_inflow_height - new_inflow_height)*self.inflow.get_area()
98
99        # set outflow
100        if old_inflow_height > 0.0 :
101                timestep_star = timestep*new_inflow_height/old_inflow_height
102        else:
103            timestep_star = 0.0
104
105        outflow_extra_height = Q*timestep_star/self.outflow.get_area()
106        outflow_direction = - self.outflow.outward_culvert_vector
107        outflow_extra_momentum = outflow_extra_height*barrel_speed*outflow_direction
108           
109        gain = outflow_extra_height*self.outflow.get_area()
110           
111        #print Q, Q*timestep, barrel_speed, outlet_depth, Qstar, factor, timestep_star
112        #print '  ', loss, gain
113
114        # Stats
115        self.discharge  = Q#outflow_extra_height*self.outflow.get_area()/timestep
116        self.velocity = barrel_speed#self.discharge/outlet_depth/self.width
117
118        new_outflow_height = self.outflow.get_average_height() + outflow_extra_height
119
120        if self.use_momentum_jet :
121            # FIXME (SR) Review momentum to account for possible hydraulic jumps at outlet
122            #new_outflow_xmom = outflow.get_average_xmom() + outflow_extra_momentum[0]
123            #new_outflow_ymom = outflow.get_average_ymom() + outflow_extra_momentum[1]
124
125            new_outflow_xmom = barrel_speed*new_outflow_height*outflow_direction[0]
126            new_outflow_ymom = barrel_speed*new_outflow_height*outflow_direction[1]
127
128        else:
129            #new_outflow_xmom = outflow.get_average_xmom()
130            #new_outflow_ymom = outflow.get_average_ymom()
131
132            new_outflow_xmom = 0.0
133            new_outflow_ymom = 0.0
134
135        self.outflow.set_heights(new_outflow_height)
136        self.outflow.set_xmoms(new_outflow_xmom)
137        self.outflow.set_ymoms(new_outflow_ymom)
138
139
140    def __determine_inflow_outflow(self):
141        # Determine flow direction based on total energy difference
142
143        if self.use_velocity_head:
144            self.delta_total_energy = self.inlets[0].get_enquiry_total_energy() - self.inlets[1].get_enquiry_total_energy()
145        else:
146            self.delta_total_energy = self.inlets[0].get_enquiry_stage() - self.inlets[1].get_enquiry_stage()
147
148
149        self.inflow  = self.inlets[0]
150        self.outflow = self.inlets[1]
151       
152
153        if self.delta_total_energy < 0:
154            self.inflow  = self.inlets[1]
155            self.outflow = self.inlets[0]
156            self.delta_total_energy = -self.delta_total_energy
157
158
159    def __create_exchange_polygons(self):
160
161        """Create polygons at the end of a culvert inlet and outlet.
162        At either end two polygons will be created; one for the actual flow to pass through and one a little further away
163        for enquiring the total energy at both ends of the culvert and transferring flow.
164        """
165
166        # Calculate geometry
167        x0, y0 = self.end_points[0]
168        x1, y1 = self.end_points[1]
169
170        dx = x1 - x0
171        dy = y1 - y0
172
173        self.culvert_vector = num.array([dx, dy])
174        self.culvert_length = math.sqrt(num.sum(self.culvert_vector**2))
175        assert self.culvert_length > 0.0, 'The length of culvert is less than 0'
176
177        # Unit direction vector and normal
178        self.culvert_vector /= self.culvert_length                      # Unit vector in culvert direction
179        self.culvert_normal = num.array([-dy, dx])/self.culvert_length  # Normal vector
180
181        # Short hands
182        w = 0.5*self.width*self.culvert_normal # Perpendicular vector of 1/2 width
183        h = self.apron*self.culvert_vector    # Vector of length=height in the
184                             # direction of the culvert
185
186        gap = (1 + self.enquiry_gap)*h
187
188        self.inlet_polygons = []
189        self.inlet_equiry_points = []
190
191        # Build exchange polygon and enquiry point
192        for i in [0, 1]:
193            i0 = (2*i-1)
194            p0 = self.end_points[i] + w
195            p1 = self.end_points[i] - w
196            p2 = p1 + i0*h
197            p3 = p0 + i0*h
198            ep = self.end_points[i] + i0*gap
199
200            self.inlet_polygons.append(num.array([p0, p1, p2, p3]))
201            self.inlet_equiry_points.append(ep)
202
203        # Check that enquiry points are outside inlet polygons
204        for i in [0,1]:
205            polygon = self.inlet_polygons[i]
206            ep = self.inlet_equiry_points[i]
207           
208            area = anuga.polygon_area(polygon)
209           
210            msg = 'Polygon %s ' %(polygon)
211            msg += ' has area = %f' % area
212            assert area > 0.0, msg
213
214            msg = 'Enquiry point falls inside an exchange polygon.'
215            assert not anuga.inside_polygon(ep, polygon), msg
216           
217   
218    def discharge_routine(self):
219       
220        pass
221           
222
223    def print_stats(self):
224
225        print '====================================='
226        print 'Generic Culvert Operator'
227        print '====================================='
228
229        print 'Culvert'
230        print self.culvert
231
232        print 'Culvert Routine'
233        print self.routine
234       
235        for i, inlet in enumerate(self.inlets):
236            print '-------------------------------------'
237            print 'Inlet %i' % i
238            print '-------------------------------------'
239
240            print 'inlet triangle indices and centres'
241            print inlet.triangle_indices
242            print self.domain.get_centroid_coordinates()[inlet.triangle_indices]
243       
244            print 'polygon'
245            print inlet.polygon
246
247        print '====================================='
248
249
250    def structure_statistics(self):
251
252        message = '---------------------------\n'
253        message += 'Structure report for structure %s:\n' % self.description
254        message += '--------------------------\n'
255        message += 'Discharge [m^3/s]: %.2f\n' % self.discharge
256        message += 'Velocity  [m/s]: %.2f\n' % self.velocity
257        message += 'Inlet Driving Energy %.2f\n' % self.driving_energy
258        message += 'delta total energy %.2f\n' % self.delta_total_energy
259
260        return message
261
262
263    def get_inlets(self):
264       
265        return self.inlets
266       
267       
268    def get_culvert_length(self):
269       
270        return self.culvert_length
271       
272       
273    def get_culvert_width(self):
274       
275        return self.width
276       
277       
278    def get_culvert_diameter(self):
279   
280            return self.width
281       
282       
283    def get_culvert_height(self):
284   
285        return self.height
286
287
288    def get_culvert_apron(self):
289
290        return self.apron
Note: See TracBrowser for help on using the repository browser.