1 | import os.path |
---|
2 | import sys |
---|
3 | |
---|
4 | from anuga.utilities.system_tools import get_pathname_from_package |
---|
5 | from anuga.geometry.polygon_function import Polygon_function |
---|
6 | |
---|
7 | from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross |
---|
8 | from anuga.abstract_2d_finite_volumes.quantity import Quantity |
---|
9 | |
---|
10 | import anuga |
---|
11 | |
---|
12 | from anuga.structures.boyd_pipe_operator import Boyd_pipe_operator |
---|
13 | |
---|
14 | #from anuga.culvert_flows.culvert_routines import boyd_generalised_culvert_model |
---|
15 | |
---|
16 | from math import pi, pow, sqrt |
---|
17 | |
---|
18 | import numpy as num |
---|
19 | |
---|
20 | |
---|
21 | |
---|
22 | |
---|
23 | |
---|
24 | """test_that_culvert_runs_rating |
---|
25 | |
---|
26 | This test exercises the culvert and checks values outside rating curve |
---|
27 | are dealt with |
---|
28 | """ |
---|
29 | |
---|
30 | path = get_pathname_from_package('anuga.culvert_flows') |
---|
31 | |
---|
32 | length = 40. |
---|
33 | width = 15. |
---|
34 | |
---|
35 | dx = dy = 0.5 # Resolution: Length of subdivisions on both axes |
---|
36 | |
---|
37 | points, vertices, boundary = rectangular_cross(int(length/dx), |
---|
38 | int(width/dy), |
---|
39 | len1=length, |
---|
40 | len2=width) |
---|
41 | domain = anuga.Domain(points, vertices, boundary) |
---|
42 | domain.set_name('Test_culvert') # Output name |
---|
43 | domain.set_default_order(2) |
---|
44 | #domain.set_beta(1.5) |
---|
45 | |
---|
46 | |
---|
47 | #---------------------------------------------------------------------- |
---|
48 | # Setup initial conditions |
---|
49 | #---------------------------------------------------------------------- |
---|
50 | |
---|
51 | def topography(x, y): |
---|
52 | """Set up a weir |
---|
53 | |
---|
54 | A culvert will connect either side |
---|
55 | """ |
---|
56 | # General Slope of Topography |
---|
57 | z=-x/1000 |
---|
58 | |
---|
59 | N = len(x) |
---|
60 | for i in range(N): |
---|
61 | |
---|
62 | # Sloping Embankment Across Channel |
---|
63 | if 5.0 < x[i] < 10.1: |
---|
64 | # Cut Out Segment for Culvert face |
---|
65 | if 1.0+(x[i]-5.0)/5.0 < y[i] < 4.0 - (x[i]-5.0)/5.0: |
---|
66 | z[i]=z[i] |
---|
67 | else: |
---|
68 | z[i] += 0.5*(x[i] -5.0) # Sloping Segment U/S Face |
---|
69 | if 10.0 < x[i] < 12.1: |
---|
70 | z[i] += 2.5 # Flat Crest of Embankment |
---|
71 | if 12.0 < x[i] < 14.5: |
---|
72 | # Cut Out Segment for Culvert face |
---|
73 | if 2.0-(x[i]-12.0)/2.5 < y[i] < 3.0 + (x[i]-12.0)/2.5: |
---|
74 | z[i]=z[i] |
---|
75 | else: |
---|
76 | z[i] += 2.5-1.0*(x[i] -12.0) # Sloping D/S Face |
---|
77 | |
---|
78 | |
---|
79 | return z |
---|
80 | |
---|
81 | |
---|
82 | domain.set_quantity('elevation', topography) |
---|
83 | domain.set_quantity('friction', 0.01) # Constant friction |
---|
84 | domain.set_quantity('stage', |
---|
85 | expression='elevation') # Dry initial condition |
---|
86 | |
---|
87 | filename=os.path.join(path, 'example_rating_curve.csv') |
---|
88 | |
---|
89 | |
---|
90 | |
---|
91 | Boyd_pipe_operator(domain, |
---|
92 | end_point0=[9.0, 2.5], |
---|
93 | end_point1=[13.0, 2.5], |
---|
94 | losses=1.5, |
---|
95 | diameter=1.5, |
---|
96 | apron=5.0, |
---|
97 | use_momentum_jet=True, |
---|
98 | use_velocity_head=False, |
---|
99 | manning=0.013, |
---|
100 | verbose=False) |
---|
101 | |
---|
102 | |
---|
103 | #culvert2 = Generic_box_culvert(domain, |
---|
104 | # end_point0=[19.0, 2.5], |
---|
105 | # end_point1=[25.0, 2.5], |
---|
106 | # width=1.00, |
---|
107 | # verbose=False) |
---|
108 | |
---|
109 | |
---|
110 | |
---|
111 | |
---|
112 | #print domain.fractional_step_operators |
---|
113 | |
---|
114 | #domain.apply_fractional_steps() |
---|
115 | |
---|
116 | ##----------------------------------------------------------------------- |
---|
117 | ## Setup boundary conditions |
---|
118 | ##----------------------------------------------------------------------- |
---|
119 | |
---|
120 | ## Inflow based on Flow Depth and Approaching Momentum |
---|
121 | Bi = anuga.Dirichlet_boundary([2.0, 0.0, 0.0]) |
---|
122 | Br = anuga.Reflective_boundary(domain) # Solid reflective wall |
---|
123 | #Bo = anuga.Dirichlet_boundary([-5, 0, 0]) # Outflow |
---|
124 | |
---|
125 | ## Upstream and downstream conditions that will exceed the rating curve |
---|
126 | ## I.e produce delta_h outside the range [0, 10] specified in the the |
---|
127 | ## file example_rating_curve.csv |
---|
128 | #Btus = anuga.Time_boundary(domain, \ |
---|
129 | #lambda t: [100*num.sin(2*pi*(t-4)/10), 0.0, 0.0]) |
---|
130 | #Btds = anuga.Time_boundary(domain, \ |
---|
131 | #lambda t: [-5*(num.cos(2*pi*(t-4)/20)), 0.0, 0.0]) |
---|
132 | #domain.set_boundary({'left': Btus, 'right': Btds, 'top': Br, 'bottom': Br}) |
---|
133 | domain.set_boundary({'left': Bi, 'right': Br, 'top': Br, 'bottom': Br}) |
---|
134 | |
---|
135 | |
---|
136 | ##----------------------------------------------------------------------- |
---|
137 | ## Evolve system through time |
---|
138 | ##----------------------------------------------------------------------- |
---|
139 | |
---|
140 | #min_delta_w = sys.maxint |
---|
141 | #max_delta_w = -min_delta_w |
---|
142 | for t in domain.evolve(yieldstep = 1.0, finaltime = 200): |
---|
143 | domain.write_time() |
---|
144 | |
---|
145 | if domain.get_time() > 150.5 and domain.get_time() < 151.5 : |
---|
146 | Bi = anuga.Dirichlet_boundary([0.0, 0.0, 0.0]) |
---|
147 | domain.set_boundary({'left': Bi, 'right': Br, 'top': Br, 'bottom': Br}) |
---|
148 | |
---|
149 | #delta_w = culvert.inlet.stage - culvert.outlet.stage |
---|
150 | |
---|
151 | #if delta_w > max_delta_w: max_delta_w = delta_w |
---|
152 | #if delta_w < min_delta_w: min_delta_w = delta_w |
---|
153 | |
---|
154 | pass |
---|
155 | |
---|
156 | ## Check that extreme values in rating curve have been exceeded |
---|
157 | ## so that we know that condition has been exercised |
---|
158 | #assert min_delta_w < 0 |
---|
159 | #assert max_delta_w > 10 |
---|
160 | |
---|
161 | |
---|
162 | #os.remove('Test_culvert.sww') |
---|