source: trunk/anuga_core/source/anuga_parallel/parallel_boyd_box_operator.py @ 8283

Last change on this file since 8283 was 8224, checked in by janes, 12 years ago

PETE: Adding parallel fractional operator code

File size: 21.2 KB
Line 
1import anuga
2import math
3import pypar
4
5from parallel_inlet_operator import Parallel_Inlet_operator
6from parallel_structure_operator import Parallel_Structure_operator
7
8class Parallel_Boyd_box_operator(Parallel_Structure_operator):
9    """Culvert flow - transfer water from one rectangular box to another.
10    Sets up the geometry of problem
11   
12    This is the base class for culverts. Inherit from this class (and overwrite
13    compute_discharge method for specific subclasses)
14   
15    Input: Two points, pipe_size (either diameter or width, height),
16    mannings_rougness,
17    """
18
19    def __init__(self,
20                 domain,
21                 losses,
22                 width,
23                 height=None,
24                 end_points=None,
25                 exchange_lines=None,
26                 enquiry_points=None,
27                 apron=0.1,
28                 manning=0.013,
29                 enquiry_gap=0.0,
30                 use_momentum_jet=True,
31                 use_velocity_head=True,
32                 description=None,
33                 label=None,
34                 structure_type='boyd_box',
35                 logging=False,
36                 verbose=False,
37                 master_proc = 0,
38                 procs = None,
39                 inlet_master_proc = [0,0],
40                 inlet_procs = None,
41                 enquiry_proc = [0,0]):
42                     
43        Parallel_Structure_operator.__init__(self,
44                                          domain,
45                                          end_points,
46                                          exchange_lines,
47                                          enquiry_points,
48                                          width,
49                                          height,
50                                          apron,
51                                          manning,
52                                          enquiry_gap,                                                       
53                                          description,
54                                          label,
55                                          structure_type,
56                                          logging,
57                                          verbose,
58                                          master_proc,
59                                          procs,
60                                          inlet_master_proc,
61                                          inlet_procs,
62                                          enquiry_proc)
63       
64        if isinstance(losses, dict):
65            self.sum_loss = sum(losses.values())
66        elif isinstance(losses, list):
67            self.sum_loss = sum(losses)
68        else:
69            self.sum_loss = losses
70       
71        self.use_momentum_jet = use_momentum_jet
72        self.use_velocity_head = use_velocity_head
73       
74        self.culvert_length = self.get_culvert_length()
75        self.culvert_width = self.get_culvert_width()
76        self.culvert_height = self.get_culvert_height()
77
78        self.max_velocity = 10.0
79
80        self.inlets = self.get_inlets()
81
82
83        # Stats
84       
85        self.discharge = 0.0
86        self.velocity = 0.0
87       
88        self.case = 'N/A'
89
90        '''
91        print "ATTRIBUTES OF PARALLEL BOYD BOX::"
92        for attr in dir(self):
93            print "obj.%s = %s" % (attr, getattr(self, attr))
94        '''
95
96
97    def debug_discharge_routine(self):
98        local_debug ='false'
99
100        if self.use_velocity_head:
101            self.delta_total_energy = self.inlets[0].get_enquiry_total_energy() - self.inlets[1].get_enquiry_total_energy()
102        else:
103            self.delta_total_energy = self.inlets[0].get_enquiry_stage() - self.inlets[1].get_enquiry_stage()
104
105        self.inflow  = self.inlets[0]
106        self.outflow = self.inlets[1]
107
108        self.inflow_index = 0
109        self.outflow_index = 1
110
111        if self.delta_total_energy < 0:
112            self.inflow  = self.inlets[1]
113            self.outflow = self.inlets[0]
114            self.delta_total_energy = -self.delta_total_energy
115            self.inflow_index = 1
116            self.outflow_index = 0
117
118
119        if self.inflow.get_enquiry_depth() > 0.01: #this value was 0.01:
120            if local_debug =='true':
121                anuga.log.critical('Specific E & Deltat Tot E = %s, %s'
122                             % (str(self.inflow.get_enquiry_specific_energy()),
123                                str(self.delta_total_energy)))
124                anuga.log.critical('culvert type = %s' % str(culvert_type))
125            # Water has risen above inlet
126
127
128            msg = 'Specific energy at inlet is negative'
129            assert self.inflow.get_enquiry_specific_energy() >= 0.0, msg
130
131            if self.use_velocity_head :
132                self.driving_energy = self.inflow.get_enquiry_specific_energy()
133            else:
134                self.driving_energy = self.inflow.get_enquiry_depth()
135
136            depth = self.culvert_height
137            width = self.culvert_width
138            flow_width = self.culvert_width
139            # intially assume the culvert flow is controlled by the inlet
140            # check unsubmerged and submerged condition and use Min Q
141            # but ensure the correct flow area and wetted perimeter are used
142            Q_inlet_unsubmerged = 0.544*anuga.g**0.5*width*self.driving_energy**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
143            Q_inlet_submerged = 0.702*anuga.g**0.5*width*depth**0.89*self.driving_energy**0.61  # Flow based on Inlet Ctrl Inlet Submerged
144
145            # FIXME(Ole): Are these functions really for inlet control?
146            if Q_inlet_unsubmerged < Q_inlet_submerged:
147                Q = Q_inlet_unsubmerged
148                dcrit = (Q**2/anuga.g/width**2)**0.333333
149                if dcrit > depth:
150                    dcrit = depth
151                    flow_area = width*dcrit
152                    perimeter= 2.0*(width+dcrit)
153                else: # dcrit < depth
154                    flow_area = width*dcrit
155                    perimeter= 2.0*dcrit+width
156                outlet_culvert_depth = dcrit
157                self.case = 'Inlet unsubmerged Box Acts as Weir'
158            else: # Inlet Submerged but check internal culvert flow depth
159                Q = Q_inlet_submerged
160                dcrit = (Q**2/anuga.g/width**2)**0.333333
161                if dcrit > depth:
162                    dcrit = depth
163                    flow_area = width*dcrit
164                    perimeter= 2.0*(width+dcrit)
165                else: # dcrit < depth
166                    flow_area = width*dcrit
167                    perimeter= 2.0*dcrit+width
168                outlet_culvert_depth = dcrit
169                self.case = 'Inlet submerged Box Acts as Orifice'
170
171            dcrit = (Q**2/anuga.g/width**2)**0.333333
172            # May not need this .... check if same is done above
173            outlet_culvert_depth = dcrit
174            if outlet_culvert_depth > depth:
175                outlet_culvert_depth = depth  # Once again the pipe is flowing full not partfull
176                flow_area = width*depth  # Cross sectional area of flow in the culvert
177                perimeter = 2*(width+depth)
178                self.case = 'Inlet CTRL Outlet unsubmerged PIPE PART FULL'
179            else:
180                flow_area = width * outlet_culvert_depth
181                perimeter = width+2*outlet_culvert_depth
182                self.case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'
183            # Initial Estimate of Flow for Outlet Control using energy slope
184            #( may need to include Culvert Bed Slope Comparison)
185            hyd_rad = flow_area/perimeter
186            culvert_velocity = math.sqrt(self.delta_total_energy/((self.sum_loss/2/anuga.g)+(self.manning**2*self.culvert_length)/hyd_rad**1.33333))
187            Q_outlet_tailwater = flow_area * culvert_velocity
188           
189           
190            if self.delta_total_energy < self.driving_energy:
191                # Calculate flows for outlet control
192
193                # Determine the depth at the outlet relative to the depth of flow in the Culvert
194                if self.outflow.get_enquiry_depth() > depth:        # The Outlet is Submerged
195                    outlet_culvert_depth=depth
196                    flow_area=width*depth       # Cross sectional area of flow in the culvert
197                    perimeter=2.0*(width+depth)
198                    self.case = 'Outlet submerged'
199                else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
200                    dcrit = (Q**2/anuga.g/width**2)**0.333333
201                    outlet_culvert_depth=dcrit   # For purpose of calculation assume the outlet depth = Critical Depth
202                    if outlet_culvert_depth > depth:
203                        outlet_culvert_depth=depth
204                        flow_area=width*depth
205                        perimeter=2.0*(width+depth)
206                        self.case = 'Outlet is Flowing Full'
207                    else:
208                        flow_area=width*outlet_culvert_depth
209                        perimeter=(width+2.0*outlet_culvert_depth)
210                        self.case = 'Outlet is open channel flow'
211
212                hyd_rad = flow_area/perimeter
213
214
215
216                # Final Outlet control velocity using tail water
217                culvert_velocity = math.sqrt(self.delta_total_energy/((self.sum_loss/2/anuga.g)+(self.manning**2*self.culvert_length)/hyd_rad**1.33333))
218                Q_outlet_tailwater = flow_area * culvert_velocity
219
220                Q = min(Q, Q_outlet_tailwater)
221            else:
222               
223                pass
224                #FIXME(Ole): What about inlet control?
225
226            culv_froude=math.sqrt(Q**2*flow_width/(anuga.g*flow_area**3))
227       
228            if local_debug =='true':
229                anuga.log.critical('FLOW AREA = %s' % str(flow_area))
230                anuga.log.critical('PERIMETER = %s' % str(perimeter))
231                anuga.log.critical('Q final = %s' % str(Q))
232                anuga.log.critical('FROUDE = %s' % str(culv_froude))
233
234            # Determine momentum at the outlet
235            barrel_velocity = Q/(flow_area + anuga.velocity_protection/flow_area)
236
237        # END CODE BLOCK for DEPTH  > Required depth for CULVERT Flow
238
239        else: # self.inflow.get_enquiry_depth() < 0.01:
240            Q = barrel_velocity = outlet_culvert_depth = 0.0
241
242        # Temporary flow limit
243        if barrel_velocity > self.max_velocity:
244            barrel_velocity = self.max_velocity
245            Q = flow_area * barrel_velocity
246
247        return Q, barrel_velocity, outlet_culvert_depth
248
249    def discharge_routine(self):
250
251        local_debug ='false'
252
253        #Send attributes of both enquiry points to the master proc
254        if self.myid == self.master_proc:
255
256            if self.myid == self.enquiry_proc[0]:
257                enq_total_energy0 = self.inlets[0].get_enquiry_total_energy()
258                enq_stage0 = self.inlets[0].get_enquiry_stage()
259            else:
260                enq_total_energy0 = pypar.receive(self.enquiry_proc[0])
261                enq_stage0 = pypar.receive(self.enquiry_proc[0])
262
263
264            if self.myid == self.enquiry_proc[1]:
265                enq_total_energy1 = self.inlets[1].get_enquiry_total_energy()
266                enq_stage1 = self.inlets[1].get_enquiry_stage()
267            else:
268                enq_total_energy1 = pypar.receive(self.enquiry_proc[1])
269                enq_stage1 = pypar.receive(self.enquiry_proc[1])
270
271        else:
272            if self.myid == self.enquiry_proc[0]:
273                pypar.send(self.inlets[0].get_enquiry_total_energy(), self.master_proc)
274                pypar.send(self.inlets[0].get_enquiry_stage(), self.master_proc)
275
276            if self.myid == self.enquiry_proc[1]:
277                pypar.send(self.inlets[1].get_enquiry_total_energy(), self.master_proc)
278                pypar.send(self.inlets[1].get_enquiry_stage(), self.master_proc)
279
280
281        # Determine the direction of the flow
282        if self.myid == self.master_proc:
283            if self.use_velocity_head:
284                self.delta_total_energy = enq_total_energy0 - enq_total_energy1
285            else:
286                self.delta_total_energy = enq_stage0 - enq_stage1
287
288        self.inflow_index = 0
289        self.outflow_index = 1
290        # master proc orders reversal if applicable
291        if self.myid == self.master_proc:
292            # Reverse the inflow and outflow direction?
293            if self.delta_total_energy < 0:
294                self.inflow_index = 1
295                self.outflow_index = 0
296
297                self.delta_total_energy = -self.delta_total_energy
298
299                for i in self.procs:
300                    if i == self.master_proc: continue
301                    pypar.send(True, i)
302            else:
303                for i in self.procs:
304                    if i == self.master_proc: continue
305                    pypar.send(False, i)
306
307            #print "ZZZZ: Delta total energy = %f" %(self.delta_total_energy)
308        else:
309            reverse = pypar.receive(self.master_proc)
310
311            if reverse:
312                self.inflow_index = 1
313                self.outflow_index = 0
314
315        # Get attribute from inflow enquiry point
316        if self.myid == self.master_proc:
317
318            if self.myid == self.enquiry_proc[self.inflow_index]:
319                    inflow_enq_depth = self.inlets[self.inflow_index].get_enquiry_depth()
320                    inflow_enq_specific_energy = self.inlets[self.inflow_index].get_enquiry_specific_energy()
321            else:
322                    inflow_enq_depth = pypar.receive(self.enquiry_proc[self.inflow_index])
323                    inflow_enq_specific_energy = pypar.receive(self.enquiry_proc[self.inflow_index])
324        else:
325            if self.myid == self.enquiry_proc[self.inflow_index]:
326                pypar.send(self.inlets[self.inflow_index].get_enquiry_depth(), self.master_proc)
327                pypar.send(self.inlets[self.inflow_index].get_enquiry_specific_energy(), self.master_proc)
328
329        # Get attribute from outflow enquiry point
330        if self.myid == self.master_proc:
331            if self.myid == self.enquiry_proc[self.outflow_index]:
332                outflow_enq_depth = self.inlets[self.outflow_index].get_enquiry_depth()
333            else:
334                outflow_enq_depth = pypar.receive(self.enquiry_proc[self.outflow_index])
335
336            #print "ZZZZZ: outflow_enq_depth = %f" %(outflow_enq_depth)
337
338        else:
339            if self.myid == self.enquiry_proc[self.outflow_index]:
340                pypar.send(self.inlets[self.outflow_index].get_enquiry_depth(), self.master_proc)
341
342        # Master proc computes return values
343        if self.myid == self.master_proc:
344            if inflow_enq_depth > 0.01: #this value was 0.01:
345                if local_debug =='true':
346                    anuga.log.critical('Specific E & Deltat Tot E = %s, %s'
347                                 % (str(inflow_enq_specific_energy),
348                                    str(self.delta_total_energy)))
349
350                    anuga.log.critical('culvert type = %s' % str(culvert_type))
351
352                # Water has risen above inlet
353
354
355                msg = 'Specific energy at inlet is negative'
356                assert inflow_enq_specific_energy >= 0.0, msg
357
358                if self.use_velocity_head :
359                    self.driving_energy = inflow_enq_specific_energy
360                else:
361                    self.driving_energy = inflow_enq_depth
362                   
363                #print "ZZZZZ: driving energy = %f" %(self.driving_energy)
364
365                depth = self.culvert_height
366                width = self.culvert_width
367                flow_width = self.culvert_width
368                # intially assume the culvert flow is controlled by the inlet
369                # check unsubmerged and submerged condition and use Min Q
370                # but ensure the correct flow area and wetted perimeter are used
371                Q_inlet_unsubmerged = 0.544*anuga.g**0.5*width*self.driving_energy**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
372                Q_inlet_submerged = 0.702*anuga.g**0.5*width*depth**0.89*self.driving_energy**0.61  # Flow based on Inlet Ctrl Inlet Submerged
373
374                # FIXME(Ole): Are these functions really for inlet control?
375                if Q_inlet_unsubmerged < Q_inlet_submerged:
376                    Q = Q_inlet_unsubmerged
377                    dcrit = (Q**2/anuga.g/width**2)**0.333333
378                    if dcrit > depth:
379                        dcrit = depth
380                        flow_area = width*dcrit
381                        perimeter= 2.0*(width+dcrit)
382                    else: # dcrit < depth
383                        flow_area = width*dcrit
384                        perimeter= 2.0*dcrit+width
385                    outlet_culvert_depth = dcrit
386                    self.case = 'Inlet unsubmerged Box Acts as Weir'
387                else: # Inlet Submerged but check internal culvert flow depth
388                    Q = Q_inlet_submerged
389                    dcrit = (Q**2/anuga.g/width**2)**0.333333
390                    if dcrit > depth:
391                        dcrit = depth
392                        flow_area = width*dcrit
393                        perimeter= 2.0*(width+dcrit)
394                    else: # dcrit < depth
395                        flow_area = width*dcrit
396                        perimeter= 2.0*dcrit+width
397                    outlet_culvert_depth = dcrit
398                    self.case = 'Inlet submerged Box Acts as Orifice'
399
400                dcrit = (Q**2/anuga.g/width**2)**0.333333
401                # May not need this .... check if same is done above
402                outlet_culvert_depth = dcrit
403                if outlet_culvert_depth > depth:
404                    outlet_culvert_depth = depth  # Once again the pipe is flowing full not partfull
405                    flow_area = width*depth  # Cross sectional area of flow in the culvert
406                    perimeter = 2*(width+depth)
407                    self.case = 'Inlet CTRL Outlet unsubmerged PIPE PART FULL'
408                else:
409                    flow_area = width * outlet_culvert_depth
410                    perimeter = width+2*outlet_culvert_depth
411                    self.case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'
412                # Initial Estimate of Flow for Outlet Control using energy slope
413                #( may need to include Culvert Bed Slope Comparison)
414                hyd_rad = flow_area/perimeter
415                culvert_velocity = math.sqrt(self.delta_total_energy/((self.sum_loss/2/anuga.g)+(self.manning**2*self.culvert_length)/hyd_rad**1.33333))
416                Q_outlet_tailwater = flow_area * culvert_velocity
417
418
419                if self.delta_total_energy < self.driving_energy:
420                    # Calculate flows for outlet control
421
422                    # Determine the depth at the outlet relative to the depth of flow in the Culvert
423                    if outflow_enq_depth > depth:        # The Outlet is Submerged
424                        outlet_culvert_depth=depth
425                        flow_area=width*depth       # Cross sectional area of flow in the culvert
426                        perimeter=2.0*(width+depth)
427                        self.case = 'Outlet submerged'
428                    else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
429                        dcrit = (Q**2/anuga.g/width**2)**0.333333
430                        outlet_culvert_depth=dcrit   # For purpose of calculation assume the outlet depth = Critical Depth
431                        if outlet_culvert_depth > depth:
432                            outlet_culvert_depth=depth
433                            flow_area=width*depth
434                            perimeter=2.0*(width+depth)
435                            self.case = 'Outlet is Flowing Full'
436                        else:
437                            flow_area=width*outlet_culvert_depth
438                            perimeter=(width+2.0*outlet_culvert_depth)
439                            self.case = 'Outlet is open channel flow'
440
441                    hyd_rad = flow_area/perimeter
442
443
444
445                    # Final Outlet control velocity using tail water
446                    culvert_velocity = math.sqrt(self.delta_total_energy/((self.sum_loss/2/anuga.g)+(self.manning**2*self.culvert_length)/hyd_rad**1.33333))
447                    Q_outlet_tailwater = flow_area * culvert_velocity
448
449                    Q = min(Q, Q_outlet_tailwater)
450                else:
451
452                    pass
453                    #FIXME(Ole): What about inlet control?
454
455                culv_froude=math.sqrt(Q**2*flow_width/(anuga.g*flow_area**3))
456               
457                if local_debug =='true':
458                    anuga.log.critical('FLOW AREA = %s' % str(flow_area))
459                    anuga.log.critical('PERIMETER = %s' % str(perimeter))
460                    anuga.log.critical('Q final = %s' % str(Q))
461                    anuga.log.critical('FROUDE = %s' % str(culv_froude))
462
463                # Determine momentum at the outlet
464                barrel_velocity = Q/(flow_area + anuga.velocity_protection/flow_area)
465
466            # END CODE BLOCK for DEPTH  > Required depth for CULVERT Flow
467
468            else: # self.inflow.get_enquiry_depth() < 0.01:
469                Q = barrel_velocity = outlet_culvert_depth = 0.0
470
471            # Temporary flow limit
472            if barrel_velocity > self.max_velocity:
473                barrel_velocity = self.max_velocity
474                Q = flow_area * barrel_velocity
475
476            return Q, barrel_velocity, outlet_culvert_depth
477        else:
478            return None, None, None
479       
480       
Note: See TracBrowser for help on using the repository browser.