Changeset 7955
 Timestamp:
 Aug 19, 2010, 4:41:10 PM (9 years ago)
 File:

 1 moved
Legend:
 Unmodified
 Added
 Removed

trunk/anuga_core/source/anuga/structures/culvert_operator.py
r7952 r7955 4 4 from anuga.structures.culvert_polygons import create_culvert_polygons 5 5 from anuga.utilities.system_tools import log_to_file 6 from anuga.geometry.polygon import inside_polygon, is_inside_polygon, plot_polygons 6 from anuga.geometry.polygon import inside_polygon, is_inside_polygon 7 from anuga.geometry.polygon import plot_polygons, polygon_area 8 7 9 8 10 from anuga.utilities.numerical_tools import mean … … 111 113 112 114 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. 115 class 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) 118 121 119 122 Input: Two points, pipe_size (either diameter or width, height), … … 123 126 def __init__(self, 124 127 domain, 125 culvert_description_filename=None,126 culvert_routine=None,127 128 end_point0=None, 128 129 end_point1=None, 129 enquiry_point0=None, 130 enquiry_point1=None, 131 type='box', 130 enquiry_gap_factor=0.2, 132 131 width=None, 133 132 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 happens138 manning=None, # Mannings Roughness for Culvert139 sum_loss=None,140 use_velocity_head=False, # FIXME(Ole): Get rid of  always True141 use_momentum_jet=False, # FIXME(Ole): Not yet implemented142 label=None,143 description=None,144 update_interval=None,145 log_file=False,146 discharge_hydrograph=False,147 133 verbose=False): 148 134 149 150 151 135 # Input check 152 136 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 167 139 168 140 self.height = height 169 141 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 212 150 213 151 # 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() 241 154 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 253 156 254 157 # Establish initial values at each enquiry point … … 326 229 s = 'Culvert Direction is %s\n' %str(self.vector) 327 230 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 = x1x0 261 dy = y1y0 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 332 322 333 323 def __call__(self, domain): … … 367 357 368 358 369 def get_enquiry_indices(self):359 def compute_enquiry_indices(self): 370 360 """Get indices for nearest centroids to self.enquiry_points 371 361 """ … … 402 392 raise Exception, msg 403 393 404 returnenquiry_indices394 self.enquiry_indices = enquiry_indices 405 395 406 396 … … 806 796 807 797 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 another814 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 topdown_blockage_factor and bottom_up_blockage_factor820 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 information841 self.label = label842 self.description = description843 self.culvert_type = type844 self.rating_curve = ensure_numeric(rating_curve)845 self.number_of_barrels = number_of_barrels846 847 if label is None: label = 'culvert_flow'848 label += '_' + str(id(self))849 self.label = label850 851 # File for storing discharge_hydrograph852 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 output859 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 POLYGON867 #if self.culvert_type == 'circle':868 # # Redefine width and height for use with create_culvert_polygons869 # width = height = diameter870 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 points878 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 pass886 #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 centroid900 N = len(domain)901 points = domain.get_centroid_coordinates(absolute=True)902 903 # Calculate indices in exchange area for this forcing term904 905 triangle_id = min_dist = sys.maxint906 for k in range(N):907 x, y = points[k,:] # Centroid908 909 c = point910 distance = (xc[0])**2+(yc[1])**2911 if distance < min_dist:912 min_dist = distance913 triangle_id = k914 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)' %point924 raise Exception, msg925 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), msg938 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.quantities951 for i, opening in enumerate(self.openings):952 elevation = dq['elevation'].get_values(location='centroids',953 indices=[self.enquiry_indices[i]])954 opening.elevation = elevation955 opening.stage = elevation # Simple assumption that culvert is dry initially956 957 # Assume two openings for now: Referred to as 0 and 1958 assert len(self.openings) == 2959 960 # Determine pipe direction961 self.delta_z = delta_z = self.openings[0].elevation  self.openings[1].elevation962 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 geometry971 self.end_points = [end_point0, end_point1]972 self.vector = P['vector']973 self.length = P['length']; assert self.length > 0.0974 if not num.allclose(self.length, length, rtol=1.0e2, atol=1.0e2):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.length978 log.critical(msg)979 980 self.verbose = verbose981 self.last_update = 0.0 # For use with update_interval982 self.last_time = 0.0983 self.update_interval = update_interval984 985 986 # Print something987 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 stuff1001 time = domain.get_time()1002 1003 1004 update = False1005 if self.update_interval is None:1006 update = True1007 delta_t = domain.timestep # Next timestep has been computed in domain.py1008 else:1009 if time  self.last_update > self.update_interval or time == 0.0:1010 update = True1011 delta_t = self.update_interval1012 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 = time1020 1021 dq = domain.quantities1022 1023 # Get average water depths at each opening1024 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 the1028 # enquiry area in front of the culvert1029 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 object1034 opening.depth = stage  opening.elevation1035 opening.stage = stage1036 1037 1038 1039 ################# End of the FOR loop ################################################1040 1041 # We now need to deal with each opening individually1042 1043 # Determine flow direction based on total energy difference1044 1045 delta_w = self.inlet.stage  self.outlet.stage1046 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.01056 else:1057 # Calculate discharge for one barrel and set inlet.rate and outlet.rate1058 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: ' %time1064 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: ' %time1072 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 barrels1082 Q *= self.number_of_barrels1083 1084 1085 # Adjust Q downwards depending on available water at inlet1086 stage = self.inlet.get_quantity_values(quantity_name='stage')1087 elevation = self.inlet.get_quantity_values(quantity_name='elevation')1088 depth = stageelevation1089 1090 1091 V = 01092 for i, d in enumerate(depth):1093 V += d * domain.areas[i]1094 1095 dt = delta_t1096 if Q*dt > V:1097 1098 Q_reduced = 0.9*V/dt # Reduce with safety factor1099 1100 msg = '%.2fs: Computed extraction for this time interval (Q*dt) is ' % time1101 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_reduced1110 1111 self.inlet.rate = Q1112 self.outlet.rate = Q1113 1114 # Log timeseries to file1115 try:1116 fid = open(self.timeseries_filename, 'a')1117 except:1118 pass1119 else:1120 fid.write('%.2f, %.2f\n' %(time, Q))1121 fid.close()1122 1123 # Store value of time1124 self.last_time = time1125 1126 1127 1128 # Execute flow term for each opening1129 # This is where Inflow objects are evaluated using the last rate that has been calculated1130 #1131 # This will take place at every internal timestep and update the domain1132 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 another1142 1143 Using Momentum as Calculated by Culvert Flow !!1144 Could be Several Methods Investigated to do This !!!1145 1146 2008_May_081147 To Ole:1148 OK so here we need to get the Polygon Creating code to create1149 polygons for the culvert Based on1150 the two input Points (X0,Y0) and (X1,Y1)  need to be passed1151 to create polygon1152 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 topdown_blockage_factor and bottom_up_blockage_factor1160 1161 1162 And the Delta H enquiry should be change from Openings in line 4121163 to the enquiry Polygons infront of the culvert1164 At the moment this script uses only Depth, later we can change it to1165 Energy...1166 1167 Once we have Delta H can calculate a Flow Rate and from Flow Rate1168 an Outlet Velocity1169 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 the1172 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 Culvert1186 invert_level0=None, # Invert level at opening 01187 invert_level1=None, # Invert level at opening 11188 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 check1202 if diameter is not None:1203 self.culvert_type = 'circle'1204 self.diameter = diameter1205 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, msg1209 else:1210 if height is not None:1211 if width is None:1212 self.culvert_type = 'square'1213 width = height1214 else:1215 self.culvert_type = 'rectangle'1216 elif width is not None:1217 if height is None:1218 self.culvert_type = 'square'1219 height = width1220 else:1221 msg = 'Either diameter or width&height must be specified.'1222 raise Exception, msg1223 1224 if height == width:1225 self.culvert_type = 'square'1226 1227 self.height = height1228 self.width = width1229 1230 1231 assert self.culvert_type in ['circle', 'square', 'rectangle']1232 1233 assert number_of_barrels >= 11234 self.number_of_barrels = number_of_barrels1235 1236 1237 # Set defaults1238 if manning is None: manning = 0.012 # Default roughness for pipe1239 if loss_exit is None: loss_exit = 1.001240 if loss_entry is None: loss_entry = 0.501241 if loss_bend is None: loss_bend=0.001242 if loss_special is None: loss_special=0.001243 if blockage_topdwn is None: blockage_topdwn=0.001244 if blockage_bottup is None: blockage_bottup=0.001245 if culvert_routine is None:1246 culvert_routine=boyd_generalised_culvert_model1247 1248 if label is None: label = 'culvert_flow'1249 label += '_' + str(id(self))1250 self.label = label1251 1252 # File for storing culvert quantities1253 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 output1259 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 POLYGON1266 if self.culvert_type == 'circle':1267 # Redefine width and height for use with create_culvert_polygons1268 width = height = diameter1269 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 points1277 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 pass1285 #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 centroid1299 N = len(domain)1300 points = domain.get_centroid_coordinates(absolute=True)1301 1302 # Calculate indices in exchange area for this forcing term1303 1304 triangle_id = min_dist = sys.maxint1305 for k in range(N):1306 x, y = points[k,:] # Centroid1307 1308 c = point1309 distance = (xc[0])**2+(yc[1])**21310 if distance < min_dist:1311 min_dist = distance1312 triangle_id = k1313 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)' %point1323 raise Exception, msg1324 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), msg1341 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 11353 assert len(self.openings) == 21354 1355 # Store basic geometry1356 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.01363 self.verbose = verbose1364 self.last_time = 0.01365 self.last_update = 0.0 # For use with update_interval1366 self.update_interval = update_interval1367 1368 1369 # Store hydraulic parameters1370 self.manning = manning1371 self.loss_exit = loss_exit1372 self.loss_entry = loss_entry1373 self.loss_bend = loss_bend1374 self.loss_special = loss_special1375 self.sum_loss = loss_exit + loss_entry + loss_bend + loss_special1376 self.blockage_topdwn = blockage_topdwn1377 self.blockage_bottup = blockage_bottup1378 1379 # Store culvert routine1380 self.culvert_routine = culvert_routine1381 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 something1400 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_filename1410 1411 # Time stuff1412 time = domain.get_time()1413 1414 # Short hand1415 dq = domain.quantities1416 1417 1418 update = False1419 if self.update_interval is None:1420 update = True1421 delta_t = domain.timestep # Next timestep has been computed in domain.py1422 else:1423 if time  self.last_update > self.update_interval or time == 0.0:1424 update = True1425 delta_t = self.update_interval1426 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 = time1434 1435 msg = 'Time did not advance'1436 if time > 0.0: assert delta_t > 0.0, msg1437 1438 1439 # Get average water depths at each opening1440 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 the1444 # exchange area in front of the culvert1445 1446 stage = opening.get_quantity_values(quantity_name='stage')1447 w = mean(stage) # Average stage1448 1449 # Use invert level instead of elevation if specified1450 invert_level = self.invert_levels[i]1451 if invert_level is not None:1452 z = invert_level1453 else:1454 elevation = opening.get_quantity_values(quantity_name='elevation')1455 z = mean(elevation) # Average elevation1456 1457 # Estimated depth above the culvert inlet1458 d = w  z # Used for calculations involving depth1459 if d < 0.0:1460 # This is possible since w and z are taken at different locations1461 #msg = 'D < 0.0: %f' %d1462 #raise msg1463 d = 0.01464 1465 1466 # Ratio of depth to culvert height.1467 # If ratio > 1 then culvert is running full1468 if self.culvert_type == 'circle':1469 ratio = d/self.diameter1470 else:1471 ratio = d/self.height1472 opening.ratio = ratio1473 1474 1475 # Average measures of energy in front of this opening1476 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 = stageelevation1487 if depth > 0.0:1488 u = xmomentum/(depth + velocity_protection/depth)1489 v = ymomentum/(depth + velocity_protection/depth)1490 else:1491 u = v = 0.01492 1493 1494 opening.total_energy = 0.5*(u*u + v*v)/g + stage1495 1496 # Store current average stage and depth with each opening object1497 opening.depth = d1498 opening.depth_trigger = d # Use this for now1499 opening.stage = w1500 opening.elevation = z1501 1502 1503 ################# End of the FOR loop ################################################1504 1505 # We now need to deal with each opening individually1506 1507 # Determine flow direction based on total energy difference1508 delta_Et = openings[0].total_energy  openings[1].total_energy1509 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_Et1525 1526 self.inlet = inlet1527 self.outlet = outlet1528 1529 msg = 'Total energy difference is negative'1530 assert delta_Et > 0.0, msg1531 1532 delta_h = inlet.stage  outlet.stage1533 delta_z = inlet.elevation  outlet.elevation1534 culvert_slope = (delta_z/self.length)1535 1536 if culvert_slope < 0.0:1537 # Adverse gradient  flow is running uphill1538 # Flow will be purely controlled by uphill outlet face1539 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.rate1550 Q, barrel_velocity, culvert_outlet_depth = self.culvert_routine(self, inlet, outlet, delta_Et, g)1551 1552 # Adjust discharge for multiple barrels1553 Q *= self.number_of_barrels1554 1555 # Compute barrel momentum1556 barrel_momentum = barrel_velocity*culvert_outlet_depth1557 1558 s = 'Barrel velocity = %f' %barrel_velocity1559 log_to_file(log_filename, s)1560 1561 # Compute momentum vector at outlet1562 outlet_mom_x, outlet_mom_y = self.vector * barrel_momentum1563 1564 s = 'Directional momentum = (%f, %f)' %(outlet_mom_x, outlet_mom_y)1565 log_to_file(log_filename, s)1566 1567 # Log timeseries to file1568 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 momentum1578 if delta_t > 0.0:1579 xmomentum_rate = outlet_mom_x  outlet.momentum[0].value1580 xmomentum_rate /= delta_t1581 1582 ymomentum_rate = outlet_mom_y  outlet.momentum[1].value1583 ymomentum_rate /= delta_t1584 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.01589 1590 1591 # Set momentum rates for outlet jet1592 outlet.momentum[0].rate = xmomentum_rate1593 outlet.momentum[1].rate = ymomentum_rate1594 1595 # Remember this value for next step (IMPORTANT)1596 outlet.momentum[0].value = outlet_mom_x1597 outlet.momentum[1].value = outlet_mom_y1598 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 time1611 self.last_time = time1612 1613 1614 1615 # Execute flow term for each opening1616 # This is where Inflow objects are evaluated and update the domain1617 for opening in self.openings:1618 opening(domain)1619 1620 # Execute momentum terms1621 # This is where Inflow objects are evaluated and update the domain1622 self.outlet.momentum[0](domain)1623 self.outlet.momentum[1](domain)1624 1625 1626 1627 798 Culvert_flow = Culvert_flow_general
Note: See TracChangeset
for help on using the changeset viewer.