Changeset 7176


Ignore:
Timestamp:
Jun 10, 2009, 5:28:35 PM (13 years ago)
Author:
rwilson
Message:

Back-merge from Numeric to numpy.

Location:
branches/numpy/anuga
Files:
12 edited

Legend:

Unmodified
Added
Removed
  • branches/numpy/anuga/abstract_2d_finite_volumes/domain.py

    r7035 r7176  
    936936                k = self.k = triangle_id
    937937
    938             x, y = self.get_centroid_coordinates()[k]
     938            x, y = self.get_centroid_coordinates(absolute=True)[k]
    939939            radius = self.get_radii()[k]
    940940            area = self.get_areas()[k]
  • branches/numpy/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py

    r7037 r7176  
    103103                 verbose=False):
    104104        Boundary.__init__(self)
     105
    105106        self.default_boundary = default_boundary
    106107        self.default_boundary_invoked = False    # Flag
  • branches/numpy/anuga/abstract_2d_finite_volumes/quantity.py

    r7038 r7176  
    11781178        # Edges have already been deprecated in set_values, see changeset:5521,
    11791179        # but *might* be useful in get_values. Any thoughts anyone?
    1180         # YES (Ole): Edge values are necessary for volumetric balance check
     1180        # YES (Ole): Edge values are necessary for volumetric balance check and inflow boundary. Keep them.
    11811181
    11821182        if location not in ['vertices', 'centroids',
  • branches/numpy/anuga/config.py

    r6553 r7176  
    88
    99################################################################################
    10 # numerical constants
     10# Numerical constants
    1111################################################################################
    1212
     
    1717                                     # least squares fitting
    1818single_precision = 1.0e-6            # Smallest single precision number
    19 velocity_protection = 1.0e-6                                     
     19velocity_protection = 1.0e-6         # Used to compute velocity from momentum
     20                                     # See section 7.4 on Flux limiting
     21                                     # in the user manual
     22                           
    2023
    2124################################################################################
     
    2629version_filename = 'stored_version_info.py'
    2730default_datadir = '.'
    28 time_format = '%d/%m/%y %H:%M:%S'
    29 umask = 002  # Controls file and directory permission created by anuga
     31time_format = '%d/%m/%y %H:%M:%S'    # Used with timefile2netcdf
     32umask = 002  # Controls file and directory permission created by anuga (UNIX)
    3033default_boundary_tag = 'exterior'
    3134
  • branches/numpy/anuga/coordinate_transforms/geo_reference.py

    r7061 r7176  
    142142        self.zone = int(infile.zone[0])
    143143
    144 ##        # Fix some assertion failures
    145 ##        if isinstance(self.zone, num.ndarray) and self.zone.shape == ():
    146 ##            self.zone = self.zone[0]
    147 ##        if (isinstance(self.xllcorner, num.ndarray) and
    148 ##                self.xllcorner.shape == ()):
    149 ##            self.xllcorner = self.xllcorner[0]
    150 ##        if (isinstance(self.yllcorner, num.ndarray) and
    151 ##                self.yllcorner.shape == ()):
    152 ##            self.yllcorner = self.yllcorner[0]
    153 ##
    154 ##        assert (self.xllcorner.dtype.kind in num.typecodes['Float'] or
    155 ##                self.xllcorner.dtype.kind in num.typecodes['Integer'])
    156 ##        assert (self.yllcorner.dtype.kind in num.typecodes['Float'] or
    157 ##                self.yllcorner.dtype.kind in num.typecodes['Integer'])
    158 ##        assert (self.zone.dtype.kind in num.typecodes['Integer'])
    159 
    160144        try:
    161145            self.false_easting = int(infile.false_easting[0])
  • branches/numpy/anuga/culvert_flows/culvert_routines.py

    r7035 r7176  
    11"""Collection of culvert routines for use with Culvert_flow in culvert_class
     2
     3This module holds various routines to determine FLOW through CULVERTS and SIMPLE BRIDGES
     4
     5
     6
    27
    38Usage:
     
    3338                                   log_filename=None):
    3439
    35     """Boyd's generalisation of the US department of transportation culvert
    36     model
    37      
    38     The quantity of flow passing through a culvert is controlled by many factors
    39     It could be that the culvert is controlled by the inlet only, with it
    40     being unsubmerged this is effectively equivalent to the weir Equation
    41     Else the culvert could be controlled by the inlet, with it being
    42     submerged, this is effectively the Orifice Equation
    43     Else it may be controlled by down stream conditions where depending on
    44     the down stream depth, the momentum in the culvert etc. flow is controlled
     40    """Boyd's generalisation of the US department of transportation culvert methods
     41       
     42        WARNING THIS IS A SIMPLISTIC APPROACH and OUTLET VELOCITIES ARE LIMITED TO EITHER
     43        FULL PIPE OR CRITICAL DEPTH ONLY
     44        For Supercritical flow this is UNDERESTIMATING the Outlet Velocity
     45        The obtain the CORRECT velocity requires an iteration of Depth to Establish
     46        the Normal Depth of flow in the pipe.
     47       
     48        It is proposed to provide this in a seperate routine called
     49    boyd_generalised_culvert_model_complex
     50       
     51    The Boyd Method is based on methods described by the following:
     52        1.
     53        US Dept. Transportation Federal Highway Administration (1965)
     54        Hydraulic Chart for Selection of Highway Culverts.
     55        Hydraulic Engineering Circular No. 5 US Government Printing
     56        2.
     57        US Dept. Transportation Federal Highway Administration (1972)
     58        Capacity charts for the Hydraulic design of highway culverts.
     59        Hydraulic Engineering Circular No. 10 US Government Printing
     60    These documents provide around 60 charts for various configurations of culverts and inlets.
     61       
     62        Note these documents have been superceded by:
     63        2005 Hydraulic Design of Highway Culverts, Hydraulic Design Series No. 5 (HDS-5),
     64        Which combines culvert design information previously contained in Hydraulic Engineering Circulars
     65        (HEC) No. 5, No. 10, and No. 13 with hydrologic, storage routing, and special culvert design information.
     66        HEC-5 provides 20 Charts
     67        HEC-10 Provides an additional 36 Charts
     68        HEC-13 Discusses the Design of improved more efficient inlets
     69        HDS-5 Provides 60 sets of Charts
     70
     71        In 1985 Professor Michael Boyd Published "Head-Discharge Relations for Culverts", and in
     72        1987 published "Generalised Head Discharge Equations for Culverts".
     73        These papers reviewed the previous work by the US DOT and provided a simplistic approach for 3 configurations.
     74       
     75        It may be possible to extend the same approach for additional charts in the original work, but to date this has not been done.
     76        The additional charts cover a range of culvert shapes and inlet configurations
     77       
     78   
    4579    """
    4680
     
    4983    from anuga.utilities.numerical_tools import safe_acos as acos
    5084
    51 
     85    local_debug ='false'
    5286    if inlet_depth > 0.1: #this value was 0.01:
     87        if local_debug =='true':
     88            print 'Specific E & Deltat Tot E = ',inlet_specific_energy,delta_total_energy
     89            print 'culvert typ = ',culvert_type
    5390        # Water has risen above inlet
    5491       
     
    6198                     
    6299
    63         if culvert_type == 'circle':
    64             # Round culvert (use width as diameter)
    65             diameter = culvert_width           
    66 
     100        if culvert_type =='circle':
     101            # Round culvert (use height as diameter)
     102            diameter = culvert_height   
     103
     104            """
     105            For a circular pipe the Boyd method reviews 3 conditions
     106            1. Whether the Pipe Inlet is Unsubmerged (acting as weir flow into the inlet)
     107            2. Whether the Pipe Inlet is Fully Submerged (acting as an Orifice)
     108            3. Whether the energy loss in the pipe results in the Pipe being controlled by Channel Flow of the Pipe
     109           
     110            For these conditions we also would like to assess the pipe flow characteristics as it leaves the pipe
     111            """
     112           
    67113            # Calculate flows for inlet control           
    68114            Q_inlet_unsubmerged = 0.421*g**0.5*diameter**0.87*inlet_specific_energy**1.63 # Inlet Ctrl Inlet Unsubmerged
    69115            Q_inlet_submerged = 0.530*g**0.5*diameter**1.87*inlet_specific_energy**0.63   # Inlet Ctrl Inlet Submerged
     116            # Note for to SUBMERGED TO OCCUR inlet_specific_energy should be > 1.2 x diameter.... Should Check !!!
    70117
    71118            if log_filename is not None:                                       
    72119                s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' % (Q_inlet_unsubmerged, Q_inlet_submerged)
    73120                log_to_file(log_filename, s)
    74 
    75 
    76             # FIXME(Ole): Are these functions really for inlet control?                   
    77             if Q_inlet_unsubmerged < Q_inlet_submerged:
    78                 Q = Q_inlet_unsubmerged
    79                 alpha = acos(1 - inlet_depth/diameter)
    80                 flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))
    81                 outlet_culvert_depth = inlet_depth
    82                 case = 'Inlet unsubmerged'
    83             else:   
    84                 Q = Q_inlet_submerged
    85                 flow_area = (diameter/2)**2 * pi
    86                 outlet_culvert_depth = diameter
    87                 case = 'Inlet submerged'                   
    88 
    89 
    90 
    91             if delta_total_energy < inlet_specific_energy:
     121            Q = min(Q_inlet_unsubmerged, Q_inlet_submerged)
     122
     123            # THE LOWEST Value will Control Calcs From here
     124            # Calculate Critical Depth Based on the Adopted Flow as an Estimate
     125            dcrit1 = diameter/1.26*(Q/g**0.5*diameter**2.5)**(1/3.75)
     126            dcrit2 = diameter/0.95*(Q/g**0.5*diameter**2.5)**(1/1.95)
     127            # From Boyd Paper ESTIMATE of Dcrit has 2 criteria as
     128            if dcrit1/diameter  > 0.85:
     129                outlet_culvert_depth = dcrit2
     130            else:
     131                outlet_culvert_depth = dcrit1
     132            #outlet_culvert_depth = min(outlet_culvert_depth, diameter)
     133            # Now determine Hydraulic Radius Parameters Area & Wetted Perimeter
     134            if outlet_culvert_depth >= diameter:
     135                outlet_culvert_depth = diameter  # Once again the pipe is flowing full not partfull
     136                flow_area = (diameter/2)**2 * pi  # Cross sectional area of flow in the culvert
     137                perimeter = diameter * pi
     138                flow_width= diameter
     139                case = 'Inlet CTRL Outlet submerged Circular PIPE FULL'
     140                if local_debug =='true':
     141                    print 'Inlet CTRL Outlet submerged Circular PIPE FULL'
     142            else:
     143                #alpha = acos(1 - outlet_culvert_depth/diameter)    # Where did this Come From ????/
     144                alpha = acos(1-2*outlet_culvert_depth/diameter)*2
     145                #flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))        # Pipe is Running Partly Full at the INLET   WHRE did this Come From ?????
     146                flow_area = diameter**2/8*(alpha - sin(alpha))   # Equation from  GIECK 5th Ed. Pg. B3
     147                flow_width= diameter*sin(alpha/2.0)
     148                perimeter = alpha*diameter/2.0
     149                case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'
     150                if local_debug =='true':
     151                    print 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'
     152                    print 'Q Outlet Depth and ALPHA = ',Q,' ',outlet_culvert_depth,'  ',alpha
     153            if delta_total_energy < inlet_specific_energy:  #  OUTLET CONTROL !!!!
    92154                # Calculate flows for outlet control
    93155               
    94156                # Determine the depth at the outlet relative to the depth of flow in the Culvert
    95                 if outlet_depth > diameter:        # The Outlet is Submerged
     157                if outlet_depth > diameter:       # Outlet is submerged Assume the end of the Pipe is flowing FULL
    96158                    outlet_culvert_depth=diameter
    97159                    flow_area = (diameter/2)**2 * pi  # Cross sectional area of flow in the culvert
    98160                    perimeter = diameter * pi
     161                    flow_width= diameter
    99162                    case = 'Outlet submerged'
    100                 elif outlet_depth==0.0:
    101                     outlet_culvert_depth=inlet_depth # For purpose of calculation assume the outlet depth = the inlet depth
    102                     alpha = acos(1 - inlet_depth/diameter)
    103                     flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))
    104                     perimeter = alpha*diameter
    105 
    106                     case = 'Outlet depth is zero'                       
    107                 else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
    108                     outlet_culvert_depth=outlet_depth   # For purpose of calculation assume the outlet depth = the inlet depth
    109                     alpha = acos(1 - outlet_depth/diameter)
    110                     flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))
    111                     perimeter = alpha*diameter
    112                     case = 'Outlet is open channel flow'
    113 
    114                 hyd_rad = flow_area/perimeter
    115                
    116                 if log_filename is not None:
    117                     s = 'hydraulic radius at outlet = %f' %hyd_rad
    118                     log_to_file(log_filename, s)
    119 
    120                 # Outlet control velocity using tail water
    121                 culvert_velocity = sqrt(delta_total_energy/((sum_loss/2/g)+(manning**2*culvert_length)/hyd_rad**1.33333))
    122                 Q_outlet_tailwater = flow_area * culvert_velocity
    123                
    124                 if log_filename is not None:               
    125                     s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater
    126                     log_to_file(log_filename, s)
    127                    
    128                 Q = min(Q, Q_outlet_tailwater)
    129             else:
    130                 pass
     163                    if local_debug =='true':
     164                        print 'Outlet submerged'
     165                else:   # Culvert running PART FULL for PART OF ITS LENGTH   Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
     166                    # IF  outlet_depth < diameter
     167                    dcrit1 = diameter/1.26*(Q/g**0.5*diameter**2.5)**(1/3.75)
     168                    dcrit2 = diameter/0.95*(Q/g**0.5*diameter**2.5)**(1/1.95)
     169                    if dcrit1/diameter >0.85:
     170                        outlet_culvert_depth= dcrit2
     171                    else:
     172                        outlet_culvert_depth = dcrit1
     173                    if outlet_culvert_depth > diameter:
     174                        outlet_culvert_depth = diameter  # Once again the pipe is flowing full not partfull
     175                        flow_area = (diameter/2)**2 * pi  # Cross sectional area of flow in the culvert
     176                        perimeter = diameter * pi
     177                        flow_width= diameter
     178                        case = 'Outlet unsubmerged PIPE FULL'
     179                        if local_debug =='true':
     180                            print  'Outlet unsubmerged PIPE FULL'
     181                    else:
     182                        alpha = acos(1-2*outlet_culvert_depth/diameter)*2
     183                        flow_area = diameter**2/8*(alpha - sin(alpha))   # Equation from  GIECK 5th Ed. Pg. B3
     184                        flow_width= diameter*sin(alpha/2.0)
     185                        perimeter = alpha*diameter/2.0
     186                        case = 'Outlet is open channel flow we will for now assume critical depth'
     187                        if local_debug =='true':
     188                            print 'Q Outlet Depth and ALPHA = ',Q,' ',outlet_culvert_depth,'  ',alpha
     189                            print 'Outlet is open channel flow we will for now assume critical depth'
     190            if local_debug =='true':
     191                print 'FLOW AREA = ',flow_area
     192                print 'PERIMETER = ',perimeter
     193                print 'Q Interim = ',Q
     194            hyd_rad = flow_area/perimeter
     195
     196            if log_filename is not None:
     197                s = 'hydraulic radius at outlet = %f' %hyd_rad
     198                log_to_file(log_filename, s)
     199
     200            # Outlet control velocity using tail water
     201            if local_debug =='true':
     202                print 'GOT IT ALL CALCULATING Velocity'
     203                print 'HydRad = ',hyd_rad     
     204            culvert_velocity = sqrt(delta_total_energy/((sum_loss/2/g)+(manning**2*culvert_length)/hyd_rad**1.33333))
     205            Q_outlet_tailwater = flow_area * culvert_velocity
     206            if local_debug =='true':
     207                print 'VELOCITY = ',culvert_velocity
     208                print 'Outlet Ctrl Q = ',Q_outlet_tailwater
     209            if log_filename is not None:               
     210                s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater
     211                log_to_file(log_filename, s)
     212            Q = min(Q, Q_outlet_tailwater)
     213            if local_debug =='true':
     214                print ('%s,%.3f,%.3f' %('dcrit 1 , dcit2 =',dcrit1,dcrit2))
     215                print ('%s,%.3f,%.3f,%.3f' %('Q and Velocity and Depth=',Q,culvert_velocity,outlet_culvert_depth))
     216
     217        else:
     218            pass
    131219                #FIXME(Ole): What about inlet control?
    132 
    133 
    134         else:
    135             # Box culvert (rectangle or square)
     220        # ====  END OF CODE BLOCK FOR "IF" CIRCULAR PIPE
     221
     222        #else....
     223        if culvert_type =='box':
     224            if local_debug =='true':
     225                print 'BOX CULVERT'
     226            # Box culvert (rectangle or square)   ========================================================================================================================
    136227
    137228            # Calculate flows for inlet control
    138229            height = culvert_height
    139             width = culvert_width           
     230            width = culvert_width
     231            flow_width=culvert_width           
    140232           
    141233            Q_inlet_unsubmerged = 0.540*g**0.5*width*inlet_specific_energy**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
     
    146238                log_to_file(log_filename, s)
    147239
    148 
    149240            # FIXME(Ole): Are these functions really for inlet control?   
    150241            if Q_inlet_unsubmerged < Q_inlet_submerged:
    151242                Q = Q_inlet_unsubmerged
    152                 flow_area = width*inlet_depth
    153                 outlet_culvert_depth = inlet_depth
    154                 case = 'Inlet unsubmerged'
     243                dcrit = (Q**2/g/width**2)**0.333333
     244                if dcrit > height:
     245                    dcrit = height
     246                flow_area = width*dcrit
     247                outlet_culvert_depth = dcrit
     248                case = 'Inlet unsubmerged Box Acts as Weir'
    155249            else:   
    156250                Q = Q_inlet_submerged
    157251                flow_area = width*height
    158252                outlet_culvert_depth = height
    159                 case = 'Inlet submerged'                   
    160 
     253                case = 'Inlet submerged Box Acts as Orifice'                   
     254
     255            dcrit = (Q**2/g/width**2)**0.333333
     256
     257            outlet_culvert_depth = dcrit
     258            if outlet_culvert_depth > height:
     259                outlet_culvert_depth = height  # Once again the pipe is flowing full not partfull
     260                flow_area = width*height  # Cross sectional area of flow in the culvert
     261                perimeter = 2*(width+height)
     262                case = 'Inlet CTRL Outlet unsubmerged PIPE PART FULL'
     263            else:
     264                flow_area = width * outlet_culvert_depth
     265                perimeter = width+2*outlet_culvert_depth
     266                case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'
     267       
    161268            if delta_total_energy < inlet_specific_energy:
    162269                # Calculate flows for outlet control
     
    168275                    perimeter=2.0*(width+height)
    169276                    case = 'Outlet submerged'
    170                 elif outlet_depth==0.0:
    171                     outlet_culvert_depth=inlet_depth   # For purpose of calculation assume the outlet depth = the inlet depth
    172                     flow_area=width*inlet_depth
    173                     perimeter=(width+2.0*inlet_depth)
    174                     case = 'Outlet depth is zero'                       
    175277                else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
    176                     outlet_culvert_depth=outlet_depth
    177                     flow_area=width*outlet_depth
    178                     perimeter=(width+2.0*outlet_depth)
    179                     case = 'Outlet is open channel flow'
     278                    dcrit = (Q**2/g/width**2)**0.333333
     279                    outlet_culvert_depth=dcrit   # For purpose of calculation assume the outlet depth = Critical Depth
     280                    if outlet_culvert_depth > height:
     281                        outlet_culvert_depth=height
     282                        flow_area=width*height
     283                        perimeter=2.0*(width+height)
     284                        case = 'Outlet is Flowing Full'
     285                    else:
     286                        flow_area=width*outlet_culvert_depth
     287                        perimeter=(width+2.0*outlet_culvert_depth)
     288                        case = 'Outlet is open channel flow'
    180289
    181290                hyd_rad = flow_area/perimeter
    182                
     291
    183292                if log_filename is not None:                               
    184293                    s = 'hydraulic radius at outlet = %f' % hyd_rad
     
    196305                pass
    197306                #FIXME(Ole): What about inlet control?
     307        # ====  END OF CODE BLOCK FOR "IF" BOX
    198308
    199309               
    200         # Common code for circle and square geometries
     310        # Common code for circle and box geometries ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    201311        if log_filename is not None:
    202312            log_to_file(log_filename, 'Case: "%s"' % case)
     
    206316            log_to_file(log_filename, s)
    207317
    208            
    209         culv_froude=sqrt(Q**2*width/(g*flow_area**3))
    210        
     318        culv_froude=sqrt(Q**2*flow_width/(g*flow_area**3))
     319        if local_debug =='true':
     320            print 'FLOW AREA = ',flow_area
     321            print 'PERIMETER = ',perimeter
     322            print 'Q final = ',Q
     323            print 'FROUDE = ',culv_froude
    211324        if log_filename is not None:                           
    212325            s = 'Froude in Culvert = %f' % culv_froude
     
    216329        barrel_velocity = Q/(flow_area + velocity_protection/flow_area)
    217330
     331    # END CODE BLOCK for DEPTH  > Required depth for CULVERT Flow 
    218332
    219333    else: # inlet_depth < 0.01:
  • branches/numpy/anuga/geospatial_data/test_geospatial_data.py

    r7035 r7176  
    112112        assert num.allclose(V, [79.4, -7])
    113113
     114                # FIXME: use failUnlessRaises()
    114115        try:
    115116            V = G.get_attributes('hdnoatedu') #Invalid
  • branches/numpy/anuga/load_mesh/loadASCII.py

    r7035 r7176  
    700700                               ('num_of_segments', 'num_of_segment_ends'))
    701701        outfile.variables['segments'][:] = mesh['segments']
    702         if mesh['segment_tags'].shape[0] > 0:
     702        if (mesh['segment_tags'].shape[1] > 0):
    703703            outfile.createDimension('num_of_segment_tag_chars',
    704704                                    mesh['segment_tags'].shape[1])
  • branches/numpy/anuga/mesh_engine/test_generate_mesh.py

    r6410 r7176  
    412412                        'Failed!')
    413413       
    414         correct = num.array(segattlist)
     414        correct = num.array(segattlist, num.Int)
    415415        self.failUnless(num.allclose(data['generatedsegmentmarkerlist'].flat,
    416416                                     correct.flat),
  • branches/numpy/anuga/shallow_water/data_manager.py

    r7035 r7176  
    566566                # Define a zero vector of same size and type as A
    567567                # for use with momenta
    568                 null = num.zeros(num.size(A), A.dtype.char) #??#
     568                null = num.zeros(num.size(A), A.dtype.char)
    569569
    570570                # Get xmomentum where depth exceeds minimum_storable_height
     
    14331433    for name in infile.variables:
    14341434        var = infile.variables[name]
    1435         outfile.createVariable(name, var.dtype.char, var.dimensions)    #??#
     1435        outfile.createVariable(name, var.dtype.char, var.dimensions)
    14361436
    14371437    # Copy the static variables
     
    21932193    swwfile = basename_in + '.sww'
    21942194    demfile = basename_out + '.' + format
    2195     # Note the use of a .ers extension is optional (write_ermapper_grid will
    2196     # deal with either option
    21972195
    21982196    # Read sww file
     
    26972695    # Interpolate using quantity values
    26982696    if verbose: print 'Interpolating'
    2699     interpolated_values = interp.interpolate(q, data_points).flatten()    #????#
     2697    interpolated_values = interp.interpolate(q, data_points).flatten()
    27002698
    27012699    if verbose:
  • branches/numpy/anuga/shallow_water/shallow_water_domain.py

    r7035 r7176  
    14621462
    14631463
    1464 ##
    1465 # @brief A class for a 'field' boundary.
    1466 # @note Inherits from Boundary.
     1464class Inflow_boundary(Boundary):
     1465    """Apply given flow in m^3/s to boundary segment.
     1466    Depth and momentum is derived using Manning's formula.
     1467
     1468    Underlying domain must be specified when boundary is instantiated
     1469    """
     1470   
     1471    # FIXME (Ole): This is work in progress and definitely not finished.
     1472    # The associated test has been disabled
     1473
     1474    def __init__(self, domain=None, rate=0.0):
     1475        Boundary.__init__(self)
     1476
     1477        if domain is None:
     1478            msg = 'Domain must be specified for '
     1479            msg += 'Inflow boundary'
     1480            raise Exception, msg
     1481
     1482        self.domain = domain
     1483       
     1484        # FIXME(Ole): Allow rate to be time dependent as well
     1485        self.rate = rate
     1486        self.tag = None # Placeholder for tag associated with this object.
     1487
     1488    def __repr__(self):
     1489        return 'Inflow_boundary(%s)' %self.domain
     1490
     1491    def evaluate(self, vol_id, edge_id):
     1492        """Apply inflow rate at each edge of this boundary
     1493        """
     1494       
     1495        # First find all segments having the same tag is vol_id, edge_id
     1496        # This will be done the first time evaluate is called.
     1497        if self.tag is None:
     1498            boundary = self.domain.boundary
     1499            self.tag = boundary[(vol_id, edge_id)]       
     1500           
     1501            # Find total length of boundary with this tag
     1502            length = 0.0
     1503            for v_id, e_id in boundary:
     1504                if self.tag == boundary[(v_id, e_id)]:
     1505                    length += self.domain.mesh.get_edgelength(v_id, e_id)           
     1506
     1507            self.length = length
     1508            self.average_momentum = self.rate/length
     1509           
     1510           
     1511        # Average momentum has now been established across this boundary
     1512        # Compute momentum in the inward normal direction
     1513       
     1514        inward_normal = -self.domain.mesh.get_normal(vol_id, edge_id)       
     1515        xmomentum, ymomentum = self.average_momentum * inward_normal
     1516           
     1517        # Compute depth based on Manning's formula v = 1/n h^{2/3} sqrt(S)
     1518        # Where v is velocity, n is manning's coefficient, h is depth and S is the slope into the domain.
     1519        # Let mu be the momentum (vh), then this equation becomes: mu = 1/n h^{5/3} sqrt(S)
     1520        # from which we can isolate depth to get
     1521        # h = (mu n/sqrt(S) )^{3/5}
     1522       
     1523        slope = 0 # get gradient for this triangle dot normal
     1524       
     1525        # get manning coef from this triangle
     1526        friction = self.domain.get_quantity('friction').get_values(location='edges',
     1527                                                                   indices=[vol_id])[0]
     1528        mannings_n = friction[edge_id]
     1529
     1530        if slope > epsilon and mannings_n > epsilon:
     1531            depth = pow(self.average_momentum * mannings_n/math.sqrt(slope), 3.0/5)
     1532        else:
     1533            depth = 1.0
     1534           
     1535        # Elevation on this edge   
     1536       
     1537        z = self.domain.get_quantity('elevation').get_values(location='edges',
     1538                                                             indices=[vol_id])[0]
     1539        elevation = z[edge_id]
     1540           
     1541        # Assign conserved quantities and return
     1542        q = num.array([elevation + depth, xmomentum, ymomentum], num.Float)
     1543        return q
     1544
     1545
     1546       
     1547   
     1548           
     1549       
    14671550class Field_boundary(Boundary):
    14681551    """Set boundary from given field represented in an sww file containing
  • branches/numpy/anuga/utilities/test_polygon.py

    r7038 r7176  
    178178
    179179    def test_inside_polygon_main(self):
     180        """test_is_inside_polygon_quick
     181       
     182        Test fast version of of is_inside_polygon
     183        """
     184
    180185        # Simplest case: Polygon is the unit square
    181186        polygon = [[0,0], [1,0], [1,1], [0,1]]
Note: See TracChangeset for help on using the changeset viewer.