Changeset 6643
- Timestamp:
- Mar 27, 2009, 2:46:01 PM (16 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py
r6637 r6643 2529 2529 # Setup only one forcing term, constant rainfall 2530 2530 domain.forcing_terms = [] 2531 domain.forcing_terms.append( Rainfall(domain, rate=2.0))2531 domain.forcing_terms.append( Rainfall(domain, rate=2.0) ) 2532 2532 2533 2533 domain.compute_forcing_terms() … … 2567 2567 polygon = [[1,1], [2,1], [2,2], [1,2]]) 2568 2568 2569 assert num.allclose(R.exchange_area, 2)2569 assert num.allclose(R.exchange_area, 1) 2570 2570 2571 2571 domain.forcing_terms.append(R) … … 2734 2734 polygon=rainfall_poly) 2735 2735 2736 assert num.allclose(R.exchange_area, 2)2736 assert num.allclose(R.exchange_area, 1) 2737 2737 2738 2738 domain.forcing_terms.append(R) … … 2807 2807 default_rate=5.0) 2808 2808 2809 assert num.allclose(R.exchange_area, 2)2809 assert num.allclose(R.exchange_area, 1) 2810 2810 2811 2811 domain.forcing_terms.append(R) … … 2882 2882 default_rate=5.0) 2883 2883 2884 assert num.allclose(R.exchange_area, 2)2884 assert num.allclose(R.exchange_area, 1) 2885 2885 2886 2886 domain.forcing_terms.append(R) … … 2923 2923 2924 2924 # Setup only one forcing term, constant inflow of 2 m^3/s on a circle affecting triangles #0 and #1 (bac and bce) 2925 2926 I = Inflow(domain, rate=2.0, center=(1,1), radius=1)2927 2925 domain.forcing_terms = [] 2928 domain.forcing_terms.append( I)2926 domain.forcing_terms.append( Inflow(domain, rate=2.0, center=(1,1), radius=1) ) 2929 2927 2930 2928 domain.compute_forcing_terms() 2931 2932 ref_dw = 2.0/I.exchange_area 2933 2934 assert num.allclose(domain.quantities['stage'].explicit_update[1], ref_dw) 2935 assert num.allclose(domain.quantities['stage'].explicit_update[0], ref_dw) 2929 #print domain.quantities['stage'].explicit_update 2930 2931 assert num.allclose(domain.quantities['stage'].explicit_update[1], 2.0/pi) 2932 assert num.allclose(domain.quantities['stage'].explicit_update[0], 2.0/pi) 2936 2933 assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0) 2937 2934 … … 2964 2961 # Setup only one forcing term, time dependent inflow of 2 m^3/s on a circle affecting triangles #0 and #1 (bac and bce) 2965 2962 domain.forcing_terms = [] 2966 I = Inflow(domain, rate=lambda t: 2., center=(1,1), radius=1) 2967 domain.forcing_terms.append(I) 2963 domain.forcing_terms.append( Inflow(domain, rate=lambda t: 2., center=(1,1), radius=1) ) 2968 2964 2969 2965 domain.compute_forcing_terms() 2970 2966 2971 ref_dw = 2.0/I.exchange_area 2972 2973 assert num.allclose(domain.quantities['stage'].explicit_update[1], ref_dw) 2974 assert num.allclose(domain.quantities['stage'].explicit_update[0], ref_dw) 2967 assert num.allclose(domain.quantities['stage'].explicit_update[1], 2.0/pi) 2968 assert num.allclose(domain.quantities['stage'].explicit_update[0], 2.0/pi) 2975 2969 assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0) 2976 2970 … … 6735 6729 6736 6730 6737 def test_inflow_using_flowline(self):6731 def Xtest_inflow_using_flowline(self): 6738 6732 """test_inflow_using_flowline 6739 6733 Test the ability of a flowline to match inflow above the flowline by 6740 6734 creating constant inflow onto a circle at the head of a 20m 6741 wide by 300m long plane dipping at 1:300with a perpendicular flowline and gauge6742 downstream of the inflow and a 45 degree flowlines at 60m downstream6735 wide by 300m long plane dipping at various slopes with a perpendicular flowline and gauge 6736 downstream of the inflow and a 45 degree flowlines at 200m downstream 6743 6737 """ 6744 6738 6745 verbose = False6739 verbose = True 6746 6740 6747 6741 … … 6762 6756 #------------------------------------------------------------------------------ 6763 6757 number_of_inflows = 2 # Number of inflows on top of each other 6758 finaltime = 300.0 6764 6759 6765 6760 length = 300. 6766 6761 width = 20. 6767 dx = dy = 2 6762 dx = dy = 2 # Resolution: of grid on both axes 6768 6763 6769 6764 points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy), 6770 6765 len1=length, len2=width) 6771 6766 6772 6773 for slope in [1.0/300, 1.0/200, 1.0/100, 1.0/50]: 6774 # Loop over a range of bedslopes representing sub to super critical flows 6775 6776 #print slope 6777 domain = Domain(points, vertices, boundary) 6778 domain.set_name('Inflow_flowline_test') # Output name 6779 6780 6781 #------------------------------------------------------------------------------ 6782 # Setup initial conditions 6783 #------------------------------------------------------------------------------ 6784 6785 def topography(x, y): 6786 z=-x * slope 6787 return z 6788 6789 6790 domain.set_quantity('elevation', topography) # Use function for elevation 6791 domain.set_quantity('friction', 0.012) # Constant friction of conc surface 6792 domain.set_quantity('stage', 6793 expression='elevation') # Dry initial condition 6794 6795 6796 #------------------------------------------------------------------------------ 6797 # Setup boundary conditions 6798 #------------------------------------------------------------------------------ 6799 6800 Br = Reflective_boundary(domain) # Solid reflective wall 6801 Bo = Dirichlet_boundary([-100, 0, 0]) # Outflow stsge at -0.9m d=0.1m 6802 6803 domain.set_boundary({'left': Br, 'right': Bo, 'top': Br, 'bottom': Br}) 6804 6805 #------------------------------------------------------------------------------ 6806 # Seup Inflow 6807 #------------------------------------------------------------------------------ 6808 6809 # Fixed Flowrate onto Area 6810 fixed_inflow = Inflow(domain, 6811 center=(10.0, 10.0), 6812 radius=5.00, 6813 rate=1.00) 6814 6815 # Stack this flow 6816 for i in range(number_of_inflows): 6817 domain.forcing_terms.append(fixed_inflow) 6818 6819 ref_flow = fixed_inflow.rate*number_of_inflows 6820 6821 #------------------------------------------------------------------------------ 6822 # Evolve system through time 6823 #------------------------------------------------------------------------------ 6824 6825 for t in domain.evolve(yieldstep = 10.0, finaltime = 300): 6826 pass 6827 #print t 6828 6829 #------------------------------------------------------------------------------ 6830 # Compute flow thru flowlines ds of inflow 6831 #------------------------------------------------------------------------------ 6767 for mannings_n in [0.0, 0.012, 0.035, 0.070, 0.150]: 6768 for slope in [1.0/300, 1.0/150, 1.0/75]: 6769 # Loop over a range of bedslopes representing sub to super critical flows 6770 6771 if verbose: 6772 print 6773 print 'Slope:', slope, 'Mannings n:', mannings_n 6774 domain = Domain(points, vertices, boundary) 6775 domain.set_name('Inflow_flowline_test') # Output name 6832 6776 6833 # Square on flowline at 30m 6834 q=domain.get_flow_through_cross_section([[200.0,0.0],[200.0,20.0]]) 6835 msg = 'Predicted flow was %f, should have been %f' % (q, ref_flow) 6836 #print q, ref_flow 6837 assert num.allclose(q, ref_flow, rtol=1.0e-2), msg 6838 6839 6840 # 45 degree flowline at 60m 6841 q=domain.get_flow_through_cross_section([[200.0,0.0],[220.0,20.0]]) 6842 msg = 'Predicted flow was %f, should have been %f' % (q, ref_flow) 6843 #print q, ref_flow 6844 assert num.allclose(q, ref_flow, rtol=1.0e-2), msg 6845 6846 6777 6778 #------------------------------------------------------------------------------ 6779 # Setup initial conditions 6780 #------------------------------------------------------------------------------ 6781 6782 def topography(x, y): 6783 z=-x * slope 6784 return z 6785 6786 domain.set_quantity('elevation', topography) # Use function for elevation 6787 domain.set_quantity('friction', mannings_n) # Constant friction of conc surface 6788 domain.set_quantity('stage', 6789 expression='elevation') # Dry initial condition 6790 6791 6792 #------------------------------------------------------------------------------ 6793 # Setup boundary conditions 6794 #------------------------------------------------------------------------------ 6795 6796 Br = Reflective_boundary(domain) # Solid reflective wall 6797 Bo = Dirichlet_boundary([-100, 0, 0]) # Outflow stsge at -0.9m d=0.1m 6798 6799 domain.set_boundary({'left': Br, 'right': Bo, 'top': Br, 'bottom': Br}) 6800 6801 #------------------------------------------------------------------------------ 6802 # Seup Inflow 6803 #------------------------------------------------------------------------------ 6804 6805 # Fixed Flowrate onto Area 6806 fixed_inflow = Inflow(domain, 6807 center=(10.0, 10.0), 6808 radius=5.00, 6809 rate=10.00) 6810 6811 # Stack this flow 6812 for i in range(number_of_inflows): 6813 domain.forcing_terms.append(fixed_inflow) 6814 6815 ref_flow = fixed_inflow.rate*number_of_inflows 6816 6817 #------------------------------------------------------------------------------ 6818 # Evolve system through time 6819 #------------------------------------------------------------------------------ 6820 6821 6822 for t in domain.evolve(yieldstep=100.0, finaltime=finaltime): 6823 if verbose : 6824 print domain.timestepping_statistics() 6825 6826 6827 x=200.0 6828 y=10.00 6829 6830 6831 6832 #------------------------------------------------------------------------------ 6833 # Compute flow thru flowlines ds of inflow 6834 #------------------------------------------------------------------------------ 6835 6836 # Square on flowline at 200m 6837 q=domain.get_flow_through_cross_section([[200.0,0.0],[200.0,20.0]]) 6838 msg = 'Predicted flow was %f, should have been %f' % (q, ref_flow) 6839 if verbose: 6840 print '90 degree flowline: ANUGA = %f, Ref = %f' % (q, ref_flow) 6841 #assert num.allclose(q, ref_flow, rtol=1.0e-2), msg 6842 6843 6844 # 45 degree flowline at 200m 6845 q=domain.get_flow_through_cross_section([[200.0,0.0],[220.0,20.0]]) 6846 msg = 'Predicted flow was %f, should have been %f' % (q, ref_flow) 6847 if verbose: 6848 print '45 degree flowline: ANUGA = %f, Ref = %f' % (q, ref_flow) 6849 6850 #assert num.allclose(q, ref_flow, rtol=1.0e-2), msg 6851 6852 # Stage recorder (gauge) in middle of plane at 200m 6853 w = domain.get_quantity('stage').get_values(interpolation_points=[[x, y]])[0] 6854 z = domain.get_quantity('elevation').get_values(interpolation_points=[[x, y]])[0] 6855 domain_depth = w-z 6856 6857 6858 # Compute normal depth at gauge location using Manning equation 6859 # v=1/n*(r^2/3)*(s^0.5) or r=(Q*n/(s^0.5*W))^0.6 6860 normal_depth=(ref_flow*mannings_n/(slope**0.5*width))**0.6 6861 msg = 'Predicted depth of flow was %f, should have been %f' % (normal_depth, domain_depth) 6862 if verbose: 6863 print 'Depth: ANUGA = %f, Mannings = %f' % (domain_depth, normal_depth) 6864 6865 if slope >= 1.0/100: 6866 # Really super critical flow is not as stable. 6867 #assert num.allclose(domain_depth,normal_depth, rtol=1.0e-1), msg 6868 pass 6869 else: 6870 pass 6871 #assert num.allclose(domain_depth,normal_depth, rtol=1.0e-2), msg 6847 6872 6848 6873 … … 6851 6876 6852 6877 if __name__ == "__main__": 6853 suite = unittest.makeSuite(Test_Shallow_Water, 'test ')6878 suite = unittest.makeSuite(Test_Shallow_Water, 'test_inflow_using_flowline') 6854 6879 runner = unittest.TextTestRunner(verbosity=1) 6855 6880 runner.run(suite)
Note: See TracChangeset
for help on using the changeset viewer.
