Changeset 7496
- Timestamp:
- Sep 7, 2009, 6:51:23 PM (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/culvert_flows/test_culvert_class.py
r7492 r7496 22 22 23 23 import numpy as num 24 25 26 # Helper functions 27 def run_culvert_flow_problem(depth): 28 """Run flow with culvert given depth 29 """ 30 31 length = 40. 32 width = 5. 33 34 dx = dy = 1 # Resolution: Length of subdivisions on both axes 35 36 points, vertices, boundary = rectangular_cross(int(length/dx), 37 int(width/dy), 38 len1=length, 39 len2=width) 40 domain = Domain(points, vertices, boundary) 41 domain.set_name('Test_culvert_shallow') # Output name 42 domain.set_default_order(2) 43 44 45 #---------------------------------------------------------------------- 46 # Setup initial conditions 47 #---------------------------------------------------------------------- 48 49 def topography(x, y): 50 """Set up a weir 51 52 A culvert will connect either side 53 """ 54 # General Slope of Topography 55 z=-x/1000 56 57 N = len(x) 58 for i in range(N): 59 60 # Sloping Embankment Across Channel 61 if 5.0 < x[i] < 10.1: 62 # Cut Out Segment for Culvert face 63 if 1.0+(x[i]-5.0)/5.0 < y[i] < 4.0 - (x[i]-5.0)/5.0: 64 z[i]=z[i] 65 else: 66 z[i] += 0.5*(x[i] -5.0) # Sloping Segment U/S Face 67 if 10.0 < x[i] < 12.1: 68 z[i] += 2.5 # Flat Crest of Embankment 69 if 12.0 < x[i] < 14.5: 70 # Cut Out Segment for Culvert face 71 if 2.0-(x[i]-12.0)/2.5 < y[i] < 3.0 + (x[i]-12.0)/2.5: 72 z[i]=z[i] 73 else: 74 z[i] += 2.5-1.0*(x[i] -12.0) # Sloping D/S Face 75 76 77 return z 78 79 80 domain.set_quantity('elevation', topography) 81 domain.set_quantity('friction', 0.01) # Constant friction 82 domain.set_quantity('stage', 83 expression='elevation + %f' % depth) # Shallow initial condition 84 85 # Boyd culvert 86 culvert = Culvert_flow(domain, 87 label='Culvert No. 1', 88 description='This culvert is a test unit 1.2m Wide by 0.75m High', 89 end_point0=[9.0, 2.5], 90 end_point1=[13.0, 2.5], 91 width=1.20, height=0.75, 92 culvert_routine=boyd_generalised_culvert_model, 93 number_of_barrels=1, 94 update_interval=2, 95 verbose=False) 96 97 98 domain.forcing_terms.append(culvert) 99 100 101 #----------------------------------------------------------------------- 102 # Setup boundary conditions 103 #----------------------------------------------------------------------- 104 105 # Inflow based on Flow Depth and Approaching Momentum 106 107 Br = Reflective_boundary(domain) # Solid reflective wall 108 domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br}) 109 110 111 112 #----------------------------------------------------------------------- 113 # Evolve system through time 114 #----------------------------------------------------------------------- 115 116 #print 'depth', depth 117 ref_volume = domain.get_quantity('stage').get_integral() 118 for t in domain.evolve(yieldstep = 0.1, finaltime = 25): 119 new_volume = domain.get_quantity('stage').get_integral() 120 121 msg = ('Total volume has changed: Is %.8f m^3 should have been %.8f m^3' 122 % (new_volume, ref_volume)) 123 assert num.allclose(new_volume, ref_volume), msg 124 24 125 25 126 … … 100 201 end_point0=[9.0, 2.5], 101 202 end_point1=[13.0, 2.5], 203 width=1.00, 102 204 use_velocity_head=True, 103 205 verbose=False) … … 214 316 end_point0=[9.0, 2.5], 215 317 end_point1=[13.0, 2.5], 318 height=0.75, 216 319 verbose=False) 217 320 … … 243 346 244 347 245 def test_that_culvert_rating_limits_flow_in_shallow_inlet_condition(self): 348 349 350 def test_that_culvert_flows_conserves_volume(self): 351 """test_that_culvert_flows_conserves_volume 352 353 Test that culvert on a sloping dry bed limits flows when very little water 354 is present at inlet. 355 356 Uses helper function: run_culvert_flow_problem(depth): 357 358 """ 359 360 # Try this for a range of depths 361 for depth in [0.1, 0.2, 0.5, 1.0]: 362 run_culvert_flow_problem(depth) 363 364 365 def OBSOLETE_XXXtest_that_culvert_rating_limits_flow_in_shallow_inlet_condition(self): 246 366 """test_that_culvert_rating_limits_flow_in_shallow_inlet_condition 247 367 … … 251 371 This one is using the rating curve variant 252 372 """ 373 374 253 375 254 376 path = get_pathname_from_package('anuga.culvert_flows') … … 307 429 domain.set_quantity('stage', 308 430 expression='elevation + 0.1') # Shallow initial condition 309 310 311 filename = os.path.join(path, 'example_rating_curve.csv') 431 432 # Boyd culvert 312 433 culvert = Culvert_flow(domain, 313 culvert_description_filename=filename, 434 label='Culvert No. 1', 435 description='This culvert is a test unit 1.2m Wide by 0.75m High', 314 436 end_point0=[9.0, 2.5], 315 437 end_point1=[13.0, 2.5], 316 trigger_depth=0.01, 438 width=1.20, height=0.75, 439 culvert_routine=boyd_generalised_culvert_model, 440 number_of_barrels=1, 441 update_interval=2, 317 442 verbose=False) 443 444 # Rating curve 445 #filename = os.path.join(path, 'example_rating_curve.csv') 446 #culvert = Culvert_flow(domain, 447 # culvert_description_filename=filename, 448 # end_point0=[9.0, 2.5], 449 # end_point1=[13.0, 2.5], 450 # trigger_depth=0.01, 451 # verbose=False) 318 452 319 453 domain.forcing_terms.append(culvert) … … 335 469 #----------------------------------------------------------------------- 336 470 471 print 'depth', 0.1 337 472 ref_volume = domain.get_quantity('stage').get_integral() 338 473 for t in domain.evolve(yieldstep = 0.1, finaltime = 25): 339 474 new_volume = domain.get_quantity('stage').get_integral() 340 475 341 342 476 msg = ('Total volume has changed: Is %.8f m^3 should have been %.8f m^3' 343 477 % (new_volume, ref_volume)) … … 346 480 347 481 return 348 # Now try this for a range of other depths 349 for depth in [1.0, 0.5, 0.2, 0.1, 0.05]: 482 # Now try this again for a depth of 10 cm and for a range of other depths 483 for depth in [0.1, 0.2, 0.5, 1.0]: 484 print 'depth', depth 350 485 domain.set_time(0.0) 351 486 487 domain.set_quantity('elevation', topography) 488 domain.set_quantity('friction', 0.01) # Constant friction 352 489 domain.set_quantity('stage', 353 490 expression='elevation + %f' % depth) … … 358 495 new_volume = domain.get_quantity('stage').get_integral() 359 496 360 msg = 'Total volume has changed: Is %. 2f m^3 should have been %.2f m^3'\497 msg = 'Total volume has changed: Is %.8f m^3 should have been %.8f m^3'\ 361 498 % (new_volume, ref_volume) 362 499 … … 439 576 end_point0=[9.0, 2.5], 440 577 end_point1=[13.0, 2.5], 441 width=1.20, height=0.75,578 width=1.20, height=0.75, 442 579 culvert_routine=boyd_generalised_culvert_model, 443 580 number_of_barrels=1, 444 581 update_interval=2, 445 verbose= True)582 verbose=False) 446 583 447 584 domain.forcing_terms.append(culvert) … … 543 680 culvert_routine=boyd_generalised_culvert_model, 544 681 number_of_barrels=1, 545 verbose= True)682 verbose=False) 546 683 547 684 … … 551 688 culvert(domain) 552 689 690 691 692 693 def test_momentum_jet(self): 694 """test_momentum_jet 695 696 Test that culvert_class can accept keyword use_momentum_jet 697 This does not yet imply that the values have been tested. FIXME 698 """ 699 700 701 length = 40. 702 width = 5. 703 704 dx = dy = 1 # Resolution: Length of subdivisions on both axes 705 706 points, vertices, boundary = rectangular_cross(int(length/dx), 707 int(width/dy), 708 len1=length, 709 len2=width) 710 domain = Domain(points, vertices, boundary) 711 domain.set_name('Test_culvert_shallow') # Output name 712 domain.set_default_order(2) 713 714 715 #---------------------------------------------------------------------- 716 # Setup initial conditions 717 #---------------------------------------------------------------------- 718 719 def topography(x, y): 720 """Set up a weir 721 722 A culvert will connect either side 723 """ 724 # General Slope of Topography 725 z=-x/1000 726 727 N = len(x) 728 for i in range(N): 729 730 # Sloping Embankment Across Channel 731 if 5.0 < x[i] < 10.1: 732 # Cut Out Segment for Culvert face 733 if 1.0+(x[i]-5.0)/5.0 < y[i] < 4.0 - (x[i]-5.0)/5.0: 734 z[i]=z[i] 735 else: 736 z[i] += 0.5*(x[i] -5.0) # Sloping Segment U/S Face 737 if 10.0 < x[i] < 12.1: 738 z[i] += 2.5 # Flat Crest of Embankment 739 if 12.0 < x[i] < 14.5: 740 # Cut Out Segment for Culvert face 741 if 2.0-(x[i]-12.0)/2.5 < y[i] < 3.0 + (x[i]-12.0)/2.5: 742 z[i]=z[i] 743 else: 744 z[i] += 2.5-1.0*(x[i] -12.0) # Sloping D/S Face 745 746 747 return z 748 749 750 domain.set_quantity('elevation', topography) 751 domain.set_quantity('friction', 0.01) # Constant friction 752 domain.set_quantity('stage', 753 expression='elevation + 1.0') # Shallow initial condition 754 755 # Boyd culvert 756 culvert = Culvert_flow(domain, 757 label='Culvert No. 1', 758 description='This culvert is a test unit 1.2m Wide by 0.75m High', 759 end_point0=[9.0, 2.5], 760 end_point1=[13.0, 2.5], 761 width=1.20, height=0.75, 762 culvert_routine=boyd_generalised_culvert_model, 763 number_of_barrels=1, 764 use_momentum_jet=True, 765 update_interval=2, 766 verbose=False) 767 768 769 domain.forcing_terms.append(culvert) 770 771 772 # Call 773 culvert(domain) 774 775 776 #----------------------------------------------------------------------- 777 # Setup boundary conditions 778 #----------------------------------------------------------------------- 779 780 781 Br = Reflective_boundary(domain) # Solid reflective wall 782 domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br}) 783 784 #----------------------------------------------------------------------- 785 # Evolve system through time 786 #----------------------------------------------------------------------- 787 788 ref_volume = domain.get_quantity('stage').get_integral() 789 for t in domain.evolve(yieldstep = 0.1, finaltime = 25): 790 new_volume = domain.get_quantity('stage').get_integral() 791 792 msg = ('Total volume has changed: Is %.8f m^3 should have been %.8f m^3' 793 % (new_volume, ref_volume)) 794 assert num.allclose(new_volume, ref_volume), msg 795 796 797 798 553 799 554 800 #-------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.
