[6076] | 1 | #!/usr/bin/env python |
---|
| 2 | |
---|
| 3 | |
---|
| 4 | import unittest |
---|
| 5 | import os.path |
---|
[6104] | 6 | import sys |
---|
[6076] | 7 | |
---|
[6102] | 8 | from anuga.utilities.system_tools import get_pathname_from_package |
---|
| 9 | from anuga.utilities.polygon import Polygon_function |
---|
| 10 | |
---|
[6076] | 11 | from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross |
---|
[6102] | 12 | from anuga.abstract_2d_finite_volumes.quantity import Quantity |
---|
[6076] | 13 | |
---|
| 14 | from anuga.shallow_water import Domain, Reflective_boundary,\ |
---|
| 15 | Dirichlet_boundary,\ |
---|
| 16 | Transmissive_boundary, Time_boundary |
---|
| 17 | |
---|
[6123] | 18 | from anuga.culvert_flows.culvert_class import Culvert_flow, Culvert_flow_rating, Culvert_flow_energy |
---|
[6076] | 19 | from anuga.culvert_flows.culvert_routines import boyd_generalised_culvert_model |
---|
| 20 | |
---|
| 21 | from math import pi,pow,sqrt |
---|
| 22 | |
---|
[7276] | 23 | import numpy as num |
---|
[6076] | 24 | |
---|
| 25 | |
---|
[7496] | 26 | # Helper functions |
---|
| 27 | def run_culvert_flow_problem(depth): |
---|
| 28 | """Run flow with culvert given depth |
---|
| 29 | """ |
---|
| 30 | |
---|
| 31 | length = 40. |
---|
| 32 | width = 5. |
---|
| 33 | |
---|
| 34 | dx = dy = 1 # Resolution: Length of subdivisions on both axes |
---|
| 35 | |
---|
| 36 | points, vertices, boundary = rectangular_cross(int(length/dx), |
---|
| 37 | int(width/dy), |
---|
| 38 | len1=length, |
---|
| 39 | len2=width) |
---|
| 40 | domain = Domain(points, vertices, boundary) |
---|
| 41 | domain.set_name('Test_culvert_shallow') # Output name |
---|
| 42 | domain.set_default_order(2) |
---|
| 43 | |
---|
| 44 | |
---|
| 45 | #---------------------------------------------------------------------- |
---|
| 46 | # Setup initial conditions |
---|
| 47 | #---------------------------------------------------------------------- |
---|
| 48 | |
---|
| 49 | def topography(x, y): |
---|
| 50 | """Set up a weir |
---|
| 51 | |
---|
| 52 | A culvert will connect either side |
---|
| 53 | """ |
---|
| 54 | # General Slope of Topography |
---|
| 55 | z=-x/1000 |
---|
| 56 | |
---|
| 57 | N = len(x) |
---|
| 58 | for i in range(N): |
---|
| 59 | |
---|
| 60 | # Sloping Embankment Across Channel |
---|
| 61 | if 5.0 < x[i] < 10.1: |
---|
| 62 | # Cut Out Segment for Culvert face |
---|
| 63 | if 1.0+(x[i]-5.0)/5.0 < y[i] < 4.0 - (x[i]-5.0)/5.0: |
---|
| 64 | z[i]=z[i] |
---|
| 65 | else: |
---|
| 66 | z[i] += 0.5*(x[i] -5.0) # Sloping Segment U/S Face |
---|
| 67 | if 10.0 < x[i] < 12.1: |
---|
| 68 | z[i] += 2.5 # Flat Crest of Embankment |
---|
| 69 | if 12.0 < x[i] < 14.5: |
---|
| 70 | # Cut Out Segment for Culvert face |
---|
| 71 | if 2.0-(x[i]-12.0)/2.5 < y[i] < 3.0 + (x[i]-12.0)/2.5: |
---|
| 72 | z[i]=z[i] |
---|
| 73 | else: |
---|
| 74 | z[i] += 2.5-1.0*(x[i] -12.0) # Sloping D/S Face |
---|
| 75 | |
---|
| 76 | |
---|
| 77 | return z |
---|
| 78 | |
---|
| 79 | |
---|
| 80 | domain.set_quantity('elevation', topography) |
---|
| 81 | domain.set_quantity('friction', 0.01) # Constant friction |
---|
| 82 | domain.set_quantity('stage', |
---|
| 83 | expression='elevation + %f' % depth) # Shallow initial condition |
---|
| 84 | |
---|
| 85 | # Boyd culvert |
---|
| 86 | culvert = Culvert_flow(domain, |
---|
| 87 | label='Culvert No. 1', |
---|
| 88 | description='This culvert is a test unit 1.2m Wide by 0.75m High', |
---|
| 89 | end_point0=[9.0, 2.5], |
---|
| 90 | end_point1=[13.0, 2.5], |
---|
| 91 | width=1.20, height=0.75, |
---|
| 92 | culvert_routine=boyd_generalised_culvert_model, |
---|
| 93 | number_of_barrels=1, |
---|
| 94 | update_interval=2, |
---|
| 95 | verbose=False) |
---|
| 96 | |
---|
| 97 | |
---|
| 98 | domain.forcing_terms.append(culvert) |
---|
| 99 | |
---|
| 100 | |
---|
| 101 | #----------------------------------------------------------------------- |
---|
| 102 | # Setup boundary conditions |
---|
| 103 | #----------------------------------------------------------------------- |
---|
| 104 | |
---|
| 105 | # Inflow based on Flow Depth and Approaching Momentum |
---|
| 106 | |
---|
| 107 | Br = Reflective_boundary(domain) # Solid reflective wall |
---|
| 108 | domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br}) |
---|
| 109 | |
---|
| 110 | |
---|
| 111 | |
---|
| 112 | #----------------------------------------------------------------------- |
---|
| 113 | # Evolve system through time |
---|
| 114 | #----------------------------------------------------------------------- |
---|
| 115 | |
---|
| 116 | #print 'depth', depth |
---|
| 117 | ref_volume = domain.get_quantity('stage').get_integral() |
---|
| 118 | for t in domain.evolve(yieldstep = 0.1, finaltime = 25): |
---|
| 119 | new_volume = domain.get_quantity('stage').get_integral() |
---|
| 120 | |
---|
| 121 | msg = ('Total volume has changed: Is %.8f m^3 should have been %.8f m^3' |
---|
| 122 | % (new_volume, ref_volume)) |
---|
| 123 | assert num.allclose(new_volume, ref_volume), msg |
---|
| 124 | |
---|
| 125 | |
---|
| 126 | |
---|
[6076] | 127 | class Test_Culvert(unittest.TestCase): |
---|
| 128 | def setUp(self): |
---|
| 129 | pass |
---|
| 130 | |
---|
| 131 | def tearDown(self): |
---|
| 132 | pass |
---|
| 133 | |
---|
| 134 | |
---|
[6103] | 135 | def test_that_culvert_runs_rating(self): |
---|
| 136 | """test_that_culvert_runs_rating |
---|
[6076] | 137 | |
---|
[6103] | 138 | This test exercises the culvert and checks values outside rating curve |
---|
| 139 | are dealt with |
---|
[6076] | 140 | """ |
---|
| 141 | |
---|
| 142 | path = get_pathname_from_package('anuga.culvert_flows') |
---|
| 143 | |
---|
| 144 | length = 40. |
---|
| 145 | width = 5. |
---|
| 146 | |
---|
| 147 | dx = dy = 1 # Resolution: Length of subdivisions on both axes |
---|
| 148 | |
---|
| 149 | points, vertices, boundary = rectangular_cross(int(length/dx), |
---|
| 150 | int(width/dy), |
---|
| 151 | len1=length, |
---|
| 152 | len2=width) |
---|
| 153 | domain = Domain(points, vertices, boundary) |
---|
| 154 | domain.set_name('Test_culvert') # Output name |
---|
| 155 | domain.set_default_order(2) |
---|
| 156 | |
---|
| 157 | |
---|
| 158 | #---------------------------------------------------------------------- |
---|
| 159 | # Setup initial conditions |
---|
| 160 | #---------------------------------------------------------------------- |
---|
| 161 | |
---|
| 162 | def topography(x, y): |
---|
| 163 | """Set up a weir |
---|
| 164 | |
---|
| 165 | A culvert will connect either side |
---|
| 166 | """ |
---|
| 167 | # General Slope of Topography |
---|
| 168 | z=-x/1000 |
---|
| 169 | |
---|
| 170 | N = len(x) |
---|
| 171 | for i in range(N): |
---|
| 172 | |
---|
| 173 | # Sloping Embankment Across Channel |
---|
| 174 | if 5.0 < x[i] < 10.1: |
---|
[6144] | 175 | # Cut Out Segment for Culvert face |
---|
[6076] | 176 | if 1.0+(x[i]-5.0)/5.0 < y[i] < 4.0 - (x[i]-5.0)/5.0: |
---|
| 177 | z[i]=z[i] |
---|
| 178 | else: |
---|
| 179 | z[i] += 0.5*(x[i] -5.0) # Sloping Segment U/S Face |
---|
| 180 | if 10.0 < x[i] < 12.1: |
---|
| 181 | z[i] += 2.5 # Flat Crest of Embankment |
---|
| 182 | if 12.0 < x[i] < 14.5: |
---|
[6144] | 183 | # Cut Out Segment for Culvert face |
---|
[6076] | 184 | if 2.0-(x[i]-12.0)/2.5 < y[i] < 3.0 + (x[i]-12.0)/2.5: |
---|
| 185 | z[i]=z[i] |
---|
| 186 | else: |
---|
| 187 | z[i] += 2.5-1.0*(x[i] -12.0) # Sloping D/S Face |
---|
| 188 | |
---|
| 189 | |
---|
| 190 | return z |
---|
| 191 | |
---|
| 192 | |
---|
| 193 | domain.set_quantity('elevation', topography) |
---|
| 194 | domain.set_quantity('friction', 0.01) # Constant friction |
---|
| 195 | domain.set_quantity('stage', |
---|
| 196 | expression='elevation') # Dry initial condition |
---|
| 197 | |
---|
| 198 | filename=os.path.join(path, 'example_rating_curve.csv') |
---|
[6123] | 199 | culvert = Culvert_flow(domain, |
---|
[6076] | 200 | culvert_description_filename=filename, |
---|
| 201 | end_point0=[9.0, 2.5], |
---|
| 202 | end_point1=[13.0, 2.5], |
---|
[7496] | 203 | width=1.00, |
---|
[6128] | 204 | use_velocity_head=True, |
---|
[6076] | 205 | verbose=False) |
---|
| 206 | |
---|
| 207 | domain.forcing_terms.append(culvert) |
---|
| 208 | |
---|
| 209 | |
---|
| 210 | #----------------------------------------------------------------------- |
---|
| 211 | # Setup boundary conditions |
---|
| 212 | #----------------------------------------------------------------------- |
---|
| 213 | |
---|
| 214 | # Inflow based on Flow Depth and Approaching Momentum |
---|
| 215 | Bi = Dirichlet_boundary([0.0, 0.0, 0.0]) |
---|
| 216 | Br = Reflective_boundary(domain) # Solid reflective wall |
---|
| 217 | Bo = Dirichlet_boundary([-5, 0, 0]) # Outflow |
---|
[6103] | 218 | |
---|
| 219 | # Upstream and downstream conditions that will exceed the rating curve |
---|
| 220 | # I.e produce delta_h outside the range [0, 10] specified in the the |
---|
| 221 | # file example_rating_curve.csv |
---|
[6150] | 222 | Btus = Time_boundary(domain, lambda t: [100*num.sin(2*pi*(t-4)/10), 0.0, 0.0]) |
---|
| 223 | Btds = Time_boundary(domain, lambda t: [-5*(num.cos(2*pi*(t-4)/20)), 0.0, 0.0]) |
---|
[6076] | 224 | domain.set_boundary({'left': Btus, 'right': Btds, 'top': Br, 'bottom': Br}) |
---|
| 225 | |
---|
| 226 | |
---|
| 227 | #----------------------------------------------------------------------- |
---|
| 228 | # Evolve system through time |
---|
| 229 | #----------------------------------------------------------------------- |
---|
| 230 | |
---|
[6104] | 231 | min_delta_w = sys.maxint |
---|
| 232 | max_delta_w = -min_delta_w |
---|
[6076] | 233 | for t in domain.evolve(yieldstep = 1, finaltime = 25): |
---|
[6104] | 234 | delta_w = culvert.inlet.stage - culvert.outlet.stage |
---|
| 235 | |
---|
| 236 | if delta_w > max_delta_w: max_delta_w = delta_w |
---|
| 237 | if delta_w < min_delta_w: min_delta_w = delta_w |
---|
| 238 | |
---|
[6076] | 239 | pass |
---|
| 240 | |
---|
[6104] | 241 | # Check that extreme values in rating curve have been exceeded |
---|
| 242 | # so that we know that condition has been exercised |
---|
| 243 | assert min_delta_w < 0 |
---|
| 244 | assert max_delta_w > 10 |
---|
| 245 | |
---|
[6076] | 246 | |
---|
[6121] | 247 | def test_that_culvert_dry_bed_rating_does_not_produce_flow(self): |
---|
[6076] | 248 | """test_that_culvert_in_dry_bed_does_not_produce_flow(self): |
---|
| 249 | |
---|
| 250 | Test that culvert on a sloping dry bed doesn't produce flows |
---|
| 251 | although there will be a 'pressure' head due to delta_w > 0 |
---|
[6079] | 252 | |
---|
| 253 | This one is using the rating curve variant |
---|
[6076] | 254 | """ |
---|
| 255 | |
---|
| 256 | path = get_pathname_from_package('anuga.culvert_flows') |
---|
| 257 | |
---|
| 258 | length = 40. |
---|
| 259 | width = 5. |
---|
| 260 | |
---|
| 261 | dx = dy = 1 # Resolution: Length of subdivisions on both axes |
---|
| 262 | |
---|
| 263 | points, vertices, boundary = rectangular_cross(int(length/dx), |
---|
| 264 | int(width/dy), |
---|
| 265 | len1=length, |
---|
| 266 | len2=width) |
---|
| 267 | domain = Domain(points, vertices, boundary) |
---|
| 268 | domain.set_name('Test_culvert_dry') # Output name |
---|
| 269 | domain.set_default_order(2) |
---|
| 270 | |
---|
| 271 | |
---|
| 272 | #---------------------------------------------------------------------- |
---|
| 273 | # Setup initial conditions |
---|
| 274 | #---------------------------------------------------------------------- |
---|
| 275 | |
---|
| 276 | def topography(x, y): |
---|
| 277 | """Set up a weir |
---|
| 278 | |
---|
| 279 | A culvert will connect either side |
---|
| 280 | """ |
---|
| 281 | # General Slope of Topography |
---|
| 282 | z=-x/1000 |
---|
| 283 | |
---|
| 284 | N = len(x) |
---|
| 285 | for i in range(N): |
---|
| 286 | |
---|
| 287 | # Sloping Embankment Across Channel |
---|
| 288 | if 5.0 < x[i] < 10.1: |
---|
[6144] | 289 | # Cut Out Segment for Culvert face |
---|
[6076] | 290 | if 1.0+(x[i]-5.0)/5.0 < y[i] < 4.0 - (x[i]-5.0)/5.0: |
---|
| 291 | z[i]=z[i] |
---|
| 292 | else: |
---|
| 293 | z[i] += 0.5*(x[i] -5.0) # Sloping Segment U/S Face |
---|
| 294 | if 10.0 < x[i] < 12.1: |
---|
| 295 | z[i] += 2.5 # Flat Crest of Embankment |
---|
| 296 | if 12.0 < x[i] < 14.5: |
---|
[6144] | 297 | # Cut Out Segment for Culvert face |
---|
[6076] | 298 | if 2.0-(x[i]-12.0)/2.5 < y[i] < 3.0 + (x[i]-12.0)/2.5: |
---|
| 299 | z[i]=z[i] |
---|
| 300 | else: |
---|
| 301 | z[i] += 2.5-1.0*(x[i] -12.0) # Sloping D/S Face |
---|
| 302 | |
---|
| 303 | |
---|
| 304 | return z |
---|
| 305 | |
---|
| 306 | |
---|
| 307 | domain.set_quantity('elevation', topography) |
---|
| 308 | domain.set_quantity('friction', 0.01) # Constant friction |
---|
| 309 | domain.set_quantity('stage', |
---|
| 310 | expression='elevation') # Dry initial condition |
---|
| 311 | |
---|
| 312 | |
---|
| 313 | filename = os.path.join(path, 'example_rating_curve.csv') |
---|
[6123] | 314 | culvert = Culvert_flow(domain, |
---|
[6076] | 315 | culvert_description_filename=filename, |
---|
| 316 | end_point0=[9.0, 2.5], |
---|
| 317 | end_point1=[13.0, 2.5], |
---|
[7496] | 318 | height=0.75, |
---|
[6076] | 319 | verbose=False) |
---|
| 320 | |
---|
| 321 | domain.forcing_terms.append(culvert) |
---|
| 322 | |
---|
| 323 | |
---|
| 324 | #----------------------------------------------------------------------- |
---|
| 325 | # Setup boundary conditions |
---|
| 326 | #----------------------------------------------------------------------- |
---|
| 327 | |
---|
| 328 | # Inflow based on Flow Depth and Approaching Momentum |
---|
| 329 | |
---|
| 330 | Br = Reflective_boundary(domain) # Solid reflective wall |
---|
| 331 | domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br}) |
---|
| 332 | |
---|
| 333 | |
---|
| 334 | #----------------------------------------------------------------------- |
---|
| 335 | # Evolve system through time |
---|
| 336 | #----------------------------------------------------------------------- |
---|
| 337 | |
---|
| 338 | ref_volume = domain.get_quantity('stage').get_integral() |
---|
| 339 | for t in domain.evolve(yieldstep = 1, finaltime = 25): |
---|
| 340 | new_volume = domain.get_quantity('stage').get_integral() |
---|
| 341 | |
---|
| 342 | msg = 'Total volume has changed' |
---|
[7492] | 343 | assert num.allclose(new_volume, ref_volume, rtol=1.0e-10), msg |
---|
[6076] | 344 | pass |
---|
| 345 | |
---|
| 346 | |
---|
[6121] | 347 | |
---|
[7496] | 348 | |
---|
| 349 | |
---|
| 350 | def test_that_culvert_flows_conserves_volume(self): |
---|
| 351 | """test_that_culvert_flows_conserves_volume |
---|
| 352 | |
---|
| 353 | Test that culvert on a sloping dry bed limits flows when very little water |
---|
| 354 | is present at inlet. |
---|
| 355 | |
---|
| 356 | Uses helper function: run_culvert_flow_problem(depth): |
---|
| 357 | |
---|
| 358 | """ |
---|
| 359 | |
---|
| 360 | # Try this for a range of depths |
---|
| 361 | for depth in [0.1, 0.2, 0.5, 1.0]: |
---|
| 362 | run_culvert_flow_problem(depth) |
---|
| 363 | |
---|
| 364 | |
---|
| 365 | def OBSOLETE_XXXtest_that_culvert_rating_limits_flow_in_shallow_inlet_condition(self): |
---|
[6121] | 366 | """test_that_culvert_rating_limits_flow_in_shallow_inlet_condition |
---|
| 367 | |
---|
| 368 | Test that culvert on a sloping dry bed limits flows when very little water |
---|
| 369 | is present at inlet |
---|
[6079] | 370 | |
---|
[6121] | 371 | This one is using the rating curve variant |
---|
| 372 | """ |
---|
| 373 | |
---|
[7496] | 374 | |
---|
| 375 | |
---|
[6121] | 376 | path = get_pathname_from_package('anuga.culvert_flows') |
---|
[6079] | 377 | |
---|
[6121] | 378 | length = 40. |
---|
| 379 | width = 5. |
---|
| 380 | |
---|
| 381 | dx = dy = 1 # Resolution: Length of subdivisions on both axes |
---|
| 382 | |
---|
| 383 | points, vertices, boundary = rectangular_cross(int(length/dx), |
---|
| 384 | int(width/dy), |
---|
| 385 | len1=length, |
---|
| 386 | len2=width) |
---|
| 387 | domain = Domain(points, vertices, boundary) |
---|
| 388 | domain.set_name('Test_culvert_shallow') # Output name |
---|
| 389 | domain.set_default_order(2) |
---|
| 390 | |
---|
| 391 | |
---|
| 392 | #---------------------------------------------------------------------- |
---|
| 393 | # Setup initial conditions |
---|
| 394 | #---------------------------------------------------------------------- |
---|
| 395 | |
---|
| 396 | def topography(x, y): |
---|
| 397 | """Set up a weir |
---|
| 398 | |
---|
| 399 | A culvert will connect either side |
---|
| 400 | """ |
---|
| 401 | # General Slope of Topography |
---|
| 402 | z=-x/1000 |
---|
| 403 | |
---|
| 404 | N = len(x) |
---|
| 405 | for i in range(N): |
---|
| 406 | |
---|
| 407 | # Sloping Embankment Across Channel |
---|
| 408 | if 5.0 < x[i] < 10.1: |
---|
[6144] | 409 | # Cut Out Segment for Culvert face |
---|
[6121] | 410 | if 1.0+(x[i]-5.0)/5.0 < y[i] < 4.0 - (x[i]-5.0)/5.0: |
---|
| 411 | z[i]=z[i] |
---|
| 412 | else: |
---|
| 413 | z[i] += 0.5*(x[i] -5.0) # Sloping Segment U/S Face |
---|
| 414 | if 10.0 < x[i] < 12.1: |
---|
| 415 | z[i] += 2.5 # Flat Crest of Embankment |
---|
| 416 | if 12.0 < x[i] < 14.5: |
---|
[6144] | 417 | # Cut Out Segment for Culvert face |
---|
[6121] | 418 | if 2.0-(x[i]-12.0)/2.5 < y[i] < 3.0 + (x[i]-12.0)/2.5: |
---|
| 419 | z[i]=z[i] |
---|
| 420 | else: |
---|
| 421 | z[i] += 2.5-1.0*(x[i] -12.0) # Sloping D/S Face |
---|
| 422 | |
---|
| 423 | |
---|
| 424 | return z |
---|
| 425 | |
---|
| 426 | |
---|
| 427 | domain.set_quantity('elevation', topography) |
---|
| 428 | domain.set_quantity('friction', 0.01) # Constant friction |
---|
| 429 | domain.set_quantity('stage', |
---|
[6144] | 430 | expression='elevation + 0.1') # Shallow initial condition |
---|
[6121] | 431 | |
---|
[7496] | 432 | # Boyd culvert |
---|
[6123] | 433 | culvert = Culvert_flow(domain, |
---|
[7496] | 434 | label='Culvert No. 1', |
---|
| 435 | description='This culvert is a test unit 1.2m Wide by 0.75m High', |
---|
[6123] | 436 | end_point0=[9.0, 2.5], |
---|
| 437 | end_point1=[13.0, 2.5], |
---|
[7496] | 438 | width=1.20, height=0.75, |
---|
| 439 | culvert_routine=boyd_generalised_culvert_model, |
---|
| 440 | number_of_barrels=1, |
---|
| 441 | update_interval=2, |
---|
[6123] | 442 | verbose=False) |
---|
[7496] | 443 | |
---|
| 444 | # Rating curve |
---|
| 445 | #filename = os.path.join(path, 'example_rating_curve.csv') |
---|
| 446 | #culvert = Culvert_flow(domain, |
---|
| 447 | # culvert_description_filename=filename, |
---|
| 448 | # end_point0=[9.0, 2.5], |
---|
| 449 | # end_point1=[13.0, 2.5], |
---|
| 450 | # trigger_depth=0.01, |
---|
| 451 | # verbose=False) |
---|
[6121] | 452 | |
---|
| 453 | domain.forcing_terms.append(culvert) |
---|
| 454 | |
---|
| 455 | |
---|
| 456 | #----------------------------------------------------------------------- |
---|
| 457 | # Setup boundary conditions |
---|
| 458 | #----------------------------------------------------------------------- |
---|
| 459 | |
---|
| 460 | # Inflow based on Flow Depth and Approaching Momentum |
---|
| 461 | |
---|
| 462 | Br = Reflective_boundary(domain) # Solid reflective wall |
---|
| 463 | domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br}) |
---|
| 464 | |
---|
| 465 | |
---|
[6144] | 466 | |
---|
[6121] | 467 | #----------------------------------------------------------------------- |
---|
| 468 | # Evolve system through time |
---|
| 469 | #----------------------------------------------------------------------- |
---|
| 470 | |
---|
[7496] | 471 | print 'depth', 0.1 |
---|
[6121] | 472 | ref_volume = domain.get_quantity('stage').get_integral() |
---|
[6128] | 473 | for t in domain.evolve(yieldstep = 0.1, finaltime = 25): |
---|
[6121] | 474 | new_volume = domain.get_quantity('stage').get_integral() |
---|
[7492] | 475 | |
---|
| 476 | msg = ('Total volume has changed: Is %.8f m^3 should have been %.8f m^3' |
---|
[7317] | 477 | % (new_volume, ref_volume)) |
---|
[7492] | 478 | assert num.allclose(new_volume, ref_volume, rtol=1.0e-10), msg |
---|
[6144] | 479 | |
---|
| 480 | |
---|
[7492] | 481 | return |
---|
[7496] | 482 | # Now try this again for a depth of 10 cm and for a range of other depths |
---|
| 483 | for depth in [0.1, 0.2, 0.5, 1.0]: |
---|
| 484 | print 'depth', depth |
---|
[6144] | 485 | domain.set_time(0.0) |
---|
| 486 | |
---|
[7496] | 487 | domain.set_quantity('elevation', topography) |
---|
| 488 | domain.set_quantity('friction', 0.01) # Constant friction |
---|
[6144] | 489 | domain.set_quantity('stage', |
---|
| 490 | expression='elevation + %f' % depth) |
---|
| 491 | |
---|
| 492 | |
---|
| 493 | ref_volume = domain.get_quantity('stage').get_integral() |
---|
| 494 | for t in domain.evolve(yieldstep = 0.1, finaltime = 25): |
---|
| 495 | new_volume = domain.get_quantity('stage').get_integral() |
---|
| 496 | |
---|
[7496] | 497 | msg = 'Total volume has changed: Is %.8f m^3 should have been %.8f m^3'\ |
---|
[6144] | 498 | % (new_volume, ref_volume) |
---|
| 499 | |
---|
[7492] | 500 | assert num.allclose(new_volume, ref_volume, rtol=1.0e-10), msg |
---|
[6121] | 501 | |
---|
| 502 | |
---|
| 503 | |
---|
| 504 | def test_that_culvert_dry_bed_boyd_does_not_produce_flow(self): |
---|
| 505 | """test_that_culvert_in_dry_bed_boyd_does_not_produce_flow(self): |
---|
| 506 | |
---|
[6079] | 507 | Test that culvert on a sloping dry bed doesn't produce flows |
---|
| 508 | although there will be a 'pressure' head due to delta_w > 0 |
---|
| 509 | |
---|
| 510 | This one is using the 'Boyd' variant |
---|
| 511 | """ |
---|
| 512 | |
---|
| 513 | path = get_pathname_from_package('anuga.culvert_flows') |
---|
| 514 | |
---|
| 515 | length = 40. |
---|
| 516 | width = 5. |
---|
| 517 | |
---|
| 518 | dx = dy = 1 # Resolution: Length of subdivisions on both axes |
---|
| 519 | |
---|
| 520 | points, vertices, boundary = rectangular_cross(int(length/dx), |
---|
| 521 | int(width/dy), |
---|
| 522 | len1=length, |
---|
| 523 | len2=width) |
---|
| 524 | domain = Domain(points, vertices, boundary) |
---|
| 525 | domain.set_name('Test_culvert_dry') # Output name |
---|
| 526 | domain.set_default_order(2) |
---|
| 527 | |
---|
| 528 | |
---|
| 529 | #---------------------------------------------------------------------- |
---|
| 530 | # Setup initial conditions |
---|
| 531 | #---------------------------------------------------------------------- |
---|
| 532 | |
---|
| 533 | def topography(x, y): |
---|
| 534 | """Set up a weir |
---|
| 535 | |
---|
| 536 | A culvert will connect either side |
---|
| 537 | """ |
---|
| 538 | # General Slope of Topography |
---|
| 539 | z=-x/1000 |
---|
| 540 | |
---|
| 541 | N = len(x) |
---|
| 542 | for i in range(N): |
---|
| 543 | |
---|
| 544 | # Sloping Embankment Across Channel |
---|
| 545 | if 5.0 < x[i] < 10.1: |
---|
[6144] | 546 | # Cut Out Segment for Culvert face |
---|
[6079] | 547 | if 1.0+(x[i]-5.0)/5.0 < y[i] < 4.0 - (x[i]-5.0)/5.0: |
---|
| 548 | z[i]=z[i] |
---|
| 549 | else: |
---|
| 550 | z[i] += 0.5*(x[i] -5.0) # Sloping Segment U/S Face |
---|
| 551 | if 10.0 < x[i] < 12.1: |
---|
| 552 | z[i] += 2.5 # Flat Crest of Embankment |
---|
| 553 | if 12.0 < x[i] < 14.5: |
---|
[6144] | 554 | # Cut Out Segment for Culvert face |
---|
[6079] | 555 | if 2.0-(x[i]-12.0)/2.5 < y[i] < 3.0 + (x[i]-12.0)/2.5: |
---|
| 556 | z[i]=z[i] |
---|
| 557 | else: |
---|
| 558 | z[i] += 2.5-1.0*(x[i] -12.0) # Sloping D/S Face |
---|
| 559 | |
---|
| 560 | |
---|
| 561 | return z |
---|
| 562 | |
---|
| 563 | |
---|
| 564 | domain.set_quantity('elevation', topography) |
---|
| 565 | domain.set_quantity('friction', 0.01) # Constant friction |
---|
| 566 | domain.set_quantity('stage', |
---|
| 567 | expression='elevation') # Dry initial condition |
---|
| 568 | |
---|
| 569 | |
---|
| 570 | filename = os.path.join(path, 'example_rating_curve.csv') |
---|
| 571 | |
---|
| 572 | |
---|
[6123] | 573 | culvert = Culvert_flow(domain, |
---|
| 574 | label='Culvert No. 1', |
---|
| 575 | description='This culvert is a test unit 1.2m Wide by 0.75m High', |
---|
| 576 | end_point0=[9.0, 2.5], |
---|
| 577 | end_point1=[13.0, 2.5], |
---|
[7496] | 578 | width=1.20, height=0.75, |
---|
[6123] | 579 | culvert_routine=boyd_generalised_culvert_model, |
---|
| 580 | number_of_barrels=1, |
---|
| 581 | update_interval=2, |
---|
[7496] | 582 | verbose=False) |
---|
[6079] | 583 | |
---|
| 584 | domain.forcing_terms.append(culvert) |
---|
| 585 | |
---|
| 586 | |
---|
| 587 | #----------------------------------------------------------------------- |
---|
| 588 | # Setup boundary conditions |
---|
| 589 | #----------------------------------------------------------------------- |
---|
| 590 | |
---|
| 591 | # Inflow based on Flow Depth and Approaching Momentum |
---|
| 592 | |
---|
| 593 | Br = Reflective_boundary(domain) # Solid reflective wall |
---|
| 594 | domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br}) |
---|
| 595 | |
---|
| 596 | |
---|
| 597 | #----------------------------------------------------------------------- |
---|
| 598 | # Evolve system through time |
---|
| 599 | #----------------------------------------------------------------------- |
---|
| 600 | |
---|
| 601 | ref_volume = domain.get_quantity('stage').get_integral() |
---|
| 602 | for t in domain.evolve(yieldstep = 1, finaltime = 25): |
---|
[7276] | 603 | |
---|
[6079] | 604 | new_volume = domain.get_quantity('stage').get_integral() |
---|
[7276] | 605 | |
---|
[6079] | 606 | msg = 'Total volume has changed' |
---|
[7492] | 607 | assert num.allclose(new_volume, ref_volume, rtol=1.0e-10), msg |
---|
[6079] | 608 | pass |
---|
[6076] | 609 | |
---|
[6079] | 610 | |
---|
| 611 | |
---|
[6102] | 612 | |
---|
| 613 | |
---|
[6128] | 614 | def test_predicted_boyd_flow(self): |
---|
[6111] | 615 | """test_predicted_boyd_flow |
---|
[6102] | 616 | |
---|
[6111] | 617 | Test that flows predicted by the boyd method are consistent with what what |
---|
[6102] | 618 | is calculated in engineering codes. |
---|
| 619 | The data was supplied by Petar Milevski |
---|
| 620 | """ |
---|
| 621 | |
---|
[6128] | 622 | # FIXME(Ole) this is nowhere near finished |
---|
[6102] | 623 | path = get_pathname_from_package('anuga.culvert_flows') |
---|
| 624 | |
---|
| 625 | length = 12. |
---|
| 626 | width = 5. |
---|
| 627 | |
---|
| 628 | dx = dy = 0.5 # Resolution: Length of subdivisions on both axes |
---|
| 629 | |
---|
| 630 | points, vertices, boundary = rectangular_cross(int(length/dx), |
---|
| 631 | int(width/dy), |
---|
| 632 | len1=length, |
---|
| 633 | len2=width) |
---|
| 634 | domain = Domain(points, vertices, boundary) |
---|
| 635 | |
---|
| 636 | domain.set_name('test_culvert') # Output name |
---|
| 637 | domain.set_default_order(2) |
---|
| 638 | |
---|
| 639 | |
---|
| 640 | #---------------------------------------------------------------------- |
---|
| 641 | # Setup initial conditions |
---|
| 642 | #---------------------------------------------------------------------- |
---|
| 643 | |
---|
| 644 | def topography(x, y): |
---|
| 645 | # General Slope of Topography |
---|
| 646 | z=-x/10 |
---|
| 647 | |
---|
| 648 | return z |
---|
| 649 | |
---|
| 650 | |
---|
| 651 | domain.set_quantity('elevation', topography) |
---|
| 652 | domain.set_quantity('friction', 0.01) # Constant friction |
---|
| 653 | domain.set_quantity('stage', expression='elevation') |
---|
| 654 | |
---|
| 655 | |
---|
| 656 | Q0 = domain.get_quantity('stage') |
---|
| 657 | Q1 = Quantity(domain) |
---|
| 658 | |
---|
| 659 | # Add depths to stage |
---|
| 660 | head_water_depth = 0.169 |
---|
| 661 | tail_water_depth = 0.089 |
---|
| 662 | |
---|
| 663 | inlet_poly = [[0,0], [6,0], [6,5], [0,5]] |
---|
| 664 | outlet_poly = [[6,0], [12,0], [12,5], [6,5]] |
---|
| 665 | |
---|
| 666 | Q1.set_values(Polygon_function([(inlet_poly, head_water_depth), |
---|
| 667 | (outlet_poly, tail_water_depth)])) |
---|
| 668 | |
---|
| 669 | domain.set_quantity('stage', Q0 + Q1) |
---|
| 670 | |
---|
| 671 | |
---|
[6111] | 672 | |
---|
[6123] | 673 | culvert = Culvert_flow(domain, |
---|
| 674 | label='Test culvert', |
---|
| 675 | description='4 m test culvert', |
---|
| 676 | end_point0=[4.0, 2.5], |
---|
| 677 | end_point1=[8.0, 2.5], |
---|
| 678 | width=1.20, |
---|
| 679 | height=0.75, |
---|
| 680 | culvert_routine=boyd_generalised_culvert_model, |
---|
| 681 | number_of_barrels=1, |
---|
[7496] | 682 | verbose=False) |
---|
[6102] | 683 | |
---|
[6123] | 684 | |
---|
[6102] | 685 | domain.forcing_terms.append(culvert) |
---|
| 686 | |
---|
| 687 | # Call |
---|
| 688 | culvert(domain) |
---|
| 689 | |
---|
[7496] | 690 | |
---|
| 691 | |
---|
| 692 | |
---|
| 693 | def test_momentum_jet(self): |
---|
| 694 | """test_momentum_jet |
---|
| 695 | |
---|
| 696 | Test that culvert_class can accept keyword use_momentum_jet |
---|
| 697 | This does not yet imply that the values have been tested. FIXME |
---|
| 698 | """ |
---|
| 699 | |
---|
| 700 | |
---|
| 701 | length = 40. |
---|
| 702 | width = 5. |
---|
| 703 | |
---|
| 704 | dx = dy = 1 # Resolution: Length of subdivisions on both axes |
---|
| 705 | |
---|
| 706 | points, vertices, boundary = rectangular_cross(int(length/dx), |
---|
| 707 | int(width/dy), |
---|
| 708 | len1=length, |
---|
| 709 | len2=width) |
---|
| 710 | domain = Domain(points, vertices, boundary) |
---|
| 711 | domain.set_name('Test_culvert_shallow') # Output name |
---|
| 712 | domain.set_default_order(2) |
---|
| 713 | |
---|
| 714 | |
---|
| 715 | #---------------------------------------------------------------------- |
---|
| 716 | # Setup initial conditions |
---|
| 717 | #---------------------------------------------------------------------- |
---|
| 718 | |
---|
| 719 | def topography(x, y): |
---|
| 720 | """Set up a weir |
---|
| 721 | |
---|
| 722 | A culvert will connect either side |
---|
| 723 | """ |
---|
| 724 | # General Slope of Topography |
---|
| 725 | z=-x/1000 |
---|
| 726 | |
---|
| 727 | N = len(x) |
---|
| 728 | for i in range(N): |
---|
| 729 | |
---|
| 730 | # Sloping Embankment Across Channel |
---|
| 731 | if 5.0 < x[i] < 10.1: |
---|
| 732 | # Cut Out Segment for Culvert face |
---|
| 733 | if 1.0+(x[i]-5.0)/5.0 < y[i] < 4.0 - (x[i]-5.0)/5.0: |
---|
| 734 | z[i]=z[i] |
---|
| 735 | else: |
---|
| 736 | z[i] += 0.5*(x[i] -5.0) # Sloping Segment U/S Face |
---|
| 737 | if 10.0 < x[i] < 12.1: |
---|
| 738 | z[i] += 2.5 # Flat Crest of Embankment |
---|
| 739 | if 12.0 < x[i] < 14.5: |
---|
| 740 | # Cut Out Segment for Culvert face |
---|
| 741 | if 2.0-(x[i]-12.0)/2.5 < y[i] < 3.0 + (x[i]-12.0)/2.5: |
---|
| 742 | z[i]=z[i] |
---|
| 743 | else: |
---|
| 744 | z[i] += 2.5-1.0*(x[i] -12.0) # Sloping D/S Face |
---|
| 745 | |
---|
| 746 | |
---|
| 747 | return z |
---|
| 748 | |
---|
| 749 | |
---|
| 750 | domain.set_quantity('elevation', topography) |
---|
| 751 | domain.set_quantity('friction', 0.01) # Constant friction |
---|
| 752 | domain.set_quantity('stage', |
---|
| 753 | expression='elevation + 1.0') # Shallow initial condition |
---|
| 754 | |
---|
| 755 | # Boyd culvert |
---|
| 756 | culvert = Culvert_flow(domain, |
---|
| 757 | label='Culvert No. 1', |
---|
| 758 | description='This culvert is a test unit 1.2m Wide by 0.75m High', |
---|
| 759 | end_point0=[9.0, 2.5], |
---|
| 760 | end_point1=[13.0, 2.5], |
---|
| 761 | width=1.20, height=0.75, |
---|
| 762 | culvert_routine=boyd_generalised_culvert_model, |
---|
| 763 | number_of_barrels=1, |
---|
| 764 | use_momentum_jet=True, |
---|
| 765 | update_interval=2, |
---|
| 766 | verbose=False) |
---|
| 767 | |
---|
| 768 | |
---|
| 769 | domain.forcing_terms.append(culvert) |
---|
| 770 | |
---|
| 771 | |
---|
| 772 | # Call |
---|
| 773 | culvert(domain) |
---|
| 774 | |
---|
| 775 | |
---|
| 776 | #----------------------------------------------------------------------- |
---|
| 777 | # Setup boundary conditions |
---|
| 778 | #----------------------------------------------------------------------- |
---|
| 779 | |
---|
| 780 | |
---|
| 781 | Br = Reflective_boundary(domain) # Solid reflective wall |
---|
| 782 | domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br}) |
---|
| 783 | |
---|
| 784 | #----------------------------------------------------------------------- |
---|
| 785 | # Evolve system through time |
---|
| 786 | #----------------------------------------------------------------------- |
---|
| 787 | |
---|
| 788 | ref_volume = domain.get_quantity('stage').get_integral() |
---|
| 789 | for t in domain.evolve(yieldstep = 0.1, finaltime = 25): |
---|
| 790 | new_volume = domain.get_quantity('stage').get_integral() |
---|
| 791 | |
---|
| 792 | msg = ('Total volume has changed: Is %.8f m^3 should have been %.8f m^3' |
---|
| 793 | % (new_volume, ref_volume)) |
---|
[7497] | 794 | assert num.allclose(new_volume, ref_volume, rtol=1.0e-10), msg |
---|
[7496] | 795 | |
---|
| 796 | |
---|
| 797 | |
---|
| 798 | |
---|
[6076] | 799 | |
---|
| 800 | #------------------------------------------------------------- |
---|
[7317] | 801 | |
---|
[6076] | 802 | if __name__ == "__main__": |
---|
[6123] | 803 | suite = unittest.makeSuite(Test_Culvert, 'test') |
---|
[7276] | 804 | runner = unittest.TextTestRunner() #verbosity=2) |
---|
[6076] | 805 | runner.run(suite) |
---|
| 806 | |
---|