Changeset 7955


Ignore:
Timestamp:
Aug 19, 2010, 4:41:10 PM (15 years ago)
Author:
steve
Message:

Changing culvert forcing term class to a culvert fractional step operator

File:
1 moved

Legend:

Unmodified
Added
Removed

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

    r7952 r7955  
    44from anuga.structures.culvert_polygons import create_culvert_polygons
    55from anuga.utilities.system_tools import log_to_file
    6 from anuga.geometry.polygon import inside_polygon, is_inside_polygon, plot_polygons
     6from anuga.geometry.polygon import inside_polygon, is_inside_polygon
     7from anuga.geometry.polygon import plot_polygons, polygon_area
     8
    79
    810from anuga.utilities.numerical_tools import mean
    … …  
    111113   
    112114
    113 class Culvert_flow_general:
    114     """Culvert flow - transfer water from one hole to another
    115    
    116     This version will work with either rating curve file or with culvert
    117     routine.
     115class Generic_box_culvert:
     116    """Culvert flow - transfer water from one rectngular box to another.
     117    Sets up the geometry of problem
     118   
     119    This is the base class for culverts. Inherit from this class (and overwrite
     120    compute_discharge method for specific subclasses)
    118121   
    119122    Input: Two points, pipe_size (either diameter or width, height),
    … …  
    123126    def __init__(self,
    124127                 domain,
    125                  culvert_description_filename=None,
    126                  culvert_routine=None,
    127128                 end_point0=None,
    128129                 end_point1=None,
    129                  enquiry_point0=None,
    130                  enquiry_point1=None,
    131                  type='box',
     130                 enquiry_gap_factor=0.2,
    132131                 width=None,
    133132                 height=None,
    134                  length=None,
    135                  number_of_barrels=1,
    136                  number_of_smoothing_steps=2000,
    137                  trigger_depth=0.01, # Depth below which no flow happens
    138                  manning=None,          # Mannings Roughness for Culvert
    139                  sum_loss=None,
    140                  use_velocity_head=False, # FIXME(Ole): Get rid of - always True
    141                  use_momentum_jet=False, # FIXME(Ole): Not yet implemented
    142                  label=None,
    143                  description=None,
    144                  update_interval=None,
    145                  log_file=False,
    146                  discharge_hydrograph=False,
    147133                 verbose=False):
    148134       
    149 
    150        
    151135        # Input check
    152136       
    153         if height is None: height = width
    154        
    155         assert number_of_barrels >= 1
    156         assert use_velocity_head is True or use_velocity_head is False
    157        
    158         #msg = 'Momentum jet not yet moved to general culvert'
    159         #assert use_momentum_jet is False, msg
    160         self.use_momentum_jet = use_momentum_jet
    161  
    162         self.culvert_routine = culvert_routine       
    163         self.culvert_description_filename = culvert_description_filename
    164         if culvert_description_filename is not None:
    165             label, type, width, height, length, number_of_barrels, description, rating_curve = read_culvert_description(culvert_description_filename)
    166             self.rating_curve = ensure_numeric(rating_curve)           
     137        if height is None:
     138            height = width
    167139
    168140        self.height = height
    169141        self.width = width
    170 
    171        
    172         self.domain = domain
    173         self.trigger_depth = trigger_depth       
    174                
    175         if manning is None:
    176             self.manning = 0.012   # Default roughness for pipe
    177        
    178         if sum_loss is None:
    179             self.sum_loss = 0.0
    180            
    181        
    182                        
    183         # Store culvert information
    184         self.label = label
    185         self.description = description
    186         self.culvert_type = type
    187         self.number_of_barrels = number_of_barrels
    188        
    189         # Store options       
    190         self.use_velocity_head = use_velocity_head
    191 
    192         if label is None: label = 'culvert_flow'
    193         label += '_' + str(id(self))
    194         self.label = label
    195        
    196         # File for storing discharge_hydrograph
    197         if discharge_hydrograph is True:
    198             self.timeseries_filename = label + '_timeseries.csv'
    199             fid = open(self.timeseries_filename, 'w')
    200             fid.write('time, discharge\n')
    201             fid.close()
    202 
    203         # Log file for storing general textual output
    204         if log_file is True:       
    205             self.log_filename = label + '.log'         
    206             log_to_file(self.log_filename, self.label)       
    207             log_to_file(self.log_filename, description)
    208             log_to_file(self.log_filename, self.culvert_type)       
    209         else:
    210             self.log_filename = None
    211 
     142       
     143        self.domain = domain
     144        self.end_points= [end_point0, end_point1]
     145        self.enquiry_gap_factor=enquiry_gap_factor
     146        self.verbose=verbose
     147       
     148        self.filename = None
     149       
    212150
    213151        # Create the fundamental culvert polygons from polygon
    214         P = create_culvert_polygons(end_point0,
    215                                     end_point1,
    216                                     width=width,   
    217                                     height=height,
    218                                     number_of_barrels=number_of_barrels)
    219         self.culvert_polygons = P
    220        
    221         # Select enquiry points
    222         if enquiry_point0 is None:
    223             enquiry_point0 = P['enquiry_point0']
    224            
    225         if enquiry_point1 is None:
    226             enquiry_point1 = P['enquiry_point1']           
    227            
    228         if verbose is True:
    229             pass
    230             #plot_polygons([[end_point0, end_point1],
    231             #               P['exchange_polygon0'],
    232             #               P['exchange_polygon1'],
    233             #               [enquiry_point0, 1.005*enquiry_point0],
    234             #               [enquiry_point1, 1.005*enquiry_point1]],
    235             #              figname='culvert_polygon_output')
    236 
    237            
    238            
    239         self.enquiry_points = [enquiry_point0, enquiry_point1]
    240         self.enquiry_indices = self.get_enquiry_indices()                 
     152        self.create_culvert_polygons()
     153        self.compute_enquiry_indices()
    241154        self.check_culvert_inside_domain()
    242        
    243                    
    244         # Create inflow object at each end of the culvert.
    245         self.openings = []
    246         self.openings.append(Inflow(domain,
    247                                     polygon=P['exchange_polygon0']))
    248         self.openings.append(Inflow(domain,
    249                                     polygon=P['exchange_polygon1']))
    250                                    
    251         # Assume two openings for now: Referred to as 0 and 1
    252         assert len(self.openings) == 2
     155
    253156
    254157        # Establish initial values at each enquiry point
    … …  
    326229            s = 'Culvert Direction is %s\n' %str(self.vector)
    327230            log_to_file(self.log_filename, s)       
    328        
    329        
    330        
    331        
     231
     232
     233    def set_store_hydrograph_discharge(self,filename=None):
     234
     235        if filename is None:
     236            self.filename = 'culvert_discharge_hydrograph'
     237        else:
     238            self.filename = filename
     239
     240        self.discharge_hydrograph = True
     241       
     242        self.timeseries_filename = self.filename + '_timeseries.csv'
     243        fid = open(self.timeseries_filename, 'w')
     244        fid.write('time, discharge\n')
     245        fid.close()
     246
     247
     248
     249
     250    def create_culvert_polygons(self):
     251        """Create polygons at the end of a culvert inlet and outlet.
     252        At either end two polygons will be created; one for the actual flow to pass through and one a little further away
     253        for enquiring the total energy at both ends of the culvert and transferring flow.
     254        """
     255
     256        # Calculate geometry
     257        x0, y0 = self.end_point0
     258        x1, y1 = self.end_point1
     259
     260        dx = x1-x0
     261        dy = y1-y0
     262
     263        self.culvert_vector = num.array([dx, dy])
     264        self.culvert_length = sqrt(num.sum(self.culvert_vector**2))
     265
     266
     267        # Unit direction vector and normal
     268        self.culvert_vector /= self.culvert_length                      # Unit vector in culvert direction
     269        self.culvert_normal = num.array([-dy, dx])/self.culvert_length  # Normal vector
     270
     271        # Short hands
     272        w = 0.5*width*normal # Perpendicular vector of 1/2 width
     273        h = height*vector    # Vector of length=height in the
     274                             # direction of the culvert
     275        gap = (1 + enquiry_gap_factor)*h
     276
     277
     278        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)
     300
     301            msg = 'Polygon %s ' %(polygon)
     302            msg += ' has area = %f' % area
     303            assert area > 0.0, msg
     304
     305            for key2 in ['enquiry_point0', 'enquiry_point1']:
     306                point = culvert_polygons[key2]
     307                msg = 'Enquiry point falls inside an enquiry point.'
     308
     309                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
     321
    332322       
    333323    def __call__(self, domain):
    … …  
    367357
    368358
    369     def get_enquiry_indices(self):
     359    def compute_enquiry_indices(self):
    370360        """Get indices for nearest centroids to self.enquiry_points
    371361        """
    … …  
    402392                raise Exception, msg
    403393       
    404         return enquiry_indices
     394        self.enquiry_indices = enquiry_indices
    405395
    406396       
    … …  
    806796           
    807797
    808 
    809 
    810                            
    811 # OBSOLETE (Except for momentum jet in Culvert_flow_energy)   
    812 class Culvert_flow_rating:
    813     """Culvert flow - transfer water from one hole to another
    814    
    815 
    816     Input: Two points, pipe_size (either diameter or width, height),
    817     mannings_rougness,
    818     inlet/outlet energy_loss_coefficients, internal_bend_coefficent,
    819     top-down_blockage_factor and bottom_up_blockage_factor
    820    
    821     """
    822 
    823     def __init__(self,
    824                  domain,
    825                  culvert_description_filename=None,
    826                  end_point0=None,
    827                  end_point1=None,
    828                  enquiry_point0=None,
    829                  enquiry_point1=None,
    830                  update_interval=None,
    831                  log_file=False,
    832                  discharge_hydrograph=False,
    833                  verbose=False):
    834        
    835 
    836        
    837         label, type, width, height, length, number_of_barrels, description, rating_curve = read_culvert_description(culvert_description_filename)
    838        
    839                
    840         # Store culvert information
    841         self.label = label
    842         self.description = description
    843         self.culvert_type = type
    844         self.rating_curve = ensure_numeric(rating_curve)
    845         self.number_of_barrels = number_of_barrels
    846 
    847         if label is None: label = 'culvert_flow'
    848         label += '_' + str(id(self))
    849         self.label = label
    850        
    851         # File for storing discharge_hydrograph
    852         if discharge_hydrograph is True:
    853             self.timeseries_filename = label + '_timeseries.csv'
    854             fid = open(self.timeseries_filename, 'w')
    855             fid.write('time, discharge\n')
    856             fid.close()
    857 
    858         # Log file for storing general textual output
    859         if log_file is True:       
    860             self.log_filename = label + '.log'         
    861             log_to_file(self.log_filename, self.label)       
    862             log_to_file(self.log_filename, description)
    863             log_to_file(self.log_filename, self.culvert_type)       
    864 
    865 
    866         # Create the fundamental culvert polygons from POLYGON
    867         #if self.culvert_type == 'circle':
    868         #    # Redefine width and height for use with create_culvert_polygons
    869         #    width = height = diameter
    870        
    871         P = create_culvert_polygons(end_point0,
    872                                     end_point1,
    873                                     width=width,   
    874                                     height=height,
    875                                     number_of_barrels=number_of_barrels)
    876        
    877         # Select enquiry points
    878         if enquiry_point0 is None:
    879             enquiry_point0 = P['enquiry_point0']
    880            
    881         if enquiry_point1 is None:
    882             enquiry_point1 = P['enquiry_point1']           
    883            
    884         if verbose is True:
    885             pass
    886             #plot_polygons([[end_point0, end_point1],
    887             #               P['exchange_polygon0'],
    888             #               P['exchange_polygon1'],
    889             #               [enquiry_point0, 1.005*enquiry_point0],
    890             #               [enquiry_point1, 1.005*enquiry_point1]],                           
    891             #              figname='culvert_polygon_output')
    892 
    893            
    894            
    895         self.enquiry_points = [enquiry_point0, enquiry_point1]                           
    896 
    897         self.enquiry_indices = []                 
    898         for point in self.enquiry_points:
    899             # Find nearest centroid
    900             N = len(domain)   
    901             points = domain.get_centroid_coordinates(absolute=True)
    902 
    903             # Calculate indices in exchange area for this forcing term
    904            
    905             triangle_id = min_dist = sys.maxint
    906             for k in range(N):
    907                 x, y = points[k,:] # Centroid
    908 
    909                 c = point                               
    910                 distance = (x-c[0])**2+(y-c[1])**2
    911                 if distance < min_dist:
    912                     min_dist = distance
    913                     triangle_id = k
    914 
    915                    
    916             if triangle_id < sys.maxint:
    917                 msg = 'found triangle with centroid (%f, %f)'\
    918                     %tuple(points[triangle_id, :])
    919                 msg += ' for point (%f, %f)' %tuple(point)
    920                
    921                 self.enquiry_indices.append(triangle_id)
    922             else:
    923                 msg = 'Triangle not found for point (%f, %f)' %point
    924                 raise Exception, msg
    925        
    926                          
    927 
    928         # Check that all polygons lie within the mesh.
    929         bounding_polygon = domain.get_boundary_polygon()
    930         for key in P.keys():
    931             if key in ['exchange_polygon0',
    932                        'exchange_polygon1']:
    933                 for point in list(P[key]) + self.enquiry_points:
    934                     msg = 'Point %s in polygon %s for culvert %s did not'\
    935                         %(str(point), key, self.label)
    936                     msg += 'fall within the domain boundary.'
    937                     assert is_inside_polygon(point, bounding_polygon), msg
    938        
    939                    
    940         # Create inflow object at each end of the culvert.
    941         self.openings = []
    942         self.openings.append(Inflow(domain,
    943                                     polygon=P['exchange_polygon0']))
    944 
    945         self.openings.append(Inflow(domain,
    946                                     polygon=P['exchange_polygon1']))                                   
    947 
    948 
    949 
    950         dq = domain.quantities                                           
    951         for i, opening in enumerate(self.openings):                           
    952             elevation = dq['elevation'].get_values(location='centroids',
    953                                                    indices=[self.enquiry_indices[i]])           
    954             opening.elevation = elevation
    955             opening.stage = elevation # Simple assumption that culvert is dry initially
    956 
    957         # Assume two openings for now: Referred to as 0 and 1
    958         assert len(self.openings) == 2
    959                            
    960         # Determine pipe direction     
    961         self.delta_z = delta_z = self.openings[0].elevation - self.openings[1].elevation
    962         if delta_z > 0.0:
    963             self.inlet = self.openings[0]
    964             self.outlet = self.openings[1]
    965         else:               
    966             self.outlet = self.openings[0]
    967             self.inlet = self.openings[1]       
    968        
    969        
    970         # Store basic geometry
    971         self.end_points = [end_point0, end_point1]
    972         self.vector = P['vector']
    973         self.length = P['length']; assert self.length > 0.0
    974         if not num.allclose(self.length, length, rtol=1.0e-2, atol=1.0e-2):
    975             msg = 'WARNING: barrel length specified in "%s" (%.2f m)' %(culvert_description_filename, length)
    976             msg += ' does not match distance between specified'
    977             msg += ' end points (%.2f m)' %self.length
    978             log.critical(msg)
    979        
    980         self.verbose = verbose
    981         self.last_update = 0.0 # For use with update_interval       
    982         self.last_time = 0.0               
    983         self.update_interval = update_interval
    984        
    985 
    986         # Print something
    987         if hasattr(self, 'log_filename'):
    988             s = 'Culvert Effective Length = %.2f m' %(self.length)
    989             log_to_file(self.log_filename, s)
    990    
    991             s = 'Culvert Direction is %s\n' %str(self.vector)
    992             log_to_file(self.log_filename, s)       
    993        
    994        
    995        
    996        
    997        
    998     def __call__(self, domain):
    999 
    1000         # Time stuff
    1001         time = domain.get_time()
    1002        
    1003        
    1004         update = False
    1005         if self.update_interval is None:
    1006             update = True
    1007             delta_t = domain.timestep # Next timestep has been computed in domain.py
    1008         else:   
    1009             if time - self.last_update > self.update_interval or time == 0.0:
    1010                 update = True
    1011             delta_t = self.update_interval
    1012            
    1013         s = '\nTime = %.2f, delta_t = %f' %(time, delta_t)
    1014         if hasattr(self, 'log_filename'):           
    1015             log_to_file(self.log_filename, s)
    1016                
    1017                                
    1018         if update is True:
    1019             self.last_update = time
    1020        
    1021             dq = domain.quantities
    1022                        
    1023             # Get average water depths at each opening       
    1024             openings = self.openings   # There are two Opening [0] and [1]
    1025             for i, opening in enumerate(openings):
    1026                
    1027                 # Compute mean values of selected quantitites in the
    1028                 # enquiry area in front of the culvert
    1029                
    1030                 stage = dq['stage'].get_values(location='centroids',
    1031                                                indices=[self.enquiry_indices[i]])
    1032                
    1033                 # Store current average stage and depth with each opening object
    1034                 opening.depth = stage - opening.elevation
    1035                 opening.stage = stage
    1036 
    1037                
    1038 
    1039             #################  End of the FOR loop ################################################
    1040 
    1041             # We now need to deal with each opening individually
    1042                
    1043             # Determine flow direction based on total energy difference
    1044 
    1045             delta_w = self.inlet.stage - self.outlet.stage
    1046            
    1047             if hasattr(self, 'log_filename'):
    1048                 s = 'Time=%.2f, inlet stage = %.2f, outlet stage = %.2f' %(time,
    1049                                                                            self.inlet.stage,
    1050                                                                            self.outlet.stage)
    1051                 log_to_file(self.log_filename, s)
    1052 
    1053 
    1054             if self.inlet.depth <= 0.01:
    1055                 Q = 0.0
    1056             else:
    1057                 # Calculate discharge for one barrel and set inlet.rate and outlet.rate
    1058                
    1059                 try:
    1060                     Q = interpolate_linearly(delta_w, self.rating_curve[:,0], self.rating_curve[:,1])
    1061                 except Below_interval, e:
    1062                     Q = self.rating_curve[0,1]             
    1063                     msg = '%.2fs: Delta head smaller than rating curve minimum: ' %time
    1064                     msg += str(e)
    1065                     msg += '\n        I will use minimum discharge %.2f m^3/s for culvert "%s"'\
    1066                         %(Q, self.label)
    1067                     if hasattr(self, 'log_filename'):                   
    1068                         log_to_file(self.log_filename, msg)
    1069                 except Above_interval, e:
    1070                     Q = self.rating_curve[-1,1]         
    1071                     msg = '%.2fs: Delta head greater than rating curve maximum: ' %time
    1072                     msg += str(e)
    1073                     msg += '\n        I will use maximum discharge %.2f m^3/s for culvert "%s"'\
    1074                         %(Q, self.label)
    1075                     if hasattr(self, 'log_filename'):                   
    1076                         log_to_file(self.log_filename, msg)                   
    1077 
    1078                
    1079                
    1080            
    1081             # Adjust discharge for multiple barrels
    1082             Q *= self.number_of_barrels
    1083            
    1084 
    1085             # Adjust Q downwards depending on available water at inlet
    1086             stage = self.inlet.get_quantity_values(quantity_name='stage')
    1087             elevation = self.inlet.get_quantity_values(quantity_name='elevation')
    1088             depth = stage-elevation
    1089            
    1090            
    1091             V = 0
    1092             for i, d in enumerate(depth):
    1093                 V += d * domain.areas[i]
    1094            
    1095             dt = delta_t           
    1096             if Q*dt > V:
    1097            
    1098                 Q_reduced = 0.9*V/dt # Reduce with safety factor
    1099                
    1100                 msg = '%.2fs: Computed extraction for this time interval (Q*dt) is ' % time
    1101                 msg += 'greater than current volume (V) at inlet.\n'
    1102                 msg += ' Q will be reduced from %.2f m^3/s to %.2f m^3/s.' % (Q, Q_reduced)
    1103                
    1104                 if self.verbose is True:
    1105                     log.critical(msg)
    1106                 if hasattr(self, 'log_filename'):                   
    1107                     log_to_file(self.log_filename, msg)                                       
    1108                
    1109                 Q = Q_reduced
    1110        
    1111             self.inlet.rate = -Q
    1112             self.outlet.rate = Q
    1113 
    1114             # Log timeseries to file
    1115             try:
    1116                 fid = open(self.timeseries_filename, 'a')       
    1117             except:
    1118                 pass
    1119             else:   
    1120                 fid.write('%.2f, %.2f\n' %(time, Q))
    1121                 fid.close()
    1122 
    1123             # Store value of time
    1124             self.last_time = time
    1125 
    1126            
    1127    
    1128         # Execute flow term for each opening
    1129         # This is where Inflow objects are evaluated using the last rate that has been calculated
    1130         #
    1131         # This will take place at every internal timestep and update the domain
    1132         for opening in self.openings:
    1133             opening(domain)
    1134 
    1135 
    1136 
    1137        
    1138        
    1139 
    1140 class Culvert_flow_energy:
    1141     """Culvert flow - transfer water from one hole to another
    1142    
    1143     Using Momentum as Calculated by Culvert Flow !!
    1144     Could be Several Methods Investigated to do This !!!
    1145 
    1146     2008_May_08
    1147     To Ole:
    1148     OK so here we need to get the Polygon Creating code to create
    1149     polygons for the culvert Based on
    1150     the two input Points (X0,Y0) and (X1,Y1) - need to be passed
    1151     to create polygon
    1152 
    1153     The two centers are now passed on to create_polygon.
    1154    
    1155 
    1156     Input: Two points, pipe_size (either diameter or width, height),
    1157     mannings_rougness,
    1158     inlet/outlet energy_loss_coefficients, internal_bend_coefficent,
    1159     top-down_blockage_factor and bottom_up_blockage_factor
    1160    
    1161    
    1162     And the Delta H enquiry should be change from Openings in line 412
    1163     to the enquiry Polygons infront of the culvert
    1164     At the moment this script uses only Depth, later we can change it to
    1165     Energy...
    1166 
    1167     Once we have Delta H can calculate a Flow Rate and from Flow Rate
    1168     an Outlet Velocity
    1169     The Outlet Velocity x Outlet Depth = Momentum to be applied at the Outlet...
    1170        
    1171     Invert levels are optional. If left out they will default to the
    1172     elevation at the opening.
    1173        
    1174     """
    1175 
    1176     def __init__(self,
    1177                  domain,
    1178                  label=None,
    1179                  description=None,
    1180                  end_point0=None,
    1181                  end_point1=None,
    1182                  width=None,
    1183                  height=None,
    1184                  diameter=None,
    1185                  manning=None,          # Mannings Roughness for Culvert
    1186                  invert_level0=None,    # Invert level at opening 0
    1187                  invert_level1=None,    # Invert level at opening 1
    1188                  loss_exit=None,
    1189                  loss_entry=None,
    1190                  loss_bend=None,
    1191                  loss_special=None,
    1192                  blockage_topdwn=None,
    1193                  blockage_bottup=None,
    1194                  culvert_routine=None,
    1195                  number_of_barrels=1,
    1196                  enquiry_point0=None,
    1197                  enquiry_point1=None,
    1198                  update_interval=None,
    1199                  verbose=False):
    1200        
    1201         # Input check
    1202         if diameter is not None:
    1203             self.culvert_type = 'circle'
    1204             self.diameter = diameter
    1205             if height is not None or width is not None:
    1206                 msg = 'Either diameter or width&height must be specified, '
    1207                 msg += 'but not both.'
    1208                 raise Exception, msg
    1209         else:
    1210             if height is not None:
    1211                 if width is None:
    1212                     self.culvert_type = 'square'                               
    1213                     width = height
    1214                 else:
    1215                     self.culvert_type = 'rectangle'
    1216             elif width is not None:
    1217                 if height is None:
    1218                     self.culvert_type = 'square'                               
    1219                     height = width
    1220             else:
    1221                 msg = 'Either diameter or width&height must be specified.'
    1222                 raise Exception, msg               
    1223                
    1224             if height == width:
    1225                 self.culvert_type = 'square'                                               
    1226                
    1227             self.height = height
    1228             self.width = width
    1229 
    1230            
    1231         assert self.culvert_type in ['circle', 'square', 'rectangle']
    1232        
    1233         assert number_of_barrels >= 1
    1234         self.number_of_barrels = number_of_barrels
    1235        
    1236        
    1237         # Set defaults
    1238         if manning is None: manning = 0.012   # Default roughness for pipe
    1239         if loss_exit is None: loss_exit = 1.00
    1240         if loss_entry is None: loss_entry = 0.50
    1241         if loss_bend is None: loss_bend=0.00
    1242         if loss_special is None: loss_special=0.00
    1243         if blockage_topdwn is None: blockage_topdwn=0.00
    1244         if blockage_bottup is None: blockage_bottup=0.00
    1245         if culvert_routine is None:
    1246             culvert_routine=boyd_generalised_culvert_model
    1247            
    1248         if label is None: label = 'culvert_flow'
    1249         label += '_' + str(id(self))
    1250         self.label = label
    1251        
    1252         # File for storing culvert quantities
    1253         self.timeseries_filename = label + '_timeseries.csv'
    1254         fid = open(self.timeseries_filename, 'w')
    1255         fid.write('time, E0, E1, Velocity, Discharge\n')
    1256         fid.close()
    1257 
    1258         # Log file for storing general textual output
    1259         self.log_filename = label + '.log'         
    1260         log_to_file(self.log_filename, self.label)       
    1261         log_to_file(self.log_filename, description)
    1262         log_to_file(self.log_filename, self.culvert_type)       
    1263 
    1264 
    1265         # Create the fundamental culvert polygons from POLYGON
    1266         if self.culvert_type == 'circle':
    1267             # Redefine width and height for use with create_culvert_polygons
    1268             width = height = diameter
    1269        
    1270         P = create_culvert_polygons(end_point0,
    1271                                     end_point1,
    1272                                     width=width,   
    1273                                     height=height,
    1274                                     number_of_barrels=number_of_barrels)
    1275        
    1276         # Select enquiry points
    1277         if enquiry_point0 is None:
    1278             enquiry_point0 = P['enquiry_point0']
    1279            
    1280         if enquiry_point1 is None:
    1281             enquiry_point1 = P['enquiry_point1']           
    1282            
    1283         if verbose is True:
    1284             pass
    1285             #plot_polygons([[end_point0, end_point1],
    1286             #               P['exchange_polygon0'],
    1287             #               P['exchange_polygon1'],
    1288             #               [enquiry_point0, 1.005*enquiry_point0],
    1289             #               [enquiry_point1, 1.005*enquiry_point1]],
    1290             #              figname='culvert_polygon_output')
    1291 
    1292 
    1293         self.enquiry_points = [enquiry_point0, enquiry_point1]                           
    1294        
    1295        
    1296         self.enquiry_indices = []                 
    1297         for point in self.enquiry_points:
    1298             # Find nearest centroid
    1299             N = len(domain)   
    1300             points = domain.get_centroid_coordinates(absolute=True)
    1301 
    1302             # Calculate indices in exchange area for this forcing term
    1303            
    1304             triangle_id = min_dist = sys.maxint
    1305             for k in range(N):
    1306                 x, y = points[k,:] # Centroid
    1307 
    1308                 c = point                               
    1309                 distance = (x-c[0])**2+(y-c[1])**2
    1310                 if distance < min_dist:
    1311                     min_dist = distance
    1312                     triangle_id = k
    1313 
    1314                    
    1315             if triangle_id < sys.maxint:
    1316                 msg = 'found triangle with centroid (%f, %f)'\
    1317                     %tuple(points[triangle_id, :])
    1318                 msg += ' for point (%f, %f)' %tuple(point)
    1319                
    1320                 self.enquiry_indices.append(triangle_id)
    1321             else:
    1322                 msg = 'Triangle not found for point (%f, %f)' %point
    1323                 raise Exception, msg
    1324        
    1325                          
    1326 
    1327            
    1328            
    1329 
    1330         # Check that all polygons lie within the mesh.
    1331         bounding_polygon = domain.get_boundary_polygon()
    1332         for key in P.keys():
    1333             if key in ['exchange_polygon0',
    1334                        'exchange_polygon1']:
    1335                 for point in P[key]:
    1336                
    1337                     msg = 'Point %s in polygon %s for culvert %s did not'\
    1338                         %(str(point), key, self.label)
    1339                     msg += 'fall within the domain boundary.'
    1340                     assert is_inside_polygon(point, bounding_polygon), msg
    1341        
    1342 
    1343         # Create inflow object at each end of the culvert.
    1344         self.openings = []
    1345         self.openings.append(Inflow(domain,
    1346                                     polygon=P['exchange_polygon0']))
    1347 
    1348         self.openings.append(Inflow(domain,
    1349                                     polygon=P['exchange_polygon1']))                                   
    1350 
    1351 
    1352         # Assume two openings for now: Referred to as 0 and 1
    1353         assert len(self.openings) == 2
    1354        
    1355         # Store basic geometry
    1356         self.end_points = [end_point0, end_point1]
    1357         self.invert_levels = [invert_level0, invert_level1]               
    1358         #self.enquiry_polygons = [P['enquiry_polygon0'], P['enquiry_polygon1']]
    1359         #self.enquiry_polylines = [P['enquiry_polygon0'][:2],
    1360         #                          P['enquiry_polygon1'][:2]]
    1361         self.vector = P['vector']
    1362         self.length = P['length']; assert self.length > 0.0
    1363         self.verbose = verbose
    1364         self.last_time = 0.0       
    1365         self.last_update = 0.0 # For use with update_interval       
    1366         self.update_interval = update_interval
    1367        
    1368 
    1369         # Store hydraulic parameters
    1370         self.manning = manning
    1371         self.loss_exit = loss_exit
    1372         self.loss_entry = loss_entry
    1373         self.loss_bend = loss_bend
    1374         self.loss_special = loss_special
    1375         self.sum_loss = loss_exit + loss_entry + loss_bend + loss_special
    1376         self.blockage_topdwn = blockage_topdwn
    1377         self.blockage_bottup = blockage_bottup
    1378 
    1379         # Store culvert routine
    1380         self.culvert_routine = culvert_routine
    1381 
    1382        
    1383         # Create objects to update momentum (a bit crude at this stage)
    1384         xmom0 = General_forcing(domain, 'xmomentum',
    1385                                 polygon=P['exchange_polygon0'])
    1386 
    1387         xmom1 = General_forcing(domain, 'xmomentum',
    1388                                 polygon=P['exchange_polygon1'])
    1389 
    1390         ymom0 = General_forcing(domain, 'ymomentum',
    1391                                 polygon=P['exchange_polygon0'])
    1392 
    1393         ymom1 = General_forcing(domain, 'ymomentum',
    1394                                 polygon=P['exchange_polygon1'])
    1395 
    1396         self.opening_momentum = [ [xmom0, ymom0], [xmom1, ymom1] ]
    1397        
    1398 
    1399         # Print something
    1400         s = 'Culvert Effective Length = %.2f m' %(self.length)
    1401         log_to_file(self.log_filename, s)
    1402    
    1403         s = 'Culvert Direction is %s\n' %str(self.vector)
    1404         log_to_file(self.log_filename, s)       
    1405        
    1406        
    1407     def __call__(self, domain):
    1408 
    1409         log_filename = self.log_filename
    1410          
    1411         # Time stuff
    1412         time = domain.get_time()
    1413        
    1414         # Short hand
    1415         dq = domain.quantities
    1416                
    1417 
    1418         update = False
    1419         if self.update_interval is None:
    1420             update = True
    1421             delta_t = domain.timestep # Next timestep has been computed in domain.py
    1422         else:   
    1423             if time - self.last_update > self.update_interval or time == 0.0:
    1424                 update = True
    1425             delta_t = self.update_interval
    1426            
    1427         s = '\nTime = %.2f, delta_t = %f' %(time, delta_t)
    1428         if hasattr(self, 'log_filename'):           
    1429             log_to_file(log_filename, s)
    1430                
    1431                                
    1432         if update is True:
    1433             self.last_update = time
    1434                        
    1435             msg = 'Time did not advance'
    1436             if time > 0.0: assert delta_t > 0.0, msg
    1437 
    1438 
    1439             # Get average water depths at each opening       
    1440             openings = self.openings   # There are two Opening [0] and [1]
    1441             for i, opening in enumerate(openings):
    1442                
    1443                 # Compute mean values of selected quantitites in the
    1444                 # exchange area in front of the culvert
    1445                      
    1446                 stage = opening.get_quantity_values(quantity_name='stage')
    1447                 w = mean(stage) # Average stage
    1448 
    1449                 # Use invert level instead of elevation if specified
    1450                 invert_level = self.invert_levels[i]
    1451                 if invert_level is not None:
    1452                     z = invert_level
    1453                 else:
    1454                     elevation = opening.get_quantity_values(quantity_name='elevation')
    1455                     z = mean(elevation) # Average elevation
    1456 
    1457                 # Estimated depth above the culvert inlet
    1458                 d = w - z  # Used for calculations involving depth
    1459                 if d < 0.0:
    1460                     # This is possible since w and z are taken at different locations
    1461                     #msg = 'D < 0.0: %f' %d
    1462                     #raise msg
    1463                     d = 0.0
    1464                
    1465 
    1466                 # Ratio of depth to culvert height.
    1467                 # If ratio > 1 then culvert is running full
    1468                 if self.culvert_type == 'circle':
    1469                     ratio = d/self.diameter
    1470                 else:   
    1471                     ratio = d/self.height 
    1472                 opening.ratio = ratio
    1473                    
    1474                    
    1475                 # Average measures of energy in front of this opening
    1476                
    1477                 id = [self.enquiry_indices[i]]
    1478                 stage = dq['stage'].get_values(location='centroids',
    1479                                                indices=id)
    1480                 elevation = dq['elevation'].get_values(location='centroids',
    1481                                                        indices=id)                                               
    1482                 xmomentum = dq['xmomentum'].get_values(location='centroids',
    1483                                                        indices=id)                                               
    1484                 ymomentum = dq['xmomentum'].get_values(location='centroids',
    1485                                                        indices=id)                                                                                             
    1486                 depth = stage-elevation
    1487                 if depth > 0.0:
    1488                     u = xmomentum/(depth + velocity_protection/depth)
    1489                     v = ymomentum/(depth + velocity_protection/depth)
    1490                 else:
    1491                     u = v = 0.0
    1492 
    1493                    
    1494                 opening.total_energy = 0.5*(u*u + v*v)/g + stage
    1495 
    1496                 # Store current average stage and depth with each opening object
    1497                 opening.depth = d
    1498                 opening.depth_trigger = d # Use this for now
    1499                 opening.stage = w
    1500                 opening.elevation = z
    1501                
    1502 
    1503             #################  End of the FOR loop ################################################
    1504 
    1505             # We now need to deal with each opening individually
    1506                
    1507             # Determine flow direction based on total energy difference
    1508             delta_Et = openings[0].total_energy - openings[1].total_energy
    1509 
    1510             if delta_Et > 0:
    1511                 inlet = openings[0]
    1512                 outlet = openings[1]
    1513 
    1514                 inlet.momentum = self.opening_momentum[0]
    1515                 outlet.momentum = self.opening_momentum[1]
    1516 
    1517             else:
    1518                 inlet = openings[1]
    1519                 outlet = openings[0]
    1520 
    1521                 inlet.momentum = self.opening_momentum[1]
    1522                 outlet.momentum = self.opening_momentum[0]
    1523                
    1524                 delta_Et = -delta_Et
    1525 
    1526             self.inlet = inlet
    1527             self.outlet = outlet
    1528                
    1529             msg = 'Total energy difference is negative'
    1530             assert delta_Et > 0.0, msg
    1531 
    1532             delta_h = inlet.stage - outlet.stage
    1533             delta_z = inlet.elevation - outlet.elevation
    1534             culvert_slope = (delta_z/self.length)
    1535 
    1536             if culvert_slope < 0.0:
    1537                 # Adverse gradient - flow is running uphill
    1538                 # Flow will be purely controlled by uphill outlet face
    1539                 if self.verbose is True:
    1540                     log.critical('WARNING: Flow is running uphill. Watch Out! '
    1541                                  'inlet.elevation=%s, outlet.elevation%s'
    1542                                  % (str(inlet.elevation), str(outlet.elevation)))
    1543 
    1544 
    1545             s = 'Delta total energy = %.3f' %(delta_Et)
    1546             log_to_file(log_filename, s)
    1547 
    1548 
    1549             # Calculate discharge for one barrel and set inlet.rate and outlet.rate
    1550             Q, barrel_velocity, culvert_outlet_depth = self.culvert_routine(self, inlet, outlet, delta_Et, g)
    1551        
    1552             # Adjust discharge for multiple barrels
    1553             Q *= self.number_of_barrels
    1554 
    1555             # Compute barrel momentum
    1556             barrel_momentum = barrel_velocity*culvert_outlet_depth
    1557                    
    1558             s = 'Barrel velocity = %f' %barrel_velocity
    1559             log_to_file(log_filename, s)
    1560 
    1561             # Compute momentum vector at outlet
    1562             outlet_mom_x, outlet_mom_y = self.vector * barrel_momentum
    1563                
    1564             s = 'Directional momentum = (%f, %f)' %(outlet_mom_x, outlet_mom_y)
    1565             log_to_file(log_filename, s)
    1566 
    1567             # Log timeseries to file
    1568             fid = open(self.timeseries_filename, 'a')       
    1569             fid.write('%f, %f, %f, %f, %f\n'\
    1570                           %(time,
    1571                             openings[0].total_energy,
    1572                             openings[1].total_energy,
    1573                             barrel_velocity,
    1574                             Q))
    1575             fid.close()
    1576 
    1577             # Update momentum       
    1578             if delta_t > 0.0:
    1579                 xmomentum_rate = outlet_mom_x - outlet.momentum[0].value
    1580                 xmomentum_rate /= delta_t
    1581                    
    1582                 ymomentum_rate = outlet_mom_y - outlet.momentum[1].value
    1583                 ymomentum_rate /= delta_t
    1584                        
    1585                 s = 'X Y MOM_RATE = (%f, %f) ' %(xmomentum_rate, ymomentum_rate)
    1586                 log_to_file(log_filename, s)                   
    1587             else:
    1588                 xmomentum_rate = ymomentum_rate = 0.0
    1589 
    1590 
    1591             # Set momentum rates for outlet jet
    1592             outlet.momentum[0].rate = xmomentum_rate
    1593             outlet.momentum[1].rate = ymomentum_rate
    1594 
    1595             # Remember this value for next step (IMPORTANT)
    1596             outlet.momentum[0].value = outlet_mom_x
    1597             outlet.momentum[1].value = outlet_mom_y                   
    1598 
    1599             if int(domain.time*100) % 100 == 0:
    1600                 s = 'T=%.5f, Culvert Discharge = %.3f f'\
    1601                     %(time, Q)
    1602                 s += ' Depth= %0.3f  Momentum = (%0.3f, %0.3f)'\
    1603                      %(culvert_outlet_depth, outlet_mom_x,outlet_mom_y)
    1604                 s += ' Momentum rate: (%.4f, %.4f)'\
    1605                      %(xmomentum_rate, ymomentum_rate)                   
    1606                 s+='Outlet Vel= %.3f'\
    1607                     %(barrel_velocity)
    1608                 log_to_file(log_filename, s)
    1609            
    1610             # Store value of time
    1611             self.last_time = time
    1612                
    1613 
    1614 
    1615         # Execute flow term for each opening
    1616         # This is where Inflow objects are evaluated and update the domain
    1617         for opening in self.openings:
    1618             opening(domain)
    1619 
    1620         # Execute momentum terms
    1621         # This is where Inflow objects are evaluated and update the domain
    1622         self.outlet.momentum[0](domain)
    1623         self.outlet.momentum[1](domain)       
    1624            
    1625 
    1626 
    1627798Culvert_flow = Culvert_flow_general       
Note: See TracChangeset for help on using the changeset viewer.