1 | """ Testing CULVERT (Changing from Horizontal Abstraction to Vertical Abstraction |
---|
2 | |
---|
3 | This example includes a Model Topography that shows a TYPICAL Headwall Configuration |
---|
4 | |
---|
5 | The aim is to change the Culvert Routine to Model more precisely the abstraction |
---|
6 | from a vertical face. |
---|
7 | |
---|
8 | The inflow must include the impact of Approach velocity. |
---|
9 | Similarly the Outflow has MOMENTUM Not just Up welling as in the Horizontal Style |
---|
10 | abstraction |
---|
11 | |
---|
12 | """ |
---|
13 | print 'Starting.... Importing Modules...' |
---|
14 | #------------------------------------------------------------------------------ |
---|
15 | # Import necessary modules |
---|
16 | #------------------------------------------------------------------------------ |
---|
17 | from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross |
---|
18 | |
---|
19 | from anuga.shallow_water import Domain, Reflective_boundary,\ |
---|
20 | Dirichlet_boundary,\ |
---|
21 | Transmissive_boundary, Time_boundary |
---|
22 | |
---|
23 | from anuga.culvert_flows.culvert_class import Culvert_flow |
---|
24 | from anuga.culvert_flows.culvert_routines import boyd_generalised_culvert_model |
---|
25 | |
---|
26 | from math import pi,pow,sqrt |
---|
27 | |
---|
28 | import Numeric as num |
---|
29 | |
---|
30 | |
---|
31 | #------------------------------------------------------------------------------ |
---|
32 | # Setup computational domain |
---|
33 | #------------------------------------------------------------------------------ |
---|
34 | print 'Setting up domain' |
---|
35 | |
---|
36 | length = 40. |
---|
37 | width = 5. |
---|
38 | |
---|
39 | dx = 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 | |
---|
44 | points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy), |
---|
45 | len1=length, len2=width) |
---|
46 | domain = Domain(points, vertices, boundary) |
---|
47 | domain.set_name('Test_Culv_Flat_WL') # Output name |
---|
48 | domain.set_default_order(2) |
---|
49 | domain.H0 = 0.01 |
---|
50 | domain.tight_slope_limiters = 1 |
---|
51 | |
---|
52 | print 'Size', len(domain) |
---|
53 | |
---|
54 | #------------------------------------------------------------------------------ |
---|
55 | # Setup initial conditions |
---|
56 | #------------------------------------------------------------------------------ |
---|
57 | |
---|
58 | def 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 | |
---|
88 | print 'Setting Quantities....' |
---|
89 | domain.set_quantity('elevation', topography) # Use function for elevation |
---|
90 | domain.set_quantity('friction', 0.01) # Constant friction |
---|
91 | domain.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 | #------------------------------------------------------------------------------ |
---|
104 | print 'DEFINING any Structures if Required' |
---|
105 | |
---|
106 | # DEFINE CULVERT INLET AND OUTLETS |
---|
107 | |
---|
108 | |
---|
109 | culvert_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 | |
---|
116 | culvert_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 | |
---|
129 | domain.forcing_terms.append(culvert_energy) |
---|
130 | |
---|
131 | #------------------------------------------------------------------------------ |
---|
132 | # Setup boundary conditions |
---|
133 | #------------------------------------------------------------------------------ |
---|
134 | print 'Setting Boundary Conditions' |
---|
135 | Bi = Dirichlet_boundary([0.0, 0.0, 0.0]) # Inflow based on Flow Depth and Approaching Momentum !!! |
---|
136 | Br = Reflective_boundary(domain) # Solid reflective wall |
---|
137 | Bo = Dirichlet_boundary([-5, 0, 0]) # Outflow |
---|
138 | Btus = Time_boundary(domain, lambda t: [0.0+ 1.25*(1+num.sin(2*pi*(t-4)/10)), 0.0, 0.0]) |
---|
139 | Btds = Time_boundary(domain, lambda t: [0.0+ 0.75*(1+num.sin(2*pi*(t-4)/20)), 0.0, 0.0]) |
---|
140 | domain.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 |
---|
155 | import time |
---|
156 | t0 = time.time() |
---|
157 | |
---|
158 | s = 'for t in domain.evolve(yieldstep = 1, finaltime = 25): domain.write_time()' |
---|
159 | |
---|
160 | import profile, pstats |
---|
161 | FN = 'profile.dat' |
---|
162 | |
---|
163 | profile.run(s, FN) |
---|
164 | |
---|
165 | print 'That took %.2f seconds' %(time.time()-t0) |
---|
166 | |
---|
167 | S = pstats.Stats(FN) |
---|
168 | #S.sort_stats('time').print_stats(20) |
---|
169 | s = S.sort_stats('cumulative').print_stats(30) |
---|
170 | |
---|
171 | print s |
---|