source: anuga_work/development/ted_rigby/test_flowhydraulics_on_plane.py @ 7475

Last change on this file since 7475 was 7475, checked in by ole, 15 years ago

Added tests from Ted Rigby that cause unexpected flows.

File size: 6.2 KB
Line 
1
2
3"""test_flow_on_plane.py
4Test the ability of ANUGA to maintain flow and match normal depth on an inclined plane
5using constant inflow thru the upstream bdry at normal_depth and equivalent momentum
6The plane is 20m wide by 400m long dipping at various slopes with a range of roughnesses
7A recording gauges is located  200m downstream
8"""
9
10verbose = True
11import numpy as num
12
13#------------------------------------------------------------------------------
14# Import necessary modules
15#------------------------------------------------------------------------------
16from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
17from anuga.shallow_water import Domain
18from anuga.shallow_water.shallow_water_domain import Reflective_boundary
19from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
20from anuga.shallow_water.shallow_water_domain import Transmissive_Momentum_Set_Stage_boundary
21from anuga.shallow_water.shallow_water_domain import Inflow
22from anuga.shallow_water.data_manager import get_flow_through_cross_section
23from anuga.abstract_2d_finite_volumes.util import sww2csv_gauges, csv2timeseries_graphs
24
25
26#------------------------------------------------------------------------------
27# Setup computational domain
28#------------------------------------------------------------------------------
29
30finaltime = 1800.0    # 1/2 hour
31
32length   = 400.
33width    = 20.
34dx = dy  = 2          # Resolution: of grid on both axes
35ref_flow = 20.0       # Inflow at head of plane
36
37points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy),
38                                               len1=length, len2=width)
39
40for 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         
Note: See TracBrowser for help on using the repository browser.