source: trunk/anuga_core/source/anuga/structures/testing_wide_bridge.py @ 8007

Last change on this file since 8007 was 8005, checked in by habili, 14 years ago

using losses

  • Property svn:executable set to *
File size: 10.5 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
27#from anuga.culvert_flows.culvert_routines import weir_orifice_channel_culvert_model
28from math import pi,pow,sqrt
29
30import numpy as num
31
32
33#------------------------------------------------------------------------------
34# Setup computational domain
35#------------------------------------------------------------------------------
36print 'Setting up domain'
37
38length = 200. #x-Dir
39width = 200.  #y-dir
40
41dx = dy = 2.0          # Resolution: Length of subdivisions on both axes
42#dx = dy = .5           # Resolution: Length of subdivisions on both axes
43#dx = dy = .5           # Resolution: Length of subdivisions on both axes
44#dx = dy = .1           # Resolution: Length of subdivisions on both axes
45
46points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy),
47                                                    len1=length, len2=width)
48domain = Domain(points, vertices, boundary)   
49domain.set_name('Test_WIDE_BRIDGE')                 # Output name
50domain.set_default_order(2)
51domain.H0 = 0.01
52domain.tight_slope_limiters = 1
53
54print 'Size', len(domain)
55
56#------------------------------------------------------------------------------
57# Setup initial conditions
58#------------------------------------------------------------------------------
59
60def topography(x, y):
61    """Set up a weir
62   
63    A culvert will connect either side
64    """
65    # General Slope of Topography
66    z=10.0-x/100.0  # % Longitudinal Slope
67   
68    #       NOW Add bits and Pieces to topography
69    bank_hgt=10.0
70    bridge_width = 50.0
71    bank_width = 10.0
72   
73    us_apron_skew = 1.0 # 1.0 = 1 Length: 1 Width, 2.0 = 2 Length : 1 Width
74    us_start_x = 10.0
75    top_start_y = 50.0
76    us_slope = 3.0  #Horiz : Vertic
77    ds_slope = 3.0
78    ds_apron_skew = 1.0 # 1.0 = 1 Length: 1 Width, 2.0 = 2 Length : 1 Width
79    centre_line_y= top_start_y+bridge_width/2.0
80
81    # CALCULATE PARAMETERS TO FORM THE EMBANKMENT
82    us_slope_length = bank_hgt*us_slope
83    us_end_x =us_start_x + us_slope_length
84    us_toe_start_y =top_start_y - us_slope_length / us_apron_skew
85    us_toe_end_y = top_start_y + bridge_width + us_slope_length / us_apron_skew
86
87    top_end_y = top_start_y + bridge_width
88    ds_slope_length = bank_hgt*ds_slope
89    ds_start_x = us_end_x + bank_width
90    ds_end_x = ds_start_x + ds_slope_length
91
92    ds_toe_start_y =top_start_y - ds_slope_length / ds_apron_skew
93    ds_toe_end_y = top_start_y + bridge_width + ds_slope_length / ds_apron_skew
94
95
96    N = len(x)
97    for i in range(N):
98
99       # Sloping Embankment Across Channel
100        if us_start_x < x[i] < us_end_x +0.1:   # For UPSLOPE on the Upstream FACE
101        #if 5.0 < x[i] < 10.1: # For a Range of X, and over a Range of Y based on X adjust Z
102            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:
103                #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
104                 z[i]=z[i] # Flat Apron
105                #z[i] += z[i] + (x[i] - us_start_x)/us_slope
106                #pass
107            else:
108               z[i] += z[i] + (x[i] - us_start_x)/us_slope    # Sloping Segment  U/S Face
109        if us_end_x < x[i] < ds_start_x + 0.1:
110           z[i] +=  z[i]+bank_hgt        # Flat Crest of Embankment
111        if ds_start_x < x[i] < ds_end_x: # DOWN SDLOPE Segment on Downstream face
112            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
113                 z[i]=z[i] # Flat Apron
114                #z[i] += z[i]+bank_hgt-(x[i] -ds_start_x)/ds_slope
115                #pass
116            else:
117               z[i] += z[i]+bank_hgt-(x[i] -ds_start_x)/ds_slope       # Sloping D/S Face
118           
119       
120
121    return z
122
123print 'Setting Quantities....'
124domain.set_quantity('elevation', topography)  # Use function for elevation
125domain.set_quantity('friction', 0.01)         # Constant friction
126domain.set_quantity('stage',
127                    expression='elevation')   # Dry initial condition
128
129
130
131
132#------------------------------------------------------------------------------
133# Setup specialised forcing terms
134#------------------------------------------------------------------------------
135
136#------------------------------------------------------------------------------
137# Setup CULVERT INLETS and OUTLETS in Current Topography
138#------------------------------------------------------------------------------
139print 'DEFINING any Structures if Required'
140
141#  DEFINE CULVERT INLET AND OUTLETS
142
143"""
144culvert_rating = Culvert_flow(domain,
145    culvert_description_filename='example_rating_curve.csv',
146    end_point0=[0.0, 75.0],
147    end_point1=[50.0, 75.0],
148    verbose=True)
149
150culvert_energy = Culvert_flow(domain,
151    label='Culvert No. 1',
152    description='This culvert is a test unit 4m diameter',   
153    end_point0=[40.0, 75.0],
154    end_point1=[50.0, 75.0],
155    width=4.0,
156    culvert_routine=boyd_generalised_culvert_model,  #culvert_routine=weir_orifice_channel_culvert_model,     
157    number_of_barrels=1,
158    number_of_smoothing_steps=10,
159    #update_interval=0.25,
160    log_file=True,
161    discharge_hydrograph=True,
162    use_velocity_head=False,
163    use_momentum_jet=False,
164    verbose=True)
165
166domain.forcing_terms.append(culvert_energy)
167"""
168#from anuga.structures.boyd_box_operator import Boyd_box_operator
169#culvert0 = Culvert_operator(domain,
170#                            end_point0=[40.0, 75.0],
171#                            end_point1=[50.0, 75.0],
172#                            width=50.0,
173#                            height=10.0,
174#                            apron=5.0,
175#                            verbose=False)
176
177
178#------------------------------------------------------------------------------
179# Setup culverts
180#------------------------------------------------------------------------------
181culverts = []
182number_of_culverts = 1 
183for i in range(number_of_culverts):
184    culvert_width = 50.0/number_of_culverts
185    y = 100-i*culvert_width - culvert_width/2.0
186    ep0 = [40.0, y]
187    ep1 = [50.0, y]
188    culverts.append(anuga.Boyd_box_operator(domain,
189                            end_point0=ep0,
190                            end_point1=ep1,
191                            losses=1.5,
192                            width=3.658,
193                            height=3.658,
194                            apron=5.0,
195                            use_momentum_jet=True,
196                            use_velocity_head=True,
197                            manning=0.013,
198                            description='bridge culvert',
199                            verbose=False))
200
201
202
203
204#culvert2 = Culvert_operator(domain,
205#                            end_point0=[40.0, 62.5],
206#                            end_point1=[50.0, 62.5],
207#                            width=25.0,
208#                            height=10.0,
209#                            apron=5.0,
210#                            manning=0.013,
211#                            verbose=False)
212
213
214
215"""
216culvert_energy = Culvert_flow(domain,
217    label='Culvert No. 1',
218    description='This culvert is a test unit 50m Wide by 5m High',   
219    end_point0=[40.0, 75.0],
220    end_point1=[50.0, 75.0],
221    width=50.0,height=5.0,
222    culvert_routine=boyd_generalised_culvert_model,  #culvert_routine=weir_orifice_channel_culvert_model,     
223    number_of_barrels=1,
224    number_of_smoothing_steps=10,
225    #update_interval=0.25,
226    log_file=True,
227    discharge_hydrograph=True,
228    use_velocity_head=False,
229    use_momentum_jet=False,
230    verbose=True)
231
232domain.forcing_terms.append(culvert_energy)
233"""
234
235#------------------------------------------------------------------------------
236# Setup boundary conditions
237#------------------------------------------------------------------------------
238print 'Setting Boundary Conditions'
239Br = anuga.Reflective_boundary(domain)              # Solid reflective wall
240Bi = anuga.Dirichlet_boundary([0.0, 0.0, 0.0])          # Inflow based on Flow Depth and Approaching Momentum !!!
241
242Bo = anuga.Dirichlet_boundary([-5.0, 0, 0])           # Outflow water at -5.0
243Bd = anuga.Dirichlet_boundary([0,0,0])                # Outflow water at 0.0
244
245#Btus = Time_boundary(domain, lambda t: [0.0+ 1.025*(1+num.sin(2*pi*(t-4)/10)), 0.0, 0.0])
246#Btds = Time_boundary(domain, lambda t: [0.0+ 0.0075*(1+num.sin(2*pi*(t-4)/20)), 0.0, 0.0])
247
248Btus = anuga.Dirichlet_boundary([20.0, 0, 0])           # Outflow water at 5.0
249Btds = anuga.Dirichlet_boundary([19.0, 0, 0])           # Outflow water at 1.0
250domain.set_boundary({'left': Btus, 'right': Btds, 'top': Br, 'bottom': Br})
251
252
253#------------------------------------------------------------------------------
254# Evolve system through time
255#------------------------------------------------------------------------------
256
257for t in domain.evolve(yieldstep = 1, finaltime = 10):
258    print domain.timestepping_statistics()
259    print domain.volumetric_balance_statistics()
260    for i, culvert in enumerate(culverts):
261        print 'culvert: ', i
262        print culvert.structure_statistics()
263   
264
265   
266"""   
267#import sys; sys.exit()
268# Profiling code
269import time
270t0 = time.time()
271   
272s = 'for t in domain.evolve(yieldstep = 0.1, finaltime = 300): domain.write_time()'
273
274import profile, pstats
275FN = 'profile.dat'
276
277profile.run(s, FN)
278   
279print 'That took %.2f seconds' %(time.time()-t0)
280
281S = pstats.Stats(FN)
282#S.sort_stats('time').print_stats(20)
283s = S.sort_stats('cumulative').print_stats(30)
284
285print s
286"""
Note: See TracBrowser for help on using the repository browser.