Changeset 6517
- Timestamp:
- Mar 16, 2009, 11:06:22 AM (16 years ago)
- Location:
- branches/numpy/anuga
- Files:
-
- 22 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/numpy/anuga/abstract_2d_finite_volumes/domain.py
r6481 r6517 319 319 return self.mesh.get_boundary_polygon(*args, **kwargs) 320 320 321 # FIXME(Ole): This doesn't seem to be required 321 322 def get_number_of_triangles_per_node(self, *args, **kwargs): 322 323 return self.mesh.get_number_of_triangles_per_node(*args, **kwargs) … … 348 349 def statistics(self, *args, **kwargs): 349 350 return self.mesh.statistics(*args, **kwargs) 351 352 def get_extent(self, *args, **kwargs): 353 return self.mesh.get_extent(*args, **kwargs) 350 354 351 355 ## -
branches/numpy/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py
r6304 r6517 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, f = None): 90 def __init__(self, domain=None, 91 f=None, 92 default_boundary=None, 93 verbose=False): 91 94 Boundary.__init__(self) 95 self.default_boundary = default_boundary 96 self.default_boundary_invoked = False # Flag 97 self.domain = domain 98 self.verbose = verbose 92 99 93 100 try: … … 122 129 def evaluate(self, vol_id=None, edge_id=None): 123 130 # FIXME (Ole): I think this should be get_time(), see ticket:306 124 return self.f(self.domain.time) 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 125 163 126 164 … … 299 337 if self.default_boundary_invoked is False: 300 338 # Issue warning the first time 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) 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 306 345 307 346 # FIXME (Ole): Replace this crude flag with -
branches/numpy/anuga/abstract_2d_finite_volumes/quantity.py
r6441 r6517 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 False: 888 #if True: # Test will fail (31 Jan 2009) 889 if False: # Test will pass (31 Jan 2009) 889 890 # Use mesh as defined by domain 890 891 # This used to cause problems for caching due to quantities -
branches/numpy/anuga/abstract_2d_finite_volumes/test_domain.py
r6410 r6517 542 542 conserved_quantities =\ 543 543 ['stage', 'xmomentum', 'ymomentum']) 544 domain.set_default_order(1) 544 545 domain.check_integrity() 545 546 … … 650 651 conserved_quantities =\ 651 652 ['stage', 'xmomentum', 'ymomentum']) 653 domain.set_default_order(1) 652 654 domain.check_integrity() 653 655 … … 800 802 [ 11.0, 11.0, 11.0], 801 803 [ 11.0, 11.0, 11.0]]) 804 805 def test_that_mesh_methods_exist(self): 806 """test_that_mesh_methods_exist 807 808 Test that relavent mesh methods are made available in 809 domain through composition 810 """ 811 from mesh_factory import rectangular 812 from shallow_water import Domain 813 814 # Create basic mesh 815 points, vertices, boundary = rectangular(1, 3) 816 817 # Create shallow water domain 818 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 802 848 803 849 #------------------------------------------------------------- -
branches/numpy/anuga/abstract_2d_finite_volumes/test_general_mesh.py
r6463 r6517 91 91 str(verts))) 92 92 self.assert_(num.allclose(num.array([nodes_absolute[1], 93 nodes_absolute[0], 94 nodes_absolute[2]]), 95 verts), msg) 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 96 102 97 103 def test_get_vertex_coordinates_triangle_id(self): -
branches/numpy/anuga/abstract_2d_finite_volumes/util.py
r6441 r6517 925 925 verbose = False): 926 926 927 # FIXME(Ole): Shouldn't print statements here be governed by verbose? 927 928 assert type(gauge_filename) == type(''), 'Gauge filename must be a string' 928 929 … … 971 972 raise msg 972 973 973 print 'swwfile', swwfile 974 if verbose: 975 print 'swwfile', swwfile 974 976 975 977 # Extract parent dir name and use as label … … 2489 2491 base_name=base, 2490 2492 verbose=verbose) 2493 #print 'sww files just after get_all_swwfiles()', sww_files 2494 # fudge to get SWW files in 'correct' order, oldest on the left 2495 sww_files.sort() 2496 2497 if verbose: 2498 print 'sww files', sww_files 2491 2499 2492 2500 #to make all the quantities lower case for file_function … … 2497 2505 2498 2506 core_quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum'] 2499 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 2500 2524 for sww_file in sww_files: 2501 2525 sww_file = join(dir_name, sww_file+'.sww') -
branches/numpy/anuga/alpha_shape/alpha_shape.py
r6304 r6517 289 289 (denom[k]< EPSILON and denom[k] > -EPSILON)] 290 290 291 if num.a ny(denom == 0.0):292 raise AlphaError293 294 dx = num.divide(y31*dist21 - y21*dist31, denom)295 dy = num.divide(x21*dist31 - x31*dist21, denom)296 291 if num.alltrue(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 AlphaError 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
r6410 r6517 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 OS 389 # and/or dependencies are upgraded 388 390 if os.name == 'posix' and os.uname()[4] in ['x86_64', 'ia64']: 389 391 # 64 bit hash values 390 f1hash = 1914027059797211698 391 f2hash = 1914027059807087171 392 f1hash = 7079146893884768701 393 f2hash = -6995306676314913340 394 395 # Prior to cluster upgrades Feb 2009 396 #f1hash = 1914027059797211698 397 #f2hash = 1914027059807087171 392 398 else: 399 # 32 bit hash values 393 400 f1hash = -758136387 394 401 f2hash = -11221564 -
branches/numpy/anuga/compile_all.py
r4978 r6517 27 27 #execfile('test_all.py') 28 28 29 if sys.platform == 'win32': 30 raw_input('Press any key') -
branches/numpy/anuga/config.py
r6361 r6517 91 91 92 92 # FIXME (Ole) Maybe get rid of order altogether and use beta_w 93 # ... and isn't it about time we make the default 2? 94 default_order = 1 93 default_order = 2 95 94 96 95 ################################################################################ -
branches/numpy/anuga/coordinate_transforms/redfearn.py
r6149 r6517 37 37 return sign*dd, mm, ss 38 38 39 def redfearn(lat, lon, false_easting=None, false_northing=None, zone=None): 39 def redfearn(lat, lon, false_easting=None, false_northing=None, 40 zone=None, central_meridian=None, scale_factor=None): 40 41 """Compute UTM projection using Redfearn's formula 41 42 … … 47 48 If zone is specified reproject lat and long to specified zone instead of 48 49 standard zone 50 If meridian is specified, reproject lat and lon to that instead of zone. In this case 51 zone will be set to -1 to indicate non-UTM projection 52 53 Note that zone and meridian cannot both be specifed 49 54 """ 50 55 … … 57 62 a = 6378137.0 #Semi major axis 58 63 inverse_flattening = 298.257222101 #1/f 59 K0 = 0.9996 #Central scale factor 64 if scale_factor is None: 65 K0 = 0.9996 #Central scale factor 66 else: 67 K0 = scale_factor 68 #print 'scale', K0 60 69 zone_width = 6 #Degrees 61 70 … … 138 147 m = term1 + term2 + term3 + term4 #OK 139 148 140 #Zone 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 141 154 if zone is None: 142 155 zone = int((lon - longitude_of_western_edge_zone0)/zone_width) 143 156 144 central_meridian = zone*zone_width+longitude_of_central_meridian_zone0 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 145 162 146 163 omega = (lon-central_meridian)*pi/180 #Relative longitude (radians) -
branches/numpy/anuga/coordinate_transforms/test_geo_reference.py
r6442 r6517 235 235 y = 3.0 236 236 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 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 248 250 g = Geo_reference() 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') 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') 257 261 258 262 # test that calling get_absolute twice does the right thing -
branches/numpy/anuga/coordinate_transforms/test_redfearn.py
r6410 r6517 122 122 123 123 def test_UTM_6_nonstandard_projection(self): 124 #Test 6 (Geraldton, WA) 125 126 #First test native projection (zone 50) 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) 127 134 zone, easting, northing = redfearn(-29.233299999,114.05) 128 135 … … 267 274 # #assert allclose(northing, 6181725.1724276) 268 275 276 277 def Xtest_nonstandard_meridian_coinciding_with_native(self): 278 """test_nonstandard_meridian_coinciding_with_native 279 280 This test will verify that redfearn can be used to project 281 points using an arbitrary central meridian that happens to 282 coincide with the standard meridian at the center of a UTM zone. 283 This is a preliminary test before testing this functionality 284 with a truly arbitrary non-standard meridian. 285 """ 286 287 # The file projection_test_points_z53.csv contains 10 points 288 # which straddle the boundary between UTM zones 53 and 54. 289 # They have been projected to zone 53 irrespective of where they 290 # 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 print 311 print 'Lat', lat 312 print 'Lon', lon 313 print 'Zone', zone 314 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 calculation 318 assert zone == forced_zone 319 print 320 #assert num.allclose(x, easting) 321 #assert num.allclose(y, northing) 322 323 324 325 326 def Xtest_nonstandard_meridian(self): 327 """test_nonstandard_meridian 328 329 This test will verify that redfearn can be used to project 330 points using an arbitrary central meridian. 331 """ 332 333 # The file projection_test_points.csv contains 10 points 334 # which straddle the boundary between UTM zones 53 and 54. 335 # They have been projected using a central meridian of 137.5 336 # degrees (the boundary is 138 so it is pretty much right 337 # 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 print 356 print 'Lat', lat 357 print 'Lon', lon 358 print 'Zone', zone 359 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 calculation 363 assert zone == -1 # Indicates non UTM projection 364 print 365 #assert num.allclose(x, easting) 366 #assert num.allclose(y, northing) 367 368 # Test that zone and meridian can't both be specified 369 try: 370 zone, easting, northing = redfearn(lat, lon, 371 zone=50, 372 central_meridian=137.5) 373 except: 374 pass 375 else: 376 msg = 'Should have raised exception' 377 raise Exception, msg 378 379 269 380 def test_convert_lats_longs(self): 270 381 -
branches/numpy/anuga/culvert_flows/Test_Culvert_Flat_Water_Lev.py
r6304 r6517 123 123 number_of_barrels=1, 124 124 update_interval=2, 125 log_file=True, 126 discharge_hydrograph=True, 125 127 verbose=True) 126 128 -
branches/numpy/anuga/culvert_flows/culvert_class.py
r6304 r6517 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 = None 206 208 207 209 … … 296 298 297 299 # Print some diagnostics to log if requested 298 if hasattr(self, 'log_filename'):300 if self.log_filename is not None: 299 301 s = 'Culvert Effective Length = %.2f m' %(self.length) 300 302 log_to_file(self.log_filename, s) … … 324 326 update = True 325 327 326 if hasattr(self, 'log_filename'): 328 329 if self.log_filename is not None: 327 330 s = '\nTime = %.2f, delta_t = %f' %(time, delta_t) 328 331 log_to_file(self.log_filename, s) … … 439 442 if self.verbose is True: 440 443 print msg 441 if hasattr(self, 'log_filename'): 444 445 if self.log_filename is not None: 442 446 log_to_file(self.log_filename, msg) 443 447 … … 524 528 print '%.2fs - WARNING: Flow is running uphill.' % time 525 529 526 if hasattr(self, 'log_filename'):530 if self.log_filename is not None: 527 531 s = 'Time=%.2f, inlet stage = %.2f, outlet stage = %.2f'\ 528 532 %(time, self.inlet.stage, self.outlet.stage) … … 573 577 msg += 'for culvert "%s"' % self.label 574 578 575 if hasattr(self, 'log_filename'):579 if self.log_filename is not None: 576 580 log_to_file(self.log_filename, msg) 577 581 else: 578 582 # User culvert routine 579 583 Q, barrel_velocity, culvert_outlet_depth =\ 580 self.culvert_routine(self, delta_total_energy, g) 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) 581 596 582 597 … … 601 616 602 617 603 618 # OBSOLETE (Except for momentum jet in Culvert_flow_energy) 604 619 class Culvert_flow_rating: 605 620 """Culvert flow - transfer water from one hole to another -
branches/numpy/anuga/culvert_flows/culvert_routines.py
r6128 r6517 8 8 9 9 #NOTE: 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 !!!! 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) 37 14 38 15 … … 40 17 41 18 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 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 49 42 """ 50 43 … … 53 46 from anuga.utilities.numerical_tools import safe_acos as acos 54 47 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 48 if inlet_depth > 0.01: 49 # Water has risen above inlet 50 51 if log_filename is not None: 52 s = 'Specific energy = %f m' % inlet_specific_energy 77 53 log_to_file(log_filename, s) 78 54 79 msg = 'Specific energy is negative'80 assert E>= 0.0, msg55 msg = 'Specific energy at inlet is negative' 56 assert inlet_specific_energy >= 0.0, msg 81 57 82 83 # Water has risen above inlet 84 if culvert.culvert_type == 'circle': 85 # Round culvert 58 if culvert_type == 'circle': 59 # Round culvert (use width as diameter) 60 diameter = culvert_width 86 61 87 62 # Calculate flows for inlet control 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) 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) 95 68 log_to_file(log_filename, s) 96 69 97 case = ''70 # FIXME(Ole): Are these functions really for inlet control? 98 71 if Q_inlet_unsubmerged < Q_inlet_submerged: 99 72 Q = Q_inlet_unsubmerged 100 73 101 alpha = acos(1 - inlet .depth/diameter)74 alpha = acos(1 - inlet_depth/diameter) 102 75 flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha)) 103 outlet_culvert_depth = inlet.depth 104 width = diameter*sin(alpha) 105 #perimeter = alpha*diameter 76 outlet_culvert_depth = inlet_depth 106 77 case = 'Inlet unsubmerged' 107 78 else: … … 109 80 flow_area = (diameter/2)**2 * pi 110 81 outlet_culvert_depth = diameter 111 width = diameter112 #perimeter = diameter113 82 case = 'Inlet submerged' 114 83 115 84 116 85 117 if delta_total_energy < E:86 if delta_total_energy < inlet_specific_energy: 118 87 # Calculate flows for outlet control 119 88 # Determine the depth at the outlet relative to the depth of flow in the Culvert 120 89 121 if outlet .depth > diameter: # The Outlet is Submerged90 if outlet_depth > diameter: # The Outlet is Submerged 122 91 outlet_culvert_depth=diameter 123 92 flow_area = (diameter/2)**2 * pi # Cross sectional area of flow in the culvert 124 93 perimeter = diameter * pi 125 width = diameter126 94 case = 'Outlet submerged' 127 elif outlet .depth==0.0:128 outlet_culvert_depth=inlet .depth# For purpose of calculation assume the outlet depth = the inlet depth129 alpha = acos(1 - inlet .depth/diameter)95 elif outlet_depth==0.0: 96 outlet_culvert_depth=inlet_depth # For purpose of calculation assume the outlet depth = the inlet depth 97 alpha = acos(1 - inlet_depth/diameter) 130 98 flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha)) 131 99 perimeter = alpha*diameter 132 width = diameter*sin(alpha)133 100 134 101 case = 'Outlet depth is zero' 135 102 else: # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity 136 outlet_culvert_depth=outlet .depth # For purpose of calculation assume the outlet depth = the inlet depth137 alpha = acos(1 - outlet .depth/diameter)103 outlet_culvert_depth=outlet_depth # For purpose of calculation assume the outlet depth = the inlet depth 104 alpha = acos(1 - outlet_depth/diameter) 138 105 flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha)) 139 106 perimeter = alpha*diameter 140 width = diameter*sin(alpha)141 107 case = 'Outlet is open channel flow' 142 108 143 109 hyd_rad = flow_area/perimeter 144 110 145 if hasattr(culvert, 'log_filename'):111 if log_filename is not None: 146 112 s = 'hydraulic radius at outlet = %f' %hyd_rad 147 113 log_to_file(log_filename, s) 148 114 149 115 # Outlet control velocity using tail water 150 culvert_velocity = sqrt(delta_total_energy/((sum_loss/2*g)+(manning**2* length)/hyd_rad**1.33333))116 culvert_velocity = sqrt(delta_total_energy/((sum_loss/2*g)+(manning**2*culvert_length)/hyd_rad**1.33333)) 151 117 Q_outlet_tailwater = flow_area * culvert_velocity 152 118 153 if hasattr(culvert, 'log_filename'):119 if log_filename is not None: 154 120 s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater 155 121 log_to_file(log_filename, s) … … 165 131 166 132 # Calculate flows for inlet control 167 height = culvert .height168 width = culvert .width133 height = culvert_height 134 width = culvert_width 169 135 170 Q_inlet_unsubmerged = 0.540*g**0.5*width* E**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged171 Q_inlet_submerged = 0.702*g**0.5*width*height**0.89* E**0.61 # Flow based on Inlet Ctrl Inlet Submerged172 173 if hasattr(culvert, 'log_filename'):136 Q_inlet_unsubmerged = 0.540*g**0.5*width*inlet_specific_energy**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged 137 Q_inlet_submerged = 0.702*g**0.5*width*height**0.89*inlet_specific_energy**0.61 # Flow based on Inlet Ctrl Inlet Submerged 138 139 if log_filename is not None: 174 140 s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' %(Q_inlet_unsubmerged, Q_inlet_submerged) 175 141 log_to_file(log_filename, s) 176 142 177 case = ''143 # FIXME(Ole): Are these functions really for inlet control? 178 144 if Q_inlet_unsubmerged < Q_inlet_submerged: 179 145 Q = Q_inlet_unsubmerged 180 flow_area = width*inlet.depth 181 outlet_culvert_depth = inlet.depth 182 #perimeter=(width+2.0*inlet.depth) 146 flow_area = width*inlet_depth 147 outlet_culvert_depth = inlet_depth 183 148 case = 'Inlet unsubmerged' 184 149 else: … … 189 154 case = 'Inlet submerged' 190 155 191 if delta_total_energy < E:156 if delta_total_energy < inlet_specific_energy: 192 157 # Calculate flows for outlet control 193 158 # Determine the depth at the outlet relative to the depth of flow in the Culvert 194 159 195 if outlet .depth > height: # The Outlet is Submerged160 if outlet_depth > height: # The Outlet is Submerged 196 161 outlet_culvert_depth=height 197 162 flow_area=width*height # Cross sectional area of flow in the culvert 198 163 perimeter=2.0*(width+height) 199 164 case = 'Outlet submerged' 200 elif outlet .depth==0.0:201 outlet_culvert_depth=inlet .depth # For purpose of calculation assume the outlet depth = the inlet depth202 flow_area=width*inlet .depth203 perimeter=(width+2.0*inlet .depth)165 elif outlet_depth==0.0: 166 outlet_culvert_depth=inlet_depth # For purpose of calculation assume the outlet depth = the inlet depth 167 flow_area=width*inlet_depth 168 perimeter=(width+2.0*inlet_depth) 204 169 case = 'Outlet depth is zero' 205 170 else: # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity 206 outlet_culvert_depth=outlet .depth207 flow_area=width*outlet .depth208 perimeter=(width+2.0*outlet .depth)171 outlet_culvert_depth=outlet_depth 172 flow_area=width*outlet_depth 173 perimeter=(width+2.0*outlet_depth) 209 174 case = 'Outlet is open channel flow' 210 175 211 176 hyd_rad = flow_area/perimeter 212 177 213 if hasattr(culvert, 'log_filename'):214 s = 'hydraulic radius at outlet = %f' % hyd_rad178 if log_filename is not None: 179 s = 'hydraulic radius at outlet = %f' % hyd_rad 215 180 log_to_file(log_filename, s) 216 181 217 182 # Outlet control velocity using tail water 218 culvert_velocity = sqrt(delta_total_energy/((sum_loss/2*g)+(manning**2* length)/hyd_rad**1.33333))183 culvert_velocity = sqrt(delta_total_energy/((sum_loss/2*g)+(manning**2*culvert_length)/hyd_rad**1.33333)) 219 184 Q_outlet_tailwater = flow_area * culvert_velocity 220 185 221 if hasattr(culvert, 'log_filename'):222 s = 'Q_outlet_tailwater = %.6f' % Q_outlet_tailwater186 if log_filename is not None: 187 s = 'Q_outlet_tailwater = %.6f' % Q_outlet_tailwater 223 188 log_to_file(log_filename, s) 224 189 Q = min(Q, Q_outlet_tailwater) … … 228 193 229 194 # Common code for circle and square geometries 230 if hasattr(culvert, 'log_filename'):231 log_to_file(log_filename, 'Case: "%s"' % case)195 if log_filename is not None: 196 log_to_file(log_filename, 'Case: "%s"' % case) 232 197 233 flow_rate_control=Q 234 235 if hasattr(culvert, 'log_filename'): 236 s = 'Flow Rate Control = %f' %flow_rate_control 198 if log_filename is not None: 199 s = 'Flow Rate Control = %f' % Q 237 200 log_to_file(log_filename, s) 238 201 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 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 246 206 log_to_file(log_filename, s) 247 207 … … 250 210 251 211 252 else: # inlet.depth < 0.01:212 else: # inlet_depth < 0.01: 253 213 Q = barrel_velocity = outlet_culvert_depth = 0.0 254 214 -
branches/numpy/anuga/fit_interpolate/fit.py
r6304 r6517 310 310 assert flag is True, msg 311 311 312 # FIXME(Ole): This is the message referred to in ticket:314 312 313 minx = min(self.mesh_boundary_polygon[:,0]) 313 314 maxx = max(self.mesh_boundary_polygon[:,0]) … … 317 318 msg += 'Mesh boundary extent is (%.f, %.f), (%.f, %.f)'\ 318 319 % (minx, maxx, miny, maxy) 320 raise RuntimeError, msg 319 321 320 322 … … 375 377 self.build_fit_subset(points, z, verbose=verbose) 376 378 379 # FIXME(Ole): I thought this test would make sense here 380 # See test_fitting_example_that_crashed_2 in test_shallow_water_domain.py 381 # Committed 11 March 2009 382 msg = 'Matrix AtA was not built' 383 assert self.AtA is not None, msg 384 385 #print 'Matrix was built OK' 386 377 387 378 388 point_coordinates = None … … 382 392 if point_coordinates is None: 383 393 if verbose: print 'Warning: no data points in fit' 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 394 msg = 'No interpolation matrix.' 395 assert self.AtA is not None, msg 396 assert self.Atz is not None 388 397 389 398 # FIXME (DSG) - do a message -
branches/numpy/anuga/geospatial_data/geospatial_data.py
r6441 r6517 1806 1806 points = G_small.get_data_points() 1807 1807 1808 #FIXME: Remove points outside boundary polygon 1809 # 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_poly 1815 # for i,point in enumerate(points): 1816 # if is_inside_polygon(point,boundary_poly, verbose=True): 1817 # new_points[i] = point 1818 # print"WOW",i,new_points[i] 1819 # points = new_points 1820 1808 1821 if verbose: print "Number of points in sample to compare: ", len(points) 1809 1822 -
branches/numpy/anuga/shallow_water/__init__.py
r6190 r6517 12 12 Transmissive_momentum_set_stage_boundary,\ 13 13 Dirichlet_discharge_boundary,\ 14 Field_boundary 14 Field_boundary,\ 15 Transmissive_stage_zero_momentum_boundary 15 16 16 17 -
branches/numpy/anuga/shallow_water/data_manager.py
r6481 r6517 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): 5865 #if type(times) is list or isinstance(times, num.ndarray): 5866 if isinstance(times, (list, num.ndarray)): 5866 5867 number_of_times = len(times) 5867 5868 times = ensure_numeric(times) … … 5930 5931 #outfile.variables[q+Write_sww.RANGE][1] = -max_float # Max 5931 5932 5932 if type(times) is list or isinstance(times, num.ndarray): 5933 outfile.variables['time'][:] = times #Store time relative 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 5934 5936 5935 5937 if verbose: … … 6265 6267 # This is being used to seperate one number from a list. 6266 6268 # what it is actually doing is sorting lists from numeric arrays. 6267 if type(times) is list or isinstance(times, num.ndarray): 6269 #if type(times) is list or isinstance(times, num.ndarray): 6270 if isinstance(times, (list, num.ndarray)): 6268 6271 number_of_times = len(times) 6269 6272 times = ensure_numeric(times) … … 6312 6315 outfile.variables[q + Write_sts.RANGE][1] = -max_float # Max 6313 6316 6314 if type(times) is list or isinstance(times, num.ndarray): 6317 #if type(times) is list or isinstance(times, num.ndarray): 6318 if isinstance(times, (list, num.ndarray)): 6315 6319 outfile.variables['time'][:] = times #Store time relative 6316 6320 -
branches/numpy/anuga/shallow_water/test_shallow_water_domain.py
r6451 r6517 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 1674 points, vertices, boundary = rectangular_cross(10, 10) # Basic mesh 1675 domain = Domain(points, vertices, boundary) # Create domain 1676 domain.set_default_order(1) 1676 1677 domain.set_quantities_to_be_stored(None) 1677 1678 domain.set_maximum_allowed_speed(100) #FIXME (Ole): try to remove this -
branches/numpy/anuga/test_all.py
r6304 r6517 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.