Changeset 5726
 Timestamp:
 Sep 3, 2008, 4:44:23 PM (15 years ago)
 Location:
 anuga_core/source/anuga/shallow_water
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

anuga_core/source/anuga/shallow_water/data_manager.py
r5723 r5726 6466 6466 6467 6467 6468 def get_energy_through_cross_section(filename, 6469 polyline, 6470 kind = 'total', 6471 verbose=False): 6472 """Obtain flow (m^3/s) perpendicular to specified cross section. 6473 6474 Inputs: 6475 filename: Name of sww file 6476 polyline: Representation of desired cross section  it may contain 6477 multiple sections allowing for complex shapes. Assume 6478 absolute UTM coordinates. 6479 Format [[x0, y0], [x1, y1], ...] 6480 kind: Select which energy to compute. 6481 Options are 'specific' and 'total' (default) 6482 6483 Output: 6484 time: All stored times in sww file 6485 Q: Average energy [m] across given segments for all stored times. 6486 6487 The average velocity is computed for each triangle intersected by the polyline 6488 and averaged weigted by segment lengths. 6489 6490 The typical usage of this function would be to get average energy of flow in a channel, 6491 and the polyline would then be a cross section perpendicular to the flow. 6492 6493 6494 #FIXME (Ole)  need name for this energy reflecting that its dimension is [m]. 6495 """ 6496 6497 from anuga.config import g, epsilon, velocity_protection as h0 6498 6499 quantity_names =['elevation', 6500 'stage', 6501 'xmomentum', 6502 'ymomentum'] 6503 6504 6505 # Get values for quantities at each midpoint of poly line from sww file 6506 X = get_interpolated_quantities_at_polyline_midpoints(filename, 6507 quantity_names=quantity_names, 6508 polyline=polyline, 6509 verbose=verbose) 6510 segments, interpolation_function = X 6511 6512 6513 # Get vectors for time and interpolation_points 6514 time = interpolation_function.time 6515 interpolation_points = interpolation_function.interpolation_points 6516 6517 if verbose: print 'Computing %s energy' %kind 6518 6519 # Compute total length of polyline for use with weighted averages 6520 total_line_length = 0.0 6521 for segment in segments: 6522 total_line_length += segment.length 6523 6524 # Compute energy 6525 E = [] 6526 for t in time: 6527 average_energy=0.0 6528 for i, p in enumerate(interpolation_points): 6529 elevation, stage, uh, vh = interpolation_function(t, point_id=i) 6530 6531 # Depth 6532 h = depth = stageelevation 6533 6534 # Average velocity across this segment 6535 if h > epsilon: 6536 # Use protection against degenerate velocities 6537 u = uh/(h + h0/h) 6538 v = vh/(h + h0/h) 6539 else: 6540 u = v = 0.0 6541 6542 speed_squared = u*u + v*v 6543 kinetic_energy = 0.5*speed_squared/g 6544 6545 if kind == 'specific': 6546 segment_energy = depth + kinetic_energy 6547 elif kind == 'total': 6548 segment_energy = stage + kinetic_energy 6549 else: 6550 msg = 'Energy kind must be either "specific" or "total".' 6551 msg += ' I got %s' %kind 6552 6553 6554 # Add to weighted average 6555 weigth = segments[i].length/total_line_length 6556 average_energy += segment_energy*weigth 6557 6558 6559 # Store energy at this timestep 6560 E.append(average_energy) 6561 6562 6563 return time, E 6564 6565 6468 6566 6469 6567 
anuga_core/source/anuga/shallow_water/test_data_manager.py
r5672 r5726 17 17 from anuga.shallow_water import * 18 18 from anuga.shallow_water.data_manager import * 19 from anuga.config import epsilon 19 from anuga.config import epsilon, g 20 20 from anuga.utilities.anuga_exceptions import ANUGAError 21 21 from anuga.utilities.numerical_tools import ensure_numeric … … 10198 10198 The specifics are 10199 10199 u = 2 m/s 10200 h = 1m10200 h = 2 m 10201 10201 w = 5 m (width of channel) 10202 10202 10203 q = u*h*w = 10 m^3/s10204 10205 10206 This run tries it with georeferencing 10203 q = u*h*w = 20 m^3/s 10204 10205 10206 This run tries it with georeferencing and with elevation = 1 10207 10207 10208 10208 """ … … 10236 10236 domain.smooth = True 10237 10237 10238 h = 1.0 10238 e = 1.0 10239 w = 1.0 10240 h = we 10239 10241 u = 2.0 10240 10242 uh = u*h 10241 10243 10242 10244 Br = Reflective_boundary(domain) # Side walls 10243 Bd = Dirichlet_boundary([ h, uh, 0]) # 2 m/s across the 5 m inlet:10244 10245 10246 10247 10248 domain.set_quantity('elevation', 0.0)10249 domain.set_quantity('stage', h)10245 Bd = Dirichlet_boundary([w, uh, 0]) # 2 m/s across the 5 m inlet: 10246 10247 10248 10249 10250 domain.set_quantity('elevation', e) 10251 domain.set_quantity('stage', w) 10250 10252 domain.set_quantity('xmomentum', uh) 10251 10253 domain.set_boundary( {'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br}) … … 10268 10270 for t in range(t_end+1): 10269 10271 for i in range(3): 10270 assert allclose(f(t, i), [1, 2, 0]) 10272 #print i, t, f(t, i) 10273 assert allclose(f(t, i), [w, uh, 0]) 10271 10274 10272 10275 … … 10284 10287 10285 10288 10289 10290 def test_get_energy_through_cross_section(self): 10291 """test_get_energy_through_cross_section(self): 10292 10293 Test that the specific and total energy through a cross section can be 10294 correctly obtained from an sww file. 10295 10296 This test creates a flat bed with a known flow through it and tests 10297 that the function correctly returns the expected energies. 10298 10299 The specifics are 10300 u = 2 m/s 10301 h = 1 m 10302 w = 5 m (width of channel) 10303 10304 q = u*h*w = 10 m^3/s 10305 Es = h + 0.5*v*v/g # Specific energy head [m] 10306 Et = w + 0.5*v*v/g # Total energy head [m] 10307 10308 10309 This test uses georeferencing 10310 10311 """ 10312 10313 import time, os 10314 from Numeric import array, zeros, allclose, Float, concatenate 10315 from Scientific.IO.NetCDF import NetCDFFile 10316 10317 # Setup 10318 from mesh_factory import rectangular 10319 10320 # Create basic mesh (20m x 3m) 10321 width = 3 10322 length = 20 10323 t_end = 1 10324 points, vertices, boundary = rectangular(length, width, 10325 length, width) 10326 10327 # Create shallow water domain 10328 domain = Domain(points, vertices, boundary, 10329 geo_reference = Geo_reference(56,308500,6189000)) 10330 10331 domain.default_order = 2 10332 domain.set_minimum_storable_height(0.01) 10333 10334 domain.set_name('flowtest') 10335 swwfile = domain.get_name() + '.sww' 10336 10337 domain.set_datadir('.') 10338 domain.format = 'sww' 10339 domain.smooth = True 10340 10341 e = 1.0 10342 w = 1.0 10343 h = we 10344 u = 2.0 10345 uh = u*h 10346 10347 Br = Reflective_boundary(domain) # Side walls 10348 Bd = Dirichlet_boundary([w, uh, 0]) # 2 m/s across the 5 m inlet: 10349 10350 10351 domain.set_quantity('elevation', e) 10352 domain.set_quantity('stage', w) 10353 domain.set_quantity('xmomentum', uh) 10354 domain.set_boundary( {'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br}) 10355 10356 for t in domain.evolve(yieldstep=1, finaltime = t_end): 10357 pass 10358 10359 # Check that momentum is as it should be in the interior 10360 10361 I = [[0, width/2.], 10362 [length/2., width/2.], 10363 [length, width/2.]] 10364 10365 I = domain.geo_reference.get_absolute(I) 10366 f = file_function(swwfile, 10367 quantities=['stage', 'xmomentum', 'ymomentum'], 10368 interpolation_points=I, 10369 verbose=False) 10370 10371 for t in range(t_end+1): 10372 for i in range(3): 10373 #print i, t, f(t, i) 10374 assert allclose(f(t, i), [w, uh, 0]) 10375 10376 10377 # Check energies through the middle 10378 for i in range(5): 10379 x = length/2. + i*0.23674563 # Arbitrary 10380 cross_section = [[x, 0], [x, width]] 10381 10382 cross_section = domain.geo_reference.get_absolute(cross_section) 10383 10384 time, Es = get_energy_through_cross_section(swwfile, 10385 cross_section, 10386 kind='specific', 10387 verbose=False) 10388 assert allclose(Es, h + 0.5*u*u/g) 10389 10390 time, Et = get_energy_through_cross_section(swwfile, 10391 cross_section, 10392 kind='total', 10393 verbose=False) 10394 assert allclose(Et, w + 0.5*u*u/g) 10395 10396 10286 10397 10287 10398 def test_get_all_swwfiles(self):
Note: See TracChangeset
for help on using the changeset viewer.