Changeset 7957
 Timestamp:
 Aug 20, 2010, 9:12:41 AM (9 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/anuga_core/source/anuga/structures/culvert_operator.py
r7955 r7957 21 21 class Below_interval(Exception): pass 22 22 class Above_interval(Exception): pass 23 24 # FIXME(Ole): Take a good hard look at logging here25 26 27 # FIXME(Ole): Write in C and reuse this function by similar code28 # in interpolate.py29 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, msg37 38 39 # Check bounds40 if x < xvec[0]:41 msg = 'Value provided = %.2f, interpolation minimum = %.2f.'\42 % (x, xvec[0])43 raise Below_interval, msg44 45 if x > xvec[1]:46 msg = 'Value provided = %.2f, interpolation maximum = %.2f.'\47 %(x, xvec[1])48 raise Above_interval, msg49 50 51 # Find appropriate slot within bounds52 i = 053 while x > xvec[i]: i += 154 55 56 x0 = xvec[i1]57 x1 = xvec[i]58 alpha = (x  x0)/(x1  x0)59 60 y0 = yvec[i1]61 y1 = yvec[i]62 y = alpha*y1 + (1alpha)*y063 64 return y65 66 67 68 def read_culvert_description(culvert_description_filename):69 70 # Read description file71 fid = open(culvert_description_filename)72 73 read_rating_curve_data = False74 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 # Header87 continue88 if i == 1:89 # Metadata90 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 losses100 description=fields[7].strip()101 102 if line.strip() == '': continue # Skip blanks103 104 if line.startswith('Rating'):105 read_rating_curve_data = True106 # Flow data follows107 108 fid.close()109 110 return label, type, width, height, length, number_of_barrels, description, rating_curve111 112 23 113 24 … … 144 55 self.end_points= [end_point0, end_point1] 145 56 self.enquiry_gap_factor=enquiry_gap_factor 57 146 58 self.verbose=verbose 147 148 59 self.filename = None 149 60 … … 156 67 157 68 # Establish initial values at each enquiry point 69 self.enquiry_quantity_values = [] 158 70 dq = domain.quantities 159 for i , opening in enumerate(self.openings):71 for i in [0,1]: 160 72 idx = self.enquiry_indices[i] 161 73 elevation = dq['elevation'].get_values(location='centroids', … … 163 75 stage = dq['stage'].get_values(location='centroids', 164 76 indices=[idx])[0] 165 opening.elevation = elevation 166 opening.stage = stage 167 opening.depth = stageelevation 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.0e2, atol=1.0e2): 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 = stageelevation 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 231 84 232 85 … … 255 108 256 109 # Calculate geometry 257 x0, y0 = self.end_point 0258 x1, y1 = self.end_point 1110 x0, y0 = self.end_points[0] 111 x1, y1 = self.end_points[1] 259 112 260 113 dx = x1x0 … … 277 130 278 131 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) 300 151 301 152 msg = 'Polygon %s ' %(polygon) … … 303 154 assert area > 0.0, msg 304 155 305 for key2 in ['enquiry_point0', 'enquiry_point1']:306 point = culvert_polygons[key2]307 msg = 'Enquiry point falls inside a n enquiry point.'156 for j in [0,1]: 157 point = self.enquiry_points[j] 158 msg = 'Enquiry point falls inside a culvert polygon.' 308 159 309 160 assert not inside_polygon(point, polygon), msg 310 311 312 313 for key1 in ['exchange_polygon0',314 'exchange_polygon1']:315 316 317 # Return results318 self.culvert_polygons = culvert_polygons319 320 321 161 322 162 … … 794 634 795 635 636 # FIXME(Ole): Write in C and reuse this function by similar code 637 # in interpolate.py 638 def 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[i1] 666 x1 = xvec[i] 667 alpha = (x  x0)/(x1  x0) 668 669 y0 = yvec[i1] 670 y1 = yvec[i] 671 y = alpha*y1 + (1alpha)*y0 672 673 return y 674 675 676 677 def 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 796 721 797 722
Note: See TracChangeset
for help on using the changeset viewer.