- Timestamp:
- Apr 1, 2009, 3:19:07 PM (15 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/numpy/anuga/shallow_water/test_shallow_water_domain.py
r6553 r6689 2 2 3 3 import unittest, os 4 import os.path 4 5 from math import pi, sqrt 5 6 import tempfile … … 14 15 from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross 15 16 17 from anuga.utilities.system_tools import get_pathname_from_package 16 18 from shallow_water_domain import * 17 19 … … 720 722 # Check that flux seen from other triangles is inverse 721 723 (ql, qr) = (qr, ql) 722 #tmp = qr723 #qr = ql724 #ql = tmp725 724 normal = domain.get_normal(2, 2) 726 725 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, … … 738 737 # Check that flux seen from other triangles is inverse 739 738 (ql, qr) = (qr, ql) 740 #tmp = qr741 #qr = ql742 #ql = tmp743 739 normal = domain.get_normal(3, 0) 744 740 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, … … 756 752 # Check that flux seen from other triangles is inverse 757 753 (ql, qr) = (qr, ql) 758 #tmp = qr759 #qr=ql760 #ql=tmp761 754 normal = domain.get_normal(0, 1) 762 755 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, … … 2472 2465 R = Rainfall(domain, rate=2.0, polygon=[[1,1], [2,1], [2,2], [1,2]]) 2473 2466 2474 assert num.allclose(R.exchange_area, 1)2467 assert num.allclose(R.exchange_area, 2) 2475 2468 2476 2469 domain.forcing_terms.append(R) … … 2508 2501 # restricted to a polygon enclosing triangle #1 (bce) 2509 2502 domain.forcing_terms = [] 2510 R = Rainfall(domain, rate=lambda t: 3*t + 7, 2511 polygon=[[1,1], [2,1], [2,2], [1,2]]) 2512 2513 assert num.allclose(R.exchange_area, 1) 2514 2503 R = Rainfall(domain, 2504 rate=lambda t: 3*t + 7, 2505 polygon = [[1,1], [2,1], [2,2], [1,2]]) 2506 2507 assert num.allclose(R.exchange_area, 2) 2508 2515 2509 domain.forcing_terms.append(R) 2516 2510 … … 2551 2545 # restricted to a polygon enclosing triangle #1 (bce) 2552 2546 domain.forcing_terms = [] 2553 R = Rainfall(domain, rate=lambda t: 3*t + 7, 2554 polygon=rainfall_poly) 2555 2556 assert num.allclose(R.exchange_area, 1) 2557 2547 R = Rainfall(domain, 2548 rate=lambda t: 3*t + 7, 2549 polygon=rainfall_poly) 2550 2551 assert num.allclose(R.exchange_area, 2) 2552 2558 2553 domain.forcing_terms.append(R) 2559 2554 … … 2616 2611 # restricted to a polygon enclosing triangle #1 (bce) 2617 2612 domain.forcing_terms = [] 2618 R = Rainfall(domain, rate=lambda t: 3*t + 7, polygon=rainfall_poly) 2619 2620 assert num.allclose(R.exchange_area, 1) 2621 2613 R = Rainfall(domain, 2614 rate=lambda t: 3*t + 7, 2615 polygon=rainfall_poly) 2616 2617 assert num.allclose(R.exchange_area, 2) 2618 2622 2619 domain.forcing_terms.append(R) 2623 2620 … … 2680 2677 2681 2678 domain.forcing_terms = [] 2682 R = Rainfall(domain, rate=main_rate, 2683 polygon=[[1,1], [2,1], [2,2], [1,2]], default_rate=5.0) 2684 2685 assert num.allclose(R.exchange_area, 1) 2686 2679 R = Rainfall(domain, 2680 rate=main_rate, 2681 polygon = [[1,1], [2,1], [2,2], [1,2]], 2682 default_rate=5.0) 2683 2684 assert num.allclose(R.exchange_area, 2) 2685 2687 2686 domain.forcing_terms.append(R) 2688 2687 … … 2746 2745 2747 2746 domain.forcing_terms = [] 2748 R = Rainfall(domain, rate=main_rate, 2749 polygon=[[1,1], [2,1], [2,2], [1,2]], default_rate=5.0) 2750 2751 assert num.allclose(R.exchange_area, 1) 2752 2747 R = Rainfall(domain, 2748 rate=main_rate, 2749 polygon=[[1,1], [2,1], [2,2], [1,2]], 2750 default_rate=5.0) 2751 2752 assert num.allclose(R.exchange_area, 2) 2753 2753 2754 domain.forcing_terms.append(R) 2754 2755 … … 2785 2786 # on a circle affecting triangles #0 and #1 (bac and bce) 2786 2787 domain.forcing_terms = [] 2787 domain.forcing_terms.append(Inflow(domain, rate=2.0,2788 center=(1,1), radius=1))2789 2788 2789 I = Inflow(domain, rate=2.0, center=(1,1), radius=1) 2790 domain.forcing_terms.append(I) 2790 2791 domain.compute_forcing_terms() 2791 2792 2792 assert num.allclose(domain.quantities['stage'].explicit_update[1], 2793 2.0/pi) 2794 assert num.allclose(domain.quantities['stage'].explicit_update[0], 2795 2.0/pi) 2796 assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0) 2793 2794 A = I.exchange_area 2795 assert num.allclose(A, 4) # Two triangles 2796 2797 assert num.allclose(domain.quantities['stage'].explicit_update[1], 2.0/A) 2798 assert num.allclose(domain.quantities['stage'].explicit_update[0], 2.0/A) 2799 assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0) 2800 2797 2801 2798 2802 def test_inflow_using_circle_function(self): … … 2823 2827 # on a circle affecting triangles #0 and #1 (bac and bce) 2824 2828 domain.forcing_terms = [] 2825 domain.forcing_terms.append(Inflow(domain, rate=lambda t: 2., 2826 center=(1,1), radius=1)) 2827 2828 domain.compute_forcing_terms() 2829 2830 assert num.allclose(domain.quantities['stage'].explicit_update[1], 2831 2.0/pi) 2832 assert num.allclose(domain.quantities['stage'].explicit_update[0], 2833 2.0/pi) 2834 assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0) 2829 I = Inflow(domain, rate=lambda t: 2., center=(1,1), radius=1) 2830 domain.forcing_terms.append(I) 2831 2832 domain.compute_forcing_terms() 2833 2834 A = I.exchange_area 2835 assert num.allclose(A, 4) # Two triangles 2836 2837 assert num.allclose(domain.quantities['stage'].explicit_update[1], 2.0/A) 2838 assert num.allclose(domain.quantities['stage'].explicit_update[0], 2.0/A) 2839 assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0) 2840 2841 2842 2835 2843 2836 2844 def test_inflow_catch_too_few_triangles(self): … … 6233 6241 from anuga.geospatial_data.geospatial_data import Geospatial_data 6234 6242 6235 #-------------------------------------------------------------------- 6243 6244 # Get path where this test is run 6245 path = get_pathname_from_package('anuga.shallow_water') 6246 6247 6248 #---------------------------------------------------------------------- 6236 6249 # Create domain 6237 6250 #-------------------------------------------------------------------- … … 6276 6289 interior_regions = [] 6277 6290 for polygon in offending_regions: 6278 interior_regions.append( [polygon, 100] ) 6279 6280 meshname = 'offending_mesh.msh'6291 interior_regions.append( [polygon, 100] ) 6292 6293 meshname = os.path.join(path, 'offending_mesh.msh') 6281 6294 create_mesh_from_regions(bounding_polygon, 6282 6295 boundary_tags={'south': [0], 'east': [1], … … 6290 6303 domain = Domain(meshname, use_cache=False, verbose=verbose) 6291 6304 6292 #-------------------------------------------------------------------- 6305 #---------------------------------------------------------------------- 6293 6306 # Fit data point to mesh 6294 #-------------------------------------------------------------------- 6295 points_file = 'offending_point.pts' 6296 6297 # This is the offending point 6307 #---------------------------------------------------------------------- 6308 6309 points_file = os.path.join(path, 'offending_point.pts') 6310 6311 # Offending point 6298 6312 G = Geospatial_data(data_points=[[306953.344, 6194461.5]], 6299 6313 attributes=[1]) 6300 6314 G.export_points_file(points_file) 6301 6302 domain.set_quantity('elevation', filename=points_file, use_cache=False, 6303 verbose=verbose, alpha=0.01) 6315 6316 try: 6317 domain.set_quantity('elevation', filename=points_file, 6318 use_cache=False, verbose=verbose, alpha=0.01) 6319 except RuntimeError, e: 6320 msg = 'Test failed: %s' % str(e) 6321 raise Exception, msg 6322 # clean up in case raise fails 6323 os.remove(meshname) 6324 os.remove(points_file) 6325 else: 6326 os.remove(meshname) 6327 os.remove(points_file) 6328 6329 #finally: 6330 # Cleanup regardless 6331 #FIXME(Ole): Finally does not work like this in python2.4 6332 #FIXME(Ole): Reinstate this when Python2.4 is out of the way 6333 #FIXME(Ole): Python 2.6 apparently introduces something called 'with' 6334 #os.remove(meshname) 6335 #os.remove(points_file) 6336 6337 6338 def test_fitting_example_that_crashed_2(self): 6339 """test_fitting_example_that_crashed_2 6340 6341 This unit test has been derived from a real world example 6342 (the JJKelly study, by Petar Milevski). 6343 6344 It shows a condition where set_quantity crashes due to AtA 6345 not being built properly 6346 6347 See ticket:314 6348 """ 6349 6350 verbose = False 6351 6352 from anuga.shallow_water import Domain 6353 from anuga.pmesh.mesh_interface import create_mesh_from_regions 6354 from anuga.geospatial_data import Geospatial_data 6355 6356 # Get path where this test is run 6357 path = get_pathname_from_package('anuga.shallow_water') 6358 6359 meshname = os.path.join(path, 'test_mesh.msh') 6360 W = 304180 6361 S = 6185270 6362 E = 307650 6363 N = 6189040 6364 maximum_triangle_area = 1000000 6365 6366 bounding_polygon = [[W, S], [E, S], [E, N], [W, N]] 6367 6368 create_mesh_from_regions(bounding_polygon, 6369 boundary_tags={'south': [0], 6370 'east': [1], 6371 'north': [2], 6372 'west': [3]}, 6373 maximum_triangle_area=maximum_triangle_area, 6374 filename=meshname, 6375 use_cache=False, 6376 verbose=verbose) 6377 6378 domain = Domain(meshname, use_cache=True, verbose=verbose) 6379 6380 # Large test set revealed one problem 6381 points_file = os.path.join(path, 'test_points_large.csv') 6382 try: 6383 domain.set_quantity('elevation', filename=points_file, 6384 use_cache=False, verbose=verbose) 6385 except AssertionError, e: 6386 msg = 'Test failed: %s' % str(e) 6387 raise Exception, msg 6388 # Cleanup in case this failed 6389 os.remove(meshname) 6390 6391 # Small test set revealed another problem 6392 points_file = os.path.join(path, 'test_points_small.csv') 6393 try: 6394 domain.set_quantity('elevation', filename=points_file, 6395 use_cache=False, verbose=verbose) 6396 except AssertionError, e: 6397 msg = 'Test failed: %s' % str(e) 6398 raise Exception, msg 6399 # Cleanup in case this failed 6400 os.remove(meshname) 6401 else: 6402 os.remove(meshname) 6403 6404 6405 def test_total_volume(self): 6406 """test_total_volume 6407 6408 Test that total volume can be computed correctly 6409 """ 6410 6411 #---------------------------------------------------------------------- 6412 # Import necessary modules 6413 #---------------------------------------------------------------------- 6414 from anuga.abstract_2d_finite_volumes.mesh_factory \ 6415 import rectangular_cross 6416 from anuga.shallow_water import Domain 6417 6418 #---------------------------------------------------------------------- 6419 # Setup computational domain 6420 #---------------------------------------------------------------------- 6421 6422 length = 100. 6423 width = 20. 6424 dx = dy = 5 # Resolution: of grid on both axes 6425 6426 points, vertices, boundary = rectangular_cross(int(length/dx), 6427 int(width/dy), 6428 len1=length, 6429 len2=width) 6430 domain = Domain(points, vertices, boundary) 6431 6432 #---------------------------------------------------------------------- 6433 # Simple flat bottom bathtub 6434 #---------------------------------------------------------------------- 6435 6436 d = 1.0 6437 domain.set_quantity('elevation', 0.0) 6438 domain.set_quantity('stage', d) 6439 6440 assert num.allclose(domain.compute_total_volume(), length*width*d) 6441 6442 #---------------------------------------------------------------------- 6443 # Slope 6444 #---------------------------------------------------------------------- 6445 6446 slope = 1.0/10 # RHS drops to -10m 6447 def topography(x, y): 6448 return -x * slope 6449 6450 domain.set_quantity('elevation', topography) 6451 domain.set_quantity('stage', 0.0) # Domain full 6452 6453 V = domain.compute_total_volume() 6454 assert num.allclose(V, length*width*10/2) 6455 6456 domain.set_quantity('stage', -5.0) # Domain 'half' full 6457 6458 # IMPORTANT: Adjust stage to match elevation 6459 domain.distribute_to_vertices_and_edges() 6460 6461 V = domain.compute_total_volume() 6462 assert num.allclose(V, width*(length/2)*5.0/2) 6463 6464 6465 def test_volumetric_balance_computation(self): 6466 """test_volumetric_balance_computation 6467 6468 Test that total in and out flows are computed correctly 6469 FIXME(Ole): This test is more about looking at the printed report 6470 """ 6471 6472 verbose = False 6473 6474 #---------------------------------------------------------------------- 6475 # Import necessary modules 6476 #---------------------------------------------------------------------- 6477 6478 from anuga.abstract_2d_finite_volumes.mesh_factory \ 6479 import rectangular_cross 6480 from anuga.shallow_water import Domain 6481 from anuga.shallow_water.shallow_water_domain import Reflective_boundary 6482 from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary 6483 from anuga.shallow_water.shallow_water_domain import Inflow 6484 from anuga.shallow_water.data_manager \ 6485 import get_flow_through_cross_section 6486 6487 #---------------------------------------------------------------------- 6488 # Setup computational domain 6489 #---------------------------------------------------------------------- 6490 6491 finaltime = 300.0 6492 length = 300. 6493 width = 20. 6494 dx = dy = 5 # Resolution: of grid on both axes 6495 6496 # Input parameters 6497 uh = 1.0 6498 vh = 0.0 6499 d = 1.0 6500 6501 # 20 m^3/s in the x direction across entire domain 6502 ref_flow = uh*d*width 6503 6504 points, vertices, boundary = rectangular_cross(int(length/dx), 6505 int(width/dy), 6506 len1=length, 6507 len2=width) 6508 6509 domain = Domain(points, vertices, boundary) 6510 domain.set_name('Inflow_flowline_test') # Output name 6511 6512 #---------------------------------------------------------------------- 6513 # Setup initial conditions 6514 #---------------------------------------------------------------------- 6515 6516 slope = 0.0 6517 def topography(x, y): 6518 z=-x * slope 6519 return z 6520 6521 # Use function for elevation 6522 domain.set_quantity('elevation', topography) 6523 domain.set_quantity('friction', 0.0) # Constant friction 6524 6525 domain.set_quantity('stage', expression='elevation') 6526 6527 #---------------------------------------------------------------------- 6528 # Setup boundary conditions 6529 #---------------------------------------------------------------------- 6530 6531 Br = Reflective_boundary(domain) # Solid reflective wall 6532 6533 # Constant flow into domain 6534 # Depth = 1m, uh=1 m/s, i.e. a flow of 20 m^3/s 6535 6536 Bi = Dirichlet_boundary([d, uh, vh]) 6537 Bo = Dirichlet_boundary([0, 0, 0]) 6538 6539 domain.set_boundary({'left': Bi, 'right': Bo, 'top': Br, 'bottom': Br}) 6540 6541 #---------------------------------------------------------------------- 6542 # Evolve system through time 6543 #---------------------------------------------------------------------- 6544 6545 for t in domain.evolve(yieldstep=100.0, finaltime=finaltime): 6546 S = domain.volumetric_balance_statistics() 6547 if verbose : 6548 print domain.timestepping_statistics() 6549 print S 6550 6551 6552 def Xtest_inflow_using_flowline(self): 6553 """test_inflow_using_flowline 6554 6555 Test the ability of a flowline to match inflow above the flowline by 6556 creating constant inflow onto a circle at the head of a 20m wide by 300m 6557 long plane dipping at various slopes with a perpendicular flowline and 6558 gauge downstream of the inflow and a 45 degree flowlines at 200m 6559 downstream. 6560 """ 6561 6562 # FIXME(Ole): Move this and the following test to validate_all.py as 6563 # they are rather time consuming 6564 6565 verbose = True 6566 6567 #---------------------------------------------------------------------- 6568 # Import necessary modules 6569 #---------------------------------------------------------------------- 6570 6571 from anuga.abstract_2d_finite_volumes.mesh_factory \ 6572 import rectangular_cross 6573 from anuga.shallow_water import Domain 6574 from anuga.shallow_water.shallow_water_domain import Reflective_boundary 6575 from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary 6576 from anuga.shallow_water.shallow_water_domain import Inflow 6577 from anuga.shallow_water.data_manager \ 6578 import get_flow_through_cross_section 6579 from anuga.abstract_2d_finite_volumes.util \ 6580 import sww2csv_gauges, csv2timeseries_graphs 6581 6582 #---------------------------------------------------------------------- 6583 # Setup computational domain 6584 #---------------------------------------------------------------------- 6585 6586 number_of_inflows = 2 # Number of inflows on top of each other 6587 finaltime = 1000.0 # 6000.0 6588 6589 length = 300. 6590 width = 20. 6591 dx = dy = 5 # Resolution: of grid on both axes 6592 6593 points, vertices, boundary = rectangular_cross(int(length/dx), 6594 int(width/dy), 6595 len1=length, 6596 len2=width) 6597 6598 for mannings_n in [0.150, 0.07, 0.035]: #[0.012, 0.035, 0.070, 0.150]: 6599 for slope in [1.0/300, 1.0/150, 1.0/75]: 6600 # Loop over a range of bedslopes representing sub to 6601 # super critical flows 6602 if verbose: 6603 print 6604 print 'Slope:', slope, 'Mannings n:', mannings_n 6605 domain = Domain(points, vertices, boundary) 6606 domain.set_name('Inflow_flowline_test') # Output name 6607 6608 #-------------------------------------------------------------- 6609 # Setup initial conditions 6610 #-------------------------------------------------------------- 6611 6612 def topography(x, y): 6613 z = -x * slope 6614 return z 6615 6616 # Use function for elevation 6617 domain.set_quantity('elevation', topography) 6618 # Constant friction of conc surface 6619 domain.set_quantity('friction', mannings_n) 6620 # Dry initial condition 6621 domain.set_quantity('stage', expression='elevation') 6622 6623 #-------------------------------------------------------------- 6624 # Setup boundary conditions 6625 #-------------------------------------------------------------- 6626 6627 Br = Reflective_boundary(domain) # Solid reflective wall 6628 # Outflow stage at -0.9m d=0.1m 6629 Bo = Dirichlet_boundary([-100, 0, 0]) 6630 6631 domain.set_boundary({'left': Br, 'right': Bo, 6632 'top': Br, 'bottom': Br}) 6633 6634 #-------------------------------------------------------------- 6635 # Setup Inflow 6636 #-------------------------------------------------------------- 6637 6638 # Fixed Flowrate onto Area 6639 fixed_inflow = Inflow(domain, center=(10.0, 10.0), 6640 radius=5.00, rate=10.00) 6641 6642 # Stack this flow 6643 for i in range(number_of_inflows): 6644 domain.forcing_terms.append(fixed_inflow) 6645 6646 ref_flow = fixed_inflow.rate*number_of_inflows 6647 6648 #-------------------------------------------------------------- 6649 # Evolve system through time 6650 #-------------------------------------------------------------- 6651 6652 for t in domain.evolve(yieldstep=100.0, finaltime=finaltime): 6653 pass 6654 #if verbose : 6655 # print domain.timestepping_statistics() 6656 6657 if verbose: 6658 print domain.volumetric_balance_statistics() 6659 #-------------------------------------------------------------- 6660 # Compute flow thru flowlines ds of inflow 6661 #-------------------------------------------------------------- 6662 6663 # Square on flowline at 200m 6664 q = domain.get_flow_through_cross_section([[200.0, 0.0], 6665 [200.0, 20.0]]) 6666 if verbose: 6667 print ('90 degree flowline: ANUGA = %f, Ref = %f' 6668 % (q, ref_flow)) 6669 6670 msg = ('Predicted flow was %f, should have been %f' 6671 % (q, ref_flow)) 6672 assert num.allclose(q, ref_flow, rtol=1.0e-2), msg 6673 6674 6675 # 45 degree flowline at 200m 6676 q = domain.get_flow_through_cross_section([[200.0, 0.0], 6677 [220.0, 20.0]]) 6678 if verbose: 6679 print ('45 degree flowline: ANUGA = %f, Ref = %f' 6680 % (q, ref_flow)) 6681 6682 msg = ('Predicted flow was %f, should have been %f' 6683 % (q, ref_flow)) 6684 assert num.allclose(q, ref_flow, rtol=1.0e-2), msg 6685 6686 # Stage recorder (gauge) in middle of plane at 200m 6687 x = 200.0 6688 y = 10.00 6689 w = domain.get_quantity('stage').\ 6690 get_values(interpolation_points=[[x, y]])[0] 6691 z = domain.get_quantity('elevation').\ 6692 get_values(interpolation_points=[[x, y]])[0] 6693 domain_depth = w-z 6694 6695 # Compute normal depth at gauge location using Manning equation 6696 # v=1/n*(r^2/3)*(s^0.5) or r=(Q*n/(s^0.5*W))^0.6 6697 normal_depth = (ref_flow*mannings_n/(slope**0.5*width))**0.6 6698 msg = ('Predicted depth of flow was %f, should have been %f' 6699 % (domain_depth, normal_depth)) 6700 if verbose: 6701 diff = abs(domain_depth-normal_depth) 6702 print ('Depth: ANUGA = %f, Mannings = %f, E=%f' 6703 % (domain_depth, normal_depth, diff/domain_depth)) 6704 assert diff < 0.1 6705 6706 if slope >= 1.0/100: 6707 # Really super critical flow is not as stable. 6708 #assert num.allclose(domain_depth, normal_depth, 6709 # rtol=1.0e-1), msg 6710 pass 6711 else: 6712 pass 6713 #assert num.allclose(domain_depth, normal_depth, 6714 # rtol=1.0e-2), msg 6715 6716 6717 def Xtest_friction_dependent_flow_using_flowline(self): 6718 """test_friction_dependent_flow_using_flowline 6719 6720 Test the internal flow (using flowline) as a function of 6721 different values of Mannings n and different slopes. 6722 6723 Flow is applied in the form of boundary conditions with fixed momentum. 6724 """ 6725 6726 verbose = True 6727 6728 #---------------------------------------------------------------------- 6729 # Import necessary modules 6730 #---------------------------------------------------------------------- 6731 6732 from anuga.abstract_2d_finite_volumes.mesh_factory \ 6733 import rectangular_cross 6734 from anuga.shallow_water import Domain 6735 from anuga.shallow_water.shallow_water_domain import Reflective_boundary 6736 from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary 6737 from anuga.shallow_water.shallow_water_domain import Inflow 6738 from anuga.shallow_water.data_manager \ 6739 import get_flow_through_cross_section 6740 from anuga.abstract_2d_finite_volumes.util \ 6741 import sww2csv_gauges, csv2timeseries_graphs 6742 6743 6744 #---------------------------------------------------------------------- 6745 # Setup computational domain 6746 #---------------------------------------------------------------------- 6747 6748 finaltime = 1000.0 6749 6750 length = 300. 6751 width = 20. 6752 dx = dy = 5 # Resolution: of grid on both axes 6753 6754 # Input parameters 6755 uh = 1.0 6756 vh = 0.0 6757 d = 1.0 6758 6759 ref_flow = uh*d*width # 20 m^3/s in the x direction across entire domain 6760 6761 points, vertices, boundary = rectangular_cross(int(length/dx), 6762 int(width/dy), 6763 len1=length, 6764 len2=width) 6765 6766 for mannings_n in [0.035]: #[0.0, 0.012, 0.035]: 6767 for slope in [1.0/300]: #[0.0, 1.0/300, 1.0/150]: 6768 # Loop over a range of bedslopes representing 6769 # sub to super critical flows 6770 if verbose: 6771 print 6772 print 'Slope:', slope, 'Mannings n:', mannings_n 6773 domain = Domain(points, vertices, boundary) 6774 domain.set_name('Inflow_flowline_test') # Output name 6775 6776 #-------------------------------------------------------------- 6777 # Setup initial conditions 6778 #-------------------------------------------------------------- 6779 6780 def topography(x, y): 6781 z = -x * slope 6782 return z 6783 6784 # Use function for elevation 6785 domain.set_quantity('elevation', topography) 6786 # Constant friction 6787 domain.set_quantity('friction', mannings_n) 6788 6789 #domain.set_quantity('stage', expression='elevation') 6790 6791 # Set initial flow as depth=1m, uh=1.0 m/s, vh = 0.0 6792 # making it 20 m^3/s across entire domain 6793 domain.set_quantity('stage', expression='elevation + %f' % d) 6794 domain.set_quantity('xmomentum', uh) 6795 domain.set_quantity('ymomentum', vh) 6796 6797 #-------------------------------------------------------------- 6798 # Setup boundary conditions 6799 #-------------------------------------------------------------- 6800 6801 Br = Reflective_boundary(domain) # Solid reflective wall 6802 6803 # Constant flow in and out of domain 6804 # Depth = 1m, uh=1 m/s, i.e. a flow of 20 m^3/s 6805 # across boundaries 6806 Bi = Dirichlet_boundary([d, uh, vh]) 6807 Bo = Dirichlet_boundary([-length*slope+d, uh, vh]) 6808 #Bo = Dirichlet_boundary([-100, 0, 0]) 6809 6810 domain.set_boundary({'left': Bi, 'right': Bo, 6811 'top': Br, 'bottom': Br}) 6812 6813 #-------------------------------------------------------------- 6814 # Evolve system through time 6815 #-------------------------------------------------------------- 6816 6817 for t in domain.evolve(yieldstep=100.0, finaltime=finaltime): 6818 if verbose : 6819 print domain.timestepping_statistics() 6820 print domain.volumetric_balance_statistics() 6821 6822 # 90 degree flowline at 200m 6823 q = domain.get_flow_through_cross_section([[200.0, 0.0], 6824 [200.0, 20.0]]) 6825 msg = ('Predicted flow was %f, should have been %f' 6826 % (q, ref_flow)) 6827 if verbose: 6828 print ('90 degree flowline: ANUGA = %f, Ref = %f' 6829 % (q, ref_flow)) 6830 6831 # 45 degree flowline at 200m 6832 q = domain.get_flow_through_cross_section([[200.0, 0.0], 6833 [220.0, 20.0]]) 6834 msg = ('Predicted flow was %f, should have been %f' 6835 % (q, ref_flow)) 6836 if verbose: 6837 print ('45 degree flowline: ANUGA = %f, Ref = %f' 6838 % (q, ref_flow)) 6839 6840 # Stage recorder (gauge) in middle of plane at 200m 6841 x = 200.0 6842 y = 10.00 6843 w = domain.get_quantity('stage').\ 6844 get_values(interpolation_points=[[x, y]])[0] 6845 z = domain.get_quantity('elevation').\ 6846 get_values(interpolation_points=[[x, y]])[0] 6847 domain_depth = w-z 6848 6849 xmom = domain.get_quantity('xmomentum').\ 6850 get_values(interpolation_points=[[x, y]])[0] 6851 ymom = domain.get_quantity('ymomentum').\ 6852 get_values(interpolation_points=[[x, y]])[0] 6853 if verbose: 6854 print ('At interpolation point (h, uh, vh): ', 6855 domain_depth, xmom, ymom) 6856 print 'uh * d * width = ', xmom*domain_depth*width 6857 6858 if slope > 0.0: 6859 # Compute normal depth at gauge location using Manning eqn 6860 # v=1/n*(r^2/3)*(s^0.5) or r=(Q*n/(s^0.5*W))^0.6 6861 normal_depth = (ref_flow*mannings_n/(slope**0.5*width))**0.6 6862 if verbose: 6863 print ('Depth: ANUGA = %f, Mannings = %f' 6864 % (domain_depth, normal_depth)) 6865 6866 ################################################################################# 6304 6867 6305 6868 if __name__ == "__main__":
Note: See TracChangeset
for help on using the changeset viewer.
