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

Last change on this file since 8027 was 8027, checked in by habili, 14 years ago

Added control case text to print_timestepping_statistics

File size: 11.5 KB
Line 
1import anuga
2import numpy as num
3import math
4import inlet
5
6from anuga.utilities.system_tools import log_to_file
7
8
9class Structure_operator:
10    """Structure Operator - transfer water from one rectangular box to another.
11    Sets up the geometry of problem
12   
13    This is the base class for structures (culverts, pipes, bridges etc). Inherit from this class (and overwrite
14    discharge_routine method for specific subclasses)
15   
16    Input: Two points, pipe_size (either diameter or width, height),
17    mannings_rougness,
18    """ 
19
20    counter = 0
21
22    def __init__(self,
23                 domain,
24                 end_point0, 
25                 end_point1,
26                 width,
27                 height,
28                 apron,
29                 manning,
30                 enquiry_gap,
31                 description,
32                 label,
33                 structure_type,
34                 logging,
35                 verbose):
36       
37        self.domain = domain
38        self.domain.set_fractional_step_operator(self)
39        self.end_points = [end_point0, end_point1]
40
41       
42       
43        if height is None:
44            height = width
45
46        if apron is None:
47            apron = width
48
49        self.width  = width
50        self.height = height
51        self.apron  = apron
52        self.manning = manning
53        self.enquiry_gap = enquiry_gap
54
55        if description == None:
56            self.description = ' '
57        else:
58            self.description = description
59       
60
61        if label == None:
62            self.label = "structure_%g" % Structure_operator.counter
63        else:
64            self.label = label + '_%g' % Structure_operator.counter
65
66
67        if structure_type == None:
68            self.structure_type = 'generic structure'
69        else:
70            self.structure_type = structure_type
71           
72        self.verbose = verbose
73
74       
75       
76        # Keep count of structures
77        Structure_operator.counter += 1
78
79        # Slots for recording current statistics
80        self.discharge = 0.0
81        self.velocity = 0.0
82        self.delta_total_energy = 0.0
83        self.driving_energy = 0.0
84       
85        self.__create_exchange_polygons()
86
87        self.inlets = []
88        polygon0 = self.inlet_polygons[0]
89        enquiry_point0 = self.inlet_equiry_points[0]
90        outward_vector0 = self.culvert_vector
91        self.inlets.append(inlet.Inlet(self.domain, polygon0, enquiry_point0, outward_vector0))
92
93        polygon1 = self.inlet_polygons[1]
94        exchange_polygon1 = self.inlet_equiry_points[1]
95        outward_vector1  = - self.culvert_vector
96        self.inlets.append(inlet.Inlet(self.domain, polygon1, exchange_polygon1, outward_vector1))
97
98        self.set_logging(logging)
99
100    def __call__(self):
101
102        timestep = self.domain.get_timestep()
103       
104        self.__determine_inflow_outflow()
105       
106        Q, barrel_speed, outlet_depth = self.discharge_routine()
107
108        old_inflow_height = self.inflow.get_average_height()
109        old_inflow_xmom = self.inflow.get_average_xmom()
110        old_inflow_ymom = self.inflow.get_average_ymom()
111
112        # Implement the update of flow over a timestep by
113        # using a semi-implict update. This ensures that
114        # the update does not create a negative height
115        if old_inflow_height > 0.0 :
116                Q_star = Q/old_inflow_height
117        else:
118                Q_star = 0.0
119
120        factor = 1.0/(1.0 + Q_star*timestep/self.inflow.get_area())
121
122        new_inflow_height = old_inflow_height*factor
123        new_inflow_xmom = old_inflow_xmom*factor
124        new_inflow_ymom = old_inflow_ymom*factor
125           
126        self.inflow.set_heights(new_inflow_height)
127
128        #inflow.set_xmoms(Q/inflow.get_area())
129        #inflow.set_ymoms(0.0)
130
131        self.inflow.set_xmoms(new_inflow_xmom)
132        self.inflow.set_ymoms(new_inflow_ymom)
133
134        loss = (old_inflow_height - new_inflow_height)*self.inflow.get_area()
135
136        # set outflow
137        if old_inflow_height > 0.0 :
138                timestep_star = timestep*new_inflow_height/old_inflow_height
139        else:
140            timestep_star = 0.0
141
142        outflow_extra_height = Q*timestep_star/self.outflow.get_area()
143        outflow_direction = - self.outflow.outward_culvert_vector
144        outflow_extra_momentum = outflow_extra_height*barrel_speed*outflow_direction
145           
146        gain = outflow_extra_height*self.outflow.get_area()
147           
148        #print Q, Q*timestep, barrel_speed, outlet_depth, Qstar, factor, timestep_star
149        #print '  ', loss, gain
150
151        # Stats
152        self.discharge  = Q#outflow_extra_height*self.outflow.get_area()/timestep
153        self.velocity = barrel_speed#self.discharge/outlet_depth/self.width
154
155        new_outflow_height = self.outflow.get_average_height() + outflow_extra_height
156
157        if self.use_momentum_jet :
158            # FIXME (SR) Review momentum to account for possible hydraulic jumps at outlet
159            #new_outflow_xmom = outflow.get_average_xmom() + outflow_extra_momentum[0]
160            #new_outflow_ymom = outflow.get_average_ymom() + outflow_extra_momentum[1]
161
162            new_outflow_xmom = barrel_speed*new_outflow_height*outflow_direction[0]
163            new_outflow_ymom = barrel_speed*new_outflow_height*outflow_direction[1]
164
165        else:
166            #new_outflow_xmom = outflow.get_average_xmom()
167            #new_outflow_ymom = outflow.get_average_ymom()
168
169            new_outflow_xmom = 0.0
170            new_outflow_ymom = 0.0
171
172        self.outflow.set_heights(new_outflow_height)
173        self.outflow.set_xmoms(new_outflow_xmom)
174        self.outflow.set_ymoms(new_outflow_ymom)
175
176
177    def __determine_inflow_outflow(self):
178        # Determine flow direction based on total energy difference
179
180        if self.use_velocity_head:
181            self.delta_total_energy = self.inlets[0].get_enquiry_total_energy() - self.inlets[1].get_enquiry_total_energy()
182        else:
183            self.delta_total_energy = self.inlets[0].get_enquiry_stage() - self.inlets[1].get_enquiry_stage()
184
185
186        self.inflow  = self.inlets[0]
187        self.outflow = self.inlets[1]
188       
189
190        if self.delta_total_energy < 0:
191            self.inflow  = self.inlets[1]
192            self.outflow = self.inlets[0]
193            self.delta_total_energy = -self.delta_total_energy
194
195
196    def __create_exchange_polygons(self):
197
198        """Create polygons at the end of a culvert inlet and outlet.
199        At either end two polygons will be created; one for the actual flow to pass through and one a little further away
200        for enquiring the total energy at both ends of the culvert and transferring flow.
201        """
202
203        # Calculate geometry
204        x0, y0 = self.end_points[0]
205        x1, y1 = self.end_points[1]
206
207        dx = x1 - x0
208        dy = y1 - y0
209
210        self.culvert_vector = num.array([dx, dy])
211        self.culvert_length = math.sqrt(num.sum(self.culvert_vector**2))
212        assert self.culvert_length > 0.0, 'The length of culvert is less than 0'
213
214        # Unit direction vector and normal
215        self.culvert_vector /= self.culvert_length                      # Unit vector in culvert direction
216        self.culvert_normal = num.array([-dy, dx])/self.culvert_length  # Normal vector
217
218        # Short hands
219        w = 0.5*self.width*self.culvert_normal # Perpendicular vector of 1/2 width
220        h = self.apron*self.culvert_vector    # Vector of length=height in the
221                             # direction of the culvert
222
223        gap = (1 + self.enquiry_gap)*h
224
225        self.inlet_polygons = []
226        self.inlet_equiry_points = []
227
228        # Build exchange polygon and enquiry point
229        for i in [0, 1]:
230            i0 = (2*i-1)
231            p0 = self.end_points[i] + w
232            p1 = self.end_points[i] - w
233            p2 = p1 + i0*h
234            p3 = p0 + i0*h
235            ep = self.end_points[i] + i0*gap
236
237            self.inlet_polygons.append(num.array([p0, p1, p2, p3]))
238            self.inlet_equiry_points.append(ep)
239
240        # Check that enquiry points are outside inlet polygons
241        for i in [0,1]:
242            polygon = self.inlet_polygons[i]
243            ep = self.inlet_equiry_points[i]
244           
245            area = anuga.polygon_area(polygon)
246           
247            msg = 'Polygon %s ' %(polygon)
248            msg += ' has area = %f' % area
249            assert area > 0.0, msg
250
251            msg = 'Enquiry point falls inside an exchange polygon.'
252            assert not anuga.inside_polygon(ep, polygon), msg
253           
254   
255    def discharge_routine(self):
256       
257        pass
258           
259
260    def statistics(self):
261
262
263        message  = '=====================================\n'
264        message += 'Structure Operator: %s\n' % self.label
265        message += '=====================================\n'
266
267        message += 'Structure Type: %s\n' % self.structure_type
268
269        message += 'Description\n'
270        message += '%s' % self.description
271        message += '\n'
272       
273        for i, inlet in enumerate(self.inlets):
274            message += '-------------------------------------\n'
275            message +=  'Inlet %i\n' % i
276            message += '-------------------------------------\n'
277
278            message += 'inlet triangle indices and centres\n'
279            message += '%s' % inlet.triangle_indices
280            message += '\n'
281           
282            message += '%s' % self.domain.get_centroid_coordinates()[inlet.triangle_indices]
283            message += '\n'
284
285            message += 'polygon\n'
286            message += '%s' % inlet.polygon
287            message += '\n'
288
289        message += '=====================================\n'
290
291        return message
292
293
294    def print_statistics(self):
295
296        print self.statistics()
297
298
299    def print_timestepping_statistics(self):
300
301        message = '---------------------------\n'
302        message += 'Structure report for %s:\n' % self.label
303        message += '--------------------------\n'
304        message += 'Type: %s\n' % self.structure_type
305        message += 'Discharge [m^3/s]: %.2f\n' % self.discharge
306        message += 'Velocity  [m/s]: %.2f\n' % self.velocity
307        message += 'Inlet Driving Energy %.2f\n' % self.driving_energy
308        message += 'Delta Total Energy %.2f\n' % self.delta_total_energy
309        message += 'Control at this instant: %s\n' % self.case
310
311        print message
312
313
314    def set_logging(self, flag=True):
315
316        self.logging = flag
317
318        # If flag is true open file with mode = "w" to form a clean file for logging
319        if self.logging:
320            self.log_filename = self.label + '.log'
321            log_to_file(self.log_filename, self.statistics(), mode='w')
322            log_to_file(self.log_filename, 'time,discharge,velocity,driving_energy,delta_total_energy')
323
324            #log_to_file(self.log_filename, self.culvert_type)
325
326
327    def timestepping_statistics(self):
328
329        message  = '%.5f, ' % self.domain.get_time()
330        message += '%.5f, ' % self.discharge
331        message += '%.5f, ' % self.velocity
332        message += '%.5f, ' % self.driving_energy
333        message += '%.5f' % self.delta_total_energy
334
335        return message
336
337    def log_timestepping_statistics(self):
338
339         if self.logging:
340             log_to_file(self.log_filename, self.timestepping_statistics())
341
342
343    def get_inlets(self):
344       
345        return self.inlets
346       
347       
348    def get_culvert_length(self):
349       
350        return self.culvert_length
351       
352       
353    def get_culvert_width(self):
354       
355        return self.width
356       
357       
358    def get_culvert_diameter(self):
359   
360            return self.width
361       
362       
363    def get_culvert_height(self):
364   
365        return self.height
366
367
368    def get_culvert_apron(self):
369
370        return self.apron
Note: See TracBrowser for help on using the repository browser.