source: trunk/anuga_core/examples/structures/run_wide_bridge.py @ 9680

Last change on this file since 9680 was 9440, checked in by steve, 10 years ago

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