source: trunk/anuga_core/source/anuga/structures/testing_inlet_operator.py @ 8093

Last change on this file since 8093 was 8093, checked in by habili, 13 years ago

Updated the inlet operator to deal with variable discharge rates.
test_hydrograph.tms is the test data

File size: 5.1 KB
Line 
1import os.path
2import sys
3
4from anuga.utilities.system_tools import get_pathname_from_package
5from anuga.geometry.polygon_function import Polygon_function
6       
7from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
8from anuga.abstract_2d_finite_volumes.quantity import Quantity
9from anuga.abstract_2d_finite_volumes.util import file_function
10
11import anuga
12
13from anuga.structures.boyd_box_operator import Boyd_box_operator
14from anuga.structures.inlet_operator import Inlet_operator
15                           
16#from anuga.culvert_flows.culvert_routines import boyd_generalised_culvert_model
17     
18from math import pi, pow, sqrt
19
20import numpy as num
21
22
23
24
25
26"""test_that_culvert_runs_rating
27
28This test exercises the culvert and checks values outside rating curve
29are dealt with       
30"""
31
32path = get_pathname_from_package('anuga.culvert_flows')   
33
34length = 40.
35width = 15.
36
37dx = dy = 0.5          # Resolution: Length of subdivisions on both axes
38
39points, vertices, boundary = rectangular_cross(int(length/dx),
40                                               int(width/dy),
41                                               len1=length, 
42                                               len2=width)
43domain = anuga.Domain(points, vertices, boundary)   
44domain.set_name('Testing_inlet_operator')                 # Output name
45domain.set_default_order(2)
46#domain.set_beta(1.5)
47
48
49#----------------------------------------------------------------------
50# Setup initial conditions
51#----------------------------------------------------------------------
52
53def topography(x, y):
54    """Set up a weir
55   
56    A culvert will connect either side
57    """
58    # General Slope of Topography
59    z=-x/1000
60   
61    N = len(x)
62    for i in range(N):
63
64       # Sloping Embankment Across Channel
65        if 5.0 < x[i] < 10.1:
66            # Cut Out Segment for Culvert face               
67            if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0: 
68               z[i]=z[i]
69            else:
70               z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
71        if 10.0 < x[i] < 12.1:
72           z[i] +=  2.5                    # Flat Crest of Embankment
73        if 12.0 < x[i] < 14.5:
74            # Cut Out Segment for Culvert face               
75            if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5:
76               z[i]=z[i]
77            else:
78               z[i] +=  2.5-1.0*(x[i] -12.0) # Sloping D/S Face
79                   
80       
81    return z
82
83
84domain.set_quantity('elevation', topography) 
85domain.set_quantity('friction', 0.01)         # Constant friction
86domain.set_quantity('stage',
87                    expression='elevation')   # Dry initial condition
88
89filename=os.path.join(path, 'example_rating_curve.csv')
90
91
92
93Boyd_box_operator(domain,
94                            end_points=[[9.0, 2.5],[13.0, 2.5]],
95                            losses=1.5,
96                            width=1.5,
97                            apron=5.0,
98                            use_momentum_jet=True,
99                            use_velocity_head=False,
100                            manning=0.013,
101                            verbose=False)
102
103
104line = [[0.0, 5.0], [0.0, 10.0]]
105
106Q = file_function('test_hydrograph.tms', quantities=['hydrograph'])
107
108Inlet_operator(domain, line, Q)
109
110
111
112
113##-----------------------------------------------------------------------
114## Setup boundary conditions
115##-----------------------------------------------------------------------
116
117## Inflow based on Flow Depth and Approaching Momentum
118Bi = anuga.Dirichlet_boundary([2.0, 0.0, 0.0])
119Br = anuga.Reflective_boundary(domain)              # Solid reflective wall
120#Bo = anuga.Dirichlet_boundary([-5, 0, 0])           # Outflow
121
122## Upstream and downstream conditions that will exceed the rating curve
123## I.e produce delta_h outside the range [0, 10] specified in the the
124## file example_rating_curve.csv
125#Btus = anuga.Time_boundary(domain, \
126            #lambda t: [100*num.sin(2*pi*(t-4)/10), 0.0, 0.0])
127#Btds = anuga.Time_boundary(domain, \
128            #lambda t: [-5*(num.cos(2*pi*(t-4)/20)), 0.0, 0.0])
129#domain.set_boundary({'left': Btus, 'right': Btds, 'top': Br, 'bottom': Br})
130domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
131
132
133##-----------------------------------------------------------------------
134## Evolve system through time
135##-----------------------------------------------------------------------
136
137#min_delta_w = sys.maxint
138#max_delta_w = -min_delta_w
139for t in domain.evolve(yieldstep = 1.0, finaltime = 38):
140    domain.write_time()
141
142    #if domain.get_time() > 150.5 and domain.get_time() < 151.5 :
143        #Bi = anuga.Dirichlet_boundary([0.0, 0.0, 0.0])
144        #domain.set_boundary({'left': Bi, 'right': Br, 'top': Br, 'bottom': Br})
145
146    #delta_w = culvert.inlet.stage - culvert.outlet.stage
147   
148    #if delta_w > max_delta_w: max_delta_w = delta_w
149    #if delta_w < min_delta_w: min_delta_w = delta_w
150
151    print domain.volumetric_balance_statistics()
152   
153    pass
154
155## Check that extreme values in rating curve have been exceeded
156## so that we know that condition has been exercised
157#assert min_delta_w < 0
158#assert max_delta_w > 10       
159
160
161#os.remove('Test_culvert.sww')
Note: See TracBrowser for help on using the repository browser.