# source:anuga_core/source/anuga/culvert_flows/Test_Culvert_Flat_Water_Lev.py@6143

Last change on this file since 6143 was 6143, checked in by ole, 15 years ago

Simplified flow protection at the inlet and wrote better diagnostics.

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