 Timestamp:
 Sep 3, 2008, 4:44:23 PM (15 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

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.