Changeset 7957


Ignore:
Timestamp:
Aug 20, 2010, 9:12:41 AM (9 years ago)
Author:
steve
Message:

update of culverts

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/anuga_core/source/anuga/structures/culvert_operator.py

    r7955 r7957  
    2121class Below_interval(Exception): pass
    2222class Above_interval(Exception): pass
    23 
    24 # FIXME(Ole): Take a good hard look at logging here
    25 
    26 
    27 # FIXME(Ole): Write in C and reuse this function by similar code
    28 # in interpolate.py
    29 def interpolate_linearly(x, xvec, yvec):
    30 
    31     msg = 'Input to function interpolate_linearly could not be converted '
    32     msg += 'to numerical scalar: x = %s' % str(x)
    33     try:
    34         x = float(x)
    35     except:
    36         raise Exception, msg
    37 
    38 
    39     # Check bounds
    40     if x < xvec[0]:
    41         msg = 'Value provided = %.2f, interpolation minimum = %.2f.'\
    42             % (x, xvec[0])
    43         raise Below_interval, msg
    44        
    45     if x > xvec[-1]:
    46         msg = 'Value provided = %.2f, interpolation maximum = %.2f.'\
    47             %(x, xvec[-1])
    48         raise Above_interval, msg       
    49        
    50        
    51     # Find appropriate slot within bounds           
    52     i = 0
    53     while x > xvec[i]: i += 1
    54 
    55    
    56     x0 = xvec[i-1]
    57     x1 = xvec[i]           
    58     alpha = (x - x0)/(x1 - x0)
    59            
    60     y0 = yvec[i-1]
    61     y1 = yvec[i]                       
    62     y = alpha*y1 + (1-alpha)*y0
    63    
    64     return y
    65            
    66 
    67            
    68 def read_culvert_description(culvert_description_filename):           
    69    
    70     # Read description file
    71     fid = open(culvert_description_filename)
    72    
    73     read_rating_curve_data = False       
    74     rating_curve = []
    75     for i, line in enumerate(fid.readlines()):
    76        
    77         if read_rating_curve_data is True:
    78             fields = line.split(',')
    79             head_difference = float(fields[0].strip())
    80             flow_rate = float(fields[1].strip())               
    81             barrel_velocity = float(fields[2].strip())
    82            
    83             rating_curve.append([head_difference, flow_rate, barrel_velocity])
    84        
    85         if i == 0:
    86             # Header
    87             continue
    88         if i == 1:
    89             # Metadata
    90             fields = line.split(',')
    91             label=fields[0].strip()
    92             type=fields[1].strip().lower()
    93             assert type in ['box', 'pipe']
    94            
    95             width=float(fields[2].strip())
    96             height=float(fields[3].strip())
    97             length=float(fields[4].strip())
    98             number_of_barrels=int(fields[5].strip())
    99             #fields[6] refers to losses
    100             description=fields[7].strip()               
    101                
    102         if line.strip() == '': continue # Skip blanks
    103 
    104         if line.startswith('Rating'):
    105             read_rating_curve_data = True
    106             # Flow data follows
    107                
    108     fid.close()
    109    
    110     return label, type, width, height, length, number_of_barrels, description, rating_curve
    111    
    11223
    11324   
     
    14455        self.end_points= [end_point0, end_point1]
    14556        self.enquiry_gap_factor=enquiry_gap_factor
     57
    14658        self.verbose=verbose
    147        
    14859        self.filename = None
    14960       
     
    15667
    15768        # Establish initial values at each enquiry point
     69        self.enquiry_quantity_values = []
    15870        dq = domain.quantities
    159         for i, opening in enumerate(self.openings):
     71        for i in [0,1]:
    16072            idx = self.enquiry_indices[i]
    16173            elevation = dq['elevation'].get_values(location='centroids',
     
    16375            stage = dq['stage'].get_values(location='centroids',
    16476                                           indices=[idx])[0]
    165             opening.elevation = elevation
    166             opening.stage = stage
    167             opening.depth = stage-elevation
    168 
    169            
    170                            
    171         # Determine initial pipe direction.
    172         # This may change dynamically based on the total energy difference     
    173         # Consequently, this may be superfluous
    174         delta_z = self.openings[0].elevation - self.openings[1].elevation
    175         if delta_z > 0.0:
    176             self.inlet = self.openings[0]
    177             self.outlet = self.openings[1]
    178         else:               
    179             self.outlet = self.openings[0]
    180             self.inlet = self.openings[1]       
    181        
    182        
    183         # Store basic geometry
    184         self.end_points = [end_point0, end_point1]
    185         self.vector = P['vector']
    186         self.length = P['length']; assert self.length > 0.0
    187         if culvert_description_filename is not None:
    188             if not num.allclose(self.length, length, rtol=1.0e-2, atol=1.0e-2):
    189                 msg = 'WARNING: barrel length specified in "%s" (%.2f m)'\
    190                     % (culvert_description_filename,
    191                        length)
    192                 msg += ' does not match distance between specified'
    193                 msg += ' end points (%.2f m)' %self.length
    194                 log.critical(msg)
    195        
    196         self.verbose = verbose
    197 
    198         # Circular index for flow averaging in culvert
    199         self.N = N = number_of_smoothing_steps
    200         self.Q_list = [0]*N
    201         self.i = i
    202        
    203         # For use with update_interval                       
    204         self.last_update = 0.0
    205         self.update_interval = update_interval
    206        
    207         # Create objects to update momentum (a bit crude at this stage). This is used with the momentum jet.
    208         xmom0 = General_forcing(domain, 'xmomentum',
    209                                 polygon=P['exchange_polygon0'])
    210 
    211         xmom1 = General_forcing(domain, 'xmomentum',
    212                                 polygon=P['exchange_polygon1'])
    213 
    214         ymom0 = General_forcing(domain, 'ymomentum',
    215                                 polygon=P['exchange_polygon0'])
    216 
    217         ymom1 = General_forcing(domain, 'ymomentum',
    218                                 polygon=P['exchange_polygon1'])
    219 
    220         self.opening_momentum = [[xmom0, ymom0], [xmom1, ymom1]]
    221 
    222 
    223 
    224         # Print some diagnostics to log if requested
    225         if self.log_filename is not None:
    226             s = 'Culvert Effective Length = %.2f m' %(self.length)
    227             log_to_file(self.log_filename, s)
    228    
    229             s = 'Culvert Direction is %s\n' %str(self.vector)
    230             log_to_file(self.log_filename, s)       
     77            depth = stage-elevation
     78           
     79            quantity_values = {'stage' : stage, 'elevation' : elevation, 'depth' : depth }
     80            self.enquiry_quantity_values.append(quantity_values)
     81       
     82       
     83        assert self.culvert_length > 0.0
    23184
    23285
     
    255108
    256109        # Calculate geometry
    257         x0, y0 = self.end_point0
    258         x1, y1 = self.end_point1
     110        x0, y0 = self.end_points[0]
     111        x1, y1 = self.end_points[1]
    259112
    260113        dx = x1-x0
     
    277130
    278131        self.exchange_polygons = []
    279 
    280         # Build exchange polygon and enquiry point for opening 0
    281         p0 = self.end_point0 + w
    282         p1 = self.end_point0 - w
    283         p2 = p1 - h
    284         p3 = p0 - h
    285         self.exchange_polygons.append(num.array([p0,p1,p2,p3]))
    286         self.enquiry_points= end_point0 - gap
    287 
    288 
    289         # Build exchange polygon and enquiry point for opening 1
    290         p0 = end_point1 + w
    291         p1 = end_point1 - w
    292         p2 = p1 + h
    293         p3 = p0 + h
    294         self.exchange_polygon1 = num.array([p0,p1,p2,p3])
    295         self.enquiry_point1 = end_point1 + gap
    296 
    297         # Check that enquiry point0 is outside exchange polygon0
    298         polygon = self.exchange_polygon0
    299         area = polygon_area(polygon)
     132        self.enquiry_points = []
     133
     134        # Build exchange polygon and enquiry points 0 and 1
     135        for i in [0,1]:
     136            p0 = self.end_points[i] + w
     137            p1 = self.end_point[i] - w
     138            p2 = p1 - h
     139            p3 = p0 - h
     140            self.exchange_polygons.append(num.array([p0,p1,p2,p3]))
     141            self.enquiry_points.append(end_point0 - gap)
     142
     143        self.polygon_areas = []
     144
     145        # Check that enquiry points are outside exchange polygons
     146        for i in [0,1]:
     147            polygon = self.exchange_polygons[i]
     148            # FIXME(SR) should use area of triangles associated with polygon
     149            area = polygon_area(polygon)
     150            self.polygon_areas.append(area)
    300151
    301152            msg = 'Polygon %s ' %(polygon)
     
    303154            assert area > 0.0, msg
    304155
    305             for key2 in ['enquiry_point0', 'enquiry_point1']:
    306                 point = culvert_polygons[key2]
    307                 msg = 'Enquiry point falls inside an enquiry point.'
     156            for j in [0,1]:
     157                point = self.enquiry_points[j]
     158                msg = 'Enquiry point falls inside a culvert polygon.'
    308159
    309160                assert not inside_polygon(point, polygon), msg
    310 
    311 
    312 
    313         for key1 in ['exchange_polygon0',
    314                      'exchange_polygon1']:
    315            
    316 
    317         # Return results
    318         self.culvert_polygons = culvert_polygons
    319 
    320 
    321161
    322162       
     
    794634
    795635
     636# FIXME(Ole): Write in C and reuse this function by similar code
     637# in interpolate.py
     638def interpolate_linearly(x, xvec, yvec):
     639
     640    msg = 'Input to function interpolate_linearly could not be converted '
     641    msg += 'to numerical scalar: x = %s' % str(x)
     642    try:
     643        x = float(x)
     644    except:
     645        raise Exception, msg
     646
     647
     648    # Check bounds
     649    if x < xvec[0]:
     650        msg = 'Value provided = %.2f, interpolation minimum = %.2f.'\
     651            % (x, xvec[0])
     652        raise Below_interval, msg
     653
     654    if x > xvec[-1]:
     655        msg = 'Value provided = %.2f, interpolation maximum = %.2f.'\
     656            %(x, xvec[-1])
     657        raise Above_interval, msg
     658
     659
     660    # Find appropriate slot within bounds
     661    i = 0
     662    while x > xvec[i]: i += 1
     663
     664
     665    x0 = xvec[i-1]
     666    x1 = xvec[i]
     667    alpha = (x - x0)/(x1 - x0)
     668
     669    y0 = yvec[i-1]
     670    y1 = yvec[i]
     671    y = alpha*y1 + (1-alpha)*y0
     672
     673    return y
     674
     675
     676
     677def read_culvert_description(culvert_description_filename):
     678
     679    # Read description file
     680    fid = open(culvert_description_filename)
     681
     682    read_rating_curve_data = False
     683    rating_curve = []
     684    for i, line in enumerate(fid.readlines()):
     685
     686        if read_rating_curve_data is True:
     687            fields = line.split(',')
     688            head_difference = float(fields[0].strip())
     689            flow_rate = float(fields[1].strip())
     690            barrel_velocity = float(fields[2].strip())
     691
     692            rating_curve.append([head_difference, flow_rate, barrel_velocity])
     693
     694        if i == 0:
     695            # Header
     696            continue
     697        if i == 1:
     698            # Metadata
     699            fields = line.split(',')
     700            label=fields[0].strip()
     701            type=fields[1].strip().lower()
     702            assert type in ['box', 'pipe']
     703
     704            width=float(fields[2].strip())
     705            height=float(fields[3].strip())
     706            length=float(fields[4].strip())
     707            number_of_barrels=int(fields[5].strip())
     708            #fields[6] refers to losses
     709            description=fields[7].strip()
     710
     711        if line.strip() == '': continue # Skip blanks
     712
     713        if line.startswith('Rating'):
     714            read_rating_curve_data = True
     715            # Flow data follows
     716
     717    fid.close()
     718
     719    return label, type, width, height, length, number_of_barrels, description, rating_curve
     720
    796721           
    797722
Note: See TracChangeset for help on using the changeset viewer.