1 | |
---|
2 | |
---|
3 | """test_flow_on_plane.py |
---|
4 | Test the ability of ANUGA to maintain flow and match normal depth on an inclined plane |
---|
5 | using constant inflow thru the upstream bdry at normal_depth and equivalent momentum |
---|
6 | The plane is 20m wide by 400m long dipping at various slopes with a range of roughnesses |
---|
7 | A recording gauges is located 200m downstream |
---|
8 | """ |
---|
9 | |
---|
10 | verbose = True |
---|
11 | import numpy as num |
---|
12 | |
---|
13 | #------------------------------------------------------------------------------ |
---|
14 | # Import necessary modules |
---|
15 | #------------------------------------------------------------------------------ |
---|
16 | from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross |
---|
17 | from anuga.shallow_water import Domain |
---|
18 | from anuga.shallow_water.shallow_water_domain import Reflective_boundary |
---|
19 | from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary |
---|
20 | from anuga.shallow_water.shallow_water_domain import Transmissive_Momentum_Set_Stage_boundary |
---|
21 | from anuga.shallow_water.shallow_water_domain import Inflow |
---|
22 | from anuga.shallow_water.data_manager import get_flow_through_cross_section |
---|
23 | from anuga.abstract_2d_finite_volumes.util import sww2csv_gauges, csv2timeseries_graphs |
---|
24 | |
---|
25 | |
---|
26 | #------------------------------------------------------------------------------ |
---|
27 | # Setup computational domain |
---|
28 | #------------------------------------------------------------------------------ |
---|
29 | |
---|
30 | finaltime = 1800.0 # 1/2 hour |
---|
31 | |
---|
32 | length = 400. |
---|
33 | width = 20. |
---|
34 | dx = dy = 2 # Resolution: of grid on both axes |
---|
35 | ref_flow = 20.0 # Inflow at head of plane |
---|
36 | |
---|
37 | points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy), |
---|
38 | len1=length, len2=width) |
---|
39 | |
---|
40 | for mannings_n in [0.150, 0.070, 0.03, 0.015]: |
---|
41 | # Loop over a range of roughnesses |
---|
42 | |
---|
43 | for slope in [1.0/1000, 1.0/300, 1.0/100]: |
---|
44 | # Loop over a range of bedslopes |
---|
45 | |
---|
46 | |
---|
47 | domain = Domain(points, vertices, boundary) |
---|
48 | domain.set_name('test_flowhydraulics_on_plane') # Output name |
---|
49 | |
---|
50 | |
---|
51 | #------------------------------------------------------------------------------ |
---|
52 | # Setup initial conditions |
---|
53 | #------------------------------------------------------------------------------ |
---|
54 | |
---|
55 | def topography(x, y): |
---|
56 | z=-x * slope |
---|
57 | return z |
---|
58 | |
---|
59 | domain.set_quantity('elevation', topography) # Use function for elevation |
---|
60 | domain.set_quantity('friction', mannings_n) # Constant friction of conc surface |
---|
61 | domain.set_quantity('stage', |
---|
62 | expression='elevation') # Dry initial condition |
---|
63 | |
---|
64 | |
---|
65 | #------------------------------------------------------------------------------ |
---|
66 | # Seup Inflow and compute flow characteristics |
---|
67 | #------------------------------------------------------------------------------ |
---|
68 | |
---|
69 | # Compute normal depth, Froude and momentum on plane using Mannings equation |
---|
70 | normal_depth =(ref_flow*mannings_n/(slope**0.5*width))**0.6 |
---|
71 | flow_velocity = ref_flow/(width*normal_depth) |
---|
72 | xmomentum = flow_velocity*normal_depth |
---|
73 | froude_no = flow_velocity/(9.8*normal_depth)**0.5 |
---|
74 | if verbose: |
---|
75 | print |
---|
76 | print 'Slope:', slope, 'Mannings n:', mannings_n, 'Velocity:',flow_velocity,'Froude:',froude_no |
---|
77 | |
---|
78 | |
---|
79 | #------------------------------------------------------------------------------ |
---|
80 | # Setup boundary conditions |
---|
81 | #------------------------------------------------------------------------------ |
---|
82 | |
---|
83 | Br = Reflective_boundary(domain) # Solid reflective wall |
---|
84 | Bo = Dirichlet_boundary([-(slope*length)+normal_depth, xmomentum, 0]) # Outflow stage at normal depth |
---|
85 | Bi = Dirichlet_boundary([normal_depth, xmomentum, 0]) # Inflow stage at normal depth |
---|
86 | |
---|
87 | domain.set_boundary({'left': Bi, 'right': Bo, 'top': Br, 'bottom': Br}) |
---|
88 | |
---|
89 | #------------------------------------------------------------------------------ |
---|
90 | # Evolve system through time |
---|
91 | #------------------------------------------------------------------------------ |
---|
92 | |
---|
93 | for t in domain.evolve(yieldstep=100.0, finaltime=finaltime): |
---|
94 | if verbose : |
---|
95 | print domain.timestepping_statistics() |
---|
96 | |
---|
97 | |
---|
98 | #------------------------------------------------------------------------------ |
---|
99 | # Compute depth and check final depth against normal depth |
---|
100 | #------------------------------------------------------------------------------ |
---|
101 | |
---|
102 | |
---|
103 | # Store temporal pattern at 200m in a csv file to check for steady state free of waves |
---|
104 | sww2csv_gauges('test_flowhydraulics_on_plane.sww','plane_gauge_locations.csv',quantities=['stage','elevation'],verbose=True) |
---|
105 | |
---|
106 | # Compute domain depth in middle of plane at 200m |
---|
107 | x=200.0 |
---|
108 | y=10.00 |
---|
109 | w = domain.get_quantity('stage').get_values(interpolation_points=[[x, y]])[0] |
---|
110 | z = domain.get_quantity('elevation').get_values(interpolation_points=[[x, y]])[0] |
---|
111 | domain_depth = w-z |
---|
112 | |
---|
113 | msg = 'Predicted depth of flow at 200m was %f, should have been %f' % (domain_depth, normal_depth) |
---|
114 | if verbose: |
---|
115 | print 'Depth at 200m: ANUGA = %f, Mannings = %f' % (domain_depth, normal_depth) |
---|
116 | assert num.allclose(domain_depth,normal_depth, rtol=5.0e-2), msg |
---|
117 | |
---|
118 | #------------------------------------------------------------------------------ |
---|
119 | # Compute flow thru flowlines at 200m - confirm equals inflow |
---|
120 | #------------------------------------------------------------------------------ |
---|
121 | |
---|
122 | # Square on flowline at 200m |
---|
123 | q=domain.get_flow_through_cross_section([[200.0,0.0],[200.0,20.0]]) |
---|
124 | msg = 'Predicted square on flow was %f, should have been %f' % (q, ref_flow) |
---|
125 | if verbose: |
---|
126 | print '90 degree flowline: ANUGA = %f, Ref = %f' % (q, ref_flow) |
---|
127 | assert num.allclose(q, ref_flow, rtol=1.0e-2), msg |
---|
128 | |
---|
129 | # 45 degree flowline at 200m |
---|
130 | q=domain.get_flow_through_cross_section([[200.0,0.0],[220.0,20.0]]) |
---|
131 | msg = 'Predicted 45 flow was %f, should have been %f' % (q, ref_flow) |
---|
132 | if verbose: |
---|
133 | print '45 degree flowline: ANUGA = %f, Ref = %f' % (q, ref_flow) |
---|
134 | assert num.allclose(q, ref_flow, rtol=1.0e-2), msg |
---|