Changeset 6533
- Timestamp:
- Mar 17, 2009, 4:02:54 PM (16 years ago)
- Location:
- branches/numpy/anuga
- Files:
-
- 22 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/numpy/anuga/abstract_2d_finite_volumes/domain.py
r6517 r6533 319 319 return self.mesh.get_boundary_polygon(*args, **kwargs) 320 320 321 # FIXME(Ole): This doesn't seem to be required322 321 def get_number_of_triangles_per_node(self, *args, **kwargs): 323 322 return self.mesh.get_number_of_triangles_per_node(*args, **kwargs) … … 349 348 def statistics(self, *args, **kwargs): 350 349 return self.mesh.statistics(*args, **kwargs) 351 352 def get_extent(self, *args, **kwargs):353 return self.mesh.get_extent(*args, **kwargs)354 350 355 351 ## -
branches/numpy/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py
r6517 r6533 88 88 # FIXME (Ole): We should rename f to function to be consistent with 89 89 # Transmissive_Momentum_Set_Stage_Boundary (cf posting by rrraman) 90 def __init__(self, domain=None, 91 f=None, 92 default_boundary=None, 93 verbose=False): 90 def __init__(self, domain = None, f = None): 94 91 Boundary.__init__(self) 95 self.default_boundary = default_boundary96 self.default_boundary_invoked = False # Flag97 self.domain = domain98 self.verbose = verbose99 92 100 93 try: … … 129 122 def evaluate(self, vol_id=None, edge_id=None): 130 123 # FIXME (Ole): I think this should be get_time(), see ticket:306 131 try: 132 res = self.f(self.domain.time) 133 except Modeltime_too_early, e: 134 raise Modeltime_too_early, e 135 except Modeltime_too_late, e: 136 if self.default_boundary is None: 137 raise Exception, e # Reraise exception 138 else: 139 # Pass control to default boundary 140 res = self.default_boundary.evaluate(vol_id, edge_id) 141 142 # Ensure that result cannot be manipulated 143 # This is a real danger in case the 144 # default_boundary is a Dirichlet type 145 # for instance. 146 res = res.copy() 147 148 if self.default_boundary_invoked is False: 149 if self.verbose: 150 # Issue warning the first time 151 msg = '%s' %str(e) 152 msg += 'Instead I will use the default boundary: %s\n'\ 153 %str(self.default_boundary) 154 msg += 'Note: Further warnings will be supressed' 155 print msg 156 157 # FIXME (Ole): Replace this crude flag with 158 # Python's ability to print warnings only once. 159 # See http://docs.python.org/lib/warning-filter.html 160 self.default_boundary_invoked = True 161 162 return res 124 return self.f(self.domain.time) 163 125 164 126 … … 337 299 if self.default_boundary_invoked is False: 338 300 # Issue warning the first time 339 if self.verbose: 340 msg = '%s' %str(e) 341 msg += 'Instead I will use the default boundary: %s\n'\ 342 %str(self.default_boundary) 343 msg += 'Note: Further warnings will be supressed' 344 print msg 301 msg = '%s' %str(e) 302 msg += 'Instead I will use the default boundary: %s\n'\ 303 %str(self.default_boundary) 304 msg += 'Note: Further warnings will be supressed' 305 warn(msg) 345 306 346 307 # FIXME (Ole): Replace this crude flag with -
branches/numpy/anuga/abstract_2d_finite_volumes/quantity.py
r6517 r6533 886 886 # a crash in fittng, so disabled it until I can investigate further 887 887 # Sorry. 23 Jan 2009. Logged as ticket:314 888 #if True: # Test will fail (31 Jan 2009) 889 if False: # Test will pass (31 Jan 2009) 888 if False: 890 889 # Use mesh as defined by domain 891 890 # This used to cause problems for caching due to quantities -
branches/numpy/anuga/abstract_2d_finite_volumes/test_domain.py
r6517 r6533 542 542 conserved_quantities =\ 543 543 ['stage', 'xmomentum', 'ymomentum']) 544 domain.set_default_order(1)545 544 domain.check_integrity() 546 545 … … 651 650 conserved_quantities =\ 652 651 ['stage', 'xmomentum', 'ymomentum']) 653 domain.set_default_order(1)654 652 domain.check_integrity() 655 653 … … 802 800 [ 11.0, 11.0, 11.0], 803 801 [ 11.0, 11.0, 11.0]]) 804 805 def test_that_mesh_methods_exist(self):806 """test_that_mesh_methods_exist807 808 Test that relavent mesh methods are made available in809 domain through composition810 """811 from mesh_factory import rectangular812 from shallow_water import Domain813 814 # Create basic mesh815 points, vertices, boundary = rectangular(1, 3)816 817 # Create shallow water domain818 domain = Domain(points, vertices, boundary)819 820 821 domain.get_centroid_coordinates()822 domain.get_radii()823 domain.get_areas()824 domain.get_area()825 domain.get_vertex_coordinates()826 domain.get_triangles()827 domain.get_nodes()828 domain.get_number_of_nodes()829 domain.get_normal(0,0)830 domain.get_intersecting_segments([[0.0, 0.0], [0.0, 1.0]])831 domain.get_disconnected_triangles()832 domain.get_boundary_tags()833 domain.get_boundary_polygon()834 #domain.get_number_of_triangles_per_node()835 domain.get_triangles_and_vertices_per_node()836 domain.get_interpolation_object()837 domain.get_tagged_elements()838 domain.get_lone_vertices()839 domain.get_unique_vertices()840 g = domain.get_georeference()841 domain.set_georeference(g)842 domain.build_tagged_elements_dictionary()843 domain.statistics()844 domain.get_extent()845 846 847 848 802 849 803 #------------------------------------------------------------- -
branches/numpy/anuga/abstract_2d_finite_volumes/test_general_mesh.py
r6517 r6533 91 91 str(verts))) 92 92 self.assert_(num.allclose(num.array([nodes_absolute[1], 93 nodes_absolute[0], 94 nodes_absolute[2]]), verts)) 95 verts = domain.get_vertex_coordinates(triangle_id=0, 96 absolute=True) 97 self.assert_(num.allclose(num.array([nodes_absolute[1], 98 nodes_absolute[0], 99 nodes_absolute[2]]), verts)) 100 101 93 nodes_absolute[0], 94 nodes_absolute[2]]), 95 verts), msg) 102 96 103 97 def test_get_vertex_coordinates_triangle_id(self): -
branches/numpy/anuga/abstract_2d_finite_volumes/util.py
r6517 r6533 925 925 verbose = False): 926 926 927 # FIXME(Ole): Shouldn't print statements here be governed by verbose?928 927 assert type(gauge_filename) == type(''), 'Gauge filename must be a string' 929 928 … … 972 971 raise msg 973 972 974 if verbose: 975 print 'swwfile', swwfile 973 print 'swwfile', swwfile 976 974 977 975 # Extract parent dir name and use as label … … 2491 2489 base_name=base, 2492 2490 verbose=verbose) 2493 #print 'sww files just after get_all_swwfiles()', sww_files2494 # fudge to get SWW files in 'correct' order, oldest on the left2495 sww_files.sort()2496 2497 if verbose:2498 print 'sww files', sww_files2499 2491 2500 2492 #to make all the quantities lower case for file_function … … 2505 2497 2506 2498 core_quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum'] 2507 gauge_file = out_name 2508 2509 heading = [quantity for quantity in quantities] 2510 heading.insert(0,'time') 2511 heading.insert(1,'hours') 2512 2513 #create a list of csv writers for all the points and write header 2514 points_writer = [] 2515 for point_i,point in enumerate(points): 2516 points_writer.append(writer(file(dir_name + sep + gauge_file 2517 + point_name[point_i] + '.csv', "wb"))) 2518 points_writer[point_i].writerow(heading) 2519 2520 if verbose: print 'Writing csv files' 2521 2522 quake_offset_time = None 2523 2499 2524 2500 for sww_file in sww_files: 2525 2501 sww_file = join(dir_name, sww_file+'.sww') -
branches/numpy/anuga/alpha_shape/alpha_shape.py
r6517 r6533 289 289 (denom[k]< EPSILON and denom[k] > -EPSILON)] 290 290 291 if num.a lltrue(denom != 0.0):292 dx = num.divide(y31*dist21 - y21*dist31,denom)293 dy = num.divide(x21*dist31 - x31*dist21,denom) 294 else:295 raise AlphaError296 291 if num.any(denom == 0.0): 292 raise AlphaError 293 294 dx = num.divide(y31*dist21 - y21*dist31, denom) 295 dy = num.divide(x21*dist31 - x31*dist21, denom) 296 297 297 self.triradius = 0.5*num.sqrt(dx*dx + dy*dy) 298 298 #print "triangle radii", self.triradius -
branches/numpy/anuga/caching/test_caching.py
r6517 r6533 386 386 387 387 # Check that hash value of callable objects don't change 388 # FIXME (Ole): The hash values do appear to change when OS389 # and/or dependencies are upgraded390 388 if os.name == 'posix' and os.uname()[4] in ['x86_64', 'ia64']: 391 389 # 64 bit hash values 392 f1hash = 7079146893884768701 393 f2hash = -6995306676314913340 394 395 # Prior to cluster upgrades Feb 2009 396 #f1hash = 1914027059797211698 397 #f2hash = 1914027059807087171 390 f1hash = 1914027059797211698 391 f2hash = 1914027059807087171 398 392 else: 399 # 32 bit hash values400 393 f1hash = -758136387 401 394 f2hash = -11221564 -
branches/numpy/anuga/compile_all.py
r6517 r6533 27 27 #execfile('test_all.py') 28 28 29 if sys.platform == 'win32':30 raw_input('Press any key') -
branches/numpy/anuga/config.py
r6517 r6533 91 91 92 92 # FIXME (Ole) Maybe get rid of order altogether and use beta_w 93 default_order = 2 93 # ... and isn't it about time we make the default 2? 94 default_order = 1 94 95 95 96 ################################################################################ -
branches/numpy/anuga/coordinate_transforms/redfearn.py
r6517 r6533 37 37 return sign*dd, mm, ss 38 38 39 def redfearn(lat, lon, false_easting=None, false_northing=None, 40 zone=None, central_meridian=None, scale_factor=None): 39 def redfearn(lat, lon, false_easting=None, false_northing=None, zone=None): 41 40 """Compute UTM projection using Redfearn's formula 42 41 … … 48 47 If zone is specified reproject lat and long to specified zone instead of 49 48 standard zone 50 If meridian is specified, reproject lat and lon to that instead of zone. In this case51 zone will be set to -1 to indicate non-UTM projection52 53 Note that zone and meridian cannot both be specifed54 49 """ 55 50 … … 62 57 a = 6378137.0 #Semi major axis 63 58 inverse_flattening = 298.257222101 #1/f 64 if scale_factor is None: 65 K0 = 0.9996 #Central scale factor 66 else: 67 K0 = scale_factor 68 #print 'scale', K0 59 K0 = 0.9996 #Central scale factor 69 60 zone_width = 6 #Degrees 70 61 … … 147 138 m = term1 + term2 + term3 + term4 #OK 148 139 149 if zone is not None and central_meridian is not None: 150 msg = 'You specified both zone and central_meridian. Provide only one of them' 151 raise Exception, msg 152 153 # Zone 140 #Zone 154 141 if zone is None: 155 142 zone = int((lon - longitude_of_western_edge_zone0)/zone_width) 156 143 157 # Central meridian 158 if central_meridian is None: 159 central_meridian = zone*zone_width+longitude_of_central_meridian_zone0 160 else: 161 zone = -1 144 central_meridian = zone*zone_width+longitude_of_central_meridian_zone0 162 145 163 146 omega = (lon-central_meridian)*pi/180 #Relative longitude (radians) -
branches/numpy/anuga/coordinate_transforms/test_geo_reference.py
r6517 r6533 235 235 y = 3.0 236 236 237 g = Geo_reference(56,x,y) 238 lofl = [[3.0,34.0], [64.0,6.0]] 239 new_lofl = g.get_absolute(lofl) 240 #print "lofl",lofl 241 #print "new_lofl",new_lofl 242 243 self.failUnless(type(new_lofl) == types.ListType, ' failed') 244 self.failUnless(type(new_lofl) == type(lofl), ' failed') 245 for point,new_point in map(None,lofl,new_lofl): 246 self.failUnless(point[0]+x==new_point[0], ' failed') 247 self.failUnless(point[1]+y==new_point[1], ' failed') 248 249 237 g = Geo_reference(56, x, y) 238 points = [[3.0,34.0], [64.0,6.0]] 239 new_points = g.get_absolute(points) 240 241 self.failUnless(isinstance(new_points, list), 'failed') 242 self.failUnless(type(new_points) == type(points), 'failed') 243 for point, new_point in map(None, points, new_points): 244 self.failUnless(point[0]+x == new_point[0], 'failed') 245 self.failUnless(point[1]+y == new_point[1], 'failed') 246 247 # test with no supplied offsets 250 248 g = Geo_reference() 251 lofl = [[3.0,34.0], [64.0,6.0]] 252 new_lofl = g.get_absolute(lofl) 253 #print "lofl",lofl 254 #print "new_lofl",new_lofl 255 256 self.failUnless(type(new_lofl) == types.ListType, ' failed') 257 self.failUnless(type(new_lofl) == type(lofl), ' failed') 258 for point,new_point in map(None,lofl,new_lofl): 259 self.failUnless(point[0]==new_point[0], ' failed') 260 self.failUnless(point[1]==new_point[1], ' failed') 249 points = [[3.0,34.0], [64.0,6.0]] 250 new_points = g.get_absolute(points) 251 252 self.failUnless(isinstance(new_points, list), 'failed') 253 self.failUnless(type(new_points) == type(points), 'failed') 254 for point, new_point in map(None, points, new_points): 255 self.failUnless(point[0] == new_point[0], 'failed') 256 self.failUnless(point[1] == new_point[1], 'failed') 261 257 262 258 # test that calling get_absolute twice does the right thing -
branches/numpy/anuga/coordinate_transforms/test_redfearn.py
r6517 r6533 122 122 123 123 def test_UTM_6_nonstandard_projection(self): 124 """test_UTM_6_nonstandard_projection 125 126 Test that projections can be forced to 127 use other than native zone. 128 129 Data is from Geraldton, WA 130 """ 131 132 133 # First test native projection (zone 50) 124 #Test 6 (Geraldton, WA) 125 126 #First test native projection (zone 50) 134 127 zone, easting, northing = redfearn(-29.233299999,114.05) 135 128 … … 274 267 # #assert allclose(northing, 6181725.1724276) 275 268 276 277 def Xtest_nonstandard_meridian_coinciding_with_native(self):278 """test_nonstandard_meridian_coinciding_with_native279 280 This test will verify that redfearn can be used to project281 points using an arbitrary central meridian that happens to282 coincide with the standard meridian at the center of a UTM zone.283 This is a preliminary test before testing this functionality284 with a truly arbitrary non-standard meridian.285 """286 287 # The file projection_test_points_z53.csv contains 10 points288 # which straddle the boundary between UTM zones 53 and 54.289 # They have been projected to zone 53 irrespective of where they290 # belong.291 292 path = get_pathname_from_package('anuga.coordinate_transforms')293 294 for forced_zone in [53, 54]:295 296 datafile = join(path, 'projection_test_points_z%d.csv' % forced_zone)297 fid = open(datafile)298 299 for line in fid.readlines()[1:]:300 fields = line.strip().split(',')301 302 lon = float(fields[1])303 lat = float(fields[2])304 x = float(fields[3])305 y = float(fields[4])306 307 zone, easting, northing = redfearn(lat, lon,308 zone=forced_zone)309 310 print311 print 'Lat', lat312 print 'Lon', lon313 print 'Zone', zone314 print 'Ref x', x, 'Computed x', easting, 'Close enough:', num.allclose(x, easting)315 print 'Ref y', y, 'Computed y', northing, 'Close enough:', num.allclose(y, northing)316 317 # Check calculation318 assert zone == forced_zone319 print320 #assert num.allclose(x, easting)321 #assert num.allclose(y, northing)322 323 324 325 326 def Xtest_nonstandard_meridian(self):327 """test_nonstandard_meridian328 329 This test will verify that redfearn can be used to project330 points using an arbitrary central meridian.331 """332 333 # The file projection_test_points.csv contains 10 points334 # which straddle the boundary between UTM zones 53 and 54.335 # They have been projected using a central meridian of 137.5336 # degrees (the boundary is 138 so it is pretty much right337 # in the middle of zones 53 and 54).338 339 path = get_pathname_from_package('anuga.coordinate_transforms')340 datafile = join(path, 'projection_test_points.csv')341 fid = open(datafile)342 343 for line in fid.readlines()[1:]:344 fields = line.strip().split(',')345 346 lon = float(fields[1])347 lat = float(fields[2])348 x = float(fields[3])349 y = float(fields[4])350 351 zone, easting, northing = redfearn(lat, lon,352 central_meridian=137.5,353 scale_factor=0.9996)354 355 print356 print 'Lat', lat357 print 'Lon', lon358 print 'Zone', zone359 print 'Ref x', x, 'Computed x', easting, 'Close enough:', num.allclose(x, easting)360 print 'Ref y', y, 'Computed y', northing, 'Close enough:', num.allclose(y, northing)361 362 # Check calculation363 assert zone == -1 # Indicates non UTM projection364 print365 #assert num.allclose(x, easting)366 #assert num.allclose(y, northing)367 368 # Test that zone and meridian can't both be specified369 try:370 zone, easting, northing = redfearn(lat, lon,371 zone=50,372 central_meridian=137.5)373 except:374 pass375 else:376 msg = 'Should have raised exception'377 raise Exception, msg378 379 380 269 def test_convert_lats_longs(self): 381 270 -
branches/numpy/anuga/culvert_flows/Test_Culvert_Flat_Water_Lev.py
r6517 r6533 123 123 number_of_barrels=1, 124 124 update_interval=2, 125 log_file=True,126 discharge_hydrograph=True,127 125 verbose=True) 128 126 -
branches/numpy/anuga/culvert_flows/culvert_class.py
r6517 r6533 204 204 log_to_file(self.log_filename, description) 205 205 log_to_file(self.log_filename, self.culvert_type) 206 else:207 self.log_filename = None208 206 209 207 … … 298 296 299 297 # Print some diagnostics to log if requested 300 if self.log_filename is not None:298 if hasattr(self, 'log_filename'): 301 299 s = 'Culvert Effective Length = %.2f m' %(self.length) 302 300 log_to_file(self.log_filename, s) … … 326 324 update = True 327 325 328 329 if self.log_filename is not None: 326 if hasattr(self, 'log_filename'): 330 327 s = '\nTime = %.2f, delta_t = %f' %(time, delta_t) 331 328 log_to_file(self.log_filename, s) … … 442 439 if self.verbose is True: 443 440 print msg 444 445 if self.log_filename is not None: 441 if hasattr(self, 'log_filename'): 446 442 log_to_file(self.log_filename, msg) 447 443 … … 528 524 print '%.2fs - WARNING: Flow is running uphill.' % time 529 525 530 if self.log_filename is not None:526 if hasattr(self, 'log_filename'): 531 527 s = 'Time=%.2f, inlet stage = %.2f, outlet stage = %.2f'\ 532 528 %(time, self.inlet.stage, self.outlet.stage) … … 577 573 msg += 'for culvert "%s"' % self.label 578 574 579 if self.log_filename is not None:575 if hasattr(self, 'log_filename'): 580 576 log_to_file(self.log_filename, msg) 581 577 else: 582 578 # User culvert routine 583 579 Q, barrel_velocity, culvert_outlet_depth =\ 584 self.culvert_routine(inlet.depth, 585 outlet.depth, 586 inlet.specific_energy, 587 delta_total_energy, 588 g, 589 culvert_length=self.length, 590 culvert_width=self.width, 591 culvert_height=self.height, 592 culvert_type=self.culvert_type, 593 manning=self.manning, 594 sum_loss=self.sum_loss, 595 log_filename=self.log_filename) 580 self.culvert_routine(self, delta_total_energy, g) 596 581 597 582 … … 616 601 617 602 618 # OBSOLETE (Except for momentum jet in Culvert_flow_energy)603 619 604 class Culvert_flow_rating: 620 605 """Culvert flow - transfer water from one hole to another -
branches/numpy/anuga/culvert_flows/culvert_routines.py
r6517 r6533 8 8 9 9 #NOTE: 10 # Inlet control: Delta_total_energy > inlet_specific_energy 11 # Outlet control: Delta_total_energy < inlet_specific_energy 12 # where total energy is (w + 0.5*v^2/g) and 13 # specific energy is (h + 0.5*v^2/g) 10 # Inlet control: Delta_Et > Es at the inlet 11 # Outlet control: Delta_Et < Es at the inlet 12 # where Et is total energy (w + 0.5*v^2/g) and 13 # Es is the specific energy (h + 0.5*v^2/g) 14 15 16 17 # NEW DEFINED CULVERT FLOW---- Flow from INLET 1 ------> INLET 2 (Outlet) 18 # 19 # The First Attempt has a Simple Horizontal Circle as a Hole on the Bed 20 # Flow Is Removed at a Rate of INFLOW 21 # Downstream there is a similar Circular Hole on the Bed where INFLOW effectively Surcharges 22 # 23 # This SHould be changed to a Vertical Opening Both BOX and Circular 24 # There will be several Culvert Routines such as: 25 # CULVERT_Boyd_Channel 26 # CULVERT_Orifice_and_Weir 27 # CULVERT_Simple_FLOOR 28 # CULVERT_Simple_WALL 29 # CULVERT_Eqn_Floor 30 # CULVERT_Eqn_Wall 31 # CULVERT_Tab_Floor 32 # CULVERT_Tab_Wall 33 # BRIDGES..... 34 # NOTE NEED TO DEVELOP CONCEPT 1D Model for Linked Pipe System !!!! 35 36 # COULD USE EPA SWMM Model !!!! 14 37 15 38 … … 17 40 18 41 19 def boyd_generalised_culvert_model(inlet_depth, 20 outlet_depth, 21 inlet_specific_energy, 22 delta_total_energy, 23 g, 24 culvert_length=0.0, 25 culvert_width=0.0, 26 culvert_height=0.0, 27 culvert_type='box', 28 manning=0.0, 29 sum_loss=0.0, 30 log_filename=None): 31 32 """Boyd's generalisation of the US department of transportation culvert 33 model 34 35 The quantity of flow passing through a culvert is controlled by many factors 36 It could be that the culvert is controlled by the inlet only, with it 37 being unsubmerged this is effectively equivalent to the weir Equation 38 Else the culvert could be controlled by the inlet, with it being 39 submerged, this is effectively the Orifice Equation 40 Else it may be controlled by down stream conditions where depending on 41 the down stream depth, the momentum in the culvert etc. flow is controlled 42 def boyd_generalised_culvert_model(culvert, delta_total_energy, g): 43 44 """Boyd's generalisation of the US department of transportation culvert model 45 # == The quantity of flow passing through a culvert is controlled by many factors 46 # == It could be that the culvert is controled by the inlet only, with it being Un submerged this is effectively equivalent to the WEIR Equation 47 # == Else the culvert could be controlled by the inlet, with it being Submerged, this is effectively the Orifice Equation 48 # == Else it may be controlled by Down stream conditions where depending on the down stream depth, the momentum in the culvert etc. flow is controlled 42 49 """ 43 50 … … 46 53 from anuga.utilities.numerical_tools import safe_acos as acos 47 54 48 if inlet_depth > 0.01: 55 inlet = culvert.inlet 56 outlet = culvert.outlet 57 Q_outlet_tailwater = 0.0 58 inlet.rate = 0.0 59 outlet.rate = 0.0 60 Q_inlet_unsubmerged = 0.0 61 Q_inlet_submerged = 0.0 62 Q_outlet_critical_depth = 0.0 63 64 if hasattr(culvert, 'log_filename'): 65 log_filename = culvert.log_filename 66 67 manning = culvert.manning 68 sum_loss = culvert.sum_loss 69 length = culvert.length 70 71 if inlet.depth > 0.01: 72 73 E = inlet.specific_energy 74 75 if hasattr(culvert, 'log_filename'): 76 s = 'Specific energy = %f m' %E 77 log_to_file(log_filename, s) 78 79 msg = 'Specific energy is negative' 80 assert E >= 0.0, msg 81 82 49 83 # Water has risen above inlet 50 51 if log_filename is not None: 52 s = 'Specific energy = %f m' % inlet_specific_energy 53 log_to_file(log_filename, s) 54 55 msg = 'Specific energy at inlet is negative' 56 assert inlet_specific_energy >= 0.0, msg 57 58 if culvert_type == 'circle': 59 # Round culvert (use width as diameter) 60 diameter = culvert_width 84 if culvert.culvert_type == 'circle': 85 # Round culvert 61 86 62 87 # Calculate flows for inlet control 63 Q_inlet_unsubmerged = 0.421*g**0.5*diameter**0.87*inlet_specific_energy**1.63 # Inlet Ctrl Inlet Unsubmerged 64 Q_inlet_submerged = 0.530*g**0.5*diameter**1.87*inlet_specific_energy**0.63 # Inlet Ctrl Inlet Submerged 65 66 if log_filename is not None: 67 s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' % (Q_inlet_unsubmerged, Q_inlet_submerged) 88 diameter = culvert.diameter 89 90 Q_inlet_unsubmerged = 0.421*g**0.5*diameter**0.87*E**1.63 # Inlet Ctrl Inlet Unsubmerged 91 Q_inlet_submerged = 0.530*g**0.5*diameter**1.87*E**0.63 # Inlet Ctrl Inlet Submerged 92 93 if hasattr(culvert, 'log_filename'): 94 s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' %(Q_inlet_unsubmerged, Q_inlet_submerged) 68 95 log_to_file(log_filename, s) 69 96 70 # FIXME(Ole): Are these functions really for inlet control?97 case = '' 71 98 if Q_inlet_unsubmerged < Q_inlet_submerged: 72 99 Q = Q_inlet_unsubmerged 73 100 74 alpha = acos(1 - inlet _depth/diameter)101 alpha = acos(1 - inlet.depth/diameter) 75 102 flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha)) 76 outlet_culvert_depth = inlet_depth 103 outlet_culvert_depth = inlet.depth 104 width = diameter*sin(alpha) 105 #perimeter = alpha*diameter 77 106 case = 'Inlet unsubmerged' 78 107 else: … … 80 109 flow_area = (diameter/2)**2 * pi 81 110 outlet_culvert_depth = diameter 111 width = diameter 112 #perimeter = diameter 82 113 case = 'Inlet submerged' 83 114 84 115 85 116 86 if delta_total_energy < inlet_specific_energy:117 if delta_total_energy < E: 87 118 # Calculate flows for outlet control 88 119 # Determine the depth at the outlet relative to the depth of flow in the Culvert 89 120 90 if outlet _depth > diameter: # The Outlet is Submerged121 if outlet.depth > diameter: # The Outlet is Submerged 91 122 outlet_culvert_depth=diameter 92 123 flow_area = (diameter/2)**2 * pi # Cross sectional area of flow in the culvert 93 124 perimeter = diameter * pi 125 width = diameter 94 126 case = 'Outlet submerged' 95 elif outlet _depth==0.0:96 outlet_culvert_depth=inlet _depth# For purpose of calculation assume the outlet depth = the inlet depth97 alpha = acos(1 - inlet _depth/diameter)127 elif outlet.depth==0.0: 128 outlet_culvert_depth=inlet.depth # For purpose of calculation assume the outlet depth = the inlet depth 129 alpha = acos(1 - inlet.depth/diameter) 98 130 flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha)) 99 131 perimeter = alpha*diameter 132 width = diameter*sin(alpha) 100 133 101 134 case = 'Outlet depth is zero' 102 135 else: # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity 103 outlet_culvert_depth=outlet _depth # For purpose of calculation assume the outlet depth = the inlet depth104 alpha = acos(1 - outlet _depth/diameter)136 outlet_culvert_depth=outlet.depth # For purpose of calculation assume the outlet depth = the inlet depth 137 alpha = acos(1 - outlet.depth/diameter) 105 138 flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha)) 106 139 perimeter = alpha*diameter 140 width = diameter*sin(alpha) 107 141 case = 'Outlet is open channel flow' 108 142 109 143 hyd_rad = flow_area/perimeter 110 144 111 if log_filename is not None:145 if hasattr(culvert, 'log_filename'): 112 146 s = 'hydraulic radius at outlet = %f' %hyd_rad 113 147 log_to_file(log_filename, s) 114 148 115 149 # Outlet control velocity using tail water 116 culvert_velocity = sqrt(delta_total_energy/((sum_loss/2*g)+(manning**2* culvert_length)/hyd_rad**1.33333))150 culvert_velocity = sqrt(delta_total_energy/((sum_loss/2*g)+(manning**2*length)/hyd_rad**1.33333)) 117 151 Q_outlet_tailwater = flow_area * culvert_velocity 118 152 119 if log_filename is not None:153 if hasattr(culvert, 'log_filename'): 120 154 s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater 121 155 log_to_file(log_filename, s) … … 131 165 132 166 # Calculate flows for inlet control 133 height = culvert _height134 width = culvert _width167 height = culvert.height 168 width = culvert.width 135 169 136 Q_inlet_unsubmerged = 0.540*g**0.5*width* inlet_specific_energy**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged137 Q_inlet_submerged = 0.702*g**0.5*width*height**0.89* inlet_specific_energy**0.61 # Flow based on Inlet Ctrl Inlet Submerged138 139 if log_filename is not None:170 Q_inlet_unsubmerged = 0.540*g**0.5*width*E**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged 171 Q_inlet_submerged = 0.702*g**0.5*width*height**0.89*E**0.61 # Flow based on Inlet Ctrl Inlet Submerged 172 173 if hasattr(culvert, 'log_filename'): 140 174 s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' %(Q_inlet_unsubmerged, Q_inlet_submerged) 141 175 log_to_file(log_filename, s) 142 176 143 # FIXME(Ole): Are these functions really for inlet control?177 case = '' 144 178 if Q_inlet_unsubmerged < Q_inlet_submerged: 145 179 Q = Q_inlet_unsubmerged 146 flow_area = width*inlet_depth 147 outlet_culvert_depth = inlet_depth 180 flow_area = width*inlet.depth 181 outlet_culvert_depth = inlet.depth 182 #perimeter=(width+2.0*inlet.depth) 148 183 case = 'Inlet unsubmerged' 149 184 else: … … 154 189 case = 'Inlet submerged' 155 190 156 if delta_total_energy < inlet_specific_energy:191 if delta_total_energy < E: 157 192 # Calculate flows for outlet control 158 193 # Determine the depth at the outlet relative to the depth of flow in the Culvert 159 194 160 if outlet _depth > height: # The Outlet is Submerged195 if outlet.depth > height: # The Outlet is Submerged 161 196 outlet_culvert_depth=height 162 197 flow_area=width*height # Cross sectional area of flow in the culvert 163 198 perimeter=2.0*(width+height) 164 199 case = 'Outlet submerged' 165 elif outlet _depth==0.0:166 outlet_culvert_depth=inlet _depth # For purpose of calculation assume the outlet depth = the inlet depth167 flow_area=width*inlet _depth168 perimeter=(width+2.0*inlet _depth)200 elif outlet.depth==0.0: 201 outlet_culvert_depth=inlet.depth # For purpose of calculation assume the outlet depth = the inlet depth 202 flow_area=width*inlet.depth 203 perimeter=(width+2.0*inlet.depth) 169 204 case = 'Outlet depth is zero' 170 205 else: # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity 171 outlet_culvert_depth=outlet _depth172 flow_area=width*outlet _depth173 perimeter=(width+2.0*outlet _depth)206 outlet_culvert_depth=outlet.depth 207 flow_area=width*outlet.depth 208 perimeter=(width+2.0*outlet.depth) 174 209 case = 'Outlet is open channel flow' 175 210 176 211 hyd_rad = flow_area/perimeter 177 212 178 if log_filename is not None:179 s = 'hydraulic radius at outlet = %f' % 213 if hasattr(culvert, 'log_filename'): 214 s = 'hydraulic radius at outlet = %f' %hyd_rad 180 215 log_to_file(log_filename, s) 181 216 182 217 # Outlet control velocity using tail water 183 culvert_velocity = sqrt(delta_total_energy/((sum_loss/2*g)+(manning**2* culvert_length)/hyd_rad**1.33333))218 culvert_velocity = sqrt(delta_total_energy/((sum_loss/2*g)+(manning**2*length)/hyd_rad**1.33333)) 184 219 Q_outlet_tailwater = flow_area * culvert_velocity 185 220 186 if log_filename is not None:187 s = 'Q_outlet_tailwater = %.6f' % 221 if hasattr(culvert, 'log_filename'): 222 s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater 188 223 log_to_file(log_filename, s) 189 224 Q = min(Q, Q_outlet_tailwater) … … 193 228 194 229 # Common code for circle and square geometries 195 if log_filename is not None:196 log_to_file(log_filename, 'Case: "%s"' % 230 if hasattr(culvert, 'log_filename'): 231 log_to_file(log_filename, 'Case: "%s"' %case) 197 232 198 if log_filename is not None: 199 s = 'Flow Rate Control = %f' % Q 233 flow_rate_control=Q 234 235 if hasattr(culvert, 'log_filename'): 236 s = 'Flow Rate Control = %f' %flow_rate_control 200 237 log_to_file(log_filename, s) 201 238 202 culv_froude=sqrt(Q**2*width/(g*flow_area**3)) 203 204 if log_filename is not None: 205 s = 'Froude in Culvert = %f' % culv_froude 239 inlet.rate = -flow_rate_control 240 outlet.rate = flow_rate_control 241 242 culv_froude=sqrt(flow_rate_control**2*width/(g*flow_area**3)) 243 244 if hasattr(culvert, 'log_filename'): 245 s = 'Froude in Culvert = %f' %culv_froude 206 246 log_to_file(log_filename, s) 207 247 … … 210 250 211 251 212 else: # inlet_depth < 0.01:252 else: #inlet.depth < 0.01: 213 253 Q = barrel_velocity = outlet_culvert_depth = 0.0 214 254 -
branches/numpy/anuga/fit_interpolate/fit.py
r6517 r6533 310 310 assert flag is True, msg 311 311 312 # FIXME(Ole): This is the message referred to in ticket:314313 312 minx = min(self.mesh_boundary_polygon[:,0]) 314 313 maxx = max(self.mesh_boundary_polygon[:,0]) … … 318 317 msg += 'Mesh boundary extent is (%.f, %.f), (%.f, %.f)'\ 319 318 % (minx, maxx, miny, maxy) 320 raise RuntimeError, msg321 319 322 320 … … 377 375 self.build_fit_subset(points, z, verbose=verbose) 378 376 379 # FIXME(Ole): I thought this test would make sense here380 # See test_fitting_example_that_crashed_2 in test_shallow_water_domain.py381 # Committed 11 March 2009382 msg = 'Matrix AtA was not built'383 assert self.AtA is not None, msg384 385 #print 'Matrix was built OK'386 387 377 388 378 point_coordinates = None … … 392 382 if point_coordinates is None: 393 383 if verbose: print 'Warning: no data points in fit' 394 msg = 'No interpolation matrix.' 395 assert self.AtA is not None, msg 396 assert self.Atz is not None 384 #assert self.AtA != None, 'no interpolation matrix' 385 #assert self.Atz != None 386 assert not self.AtA is None, 'no interpolation matrix' 387 assert not self.Atz is None 397 388 398 389 # FIXME (DSG) - do a message -
branches/numpy/anuga/geospatial_data/geospatial_data.py
r6517 r6533 1806 1806 points = G_small.get_data_points() 1807 1807 1808 #FIXME: Remove points outside boundary polygon1809 # print 'new point',len(points)1810 #1811 # new_points=[]1812 # new_points=array([],typecode=Float)1813 # new_points=resize(new_points,(len(points),2))1814 # print "BOUNDARY", boundary_poly1815 # for i,point in enumerate(points):1816 # if is_inside_polygon(point,boundary_poly, verbose=True):1817 # new_points[i] = point1818 # print"WOW",i,new_points[i]1819 # points = new_points1820 1821 1808 if verbose: print "Number of points in sample to compare: ", len(points) 1822 1809 -
branches/numpy/anuga/shallow_water/__init__.py
r6517 r6533 12 12 Transmissive_momentum_set_stage_boundary,\ 13 13 Dirichlet_discharge_boundary,\ 14 Field_boundary,\ 15 Transmissive_stage_zero_momentum_boundary 14 Field_boundary 16 15 17 16 -
branches/numpy/anuga/shallow_water/data_manager.py
r6517 r6533 5863 5863 # This is being used to seperate one number from a list. 5864 5864 # what it is actually doing is sorting lists from numeric arrays. 5865 #if type(times) is list or isinstance(times, num.ndarray): 5866 if isinstance(times, (list, num.ndarray)): 5865 if type(times) is list or isinstance(times, num.ndarray): 5867 5866 number_of_times = len(times) 5868 5867 times = ensure_numeric(times) … … 5931 5930 #outfile.variables[q+Write_sww.RANGE][1] = -max_float # Max 5932 5931 5933 #if type(times) is list or isinstance(times, num.ndarray): 5934 if isinstance(times, (list, num.ndarray)): 5935 outfile.variables['time'][:] = times # Store time relative 5932 if type(times) is list or isinstance(times, num.ndarray): 5933 outfile.variables['time'][:] = times #Store time relative 5936 5934 5937 5935 if verbose: … … 6267 6265 # This is being used to seperate one number from a list. 6268 6266 # what it is actually doing is sorting lists from numeric arrays. 6269 #if type(times) is list or isinstance(times, num.ndarray): 6270 if isinstance(times, (list, num.ndarray)): 6267 if type(times) is list or isinstance(times, num.ndarray): 6271 6268 number_of_times = len(times) 6272 6269 times = ensure_numeric(times) … … 6315 6312 outfile.variables[q + Write_sts.RANGE][1] = -max_float # Max 6316 6313 6317 #if type(times) is list or isinstance(times, num.ndarray): 6318 if isinstance(times, (list, num.ndarray)): 6314 if type(times) is list or isinstance(times, num.ndarray): 6319 6315 outfile.variables['time'][:] = times #Store time relative 6320 6316 -
branches/numpy/anuga/shallow_water/test_shallow_water_domain.py
r6517 r6533 1672 1672 # Setup computational domain 1673 1673 #----------------------------------------------------------------- 1674 points, vertices, boundary = rectangular_cross(10, 10) # Basic mesh 1675 domain = Domain(points, vertices, boundary) # Create domain 1676 domain.set_default_order(1) 1674 points, vertices, boundary = rectangular_cross(10, 10) # Basic mesh 1675 domain = Domain(points, vertices, boundary) # Create domain 1677 1676 domain.set_quantities_to_be_stored(None) 1678 1677 domain.set_maximum_allowed_speed(100) #FIXME (Ole): try to remove this -
branches/numpy/anuga/test_all.py
r6517 r6533 215 215 #os.remove(filename) 216 216 217 218 if sys.platform == 'win32':219 raw_input('Press any key')
Note: See TracChangeset
for help on using the changeset viewer.