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

Last change on this file since 8048 was 8048, checked in by habili, 13 years ago

Changed from using polygons to polylines. The inlets are now represented by triangles that intersect or contain a polyline.

File size: 11.0 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 = []
[8048]88        polyline0 = self.inlet_polylines[0]
[7993]89        enquiry_point0 = self.inlet_equiry_points[0]
90        outward_vector0 = self.culvert_vector
[8048]91        self.inlets.append(inlet.Inlet(self.domain, polyline0, enquiry_point0, outward_vector0))
[7993]92
[8048]93        polyline1 = self.inlet_polylines[1]
94        enquiry_point1 = self.inlet_equiry_points[1]
[7993]95        outward_vector1  = - self.culvert_vector
[8048]96        self.inlets.append(inlet.Inlet(self.domain, polyline1, enquiry_point1, outward_vector1))
[7993]97
[8018]98        self.set_logging(logging)
[8008]99
[7993]100    def __call__(self):
101
[8008]102        timestep = self.domain.get_timestep()
103       
[8035]104        self.determine_inflow_outflow()
[8008]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()
[8024]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
[8008]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
[8035]177    def determine_inflow_outflow(self):
[8008]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
[7993]196    def __create_exchange_polygons(self):
197
[8048]198        """Create polylines at the end of a culvert inlet and outlet.
199        At either end two polylines will be created; one for the actual flow to pass through and one a little further away
[7993]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
[8048]220        #h = self.apron*self.culvert_vector    # Vector of length=height in the
[7993]221                             # direction of the culvert
222
[8048]223        #gap = 1.5*h + self.enquiry_gap
224        gap = 1.5*self.culvert_vector + self.enquiry_gap
[7993]225
[8048]226        self.inlet_polylines = []
[7993]227        self.inlet_equiry_points = []
228
229        # Build exchange polygon and enquiry point
230        for i in [0, 1]:
[8034]231            i0 = (2*i-1) #i0 determines the sign of the points
[7993]232            p0 = self.end_points[i] + w
233            p1 = self.end_points[i] - w
234            ep = self.end_points[i] + i0*gap
235
[8048]236            self.inlet_polylines.append(num.array([p0, p1]))
237            self.inlet_equiry_points.append(ep)           
[7993]238   
[8008]239    def discharge_routine(self):
240       
241        pass
[7993]242           
243
[8021]244    def statistics(self):
[7993]245
246
[8018]247        message  = '=====================================\n'
248        message += 'Structure Operator: %s\n' % self.label
249        message += '=====================================\n'
[7993]250
[8018]251        message += 'Structure Type: %s\n' % self.structure_type
252
253        message += 'Description\n'
254        message += '%s' % self.description
255        message += '\n'
[7993]256       
257        for i, inlet in enumerate(self.inlets):
[8018]258            message += '-------------------------------------\n'
259            message +=  'Inlet %i\n' % i
260            message += '-------------------------------------\n'
[7993]261
[8018]262            message += 'inlet triangle indices and centres\n'
263            message += '%s' % inlet.triangle_indices
264            message += '\n'
265           
266            message += '%s' % self.domain.get_centroid_coordinates()[inlet.triangle_indices]
267            message += '\n'
[7993]268
[8018]269            message += 'polygon\n'
[8048]270            message += '%s' % inlet.polyline
[8018]271            message += '\n'
[7993]272
[8018]273        message += '=====================================\n'
[7995]274
[8018]275        return message
[7995]276
[8018]277
[8021]278    def print_statistics(self):
[8018]279
[8021]280        print self.statistics()
[8018]281
282
283    def print_timestepping_statistics(self):
284
[7995]285        message = '---------------------------\n'
[8018]286        message += 'Structure report for %s:\n' % self.label
[7995]287        message += '--------------------------\n'
[8018]288        message += 'Type: %s\n' % self.structure_type
[7995]289        message += 'Discharge [m^3/s]: %.2f\n' % self.discharge
290        message += 'Velocity  [m/s]: %.2f\n' % self.velocity
[7996]291        message += 'Inlet Driving Energy %.2f\n' % self.driving_energy
[8020]292        message += 'Delta Total Energy %.2f\n' % self.delta_total_energy
[8027]293        message += 'Control at this instant: %s\n' % self.case
[7995]294
[8018]295        print message
296
297
298    def set_logging(self, flag=True):
299
300        self.logging = flag
301
302        # If flag is true open file with mode = "w" to form a clean file for logging
303        if self.logging:
304            self.log_filename = self.label + '.log'
[8021]305            log_to_file(self.log_filename, self.statistics(), mode='w')
[8018]306            log_to_file(self.log_filename, 'time,discharge,velocity,driving_energy,delta_total_energy')
307
308            #log_to_file(self.log_filename, self.culvert_type)
309
310
311    def timestepping_statistics(self):
312
313        message  = '%.5f, ' % self.domain.get_time()
314        message += '%.5f, ' % self.discharge
315        message += '%.5f, ' % self.velocity
316        message += '%.5f, ' % self.driving_energy
317        message += '%.5f' % self.delta_total_energy
318
[7995]319        return message
320
[8018]321    def log_timestepping_statistics(self):
[8008]322
[8018]323         if self.logging:
324             log_to_file(self.log_filename, self.timestepping_statistics())
325
326
[7993]327    def get_inlets(self):
328       
329        return self.inlets
330       
331       
332    def get_culvert_length(self):
333       
334        return self.culvert_length
335       
336       
337    def get_culvert_width(self):
338       
339        return self.width
340       
341       
[7998]342    def get_culvert_diameter(self):
343   
344            return self.width
345       
346       
[7993]347    def get_culvert_height(self):
348   
349        return self.height
350
351
352    def get_culvert_apron(self):
353
354        return self.apron
Note: See TracBrowser for help on using the repository browser.