source: trunk/anuga_core/source/anuga/structures/testing_wide_bridge_old.py @ 7992

Last change on this file since 7992 was 7980, checked in by steve, 14 years ago

Added new tests

  • Property svn:executable set to *
File size: 8.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
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            else:
107               z[i] += z[i] + (x[i] - us_start_x)/us_slope    # Sloping Segment  U/S Face
108        if us_end_x < x[i] < ds_start_x + 0.1:
109           z[i] +=  z[i]+bank_hgt        # Flat Crest of Embankment
110        if ds_start_x < x[i] < ds_end_x: # DOWN SDLOPE Segment on Downstream face
111            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
112               z[i]=z[i] # Flat Apron
113            else:
114               z[i] += z[i]+bank_hgt-(x[i] -ds_start_x)/ds_slope       # Sloping D/S Face
115           
116       
117
118    return z
119
120print 'Setting Quantities....'
121domain.set_quantity('elevation', topography)  # Use function for elevation
122domain.set_quantity('friction', 0.01)         # Constant friction
123domain.set_quantity('stage',
124                    expression='elevation')   # Dry initial condition
125
126
127
128
129#------------------------------------------------------------------------------
130# Setup specialised forcing terms
131#------------------------------------------------------------------------------
132
133#------------------------------------------------------------------------------
134# Setup CULVERT INLETS and OUTLETS in Current Topography
135#------------------------------------------------------------------------------
136print 'DEFINING any Structures if Required'
137
138#  DEFINE CULVERT INLET AND OUTLETS
139
140"""
141culvert_rating = Culvert_flow(domain,
142    culvert_description_filename='example_rating_curve.csv',
143    end_point0=[40.0, 75.0],
144    end_point1=[50.0, 75.0],
145    verbose=True)
146
147culvert_energy = Culvert_flow(domain,
148    label='Culvert No. 1',
149    description='This culvert is a test unit 4m diameter',   
150    end_point0=[40.0, 75.0],
151    end_point1=[50.0, 75.0],
152    width=4.0,
153    culvert_routine=boyd_generalised_culvert_model,  #culvert_routine=weir_orifice_channel_culvert_model,     
154    number_of_barrels=1,
155    number_of_smoothing_steps=10,
156    #update_interval=0.25,
157    log_file=True,
158    discharge_hydrograph=True,
159    use_velocity_head=False,
160    use_momentum_jet=False,
161    verbose=True)
162
163domain.forcing_terms.append(culvert_energy)
164"""   
165culvert_energy = Culvert_flow(domain,
166    label='Culvert No. 1',
167    description='This culvert is a test unit 50m Wide by 5m High',   
168    end_point0=[40.0, 75.0], 
169    end_point1=[50.0, 75.0],
170    width=50.0,height=5.0,
171    culvert_routine=boyd_generalised_culvert_model,  #culvert_routine=weir_orifice_channel_culvert_model,     
172    number_of_barrels=1,
173    number_of_smoothing_steps=10,
174    #update_interval=0.25,
175    log_file=True,
176    discharge_hydrograph=True,
177    use_velocity_head=False,
178    use_momentum_jet=False,
179    verbose=True)
180
181domain.forcing_terms.append(culvert_energy)
182
183#------------------------------------------------------------------------------
184# Setup boundary conditions
185#------------------------------------------------------------------------------
186print 'Setting Boundary Conditions'
187Br = anuga.Reflective_boundary(domain)              # Solid reflective wall
188Bi = anuga.Dirichlet_boundary([0.0, 0.0, 0.0])          # Inflow based on Flow Depth and Approaching Momentum !!!
189
190Bo = anuga.Dirichlet_boundary([-5.0, 0, 0])           # Outflow water at -5.0
191Bd = anuga.Dirichlet_boundary([0,0,0])                # Outflow water at 0.0
192
193#Btus = Time_boundary(domain, lambda t: [0.0+ 1.025*(1+num.sin(2*pi*(t-4)/10)), 0.0, 0.0])
194#Btds = Time_boundary(domain, lambda t: [0.0+ 0.0075*(1+num.sin(2*pi*(t-4)/20)), 0.0, 0.0])
195
196Btus = anuga.Dirichlet_boundary([18.0, 0, 0])           # Outflow water at 5.0
197Btds = anuga.Dirichlet_boundary([0.0, 0, 0])           # Outflow water at 1.0
198domain.set_boundary({'left': Btus, 'right': Btds, 'top': Br, 'bottom': Br})
199
200
201#------------------------------------------------------------------------------
202# Evolve system through time
203#------------------------------------------------------------------------------
204
205for t in domain.evolve(yieldstep = 1, finaltime = 100):
206    print domain.timestepping_statistics()
207   
208
209   
210"""   
211#import sys; sys.exit()
212# Profiling code
213import time
214t0 = time.time()
215   
216s = 'for t in domain.evolve(yieldstep = 0.1, finaltime = 300): domain.write_time()'
217
218import profile, pstats
219FN = 'profile.dat'
220
221profile.run(s, FN)
222   
223print 'That took %.2f seconds' %(time.time()-t0)
224
225S = pstats.Stats(FN)
226#S.sort_stats('time').print_stats(20)
227s = S.sort_stats('cumulative').print_stats(30)
228
229print s
230"""
Note: See TracBrowser for help on using the repository browser.