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

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

Added logging to file for fractional step operators. There needs to an member function log_timestepping_statistics() implemented for each operator.

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:
64            self.label = label
65        print label
66
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.