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

Last change on this file since 7984 was 7984, checked in by steve, 12 years ago

Added enquiry points to culverts

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