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