source: trunk/anuga_core/source/anuga/structures/run_wide_bridge.py @ 8851

Last change on this file since 8851 was 8851, checked in by steve, 11 years ago

Just a change in name from testing_ run_

File size: 10.6 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#                            depth=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 = num.array([40.0, y])
187    ep1 = num.array([50.0, y])
188    culverts.append(anuga.Boyd_box_operator(domain,
189                                            losses=1.5,
190                                            width=3.658,
191                                            height=3.658,
192                                            end_points=[ep0, ep1],
193                                            apron=0.5,
194                                            manning=0.013,
195                                            enquiry_gap=1.0,
196                                            description='bridge culvert',
197                                            verbose=False))
198
199#culvert2 = Culvert_operator(domain,
200#                            end_point0=[40.0, 62.5],
201#                            end_point1=[50.0, 62.5],
202#                            width=25.0,
203#                            depth=10.0,
204#                            apron=5.0,
205#                            manning=0.013,
206#                            verbose=False)
207
208
209
210"""
211culvert_energy = Culvert_flow(domain,
212    label='Culvert No. 1',
213    description='This culvert is a test unit 50m Wide by 5m High',   
214    end_point0=[40.0, 75.0],
215    end_point1=[50.0, 75.0],
216    width=50.0,depth=5.0,
217    culvert_routine=boyd_generalised_culvert_model,  #culvert_routine=weir_orifice_channel_culvert_model,     
218    number_of_barrels=1,
219    number_of_smoothing_steps=10,
220    #update_interval=0.25,
221    log_file=True,
222    discharge_hydrograph=True,
223    use_velocity_head=False,
224    use_momentum_jet=False,
225    verbose=True)
226
227domain.forcing_terms.append(culvert_energy)
228"""
229
230#------------------------------------------------------------------------------
231# Setup boundary conditions
232#------------------------------------------------------------------------------
233print 'Setting Boundary Conditions'
234Br = anuga.Reflective_boundary(domain)              # Solid reflective wall
235Bi = anuga.Dirichlet_boundary([0.0, 0.0, 0.0])          # Inflow based on Flow Depth and Approaching Momentum !!!
236
237Bo = anuga.Dirichlet_boundary([-5.0, 0, 0])           # Outflow water at -5.0
238Bd = anuga.Dirichlet_boundary([0,0,0])                # Outflow water at 0.0
239
240#Btus = Time_boundary(domain, lambda t: [0.0+ 1.025*(1+num.sin(2*pi*(t-4)/10)), 0.0, 0.0])
241#Btds = Time_boundary(domain, lambda t: [0.0+ 0.0075*(1+num.sin(2*pi*(t-4)/20)), 0.0, 0.0])
242
243Btus = anuga.Dirichlet_boundary([20.0, 0, 0])           # Outflow water at 20
244Btds = anuga.Dirichlet_boundary([19.0, 0, 0])           # Outflow water at 19
245domain.set_boundary({'left': Btus, 'right': Br, 'top': Br, 'bottom': Br})
246
247
248#------------------------------------------------------------------------------
249# Evolve system through time
250#------------------------------------------------------------------------------
251
252for t in domain.evolve(yieldstep = 1, finaltime = 100):
253    print domain.timestepping_statistics()
254    print domain.volumetric_balance_statistics()
255    for i, culvert in enumerate(culverts):
256        print 'culvert: ', i
257        print culvert.statistics()
258   
259
260   
261"""   
262#import sys; sys.exit()
263# Profiling code
264import time
265t0 = time.time()
266   
267s = 'for t in domain.evolve(yieldstep = 0.1, finaltime = 300): domain.write_time()'
268
269import profile, pstats
270FN = 'profile.dat'
271
272profile.run(s, FN)
273   
274print 'That took %.2f seconds' %(time.time()-t0)
275
276S = pstats.Stats(FN)
277#S.sort_stats('time').print_stats(20)
278s = S.sort_stats('cumulative').print_stats(30)
279
280print s
281"""
Note: See TracBrowser for help on using the repository browser.