source: branches/numpy/anuga/culvert_flows/Test_Culvert_Flat_Water_Lev.py @ 7177

Last change on this file since 7177 was 6553, checked in by rwilson, 15 years ago

Merged trunk into numpy, all tests and validations work.

File size: 6.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#------------------------------------------------------------------------------
17from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
18
19from anuga.shallow_water import Domain, Reflective_boundary,\
20     Dirichlet_boundary,\
21     Transmissive_boundary, Time_boundary
22
23from anuga.culvert_flows.culvert_class import Culvert_flow
24from anuga.culvert_flows.culvert_routines import boyd_generalised_culvert_model
25     
26from math import pi,pow,sqrt
27
28import numpy as num
29
30
31#------------------------------------------------------------------------------
32# Setup computational domain
33#------------------------------------------------------------------------------
34print 'Setting up domain'
35
36length = 40.
37width = 5.
38
39dx = dy = 1           # Resolution: Length of subdivisions on both axes
40#dx = dy = .5           # Resolution: Length of subdivisions on both axes
41#dx = dy = .5           # Resolution: Length of subdivisions on both axes
42#dx = dy = .1           # Resolution: Length of subdivisions on both axes
43
44points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy),
45                                               len1=length, len2=width)
46domain = Domain(points, vertices, boundary)   
47domain.set_name('Test_Culv_Flat_WL')                 # Output name
48domain.set_default_order(2)
49domain.H0 = 0.01
50domain.tight_slope_limiters = 1
51
52print 'Size', len(domain)
53
54#------------------------------------------------------------------------------
55# Setup initial conditions
56#------------------------------------------------------------------------------
57
58def topography(x, y):
59    """Set up a weir
60   
61    A culvert will connect either side
62    """
63    # General Slope of Topography
64    z=-x/1000
65   
66    #       NOW Add bits and Pieces to topography
67    N = len(x)
68    for i in range(N):
69
70       # Sloping Embankment Across Channel
71        if 5.0 < x[i] < 10.1:
72            if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0: # Cut Out Segment for Culvert FACE
73               z[i]=z[i]
74            else:
75               z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
76        if 10.0 < x[i] < 12.1:
77           z[i] +=  2.5              # Flat Crest of Embankment
78        if 12.0 < x[i] < 14.5:
79            if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5: # Cut Out Segment for Culvert FACE
80               z[i]=z[i]
81            else:
82               z[i] +=  2.5-1.0*(x[i] -12.0)       # Sloping D/S Face
83                           
84       
85               
86    return z
87
88print 'Setting Quantities....'
89domain.set_quantity('elevation', topography)  # Use function for elevation
90domain.set_quantity('friction', 0.01)         # Constant friction
91domain.set_quantity('stage',
92                    expression='elevation')   # Dry initial condition
93
94
95
96
97#------------------------------------------------------------------------------
98# Setup specialised forcing terms
99#------------------------------------------------------------------------------
100
101#------------------------------------------------------------------------------
102# Setup CULVERT INLETS and OUTLETS in Current Topography
103#------------------------------------------------------------------------------
104print 'DEFINING any Structures if Required'
105
106#  DEFINE CULVERT INLET AND OUTLETS
107
108
109culvert_rating = Culvert_flow(domain,
110                       culvert_description_filename='example_rating_curve.csv',
111                       end_point0=[9.0, 2.5], 
112                       end_point1=[13.0, 2.5],
113                       verbose=True)
114
115
116culvert_energy = Culvert_flow(domain,
117                       label='Culvert No. 1',
118                       description='This culvert is a test unit 1.2m Wide by 0.75m High',   
119                       end_point0=[9.0, 2.5], 
120                       end_point1=[13.0, 2.5],
121                       width=1.20,height=0.75,
122                       culvert_routine=boyd_generalised_culvert_model,       
123                       number_of_barrels=1,
124                       update_interval=2,
125                       log_file=True,
126                       discharge_hydrograph=True,
127                       verbose=True)
128
129domain.forcing_terms.append(culvert_energy)
130
131#------------------------------------------------------------------------------
132# Setup boundary conditions
133#------------------------------------------------------------------------------
134print 'Setting Boundary Conditions'
135Bi = Dirichlet_boundary([0.0, 0.0, 0.0])          # Inflow based on Flow Depth and Approaching Momentum !!!
136Br = Reflective_boundary(domain)              # Solid reflective wall
137Bo = Dirichlet_boundary([-5, 0, 0])           # Outflow
138Btus = Time_boundary(domain, lambda t: [0.0+ 1.25*(1+num.sin(2*pi*(t-4)/10)), 0.0, 0.0])
139Btds = Time_boundary(domain, lambda t: [0.0+ 0.75*(1+num.sin(2*pi*(t-4)/20)), 0.0, 0.0])
140domain.set_boundary({'left': Btus, 'right': Btds, 'top': Br, 'bottom': Br})
141
142
143#------------------------------------------------------------------------------
144# Evolve system through time
145#------------------------------------------------------------------------------
146
147#for t in domain.evolve(yieldstep = 1, finaltime = 25):
148#    print domain.timestepping_statistics()
149   
150
151   
152   
153#import sys; sys.exit()
154# Profiling code
155import time
156t0 = time.time()
157   
158s = 'for t in domain.evolve(yieldstep = 1, finaltime = 25): domain.write_time()'
159
160import profile, pstats
161FN = 'profile.dat'
162
163profile.run(s, FN)
164   
165print 'That took %.2f seconds' %(time.time()-t0)
166
167S = pstats.Stats(FN)
168#S.sort_stats('time').print_stats(20)
169s = S.sort_stats('cumulative').print_stats(30)
170
171print s
Note: See TracBrowser for help on using the repository browser.