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

Last change on this file since 7987 was 7987, checked in by habili, 12 years ago

deal with manning

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