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

Last change on this file since 8019 was 8019, checked in by steve, 14 years ago

Small typo fixed

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