Changeset 8832


Ignore:
Timestamp:
Apr 16, 2013, 8:04:09 PM (12 years ago)
Author:
steve
Message:

Fixed some problems with dem2pts which turned up in merewether case study

Location:
trunk/anuga_core/source
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/anuga_core/source/anuga/file_conversion/dem2pts.py

    r8821 r8832  
    155155    dem_elevation_r = num.reshape(dem_elevation, (nrows, ncols))
    156156    totalnopoints = nrows*ncols
     157
     158
    157159
    158160
    … …  
    282284    # Do the preceeding with numpy
    283285    #========================================
    284 
    285     start = (nrows-1)*cellsize+yllcorner
    286     stop = yllcorner-cellsize
    287     step = -cellsize
    288     y = num.arange(start, stop,step)
    289 
    290     start = xllcorner
    291     stop = (ncols)*cellsize + xllcorner
    292     step = cellsize
    293     x = num.arange(start, stop,step)
    294 
    295     #print nrows,ncols
     286    y = num.arange(nrows,dtype=num.float)
     287    y = yllcorner + (nrows-1)*cellsize - y*cellsize
     288
     289    x = num.arange(ncols,dtype=num.float)
     290    x = xllcorner + x*cellsize
    296291
    297292    xx,yy = num.meshgrid(x,y)
    298 
    299     #print xx
    300     #print yy
    301293
    302294    xx = xx.flatten()
    303295    yy = yy.flatten()
    304296
     297   
    305298    flag = num.logical_and(num.logical_and((xx <= easting_max),(xx >= easting_min)),
    306299                           num.logical_and((yy <= northing_max),(yy >= northing_min)))
    307300
    308 
    309     #print flag
    310 
    311     #print xx
    312     #print yy
    313     #print easting_min, easting_max, northing_min, northing_max
    314 
     301   
    315302    dem = dem_elevation[:].flatten()
    316303
    … …  
    321308    yy = yy[id]
    322309    dem = dem[id]
     310
    323311
    324312    clippednopoints = len(dem)

  • trunk/anuga_core/source/anuga/file_conversion/test_dem2pts.py

    r8821 r8832  
    247247
    248248
     249        #print points
     250        #print new_ref_points
     251
    249252        assert num.allclose(points, new_ref_points)
     253
     254        #print elevation
     255        #print ref_elevation
     256       
    250257        assert num.allclose(elevation, ref_elevation)
    251258
    … …  
    254261
    255262
    256         os.remove(root + '.pts')
     263        #os.remove(root + '.pts')
    257264        os.remove(root + '.dem')
    258265        os.remove(root + '.asc')

  • trunk/anuga_core/source/anuga/structures/boyd_box_operator.py

    r8828 r8832  
    128128
    129129
    130             Q, barrel_velocity, outlet_culvert_depth, case = \
    131                               boyd_box_function(depth               =self.culvert_height,
    132                                                 width               =self.culvert_width,
     130            Q, barrel_velocity, outlet_culvert_depth, flow_area, case = \
     131                              boyd_box_function(width               =self.culvert_width,
     132                                                depth               =self.culvert_height,
    133133                                                flow_width          =self.culvert_width,
    134134                                                length              =self.culvert_length,
    … …  
    159159# define separately so that can be imported in parallel code.
    160160#=============================================================================
    161 def boyd_box_function(depth, width, flow_width, length, driving_energy, delta_total_energy, outlet_enquiry_depth, sum_loss, manning):
     161def boyd_box_function(width, depth, flow_width, length, driving_energy, delta_total_energy, outlet_enquiry_depth, sum_loss, manning):
    162162
    163163    # intially assume the culvert flow is controlled by the inlet
    … …  
    263263    # END CODE BLOCK for DEPTH  > Required depth for CULVERT Flow
    264264
    265     return Q, barrel_velocity, outlet_culvert_depth, case
    266 
    267 
    268 
    269 
    270 
    271 
    272 
    273 
     265    return Q, barrel_velocity, outlet_culvert_depth, flow_area, case
     266
     267
     268
     269
     270
     271
     272
     273

  • trunk/anuga_core/source/anuga/structures/boyd_pipe_operator.py

    r8827 r8832  
    33
    44
    5 #=============================================================================
    6 # define separately so that can be imported in parallel code.
    7 #=============================================================================
    8 def discharge_routine(operator):
    9 
    10     #operator.__determine_inflow_outflow()
    11 
    12     local_debug = False
    13 
    14     if operator.use_velocity_head:
    15         operator.delta_total_energy = operator.inlets[0].get_enquiry_total_energy() - operator.inlets[1].get_enquiry_total_energy()
    16     else:
    17         operator.delta_total_energy = operator.inlets[0].get_enquiry_stage() - operator.inlets[1].get_enquiry_stage()
    18 
    19     operator.inflow  = operator.inlets[0]
    20     operator.outflow = operator.inlets[1]
    21 
    22     if operator.delta_total_energy < 0:
    23         operator.inflow  = operator.inlets[1]
    24         operator.outflow = operator.inlets[0]
    25         operator.delta_total_energy = -operator.delta_total_energy
    26 
    27     if operator.inflow.get_enquiry_depth() > 0.01: #this value was 0.01: Remember this needs to be compared to the Invert Lvl
    28         if local_debug:
    29             anuga.log.critical('Specific E & Deltat Tot E = %s, %s'
    30                          % (str(operator.inflow.get_enquiry_specific_energy()),
    31                             str(operator.delta_total_energy)))
    32             anuga.log.critical('culvert type = %s' % str(culvert_type))
    33         # Water has risen above inlet
    34 
    35 
    36         msg = 'Specific energy at inlet is negative'
    37         assert operator.inflow.get_enquiry_specific_energy() >= 0.0, msg
    38 
    39         if operator.use_velocity_head :
    40             operator.driving_energy = operator.inflow.get_enquiry_specific_energy()
    41         else:
    42             operator.driving_energy = operator.inflow.get_enquiry_depth()
    43             """
    44     For a circular pipe the Boyd method reviews 3 conditions
    45     1. Whether the Pipe Inlet is Unsubmerged (acting as weir flow into the inlet)
    46     2. Whether the Pipe Inlet is Fully Submerged (acting as an Orifice)
    47     3. Whether the energy loss in the pipe results in the Pipe being controlled by Channel Flow of the Pipe
    48 
    49     For these conditions we also would like to assess the pipe flow characteristics as it leaves the pipe
    50     """
    51 
    52     diameter = operator.culvert_diameter
    53 
    54     if operator.inflow.get_average_depth() > 0.01: #this should test against invert
    55         if local_debug:
    56             anuga.log.critical('Specific E & Deltat Tot E = %s, %s'
    57                          % (str(operator.inflow.get_average_specific_energy()),
    58                             str(operator.delta_total_energy)))
    59             anuga.log.critical('culvert type = %s' % str(culvert_type))
    60         # Water has risen above inlet
    61 
    62 
    63         msg = 'Specific energy at inlet is negative'
    64         assert operator.inflow.get_average_specific_energy() >= 0.0, msg
    65 
    66         # Calculate flows for inlet control for circular pipe
    67         Q_inlet_unsubmerged = 0.421*anuga.g**0.5*diameter**0.87*operator.inflow.get_average_specific_energy()**1.63 # Inlet Ctrl Inlet Unsubmerged
    68         Q_inlet_submerged = 0.530*anuga.g**0.5*diameter**1.87*operator.inflow.get_average_specific_energy()**0.63   # Inlet Ctrl Inlet Submerged
    69         # Note for to SUBMERGED TO OCCUR operator.inflow.get_average_specific_energy() should be > 1.2 x diameter.... Should Check !!!
    70 
    71 
    72         Q = min(Q_inlet_unsubmerged, Q_inlet_submerged)
    73 
    74         # THE LOWEST Value will Control Calcs From here
    75         # Calculate Critical Depth Based on the Adopted Flow as an Estimate
    76         dcrit1 = diameter/1.26*(Q/anuga.g**0.5*diameter**2.5)**(1/3.75)
    77         dcrit2 = diameter/0.95*(Q/anuga.g**0.5*diameter**2.5)**(1/1.95)
    78         # From Boyd Paper ESTIMATE of Dcrit has 2 criteria as
    79         if dcrit1/diameter  > 0.85:
    80             outlet_culvert_depth = dcrit2
    81         else:
    82             outlet_culvert_depth = dcrit1
    83         #outlet_culvert_depth = min(outlet_culvert_depth, diameter)
    84         # Now determine Hydraulic Radius Parameters Area & Wetted Perimeter
    85         if outlet_culvert_depth >= diameter:
    86             outlet_culvert_depth = diameter  # Once again the pipe is flowing full not partfull
    87             flow_area = (diameter/2)**2 * math.pi  # Cross sectional area of flow in the culvert
    88             perimeter = diameter * math.pi
    89             flow_width= diameter
    90             operator.case = 'Inlet CTRL Outlet submerged Circular PIPE FULL'
    91             if local_debug:
    92                 anuga.log.critical('Inlet CTRL Outlet submerged Circular '
    93                              'PIPE FULL')
    94         else:
    95             #alpha = anuga.acos(1 - outlet_culvert_depth/diameter)    # Where did this Come From ????/
    96             alpha = anuga.acos(1-2*outlet_culvert_depth/diameter)*2
    97             #flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))        # Pipe is Running Partly Full at the INLET   WHRE did this Come From ?????
    98             flow_area = diameter**2/8*(alpha - math.sin(alpha))   # Equation from  GIECK 5th Ed. Pg. B3
    99             flow_width= diameter*math.sin(alpha/2.0)
    100             perimeter = alpha*diameter/2.0
    101             operator.case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'
    102             if local_debug:
    103                 anuga.log.critical('INLET CTRL Culvert is open channel flow '
    104                              'we will for now assume critical depth')
    105                 anuga.log.critical('Q Outlet Depth and ALPHA = %s, %s, %s'
    106                              % (str(Q), str(outlet_culvert_depth),
    107                                 str(alpha)))
    108         if operator.delta_total_energy < operator.inflow.get_average_specific_energy():  #  OUTLET CONTROL !!!!
    109             # Calculate flows for outlet control
    110 
    111             # Determine the depth at the outlet relative to the depth of flow in the Culvert
    112             if operator.outflow.get_average_depth() > diameter:       # Outlet is submerged Assume the end of the Pipe is flowing FULL
    113                 outlet_culvert_depth=diameter
    114                 flow_area = (diameter/2)**2 * math.pi  # Cross sectional area of flow in the culvert
    115                 perimeter = diameter * math.pi
    116                 flow_width= diameter
    117                 operator.case = 'Outlet submerged'
    118                 if local_debug:
    119                     anuga.log.critical('Outlet submerged')
    120             else:   # Culvert running PART FULL for PART OF ITS LENGTH   Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
    121                 # IF  operator.outflow.get_average_depth() < diameter
    122                 dcrit1 = diameter/1.26*(Q/anuga.g**0.5*diameter**2.5)**(1/3.75)
    123                 dcrit2 = diameter/0.95*(Q/anuga.g**0.5*diameter**2.5)**(1/1.95)
    124                 if dcrit1/diameter >0.85:
    125                     outlet_culvert_depth= dcrit2
    126                 else:
    127                     outlet_culvert_depth = dcrit1
    128                 if outlet_culvert_depth > diameter:
    129                     outlet_culvert_depth = diameter  # Once again the pipe is flowing full not partfull
    130                     flow_area = (diameter/2)**2 * math.pi  # Cross sectional area of flow in the culvert
    131                     perimeter = diameter * math.pi
    132                     flow_width= diameter
    133                     operator.case = 'Outlet unsubmerged PIPE FULL'
    134                     if local_debug:
    135                         anuga.log.critical('Outlet unsubmerged PIPE FULL')
    136                 else:
    137                     alpha = anuga.acos(1-2*outlet_culvert_depth/diameter)*2
    138                     flow_area = diameter**2/8*(alpha - math.sin(alpha))   # Equation from  GIECK 5th Ed. Pg. B3
    139                     flow_width= diameter*math.sin(alpha/2.0)
    140                     perimeter = alpha*diameter/2.0
    141                     operator.case = 'Outlet is open channel flow we will for now assume critical depth'
    142                     if local_debug:
    143                         anuga.log.critical('Q Outlet Depth and ALPHA = %s, %s, %s'
    144                                      % (str(Q), str(outlet_culvert_depth),
    145                                         str(alpha)))
    146                         anuga.log.critical('Outlet is open channel flow we '
    147                                      'will for now assume critical depth')
    148         if local_debug:
    149             anuga.log.critical('FLOW AREA = %s' % str(flow_area))
    150             anuga.log.critical('PERIMETER = %s' % str(perimeter))
    151             anuga.log.critical('Q Interim = %s' % str(Q))
    152         hyd_rad = flow_area/perimeter
    153 
    154 
    155 
    156         # Outlet control velocity using tail water
    157         if local_debug:
    158             anuga.log.critical('GOT IT ALL CALCULATING Velocity')
    159             anuga.log.critical('HydRad = %s' % str(hyd_rad))
    160         # Calculate Pipe Culvert Outlet Control Velocity.... May need initial Estimate First ??
    161 
    162         culvert_velocity = math.sqrt(operator.delta_total_energy/((operator.sum_loss/2/anuga.g)+\
    163                                                                   (operator.manning**2*operator.culvert_length)/hyd_rad**1.33333))
    164         Q_outlet_tailwater = flow_area * culvert_velocity
    165 
    166 
    167         if local_debug:
    168             anuga.log.critical('VELOCITY = %s' % str(culvert_velocity))
    169             anuga.log.critical('Outlet Ctrl Q = %s' % str(Q_outlet_tailwater))
    170 
    171         Q = min(Q, Q_outlet_tailwater)
    172         if local_debug:
    173             anuga.log.critical('%s,%.3f,%.3f'
    174                          % ('dcrit 1 , dcit2 =',dcrit1,dcrit2))
    175             anuga.log.critical('%s,%.3f,%.3f,%.3f'
    176                          % ('Q and Velocity and Depth=', Q,
    177                             culvert_velocity, outlet_culvert_depth))
    178 
    179         culv_froude=math.sqrt(Q**2*flow_width/(anuga.g*flow_area**3))
    180         if local_debug:
    181             anuga.log.critical('FLOW AREA = %s' % str(flow_area))
    182             anuga.log.critical('PERIMETER = %s' % str(perimeter))
    183             anuga.log.critical('Q final = %s' % str(Q))
    184             anuga.log.critical('FROUDE = %s' % str(culv_froude))
    185 
    186         # Determine momentum at the outlet
    187         barrel_velocity = Q/(flow_area + anuga.velocity_protection/flow_area)
    188 
    189     else: # operator.inflow.get_average_depth() < 0.01:
    190         Q = barrel_velocity = outlet_culvert_depth = 0.0
    191 
    192     # Temporary flow limit
    193     if barrel_velocity > operator.max_velocity:
    194         barrel_velocity = operator.max_velocity
    195         Q = flow_area * barrel_velocity
    196 
    197     return Q, barrel_velocity, outlet_culvert_depth
    198 
    199        
    200        
    201        
    2025
    2036
    … …  
    21720    mannings_rougness,
    21821    """
    219 
    220     discharge_routine = discharge_routine
    22122
    22223
    … …  
    28687       
    28788
     89
     90    def discharge_routine(self):
     91
     92        local_debug = False
     93
     94        if self.use_velocity_head:
     95            self.delta_total_energy = \
     96                 self.inlets[0].get_enquiry_total_energy() - self.inlets[1].get_enquiry_total_energy()
     97        else:
     98            self.delta_total_energy = \
     99                 self.inlets[0].get_enquiry_stage() - self.inlets[1].get_enquiry_stage()
     100
     101        self.inflow  = self.inlets[0]
     102        self.outflow = self.inlets[1]
     103
     104        if self.delta_total_energy < 0:
     105            self.inflow  = self.inlets[1]
     106            self.outflow = self.inlets[0]
     107            self.delta_total_energy = -self.delta_total_energy
     108
     109
     110        if self.inflow.get_enquiry_depth() > 0.01: #this value was 0.01:
     111
     112            if local_debug:
     113                anuga.log.critical('Specific E & Deltat Tot E = %s, %s'
     114                             % (str(self.inflow.get_enquiry_specific_energy()),
     115                                str(self.delta_total_energy)))
     116                anuga.log.critical('culvert type = %s' % str(culvert_type))
     117            # Water has risen above inlet
     118
     119
     120            msg = 'Specific energy at inlet is negative'
     121            assert self.inflow.get_enquiry_specific_energy() >= 0.0, msg
     122
     123            if self.use_velocity_head :
     124                self.driving_energy = self.inflow.get_enquiry_specific_energy()
     125            else:
     126                self.driving_energy = self.inflow.get_enquiry_depth()
     127
     128
     129
     130            Q, barrel_velocity, outlet_culvert_depth, flow_area, case = \
     131                              boyd_pipe_function(depth               =self.inflow.get_enquiry_depth(),
     132                                                 diameter            =self.culvert_diameter,
     133                                                 length              =self.culvert_length,
     134                                                 driving_energy      =self.driving_energy,
     135                                                 delta_total_energy  =self.delta_total_energy,
     136                                                 outlet_enquiry_depth=self.outflow.get_enquiry_depth(),
     137                                                 sum_loss            =self.sum_loss,
     138                                                 manning             =self.manning)
     139        else:
     140             Q = barrel_velocity = outlet_culvert_depth = 0.0
     141             case = 'Inlet dry'
     142
     143
     144        self.case = case
     145
     146
     147        # Temporary flow limit
     148        if barrel_velocity > self.max_velocity:
     149            barrel_velocity = self.max_velocity
     150            Q = flow_area * barrel_velocity
     151
     152        return Q, barrel_velocity, outlet_culvert_depth
     153
     154#=============================================================================
     155# define separately so that can be imported in parallel code.
     156#=============================================================================
     157def boyd_pipe_function(depth, diameter, length, driving_energy, delta_total_energy, outlet_enquiry_depth, sum_loss, manning):
     158
     159
     160    local_debug = False
     161
     162
     163    """
     164    For a circular pipe the Boyd method reviews 3 conditions
     165    1. Whether the Pipe Inlet is Unsubmerged (acting as weir flow into the inlet)
     166    2. Whether the Pipe Inlet is Fully Submerged (acting as an Orifice)
     167    3. Whether the energy loss in the pipe results in the Pipe being controlled by Channel Flow of the Pipe
     168
     169    For these conditions we also would like to assess the pipe flow characteristics as it leaves the pipe
     170    """
     171
     172    # Calculate flows for inlet control for circular pipe
     173    Q_inlet_unsubmerged = 0.421*anuga.g**0.5*diameter**0.87*driving_energy**1.63 # Inlet Ctrl Inlet Unsubmerged
     174    Q_inlet_submerged = 0.530*anuga.g**0.5*diameter**1.87*driving_energy**0.63   # Inlet Ctrl Inlet Submerged
     175    # Note for to SUBMERGED TO OCCUR operator.inflow.get_average_specific_energy() should be > 1.2 x diameter.... Should Check !!!
     176
     177
     178    Q = min(Q_inlet_unsubmerged, Q_inlet_submerged)
     179
     180    # THE LOWEST Value will Control Calcs From here
     181    # Calculate Critical Depth Based on the Adopted Flow as an Estimate
     182    dcrit1 = diameter/1.26*(Q/anuga.g**0.5*diameter**2.5)**(1/3.75)
     183    dcrit2 = diameter/0.95*(Q/anuga.g**0.5*diameter**2.5)**(1/1.95)
     184    # From Boyd Paper ESTIMATE of Dcrit has 2 criteria as
     185    if dcrit1/diameter  > 0.85:
     186        outlet_culvert_depth = dcrit2
     187    else:
     188        outlet_culvert_depth = dcrit1
     189    #outlet_culvert_depth = min(outlet_culvert_depth, diameter)
     190    # Now determine Hydraulic Radius Parameters Area & Wetted Perimeter
     191    if outlet_culvert_depth >= diameter:
     192        outlet_culvert_depth = diameter  # Once again the pipe is flowing full not partfull
     193        flow_area = (diameter/2)**2 * math.pi  # Cross sectional area of flow in the culvert
     194        perimeter = diameter * math.pi
     195        flow_width= diameter
     196        case = 'Inlet CTRL Outlet submerged Circular PIPE FULL'
     197        if local_debug:
     198            anuga.log.critical('Inlet CTRL Outlet submerged Circular '
     199                         'PIPE FULL')
     200    else:
     201        #alpha = anuga.acos(1 - outlet_culvert_depth/diameter)    # Where did this Come From ????/
     202        alpha = anuga.acos(1-2*outlet_culvert_depth/diameter)*2
     203        #flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))        # Pipe is Running Partly Full at the INLET   WHRE did this Come From ?????
     204        flow_area = diameter**2/8*(alpha - math.sin(alpha))   # Equation from  GIECK 5th Ed. Pg. B3
     205        flow_width= diameter*math.sin(alpha/2.0)
     206        perimeter = alpha*diameter/2.0
     207        case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'
     208        if local_debug:
     209            anuga.log.critical('INLET CTRL Culvert is open channel flow '
     210                         'we will for now assume critical depth')
     211            anuga.log.critical('Q Outlet Depth and ALPHA = %s, %s, %s'
     212                         % (str(Q), str(outlet_culvert_depth),
     213                            str(alpha)))
     214    if delta_total_energy < driving_energy:  #  OUTLET CONTROL !!!!
     215        # Calculate flows for outlet control
     216
     217        # Determine the depth at the outlet relative to the depth of flow in the Culvert
     218        if outlet_enquiry_depth > diameter:       # Outlet is submerged Assume the end of the Pipe is flowing FULL
     219            outlet_culvert_depth=diameter
     220            flow_area = (diameter/2)**2 * math.pi  # Cross sectional area of flow in the culvert
     221            perimeter = diameter * math.pi
     222            flow_width= diameter
     223            case = 'Outlet submerged'
     224            if local_debug:
     225                anuga.log.critical('Outlet submerged')
     226        else:   # Culvert running PART FULL for PART OF ITS LENGTH   Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
     227            # IF  operator.outflow.get_average_depth() < diameter
     228            dcrit1 = diameter/1.26*(Q/anuga.g**0.5*diameter**2.5)**(1/3.75)
     229            dcrit2 = diameter/0.95*(Q/anuga.g**0.5*diameter**2.5)**(1/1.95)
     230            if dcrit1/diameter >0.85:
     231                outlet_culvert_depth= dcrit2
     232            else:
     233                outlet_culvert_depth = dcrit1
     234            if outlet_culvert_depth > diameter:
     235                outlet_culvert_depth = diameter  # Once again the pipe is flowing full not partfull
     236                flow_area = (diameter/2)**2 * math.pi  # Cross sectional area of flow in the culvert
     237                perimeter = diameter * math.pi
     238                flow_width= diameter
     239                case = 'Outlet unsubmerged PIPE FULL'
     240                if local_debug:
     241                    anuga.log.critical('Outlet unsubmerged PIPE FULL')
     242            else:
     243                alpha = anuga.acos(1-2*outlet_culvert_depth/diameter)*2
     244                flow_area = diameter**2/8*(alpha - math.sin(alpha))   # Equation from  GIECK 5th Ed. Pg. B3
     245                flow_width= diameter*math.sin(alpha/2.0)
     246                perimeter = alpha*diameter/2.0
     247                case = 'Outlet is open channel flow we will for now assume critical depth'
     248                if local_debug:
     249                    anuga.log.critical('Q Outlet Depth and ALPHA = %s, %s, %s'
     250                                 % (str(Q), str(outlet_culvert_depth),
     251                                    str(alpha)))
     252                    anuga.log.critical('Outlet is open channel flow we '
     253                                 'will for now assume critical depth')
     254    if local_debug:
     255        anuga.log.critical('FLOW AREA = %s' % str(flow_area))
     256        anuga.log.critical('PERIMETER = %s' % str(perimeter))
     257        anuga.log.critical('Q Interim = %s' % str(Q))
     258    hyd_rad = flow_area/perimeter
     259
     260
     261
     262    # Outlet control velocity using tail water
     263    if local_debug:
     264        anuga.log.critical('GOT IT ALL CALCULATING Velocity')
     265        anuga.log.critical('HydRad = %s' % str(hyd_rad))
     266    # Calculate Pipe Culvert Outlet Control Velocity.... May need initial Estimate First ??
     267
     268    culvert_velocity = math.sqrt(delta_total_energy/((sum_loss/2/anuga.g)+\
     269                                                              (manning**2*length)/hyd_rad**1.33333))
     270    Q_outlet_tailwater = flow_area * culvert_velocity
     271
     272
     273    if local_debug:
     274        anuga.log.critical('VELOCITY = %s' % str(culvert_velocity))
     275        anuga.log.critical('Outlet Ctrl Q = %s' % str(Q_outlet_tailwater))
     276
     277    Q = min(Q, Q_outlet_tailwater)
     278    if local_debug:
     279        anuga.log.critical('%s,%.3f,%.3f'
     280                     % ('dcrit 1 , dcit2 =',dcrit1,dcrit2))
     281        anuga.log.critical('%s,%.3f,%.3f,%.3f'
     282                     % ('Q and Velocity and Depth=', Q,
     283                        culvert_velocity, outlet_culvert_depth))
     284
     285    culv_froude=math.sqrt(Q**2*flow_width/(anuga.g*flow_area**3))
     286    if local_debug:
     287        anuga.log.critical('FLOW AREA = %s' % str(flow_area))
     288        anuga.log.critical('PERIMETER = %s' % str(perimeter))
     289        anuga.log.critical('Q final = %s' % str(Q))
     290        anuga.log.critical('FROUDE = %s' % str(culv_froude))
     291
     292    # Determine momentum at the outlet
     293    barrel_velocity = Q/(flow_area + anuga.velocity_protection/flow_area)
     294
     295
     296    return Q, barrel_velocity, outlet_culvert_depth, flow_area, case
     297
     298
     299
     300

  • trunk/anuga_core/source/anuga_parallel/parallel_boyd_box_operator.py

    r8828 r8832  
    393393
    394394           
    395                 Q, barrel_velocity, outlet_culvert_depth, case = \
     395                Q, barrel_velocity, outlet_culvert_depth, flow_area, case = \
    396396                              boyd_box_function(depth               =self.culvert_height,
    397397                                                width               =self.culvert_width,

  • trunk/anuga_core/source/anuga_parallel/parallel_boyd_pipe_operator.py

    r8826 r8832  
    22import math
    33import pypar
     4
     5from anuga.structures.boyd_pipe_operator import boyd_pipe_function
    46
    57from parallel_inlet_operator import Parallel_Inlet_operator
    … …  
    7779        self.culvert_width = self.get_culvert_width()
    7880        self.culvert_height = self.get_culvert_height()
     81        self.culvert_diameter = self.get_culvert_diameter()
    7982
    8083        self.max_velocity = 10.0
    … …  
    307310        # master proc orders reversal if applicable
    308311        if self.myid == self.master_proc:
     312
     313
    309314            # Reverse the inflow and outflow direction?
    310315            if self.delta_total_energy < 0:
    … …  
    357362                pypar.send(self.inlets[self.outflow_index].get_enquiry_depth(), self.master_proc)
    358363
     364
     365
     366
    359367        # Master proc computes return values
    360368        if self.myid == self.master_proc:
     369
     370            #inflow_enq_specific_energy
     371
     372
    361373            if inflow_enq_depth > 0.01: #this value was 0.01:
    362374                if local_debug:
    … …  
    377389                else:
    378390                    self.driving_energy = inflow_enq_depth
    379                    
    380                 #print "ZZZZZ: driving energy = %f" %(self.driving_energy)
    381 
    382                 depth = self.culvert_height
    383                 width = self.culvert_width
    384                 flow_width = self.culvert_width
    385                 # intially assume the culvert flow is controlled by the inlet
    386                 # check unsubmerged and submerged condition and use Min Q
    387                 # but ensure the correct flow area and wetted perimeter are used
    388                 Q_inlet_unsubmerged = 0.544*anuga.g**0.5*width*self.driving_energy**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
    389                 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
    390 
    391                 # FIXME(Ole): Are these functions really for inlet control?
    392                 if Q_inlet_unsubmerged < Q_inlet_submerged:
    393                     Q = Q_inlet_unsubmerged
    394                     dcrit = (Q**2/anuga.g/width**2)**0.333333
    395                     if dcrit > depth:
    396                         dcrit = depth
    397                         flow_area = width*dcrit
    398                         perimeter= 2.0*(width+dcrit)
    399                     else: # dcrit < depth
    400                         flow_area = width*dcrit
    401                         perimeter= 2.0*dcrit+width
    402                     outlet_culvert_depth = dcrit
    403                     self.case = 'Inlet unsubmerged Box Acts as Weir'
    404                 else: # Inlet Submerged but check internal culvert flow depth
    405                     Q = Q_inlet_submerged
    406                     dcrit = (Q**2/anuga.g/width**2)**0.333333
    407                     if dcrit > depth:
    408                         dcrit = depth
    409                         flow_area = width*dcrit
    410                         perimeter= 2.0*(width+dcrit)
    411                     else: # dcrit < depth
    412                         flow_area = width*dcrit
    413                         perimeter= 2.0*dcrit+width
    414                     outlet_culvert_depth = dcrit
    415                     self.case = 'Inlet submerged Box Acts as Orifice'
    416 
    417                 dcrit = (Q**2/anuga.g/width**2)**0.333333
    418                 # May not need this .... check if same is done above
    419                 outlet_culvert_depth = dcrit
    420                 if outlet_culvert_depth > depth:
    421                     outlet_culvert_depth = depth  # Once again the pipe is flowing full not partfull
    422                     flow_area = width*depth  # Cross sectional area of flow in the culvert
    423                     perimeter = 2*(width+depth)
    424                     self.case = 'Inlet CTRL Outlet unsubmerged PIPE PART FULL'
    425                 else:
    426                     flow_area = width * outlet_culvert_depth
    427                     perimeter = width+2*outlet_culvert_depth
    428                     self.case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'
    429                 # Initial Estimate of Flow for Outlet Control using energy slope
    430                 #( may need to include Culvert Bed Slope Comparison)
    431                 hyd_rad = flow_area/perimeter
    432                 culvert_velocity = math.sqrt(self.delta_total_energy/((self.sum_loss/2/anuga.g)+(self.manning**2*self.culvert_length)/hyd_rad**1.33333))
    433                 Q_outlet_tailwater = flow_area * culvert_velocity
    434 
    435 
    436                 if self.delta_total_energy < self.driving_energy:
    437                     # Calculate flows for outlet control
    438 
    439                     # Determine the depth at the outlet relative to the depth of flow in the Culvert
    440                     if outflow_enq_depth > depth:        # The Outlet is Submerged
    441                         outlet_culvert_depth=depth
    442                         flow_area=width*depth       # Cross sectional area of flow in the culvert
    443                         perimeter=2.0*(width+depth)
    444                         self.case = 'Outlet submerged'
    445                     else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
    446                         dcrit = (Q**2/anuga.g/width**2)**0.333333
    447                         outlet_culvert_depth=dcrit   # For purpose of calculation assume the outlet depth = Critical Depth
    448                         if outlet_culvert_depth > depth:
    449                             outlet_culvert_depth=depth
    450                             flow_area=width*depth
    451                             perimeter=2.0*(width+depth)
    452                             self.case = 'Outlet is Flowing Full'
    453                         else:
    454                             flow_area=width*outlet_culvert_depth
    455                             perimeter=(width+2.0*outlet_culvert_depth)
    456                             self.case = 'Outlet is open channel flow'
    457 
    458                     hyd_rad = flow_area/perimeter
    459 
    460 
    461 
    462                     # Final Outlet control velocity using tail water
    463                     culvert_velocity = math.sqrt(self.delta_total_energy/((self.sum_loss/2/anuga.g)+(self.manning**2*self.culvert_length)/hyd_rad**1.33333))
    464                     Q_outlet_tailwater = flow_area * culvert_velocity
    465 
    466                     Q = min(Q, Q_outlet_tailwater)
    467                 else:
    468 
    469                     pass
    470                     #FIXME(Ole): What about inlet control?
    471 
    472                 culv_froude=math.sqrt(Q**2*flow_width/(anuga.g*flow_area**3))
    473                
    474                 if local_debug:
    475                     anuga.log.critical('FLOW AREA = %s' % str(flow_area))
    476                     anuga.log.critical('PERIMETER = %s' % str(perimeter))
    477                     anuga.log.critical('Q final = %s' % str(Q))
    478                     anuga.log.critical('FROUDE = %s' % str(culv_froude))
    479 
    480                 # Determine momentum at the outlet
    481                 barrel_velocity = Q/(flow_area + anuga.velocity_protection/flow_area)
     391
     392
     393                Q, barrel_velocity, outlet_culvert_depth, flow_area, case = \
     394                              boyd_pipe_function(depth               =inflow_enq_depth,
     395                                                diameter             =self.culvert_diameter,
     396                                                length               =self.culvert_length,
     397                                                driving_energy       =self.driving_energy,
     398                                                delta_total_energy   =self.delta_total_energy,
     399                                                outlet_enquiry_depth =outflow_enq_depth,
     400                                                sum_loss             =self.sum_loss,
     401                                                manning              =self.manning)
     402
    482403
    483404            # END CODE BLOCK for DEPTH  > Required depth for CULVERT Flow
    … …  
    485406            else: # self.inflow.get_enquiry_depth() < 0.01:
    486407                Q = barrel_velocity = outlet_culvert_depth = 0.0
     408                case = 'Inlet dry'
     409
     410
     411            self.case = case
    487412
    488413            # Temporary flow limit
    … …  
    491416                Q = flow_area * barrel_velocity
    492417
     418
     419
    493420            return Q, barrel_velocity, outlet_culvert_depth
    494421        else:
    495422            return None, None, None
    496        
    497        

  • trunk/anuga_core/source/anuga_parallel/parallel_structure_operator.py

    r8825 r8832  
    510510       
    511511    def get_culvert_diameter(self):
    512         return self.diamter
     512        return self.diameter
    513513       
    514514       

  • trunk/anuga_core/source/anuga_parallel/test_parallel_boyd_pipe_operator.py

    r8825 r8832  
    353353
    354354
    355     finalize()
    356    
    357 
     355   
     356
Note: See TracChangeset for help on using the changeset viewer.