source: trunk/anuga_core/source/anuga/structures/testing_outlet_control.py @ 8073

Last change on this file since 8073 was 8073, checked in by steve, 13 years ago

Added in a unit test for inlet operator

  • Property svn:executable set to *
File size: 9.2 KB
Line 
1""" Testing CULVERT (Changing from Horizontal Abstraction to Vertical Abstraction
2
3This example includes a Model Topography that shows a TYPICAL Headwall Configuration
4
5The aim is to change the Culvert Routine to Model more precisely the abstraction
6from a vertical face.
7
8The inflow must include the impact of Approach velocity.
9Similarly the Outflow has MOMENTUM Not just Up welling as in the Horizontal Style
10abstraction
11
12"""
13print 'Starting.... Importing Modules...'
14#------------------------------------------------------------------------------
15# Import necessary modules
16#------------------------------------------------------------------------------
17import anuga
18
19from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
20from anuga.shallow_water.shallow_water_domain import Domain
21from anuga.shallow_water.forcing import Rainfall, Inflow
22#from anuga.shallow_water.forcing import Reflective_boundary
23#from anuga.shallow_water.forcing import Dirichlet_boundary
24#from anuga.shallow_water.forcing import Transmissive_boundary, Time_boundary
25
26#from anuga.culvert_flows.culvert_class import Culvert_flow
27from anuga.structures.boyd_pipe_operator import Boyd_pipe_operator
28from anuga.structures.boyd_box_operator import Boyd_box_operator
29#from anuga.culvert_flows.culvert_routines import weir_orifice_channel_culvert_model
30from math import pi,pow,sqrt
31
32import numpy as num
33
34
35
36
37#------------------------------------------------------------------------------
38# Setup computational domain
39#------------------------------------------------------------------------------
40print 'Setting up domain'
41
42length = 120. #x-Dir
43width = 200.  #y-dir
44
45dx = dy = 2.0          # Resolution: Length of subdivisions on both axes
46#dx = dy = .5           # Resolution: Length of subdivisions on both axes
47#dx = dy = .5           # Resolution: Length of subdivisions on both axes
48#dx = dy = .1           # Resolution: Length of subdivisions on both axes
49
50points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy),
51                                                    len1=length, len2=width)
52domain = Domain(points, vertices, boundary)   
53domain.set_name('Test_Outlet_Ctrl')                 # Output name
54domain.set_default_order(2)
55domain.H0 = 0.01
56domain.tight_slope_limiters = 1
57
58print 'Size', len(domain)
59
60#------------------------------------------------------------------------------
61# Setup initial conditions
62#------------------------------------------------------------------------------
63
64def topography(x, y):
65    """Set up a weir
66   
67    A culvert will connect either side
68    """
69    # General Slope of Topography
70    z=10.0-x/100.0  # % Longitudinal Slope
71   
72    #       NOW Add bits and Pieces to topography
73    bank_hgt=10.0
74    bridge_width = 50.0
75    bank_width = 10.0
76   
77    us_apron_skew = 1.0 # 1.0 = 1 Length: 1 Width, 2.0 = 2 Length : 1 Width
78    us_start_x = 10.0
79    top_start_y = 50.0
80    us_slope = 3.0  #Horiz : Vertic
81    ds_slope = 3.0
82    ds_apron_skew = 1.0 # 1.0 = 1 Length: 1 Width, 2.0 = 2 Length : 1 Width
83    centre_line_y= top_start_y+bridge_width/2.0
84
85    # CALCULATE PARAMETERS TO FORM THE EMBANKMENT
86    us_slope_length = bank_hgt*us_slope
87    us_end_x =us_start_x + us_slope_length
88    us_toe_start_y =top_start_y - us_slope_length / us_apron_skew
89    us_toe_end_y = top_start_y + bridge_width + us_slope_length / us_apron_skew
90
91    top_end_y = top_start_y + bridge_width
92    ds_slope_length = bank_hgt*ds_slope
93    ds_start_x = us_end_x + bank_width
94    ds_end_x = ds_start_x + ds_slope_length
95
96    ds_toe_start_y =top_start_y - ds_slope_length / ds_apron_skew
97    ds_toe_end_y = top_start_y + bridge_width + ds_slope_length / ds_apron_skew
98
99
100    N = len(x)
101    for i in range(N):
102
103       # Sloping Embankment Across Channel
104        if us_start_x < x[i] < us_end_x +0.1:   # For UPSLOPE on the Upstream FACE
105        #if 5.0 < x[i] < 10.1: # For a Range of X, and over a Range of Y based on X adjust Z
106            if us_toe_start_y +(x[i] - us_start_x)/us_apron_skew < y[i] < us_toe_end_y - (x[i] - us_start_x)/us_apron_skew:
107                #if  49.0+(x[i]-5.0)/5.0 <  y[i]  < 151.0 - (x[i]-5.0)/5.0: # Cut Out Base Segment for Culvert FACE
108                 z[i]=z[i] # Flat Apron
109                #z[i] += z[i] + (x[i] - us_start_x)/us_slope
110                #pass
111            else:
112               z[i] += z[i] + (x[i] - us_start_x)/us_slope    # Sloping Segment  U/S Face
113        if us_end_x < x[i] < ds_start_x + 0.1:
114           z[i] +=  z[i]+bank_hgt        # Flat Crest of Embankment
115        if ds_start_x < x[i] < ds_end_x: # DOWN SDLOPE Segment on Downstream face
116            if  top_start_y-(x[i]-ds_start_x)/ds_apron_skew <  y[i]  < top_end_y + (x[i]-ds_start_x)/ds_apron_skew: # Cut Out Segment for Culvert FACE
117                 z[i]=z[i] # Flat Apron
118                #z[i] += z[i]+bank_hgt-(x[i] -ds_start_x)/ds_slope
119                #pass
120            else:
121               z[i] += z[i]+bank_hgt-(x[i] -ds_start_x)/ds_slope       # Sloping D/S Face
122           
123       
124
125    return z
126
127print 'Setting Quantities....'
128domain.set_quantity('elevation', topography)  # Use function for elevation
129domain.set_quantity('friction', 0.01)         # Constant friction
130domain.set_quantity('stage',
131                    expression='elevation')   # Dry initial condition
132
133
134
135
136#------------------------------------------------------------------------------
137# Setup specialised forcing terms
138#------------------------------------------------------------------------------
139
140#------------------------------------------------------------------------------
141# Setup CULVERT INLETS and OUTLETS in Current Topography
142#------------------------------------------------------------------------------
143print 'Defining Structures'
144
145#  DEFINE CULVERT INLET AND OUTLETS
146
147
148#culvert0 = Culvert_operator(domain,
149#                            end_point0=[40.0, 75.0],
150#                            end_point1=[50.0, 75.0],
151#                            width=50.0,
152#                            depth=10.0,
153#                            apron=5.0,
154#                            verbose=False)
155
156
157#------------------------------------------------------------------------------
158# Setup culverts
159#------------------------------------------------------------------------------
160
161culverts = []
162number_of_culverts = 2
163for i in range(number_of_culverts):
164    culvert_width = 50.0/number_of_culverts
165    y = 100-i*culvert_width - culvert_width/2.0
166    ep0 = num.array([40.0, y])
167    ep1 = num.array([50.0, y])
168   
169    losses = {'inlet':0.5, 'outlet':1, 'bend':0, 'grate':0, 'pier': 0, 'other': 0}
170#    culverts.append(Boyd_pipe_operator(domain,
171#                            end_point0=ep0,
172#                            end_point1=ep1,
173#                            losses=losses,
174#                            diameter=1.5, #culvert_width, #3.658,
175#                            apron=6.0,
176#                            use_momentum_jet=True,
177#                            use_velocity_head=True,
178#                            manning=0.013,
179#                            logging=True,
180#                            label='culvert',
181#                            verbose=False))
182
183
184
185
186
187    Boyd_box_operator(domain,
188                            end_points=[ep0,ep1],
189                            losses=losses,
190                            width=culvert_width,
191                            depth=10.0,
192                            apron=6.0,
193                            use_momentum_jet=True,
194                            use_velocity_head=True,
195                            manning=0.013,
196                            logging=True,
197                            label='box_culvert',
198                            verbose=False)
199
200                       
201
202#losses = {'inlet':1, 'outlet':1, 'bend':1, 'grate':1, 'pier': 1, 'other': 1}
203#culvert2 = Culvert_operator(domain,
204                            #end_point0=[40.0, 62.5],
205                            #end_point1=[50.0, 62.5],
206                            #losses,
207                            #width=25.0,
208                            #depth=10.0,
209                            #apron=5.0,
210                            #manning=0.013,
211                            #verbose=False)
212
213
214
215#------------------------------------------------------------------------------
216# Setup boundary conditions
217#------------------------------------------------------------------------------
218print 'Setting Boundary Conditions'
219Br = anuga.Reflective_boundary(domain)              # Solid reflective wall
220Bi = anuga.Dirichlet_boundary([0.0, 0.0, 0.0])          # Inflow based on Flow Depth and Approaching Momentum !!!
221
222Btus = anuga.Dirichlet_boundary([20.0, 0, 0])           # Outflow water at 10.0
223Btds = anuga.Dirichlet_boundary([19.0, 0, 0])           # Outflow water at 9.0
224domain.set_boundary({'left': Btus, 'right': Btds, 'top': Br, 'bottom': Br})
225
226
227
228#------------------------------------------------------------------------------
229# Evolve system through time
230#------------------------------------------------------------------------------
231
232for t in domain.evolve(yieldstep = 1, finaltime = 100):
233    print domain.timestepping_statistics()
234
235    domain.print_operator_timestepping_statistics()
236
Note: See TracBrowser for help on using the repository browser.