Changeset 7176
- Timestamp:
- Jun 10, 2009, 5:28:35 PM (16 years ago)
- Location:
- branches/numpy/anuga
- Files:
-
- 12 edited
-
abstract_2d_finite_volumes/domain.py (modified) (1 diff)
-
abstract_2d_finite_volumes/generic_boundary_conditions.py (modified) (1 diff)
-
abstract_2d_finite_volumes/quantity.py (modified) (1 diff)
-
config.py (modified) (3 diffs)
-
coordinate_transforms/geo_reference.py (modified) (1 diff)
-
culvert_flows/culvert_routines.py (modified) (9 diffs)
-
geospatial_data/test_geospatial_data.py (modified) (1 diff)
-
load_mesh/loadASCII.py (modified) (1 diff)
-
mesh_engine/test_generate_mesh.py (modified) (1 diff)
-
shallow_water/data_manager.py (modified) (4 diffs)
-
shallow_water/shallow_water_domain.py (modified) (1 diff)
-
utilities/test_polygon.py (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
branches/numpy/anuga/abstract_2d_finite_volumes/domain.py
r7035 r7176 936 936 k = self.k = triangle_id 937 937 938 x, y = self.get_centroid_coordinates( )[k]938 x, y = self.get_centroid_coordinates(absolute=True)[k] 939 939 radius = self.get_radii()[k] 940 940 area = self.get_areas()[k] -
branches/numpy/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py
r7037 r7176 103 103 verbose=False): 104 104 Boundary.__init__(self) 105 105 106 self.default_boundary = default_boundary 106 107 self.default_boundary_invoked = False # Flag -
branches/numpy/anuga/abstract_2d_finite_volumes/quantity.py
r7038 r7176 1178 1178 # Edges have already been deprecated in set_values, see changeset:5521, 1179 1179 # but *might* be useful in get_values. Any thoughts anyone? 1180 # YES (Ole): Edge values are necessary for volumetric balance check 1180 # YES (Ole): Edge values are necessary for volumetric balance check and inflow boundary. Keep them. 1181 1181 1182 1182 if location not in ['vertices', 'centroids', -
branches/numpy/anuga/config.py
r6553 r7176 8 8 9 9 ################################################################################ 10 # numerical constants10 # Numerical constants 11 11 ################################################################################ 12 12 … … 17 17 # least squares fitting 18 18 single_precision = 1.0e-6 # Smallest single precision number 19 velocity_protection = 1.0e-6 19 velocity_protection = 1.0e-6 # Used to compute velocity from momentum 20 # See section 7.4 on Flux limiting 21 # in the user manual 22 20 23 21 24 ################################################################################ … … 26 29 version_filename = 'stored_version_info.py' 27 30 default_datadir = '.' 28 time_format = '%d/%m/%y %H:%M:%S' 29 umask = 002 # Controls file and directory permission created by anuga 31 time_format = '%d/%m/%y %H:%M:%S' # Used with timefile2netcdf 32 umask = 002 # Controls file and directory permission created by anuga (UNIX) 30 33 default_boundary_tag = 'exterior' 31 34 -
branches/numpy/anuga/coordinate_transforms/geo_reference.py
r7061 r7176 142 142 self.zone = int(infile.zone[0]) 143 143 144 ## # Fix some assertion failures145 ## if isinstance(self.zone, num.ndarray) and self.zone.shape == ():146 ## self.zone = self.zone[0]147 ## if (isinstance(self.xllcorner, num.ndarray) and148 ## self.xllcorner.shape == ()):149 ## self.xllcorner = self.xllcorner[0]150 ## if (isinstance(self.yllcorner, num.ndarray) and151 ## self.yllcorner.shape == ()):152 ## self.yllcorner = self.yllcorner[0]153 ##154 ## assert (self.xllcorner.dtype.kind in num.typecodes['Float'] or155 ## self.xllcorner.dtype.kind in num.typecodes['Integer'])156 ## assert (self.yllcorner.dtype.kind in num.typecodes['Float'] or157 ## self.yllcorner.dtype.kind in num.typecodes['Integer'])158 ## assert (self.zone.dtype.kind in num.typecodes['Integer'])159 160 144 try: 161 145 self.false_easting = int(infile.false_easting[0]) -
branches/numpy/anuga/culvert_flows/culvert_routines.py
r7035 r7176 1 1 """Collection of culvert routines for use with Culvert_flow in culvert_class 2 3 This module holds various routines to determine FLOW through CULVERTS and SIMPLE BRIDGES 4 5 6 2 7 3 8 Usage: … … 33 38 log_filename=None): 34 39 35 """Boyd's generalisation of the US department of transportation culvert 36 model 37 38 The quantity of flow passing through a culvert is controlled by many factors 39 It could be that the culvert is controlled by the inlet only, with it 40 being unsubmerged this is effectively equivalent to the weir Equation 41 Else the culvert could be controlled by the inlet, with it being 42 submerged, this is effectively the Orifice Equation 43 Else it may be controlled by down stream conditions where depending on 44 the down stream depth, the momentum in the culvert etc. flow is controlled 40 """Boyd's generalisation of the US department of transportation culvert methods 41 42 WARNING THIS IS A SIMPLISTIC APPROACH and OUTLET VELOCITIES ARE LIMITED TO EITHER 43 FULL PIPE OR CRITICAL DEPTH ONLY 44 For Supercritical flow this is UNDERESTIMATING the Outlet Velocity 45 The obtain the CORRECT velocity requires an iteration of Depth to Establish 46 the Normal Depth of flow in the pipe. 47 48 It is proposed to provide this in a seperate routine called 49 boyd_generalised_culvert_model_complex 50 51 The Boyd Method is based on methods described by the following: 52 1. 53 US Dept. Transportation Federal Highway Administration (1965) 54 Hydraulic Chart for Selection of Highway Culverts. 55 Hydraulic Engineering Circular No. 5 US Government Printing 56 2. 57 US Dept. Transportation Federal Highway Administration (1972) 58 Capacity charts for the Hydraulic design of highway culverts. 59 Hydraulic Engineering Circular No. 10 US Government Printing 60 These documents provide around 60 charts for various configurations of culverts and inlets. 61 62 Note these documents have been superceded by: 63 2005 Hydraulic Design of Highway Culverts, Hydraulic Design Series No. 5 (HDS-5), 64 Which combines culvert design information previously contained in Hydraulic Engineering Circulars 65 (HEC) No. 5, No. 10, and No. 13 with hydrologic, storage routing, and special culvert design information. 66 HEC-5 provides 20 Charts 67 HEC-10 Provides an additional 36 Charts 68 HEC-13 Discusses the Design of improved more efficient inlets 69 HDS-5 Provides 60 sets of Charts 70 71 In 1985 Professor Michael Boyd Published "Head-Discharge Relations for Culverts", and in 72 1987 published "Generalised Head Discharge Equations for Culverts". 73 These papers reviewed the previous work by the US DOT and provided a simplistic approach for 3 configurations. 74 75 It may be possible to extend the same approach for additional charts in the original work, but to date this has not been done. 76 The additional charts cover a range of culvert shapes and inlet configurations 77 78 45 79 """ 46 80 … … 49 83 from anuga.utilities.numerical_tools import safe_acos as acos 50 84 51 85 local_debug ='false' 52 86 if inlet_depth > 0.1: #this value was 0.01: 87 if local_debug =='true': 88 print 'Specific E & Deltat Tot E = ',inlet_specific_energy,delta_total_energy 89 print 'culvert typ = ',culvert_type 53 90 # Water has risen above inlet 54 91 … … 61 98 62 99 63 if culvert_type == 'circle': 64 # Round culvert (use width as diameter) 65 diameter = culvert_width 66 100 if culvert_type =='circle': 101 # Round culvert (use height as diameter) 102 diameter = culvert_height 103 104 """ 105 For a circular pipe the Boyd method reviews 3 conditions 106 1. Whether the Pipe Inlet is Unsubmerged (acting as weir flow into the inlet) 107 2. Whether the Pipe Inlet is Fully Submerged (acting as an Orifice) 108 3. Whether the energy loss in the pipe results in the Pipe being controlled by Channel Flow of the Pipe 109 110 For these conditions we also would like to assess the pipe flow characteristics as it leaves the pipe 111 """ 112 67 113 # Calculate flows for inlet control 68 114 Q_inlet_unsubmerged = 0.421*g**0.5*diameter**0.87*inlet_specific_energy**1.63 # Inlet Ctrl Inlet Unsubmerged 69 115 Q_inlet_submerged = 0.530*g**0.5*diameter**1.87*inlet_specific_energy**0.63 # Inlet Ctrl Inlet Submerged 116 # Note for to SUBMERGED TO OCCUR inlet_specific_energy should be > 1.2 x diameter.... Should Check !!! 70 117 71 118 if log_filename is not None: 72 119 s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' % (Q_inlet_unsubmerged, Q_inlet_submerged) 73 120 log_to_file(log_filename, s) 74 75 76 # FIXME(Ole): Are these functions really for inlet control? 77 if Q_inlet_unsubmerged < Q_inlet_submerged: 78 Q = Q_inlet_unsubmerged 79 alpha = acos(1 - inlet_depth/diameter) 80 flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha)) 81 outlet_culvert_depth = inlet_depth 82 case = 'Inlet unsubmerged' 83 else: 84 Q = Q_inlet_submerged 85 flow_area = (diameter/2)**2 * pi 86 outlet_culvert_depth = diameter 87 case = 'Inlet submerged' 88 89 90 91 if delta_total_energy < inlet_specific_energy: 121 Q = min(Q_inlet_unsubmerged, Q_inlet_submerged) 122 123 # THE LOWEST Value will Control Calcs From here 124 # Calculate Critical Depth Based on the Adopted Flow as an Estimate 125 dcrit1 = diameter/1.26*(Q/g**0.5*diameter**2.5)**(1/3.75) 126 dcrit2 = diameter/0.95*(Q/g**0.5*diameter**2.5)**(1/1.95) 127 # From Boyd Paper ESTIMATE of Dcrit has 2 criteria as 128 if dcrit1/diameter > 0.85: 129 outlet_culvert_depth = dcrit2 130 else: 131 outlet_culvert_depth = dcrit1 132 #outlet_culvert_depth = min(outlet_culvert_depth, diameter) 133 # Now determine Hydraulic Radius Parameters Area & Wetted Perimeter 134 if outlet_culvert_depth >= diameter: 135 outlet_culvert_depth = diameter # Once again the pipe is flowing full not partfull 136 flow_area = (diameter/2)**2 * pi # Cross sectional area of flow in the culvert 137 perimeter = diameter * pi 138 flow_width= diameter 139 case = 'Inlet CTRL Outlet submerged Circular PIPE FULL' 140 if local_debug =='true': 141 print 'Inlet CTRL Outlet submerged Circular PIPE FULL' 142 else: 143 #alpha = acos(1 - outlet_culvert_depth/diameter) # Where did this Come From ????/ 144 alpha = acos(1-2*outlet_culvert_depth/diameter)*2 145 #flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha)) # Pipe is Running Partly Full at the INLET WHRE did this Come From ????? 146 flow_area = diameter**2/8*(alpha - sin(alpha)) # Equation from GIECK 5th Ed. Pg. B3 147 flow_width= diameter*sin(alpha/2.0) 148 perimeter = alpha*diameter/2.0 149 case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth' 150 if local_debug =='true': 151 print 'INLET CTRL Culvert is open channel flow we will for now assume critical depth' 152 print 'Q Outlet Depth and ALPHA = ',Q,' ',outlet_culvert_depth,' ',alpha 153 if delta_total_energy < inlet_specific_energy: # OUTLET CONTROL !!!! 92 154 # Calculate flows for outlet control 93 155 94 156 # Determine the depth at the outlet relative to the depth of flow in the Culvert 95 if outlet_depth > diameter: # The Outlet is Submerged157 if outlet_depth > diameter: # Outlet is submerged Assume the end of the Pipe is flowing FULL 96 158 outlet_culvert_depth=diameter 97 159 flow_area = (diameter/2)**2 * pi # Cross sectional area of flow in the culvert 98 160 perimeter = diameter * pi 161 flow_width= diameter 99 162 case = 'Outlet submerged' 100 elif outlet_depth==0.0: 101 outlet_culvert_depth=inlet_depth # For purpose of calculation assume the outlet depth = the inlet depth 102 alpha = acos(1 - inlet_depth/diameter) 103 flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha)) 104 perimeter = alpha*diameter 105 106 case = 'Outlet depth is zero' 107 else: # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity 108 outlet_culvert_depth=outlet_depth # For purpose of calculation assume the outlet depth = the inlet depth 109 alpha = acos(1 - outlet_depth/diameter) 110 flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha)) 111 perimeter = alpha*diameter 112 case = 'Outlet is open channel flow' 113 114 hyd_rad = flow_area/perimeter 115 116 if log_filename is not None: 117 s = 'hydraulic radius at outlet = %f' %hyd_rad 118 log_to_file(log_filename, s) 119 120 # Outlet control velocity using tail water 121 culvert_velocity = sqrt(delta_total_energy/((sum_loss/2/g)+(manning**2*culvert_length)/hyd_rad**1.33333)) 122 Q_outlet_tailwater = flow_area * culvert_velocity 123 124 if log_filename is not None: 125 s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater 126 log_to_file(log_filename, s) 127 128 Q = min(Q, Q_outlet_tailwater) 129 else: 130 pass 163 if local_debug =='true': 164 print 'Outlet submerged' 165 else: # Culvert running PART FULL for PART OF ITS LENGTH Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity 166 # IF outlet_depth < diameter 167 dcrit1 = diameter/1.26*(Q/g**0.5*diameter**2.5)**(1/3.75) 168 dcrit2 = diameter/0.95*(Q/g**0.5*diameter**2.5)**(1/1.95) 169 if dcrit1/diameter >0.85: 170 outlet_culvert_depth= dcrit2 171 else: 172 outlet_culvert_depth = dcrit1 173 if outlet_culvert_depth > diameter: 174 outlet_culvert_depth = diameter # Once again the pipe is flowing full not partfull 175 flow_area = (diameter/2)**2 * pi # Cross sectional area of flow in the culvert 176 perimeter = diameter * pi 177 flow_width= diameter 178 case = 'Outlet unsubmerged PIPE FULL' 179 if local_debug =='true': 180 print 'Outlet unsubmerged PIPE FULL' 181 else: 182 alpha = acos(1-2*outlet_culvert_depth/diameter)*2 183 flow_area = diameter**2/8*(alpha - sin(alpha)) # Equation from GIECK 5th Ed. Pg. B3 184 flow_width= diameter*sin(alpha/2.0) 185 perimeter = alpha*diameter/2.0 186 case = 'Outlet is open channel flow we will for now assume critical depth' 187 if local_debug =='true': 188 print 'Q Outlet Depth and ALPHA = ',Q,' ',outlet_culvert_depth,' ',alpha 189 print 'Outlet is open channel flow we will for now assume critical depth' 190 if local_debug =='true': 191 print 'FLOW AREA = ',flow_area 192 print 'PERIMETER = ',perimeter 193 print 'Q Interim = ',Q 194 hyd_rad = flow_area/perimeter 195 196 if log_filename is not None: 197 s = 'hydraulic radius at outlet = %f' %hyd_rad 198 log_to_file(log_filename, s) 199 200 # Outlet control velocity using tail water 201 if local_debug =='true': 202 print 'GOT IT ALL CALCULATING Velocity' 203 print 'HydRad = ',hyd_rad 204 culvert_velocity = sqrt(delta_total_energy/((sum_loss/2/g)+(manning**2*culvert_length)/hyd_rad**1.33333)) 205 Q_outlet_tailwater = flow_area * culvert_velocity 206 if local_debug =='true': 207 print 'VELOCITY = ',culvert_velocity 208 print 'Outlet Ctrl Q = ',Q_outlet_tailwater 209 if log_filename is not None: 210 s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater 211 log_to_file(log_filename, s) 212 Q = min(Q, Q_outlet_tailwater) 213 if local_debug =='true': 214 print ('%s,%.3f,%.3f' %('dcrit 1 , dcit2 =',dcrit1,dcrit2)) 215 print ('%s,%.3f,%.3f,%.3f' %('Q and Velocity and Depth=',Q,culvert_velocity,outlet_culvert_depth)) 216 217 else: 218 pass 131 219 #FIXME(Ole): What about inlet control? 132 133 134 else: 135 # Box culvert (rectangle or square) 220 # ==== END OF CODE BLOCK FOR "IF" CIRCULAR PIPE 221 222 #else.... 223 if culvert_type =='box': 224 if local_debug =='true': 225 print 'BOX CULVERT' 226 # Box culvert (rectangle or square) ======================================================================================================================== 136 227 137 228 # Calculate flows for inlet control 138 229 height = culvert_height 139 width = culvert_width 230 width = culvert_width 231 flow_width=culvert_width 140 232 141 233 Q_inlet_unsubmerged = 0.540*g**0.5*width*inlet_specific_energy**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged … … 146 238 log_to_file(log_filename, s) 147 239 148 149 240 # FIXME(Ole): Are these functions really for inlet control? 150 241 if Q_inlet_unsubmerged < Q_inlet_submerged: 151 242 Q = Q_inlet_unsubmerged 152 flow_area = width*inlet_depth 153 outlet_culvert_depth = inlet_depth 154 case = 'Inlet unsubmerged' 243 dcrit = (Q**2/g/width**2)**0.333333 244 if dcrit > height: 245 dcrit = height 246 flow_area = width*dcrit 247 outlet_culvert_depth = dcrit 248 case = 'Inlet unsubmerged Box Acts as Weir' 155 249 else: 156 250 Q = Q_inlet_submerged 157 251 flow_area = width*height 158 252 outlet_culvert_depth = height 159 case = 'Inlet submerged' 160 253 case = 'Inlet submerged Box Acts as Orifice' 254 255 dcrit = (Q**2/g/width**2)**0.333333 256 257 outlet_culvert_depth = dcrit 258 if outlet_culvert_depth > height: 259 outlet_culvert_depth = height # Once again the pipe is flowing full not partfull 260 flow_area = width*height # Cross sectional area of flow in the culvert 261 perimeter = 2*(width+height) 262 case = 'Inlet CTRL Outlet unsubmerged PIPE PART FULL' 263 else: 264 flow_area = width * outlet_culvert_depth 265 perimeter = width+2*outlet_culvert_depth 266 case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth' 267 161 268 if delta_total_energy < inlet_specific_energy: 162 269 # Calculate flows for outlet control … … 168 275 perimeter=2.0*(width+height) 169 276 case = 'Outlet submerged' 170 elif outlet_depth==0.0:171 outlet_culvert_depth=inlet_depth # For purpose of calculation assume the outlet depth = the inlet depth172 flow_area=width*inlet_depth173 perimeter=(width+2.0*inlet_depth)174 case = 'Outlet depth is zero'175 277 else: # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity 176 outlet_culvert_depth=outlet_depth 177 flow_area=width*outlet_depth 178 perimeter=(width+2.0*outlet_depth) 179 case = 'Outlet is open channel flow' 278 dcrit = (Q**2/g/width**2)**0.333333 279 outlet_culvert_depth=dcrit # For purpose of calculation assume the outlet depth = Critical Depth 280 if outlet_culvert_depth > height: 281 outlet_culvert_depth=height 282 flow_area=width*height 283 perimeter=2.0*(width+height) 284 case = 'Outlet is Flowing Full' 285 else: 286 flow_area=width*outlet_culvert_depth 287 perimeter=(width+2.0*outlet_culvert_depth) 288 case = 'Outlet is open channel flow' 180 289 181 290 hyd_rad = flow_area/perimeter 182 291 183 292 if log_filename is not None: 184 293 s = 'hydraulic radius at outlet = %f' % hyd_rad … … 196 305 pass 197 306 #FIXME(Ole): What about inlet control? 307 # ==== END OF CODE BLOCK FOR "IF" BOX 198 308 199 309 200 # Common code for circle and square geometries310 # Common code for circle and box geometries ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 201 311 if log_filename is not None: 202 312 log_to_file(log_filename, 'Case: "%s"' % case) … … 206 316 log_to_file(log_filename, s) 207 317 208 209 culv_froude=sqrt(Q**2*width/(g*flow_area**3)) 210 318 culv_froude=sqrt(Q**2*flow_width/(g*flow_area**3)) 319 if local_debug =='true': 320 print 'FLOW AREA = ',flow_area 321 print 'PERIMETER = ',perimeter 322 print 'Q final = ',Q 323 print 'FROUDE = ',culv_froude 211 324 if log_filename is not None: 212 325 s = 'Froude in Culvert = %f' % culv_froude … … 216 329 barrel_velocity = Q/(flow_area + velocity_protection/flow_area) 217 330 331 # END CODE BLOCK for DEPTH > Required depth for CULVERT Flow 218 332 219 333 else: # inlet_depth < 0.01: -
branches/numpy/anuga/geospatial_data/test_geospatial_data.py
r7035 r7176 112 112 assert num.allclose(V, [79.4, -7]) 113 113 114 # FIXME: use failUnlessRaises() 114 115 try: 115 116 V = G.get_attributes('hdnoatedu') #Invalid -
branches/numpy/anuga/load_mesh/loadASCII.py
r7035 r7176 700 700 ('num_of_segments', 'num_of_segment_ends')) 701 701 outfile.variables['segments'][:] = mesh['segments'] 702 if mesh['segment_tags'].shape[0] > 0:702 if (mesh['segment_tags'].shape[1] > 0): 703 703 outfile.createDimension('num_of_segment_tag_chars', 704 704 mesh['segment_tags'].shape[1]) -
branches/numpy/anuga/mesh_engine/test_generate_mesh.py
r6410 r7176 412 412 'Failed!') 413 413 414 correct = num.array(segattlist )414 correct = num.array(segattlist, num.Int) 415 415 self.failUnless(num.allclose(data['generatedsegmentmarkerlist'].flat, 416 416 correct.flat), -
branches/numpy/anuga/shallow_water/data_manager.py
r7035 r7176 566 566 # Define a zero vector of same size and type as A 567 567 # for use with momenta 568 null = num.zeros(num.size(A), A.dtype.char) #??#568 null = num.zeros(num.size(A), A.dtype.char) 569 569 570 570 # Get xmomentum where depth exceeds minimum_storable_height … … 1433 1433 for name in infile.variables: 1434 1434 var = infile.variables[name] 1435 outfile.createVariable(name, var.dtype.char, var.dimensions) #??#1435 outfile.createVariable(name, var.dtype.char, var.dimensions) 1436 1436 1437 1437 # Copy the static variables … … 2193 2193 swwfile = basename_in + '.sww' 2194 2194 demfile = basename_out + '.' + format 2195 # Note the use of a .ers extension is optional (write_ermapper_grid will2196 # deal with either option2197 2195 2198 2196 # Read sww file … … 2697 2695 # Interpolate using quantity values 2698 2696 if verbose: print 'Interpolating' 2699 interpolated_values = interp.interpolate(q, data_points).flatten() #????#2697 interpolated_values = interp.interpolate(q, data_points).flatten() 2700 2698 2701 2699 if verbose: -
branches/numpy/anuga/shallow_water/shallow_water_domain.py
r7035 r7176 1462 1462 1463 1463 1464 ## 1465 # @brief A class for a 'field' boundary. 1466 # @note Inherits from Boundary. 1464 class Inflow_boundary(Boundary): 1465 """Apply given flow in m^3/s to boundary segment. 1466 Depth and momentum is derived using Manning's formula. 1467 1468 Underlying domain must be specified when boundary is instantiated 1469 """ 1470 1471 # FIXME (Ole): This is work in progress and definitely not finished. 1472 # The associated test has been disabled 1473 1474 def __init__(self, domain=None, rate=0.0): 1475 Boundary.__init__(self) 1476 1477 if domain is None: 1478 msg = 'Domain must be specified for ' 1479 msg += 'Inflow boundary' 1480 raise Exception, msg 1481 1482 self.domain = domain 1483 1484 # FIXME(Ole): Allow rate to be time dependent as well 1485 self.rate = rate 1486 self.tag = None # Placeholder for tag associated with this object. 1487 1488 def __repr__(self): 1489 return 'Inflow_boundary(%s)' %self.domain 1490 1491 def evaluate(self, vol_id, edge_id): 1492 """Apply inflow rate at each edge of this boundary 1493 """ 1494 1495 # First find all segments having the same tag is vol_id, edge_id 1496 # This will be done the first time evaluate is called. 1497 if self.tag is None: 1498 boundary = self.domain.boundary 1499 self.tag = boundary[(vol_id, edge_id)] 1500 1501 # Find total length of boundary with this tag 1502 length = 0.0 1503 for v_id, e_id in boundary: 1504 if self.tag == boundary[(v_id, e_id)]: 1505 length += self.domain.mesh.get_edgelength(v_id, e_id) 1506 1507 self.length = length 1508 self.average_momentum = self.rate/length 1509 1510 1511 # Average momentum has now been established across this boundary 1512 # Compute momentum in the inward normal direction 1513 1514 inward_normal = -self.domain.mesh.get_normal(vol_id, edge_id) 1515 xmomentum, ymomentum = self.average_momentum * inward_normal 1516 1517 # Compute depth based on Manning's formula v = 1/n h^{2/3} sqrt(S) 1518 # Where v is velocity, n is manning's coefficient, h is depth and S is the slope into the domain. 1519 # Let mu be the momentum (vh), then this equation becomes: mu = 1/n h^{5/3} sqrt(S) 1520 # from which we can isolate depth to get 1521 # h = (mu n/sqrt(S) )^{3/5} 1522 1523 slope = 0 # get gradient for this triangle dot normal 1524 1525 # get manning coef from this triangle 1526 friction = self.domain.get_quantity('friction').get_values(location='edges', 1527 indices=[vol_id])[0] 1528 mannings_n = friction[edge_id] 1529 1530 if slope > epsilon and mannings_n > epsilon: 1531 depth = pow(self.average_momentum * mannings_n/math.sqrt(slope), 3.0/5) 1532 else: 1533 depth = 1.0 1534 1535 # Elevation on this edge 1536 1537 z = self.domain.get_quantity('elevation').get_values(location='edges', 1538 indices=[vol_id])[0] 1539 elevation = z[edge_id] 1540 1541 # Assign conserved quantities and return 1542 q = num.array([elevation + depth, xmomentum, ymomentum], num.Float) 1543 return q 1544 1545 1546 1547 1548 1549 1467 1550 class Field_boundary(Boundary): 1468 1551 """Set boundary from given field represented in an sww file containing -
branches/numpy/anuga/utilities/test_polygon.py
r7038 r7176 178 178 179 179 def test_inside_polygon_main(self): 180 """test_is_inside_polygon_quick 181 182 Test fast version of of is_inside_polygon 183 """ 184 180 185 # Simplest case: Polygon is the unit square 181 186 polygon = [[0,0], [1,0], [1,1], [0,1]]
Note: See TracChangeset
for help on using the changeset viewer.
