Ignore:
Timestamp:
Sep 20, 2008, 12:18:52 AM (14 years ago)
Author:
ole
Message:

Worked on energy line enquiry for culverts - this is rather slow at the moment

File:
1 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/culvert_flows/culvert_class.py

    r5586 r5777  
    1616    2008_May_08
    1717    To Ole:
    18     OK so here we need to get the Polygon Creating code to create polygons for the culvert Based on
    19     the two input Points (X0,Y0) and (X1,Y1) - need to be passed to create polygon
     18    OK so here we need to get the Polygon Creating code to create
     19    polygons for the culvert Based on
     20    the two input Points (X0,Y0) and (X1,Y1) - need to be passed
     21    to create polygon
    2022
    2123    The two centers are now passed on to create_polygon.
    2224   
    2325
    24     Input: Two points, pipe_size (either diameter or width, height), mannings_rougness,
     26    Input: Two points, pipe_size (either diameter or width, height),
     27    mannings_rougness,
    2528    inlet/outlet energy_loss_coefficients, internal_bend_coefficent,
    2629    top-down_blockage_factor and bottom_up_blockage_factor
    2730   
    2831   
    29     And the Delta H enquiry should be change from Openings in line 412 to the enquiry Polygons infront
    30     of the culvert
    31     At the moment this script uses only Depth, later we can change it to Energy...
    32 
    33     Once we have Delta H can calculate a Flow Rate and from Flow Rate an Outlet Velocity
     32    And the Delta H enquiry should be change from Openings in line 412
     33    to the enquiry Polygons infront of the culvert
     34    At the moment this script uses only Depth, later we can change it to
     35    Energy...
     36
     37    Once we have Delta H can calculate a Flow Rate and from Flow Rate
     38    an Outlet Velocity
    3439    The Outlet Velocity x Outlet Depth = Momentum to be applied at the Outlet...
    3540       
     41    Invert levels are optional. If left out they will default to the
     42    elevation at the opening.
     43       
    3644    """
    3745
     
    4654                 diameter=None,
    4755                 manning=None,          # Mannings Roughness for Culvert
    48                  invert_level0=None,    # Invert level if not the same as the Elevation on the Domain
    49                  invert_level1=None,    # Invert level if not the same as the Elevation on the Domain
     56                 invert_level0=None,    # Invert level at opening 0
     57                 invert_level1=None,    # Invert level at opening 1
    5058                 loss_exit=None,
    5159                 loss_entry=None,
     
    6573            self.diameter = diameter
    6674            if height is not None or width is not None:
    67                 msg = 'Either diameter or width&height must be specified, but not both.'
     75                msg = 'Either diameter or width&height must be specified, '
     76                msg += 'but not both.'
    6877                raise Exception, msg
    6978        else:
     
    96105       
    97106        # Set defaults
    98         if manning is None: manning = 0.012   # Set a Default Mannings Roughness for Pipe
     107        if manning is None: manning = 0.012   # Default roughness for pipe
    99108        if loss_exit is None: loss_exit = 1.00
    100109        if loss_entry is None: loss_entry = 0.50
     
    103112        if blockage_topdwn is None: blockage_topdwn=0.00
    104113        if blockage_bottup is None: blockage_bottup=0.00
    105         if culvert_routine is None: culvert_routine=boyd_generalised_culvert_model
     114        if culvert_routine is None:
     115            culvert_routine=boyd_generalised_culvert_model
     116           
    106117        if label is None: label = 'culvert_flow'
    107118        label += '_' + str(id(self))
    108        
    109119        self.label = label
    110120       
     
    174184        self.end_points = [end_point0, end_point1]
    175185        self.invert_levels = [invert_level0, invert_level1]               
    176         self.enquiry_polygons = [P['enquiry_polygon0'], P['enquiry_polygon1']]
     186        #self.enquiry_polygons = [P['enquiry_polygon0'], P['enquiry_polygon1']]
     187        self.enquiry_polylines = [P['enquiry_polygon0'][:2],
     188                                  P['enquiry_polygon1'][:2]]
    177189        self.vector = P['vector']
    178190        self.length = P['length']; assert self.length > 0.0
     
    245257            dq = domain.quantities
    246258           
    247             stage = dq['stage'].get_values(location='centroids',
    248                                            indices=opening.exchange_indices)
    249             elevation = dq['elevation'].get_values(location='centroids',
    250                                                    indices=opening.exchange_indices)
    251 
    252             # Indices corresponding to energy enquiry field for this opening
    253             coordinates = domain.get_centroid_coordinates(absolute=True) # Get all centroid points (x,y)
    254             enquiry_indices = inside_polygon(coordinates,
    255                                              self.enquiry_polygons[i])
    256 
    257             if len(enquiry_indices) == 0:
    258                 msg = 'No triangles have been identified in specified region: %s' %str(self.enquiry_polygons[i])
    259                 raise Exception, msg               
    260            
    261             # Get model values for points in enquiry polygon for this opening
    262             stage = dq['stage'].get_values(location='centroids',
    263                                            indices=enquiry_indices)
    264             xmomentum = dq['xmomentum'].get_values(location='centroids',
    265                                                    indices=enquiry_indices)
    266             ymomentum = dq['ymomentum'].get_values(location='centroids',
    267                                                    indices=enquiry_indices)
    268             elevation = dq['elevation'].get_values(location='centroids',
    269                                                    indices=enquiry_indices)
    270             depth = stage - elevation
    271 
    272             # Compute mean values of selected quantitites in the enquiry area in front of the culvert
    273             # Epsilon handles a dry cell case
    274             ux = xmomentum/(depth+velocity_protection/depth)   # Velocity (x-direction)
    275             uy = ymomentum/(depth+velocity_protection/depth)   # Velocity (y-direction)
    276             #print 'Velocity in culvert:', ux, uy, depth, xmomentum, ymomentum
    277             v = mean(sqrt(ux**2+uy**2))      # Average velocity
    278             w = mean(stage)                  # Average stage
    279 
    280             # Store values at enquiry field
    281             opening.velocity = v
    282 
    283 
    284259            # Compute mean values of selected quantitites in the exchange area in front of the culvert
    285260            # Stage and velocity comes from enquiry area and elevation from exchange area
    286261           
     262            stage = dq['stage'].get_values(location='centroids',
     263                                           indices=opening.exchange_indices)           
     264            w = mean(stage) # Average stage
     265
    287266            # Use invert level instead of elevation if specified
    288267            invert_level = self.invert_levels[i]
     
    291270            else:
    292271                elevation = dq['elevation'].get_values(location='centroids', indices=opening.exchange_indices)
    293                 z = mean(elevation)                   # Average elevation
     272                z = mean(elevation) # Average elevation
    294273
    295274            # Estimated depth above the culvert inlet
     
    302281           
    303282
    304 
    305             # Depth at exchange area used to trigger calculations
    306             stage = dq['stage'].get_values(location='centroids', indices=enquiry_indices)
    307             elevation = dq['elevation'].get_values(location='centroids', indices=enquiry_indices)
    308             depth = stage - elevation
    309             d_trigger = mean(depth)
    310 
    311 
    312 
    313283            # Ratio of depth to culvert height.
    314284            # If ratio > 1 then culvert is running full
     
    319289            opening.ratio = ratio
    320290               
     291               
    321292            # Average measures of energy in front of this opening
    322             Es = d + 0.5*v**2/g  #  Specific energy in exchange area
    323             Et = w + 0.5*v**2/g  #  Total energy in the enquiry area
    324             opening.total_energy = Et
    325             opening.specific_energy = Es           
    326            
     293            polyline = self.enquiry_polylines[i]
     294            #print 't = %.4f, opening=%d,' %(domain.time, i),
     295            opening.total_energy = domain.get_energy_through_cross_section(polyline,
     296                                                                           kind='total')           
     297            #print 'Et = %.3f m' %opening.total_energy
     298
    327299            # Store current average stage and depth with each opening object
    328300            opening.depth = d
    329             opening.depth_trigger = d_trigger           
     301            opening.depth_trigger = d # Use this for now
    330302            opening.stage = w
    331303            opening.elevation = z
     
    342314        if delta_Et > 0:
    343315            #print 'Flow U/S ---> D/S'
    344             inlet=openings[0]
    345             outlet=openings[1]
     316            inlet = openings[0]
     317            outlet = openings[1]
    346318
    347319            inlet.momentum = self.opening_momentum[0]
     
    350322        else:
    351323            #print 'Flow D/S ---> U/S'
    352             inlet=openings[1]
    353             outlet=openings[0]
     324            inlet = openings[1]
     325            outlet = openings[0]
    354326
    355327            inlet.momentum = self.opening_momentum[1]
Note: See TracChangeset for help on using the changeset viewer.