Changeset 6055
- Timestamp:
- Dec 10, 2008, 5:07:51 PM (16 years ago)
- Location:
- anuga_core
- Files:
-
- 5 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/documentation/user_manual/version.tex
r4788 r6055 7 7 % release version; this is used to define the 8 8 % \version macro 9 \release{1.0beta\_ 0000}9 \release{1.0beta\_6051} -
anuga_core/source/anuga/culvert_flows/Test_Culvert_Flat_Water_Lev.py
r6051 r6055 61 61 z=-x/1000 62 62 63 # Changing Slopes in the X- Direction64 # N = len(x)65 # for i in range(N):66 # if 0 <x[i] < 5.1:67 # z[i] -= -x[i]/1068 # if 5 <x[i] < 10.1:69 # z[i] -= -x[i]/10070 # if 10 <x[i]:71 # z[i] -= -x[i]/2072 # return z73 74 63 # NOW Add bits and Pieces to topography 75 64 N = len(x) … … 91 80 92 81 93 # Constriction94 #if 27 < x[i] < 29 and y[i] > 3:95 # z[i] += 296 97 # Pole98 #if (x[i] - 34)**2 + (y[i] - 2)**2 < 0.4**2:99 # z[i] += 2100 101 # HOLE For Pit at Opening[0]102 #if (x[i]-4)**2 + (y[i]-1.5)**2<0.75**2:103 # z[i]-=1104 105 # HOLE For Pit at Opening[1]106 #if (x[i]-20)**2 + (y[i]-3.5)**2<0.5**2:107 # z[i]-=1108 82 109 83 return z … … 131 105 132 106 culvert = Culvert_flow(domain, 133 label='Culvert No. 1', 134 description='This culvert is a test unit 1.2m Wide by 0.75m High', 107 culvert_description_filename='example_rating_curve.csv', 135 108 end_point0=[9.0, 2.5], 136 109 end_point1=[13.0, 2.5], 137 width=1.20,height=0.75,138 culvert_routine=boyd_generalised_culvert_model,139 number_of_barrels=1,140 update_interval=2,141 110 verbose=True) 111 112 113 #culvert = 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) 142 123 143 124 domain.forcing_terms.append(culvert) -
anuga_core/source/anuga/culvert_flows/culvert_class.py
r6054 r6055 6 6 from anuga.utilities.polygon import plot_polygons 7 7 8 8 from anuga.utilities.numerical_tools import ensure_numeric 9 from Numeric import allclose 10 11 import sys 12 13 class Below_interval(Exception): pass 14 class Above_interval(Exception): pass 15 16 17 18 # FIXME(Ole): Write in C and reuse this function by similar code in interpolate.py 19 def interpolate_linearly(x, xvec, yvec): 20 21 # Find appropriate slot 22 i = 0 23 while x > xvec[i]: i += 1 24 25 if i == 0: 26 msg = 'Value provided = %.2f, interpolation minimum = %.2f.' %(x, xvec[0]) 27 raise Below_interval, msg 28 29 if i == len(xvec): 30 msg = 'Value provided = %.2f, interpolation maximum = %.2f.' %(x, xvec[-1]) 31 raise Above_interval, msg 32 33 34 x0 = xvec[i-1] 35 x1 = xvec[i] 36 alpha = (x - x0)/(x1 - x0) 37 38 y0 = yvec[i-1] 39 y1 = yvec[i] 40 y = alpha*y1 + (1-alpha)*y0 41 42 return y 43 44 45 46 def read_culvert_description(culvert_description_filename): 47 48 # Read description file 49 fid = open(culvert_description_filename) 50 51 read_rating_curve_data = False 52 rating_curve = [] 53 for i, line in enumerate(fid.readlines()): 54 55 if read_rating_curve_data is True: 56 fields = line.split(',') 57 head_difference = float(fields[0].strip()) 58 flow_rate = float(fields[1].strip()) 59 barrel_velocity = float(fields[2].strip()) 60 61 rating_curve.append( [head_difference, flow_rate, barrel_velocity] ) 62 63 if i == 0: 64 # Header 65 continue 66 if i == 1: 67 # Metadata 68 fields = line.split(',') 69 label=fields[0].strip() 70 type=fields[1].strip().lower() 71 assert type in ['box', 'pipe'] 72 73 width=float(fields[2].strip()) 74 height=float(fields[3].strip()) 75 length=float(fields[4].strip()) 76 number_of_barrels=int(fields[5].strip()) 77 #fields[6] refers to losses 78 description=fields[7].strip() 79 80 if line.strip() == '': continue # Skip blanks 81 82 if line.startswith('Rating'): 83 read_rating_curve_data = True 84 # Flow data follows 85 86 fid.close() 87 88 return label, type, width, height, length, number_of_barrels, description, rating_curve 89 9 90 10 91 class Culvert_flow: … … 21 102 def __init__(self, 22 103 domain, 23 104 culvert_description_filename=None, 105 end_point0=None, 106 end_point1=None, 107 update_interval=None, 108 log_file=False, 109 discharge_hydrograph=False, 110 verbose=False): 111 112 from Numeric import sqrt, sum 113 114 115 label, type, width, height, length, number_of_barrels, description, rating_curve = read_culvert_description(culvert_description_filename) 116 117 118 # Store culvert information 119 self.label = label 120 self.description = description 121 self.culvert_type = type 122 self.rating_curve = ensure_numeric(rating_curve) 123 self.number_of_barrels = number_of_barrels 124 125 if label is None: label = 'culvert_flow' 126 label += '_' + str(id(self)) 127 self.label = label 128 129 # File for storing discharge_hydrograph 130 if discharge_hydrograph is True: 131 self.timeseries_filename = label + '_timeseries.csv' 132 fid = open(self.timeseries_filename, 'w') 133 fid.write('time, discharge\n') 134 fid.close() 135 136 # Log file for storing general textual output 137 if log_file is True: 138 self.log_filename = label + '.log' 139 log_to_file(self.log_filename, self.label) 140 log_to_file(self.log_filename, description) 141 log_to_file(self.log_filename, self.culvert_type) 142 143 144 # Create the fundamental culvert polygons from POLYGON 145 #if self.culvert_type == 'circle': 146 # # Redefine width and height for use with create_culvert_polygons 147 # width = height = diameter 148 149 P = create_culvert_polygons(end_point0, 150 end_point1, 151 width=width, 152 height=height, 153 number_of_barrels=number_of_barrels) 154 155 if verbose is True: 156 pass 157 #plot_polygons([[end_point0, end_point1], 158 # P['exchange_polygon0'], 159 # P['exchange_polygon1'], 160 # P['enquiry_polygon0'], 161 # P['enquiry_polygon1']], 162 # figname='culvert_polygon_output') 163 164 165 # Compute the average point for enquiry 166 enquiry_point0 = sum(P['enquiry_polygon0'][:2])/2 167 enquiry_point1 = sum(P['enquiry_polygon1'][:2])/2 168 169 self.enquiry_points = [enquiry_point0, enquiry_point1] 170 self.enquiry_indices = [] 171 for point in self.enquiry_points: 172 # Find nearest centroid 173 N = len(domain) 174 points = domain.get_centroid_coordinates(absolute=True) 175 176 # Calculate indices in exchange area for this forcing term 177 178 triangle_id = min_dist = sys.maxint 179 for k in range(N): 180 x, y = points[k,:] # Centroid 181 182 c = point 183 distance = (x-c[0])**2+(y-c[1])**2 184 if distance < min_dist: 185 min_dist = distance 186 triangle_id = k 187 188 189 if triangle_id < sys.maxint: 190 msg = 'found triangle with centroid (%f, %f)'\ 191 %tuple(points[triangle_id, :]) 192 msg += ' for point (%f, %f)' %tuple(point) 193 194 self.enquiry_indices.append(triangle_id) 195 else: 196 msg = 'Triangle not found for point (%f, %f)' %point 197 raise Exception, msg 198 199 200 201 # Check that all polygons lie within the mesh. 202 bounding_polygon = domain.get_boundary_polygon() 203 for key in P.keys(): 204 if key in ['exchange_polygon0', 205 'exchange_polygon1']: 206 for point in list(P[key]) + self.enquiry_points: 207 msg = 'Point %s in polygon %s for culvert %s did not'\ 208 %(str(point), key, self.label) 209 msg += 'fall within the domain boundary.' 210 assert is_inside_polygon(point, bounding_polygon), msg 211 212 213 # Create inflow object at each end of the culvert. 214 self.openings = [] 215 self.openings.append(Inflow(domain, 216 polygon=P['exchange_polygon0'])) 217 218 self.openings.append(Inflow(domain, 219 polygon=P['exchange_polygon1'])) 220 221 222 223 dq = domain.quantities 224 for i, opening in enumerate(self.openings): 225 #elevation = dq['elevation'].get_values(location='centroids', 226 # interpolation_points=[self.enquiry_points[i]]) 227 228 elevation = dq['elevation'].get_values(location='centroids', 229 indices=[self.enquiry_indices[i]]) 230 opening.elevation = elevation 231 232 # Assume two openings for now: Referred to as 0 and 1 233 assert len(self.openings) == 2 234 235 # Determine pipe direction 236 self.delta_z = delta_z = self.openings[0].elevation - self.openings[1].elevation 237 if delta_z > 0.0: 238 self.inlet = self.openings[0] 239 self.outlet = self.openings[1] 240 else: 241 self.outlet = self.openings[0] 242 self.inlet = self.openings[1] 243 244 245 # Store basic geometry 246 self.end_points = [end_point0, end_point1] 247 self.vector = P['vector'] 248 self.length = P['length']; assert self.length > 0.0 249 if not allclose(self.length, length, rtol=1.0e-2, atol=1.0e-2): 250 msg = 'WARNING: barrel length specified in "%s" (%.2f m)' %(culvert_description_filename, length) 251 msg += ' does not match distance between specified' 252 msg += ' end points (%.2f m)' %self.length 253 print msg 254 255 self.verbose = verbose 256 self.last_update = 0.0 # For use with update_interval 257 self.update_interval = update_interval 258 259 260 # Print something 261 if hasattr(self, 'log_filename'): 262 s = 'Culvert Effective Length = %.2f m' %(self.length) 263 log_to_file(self.log_filename, s) 264 265 s = 'Culvert Direction is %s\n' %str(self.vector) 266 log_to_file(self.log_filename, s) 267 268 269 270 271 272 def __call__(self, domain): 273 from anuga.utilities.numerical_tools import mean 274 275 from anuga.config import g, epsilon 276 from Numeric import take, sqrt 277 from anuga.config import velocity_protection 278 279 280 # Time stuff 281 time = domain.get_time() 282 283 284 update = False 285 if self.update_interval is None: 286 update = True 287 else: 288 if time - self.last_update > self.update_interval or time == 0.0: 289 update = True 290 291 292 if update is True: 293 self.last_update = time 294 dq = domain.quantities 295 296 # Get average water depths at each opening 297 openings = self.openings # There are two Opening [0] and [1] 298 for i, opening in enumerate(openings): 299 300 # Compute mean values of selected quantitites in the 301 # enquiry area in front of the culvert 302 # Stage and velocity comes from enquiry area 303 # and elevation from exchange area 304 305 stage = dq['stage'].get_values(location='centroids', 306 indices=[self.enquiry_indices[i]]) 307 308 309 # Store current average stage and depth with each opening object 310 opening.stage = stage 311 312 313 314 ################# End of the FOR loop ################################################ 315 316 # We now need to deal with each opening individually 317 318 # Determine flow direction based on total energy difference 319 320 delta_w = self.inlet.stage - self.outlet.stage 321 322 if hasattr(self, 'log_filename'): 323 s = 'Time=%.2f, inlet stage = %.2f, outlet stage = %.2f' %(time, 324 self.inlet.stage, 325 self.outlet.stage) 326 log_to_file(self.log_filename, s) 327 328 329 # Calculate discharge for one barrel and set inlet.rate and outlet.rate 330 try: 331 Q = interpolate_linearly(delta_w, self.rating_curve[:,0], self.rating_curve[:,1]) 332 except Below_interval, e: 333 Q = self.rating_curve[0,1] 334 msg = '%.2fs: Delta head smaller than rating curve minimum: ' %time 335 msg += str(e) 336 msg += '\n I will use minimum discharge %.2f m^3/s for culvert "%s"'\ 337 %(Q, self.label) 338 if hasattr(self, 'log_filename'): 339 log_to_file(self.log_filename, msg) 340 except Above_interval, e: 341 Q = self.rating_curve[-1,1] 342 msg = '%.2fs: Delta head greater than rating curve maximum: ' %time 343 msg += str(e) 344 msg += '\n I will use maximum discharge %.2f m^3/s for culvert "%s"'\ 345 %(Q, self.label) 346 if hasattr(self, 'log_filename'): 347 log_to_file(self.log_filename, msg) 348 349 350 351 352 # Adjust discharge for multiple barrels 353 Q *= self.number_of_barrels 354 355 356 self.inlet.rate = -Q 357 self.outlet.rate = Q 358 359 # Log timeseries to file 360 try: 361 fid = open(self.timeseries_filename, 'a') 362 except: 363 pass 364 else: 365 fid.write('%.2f, %.2f\n' %(time, Q)) 366 fid.close() 367 368 369 370 371 372 # Execute flow term for each opening 373 # This is where Inflow objects are evaluated using the last rate that has been calculated 374 # 375 # This will take place at every internal timestep and update the domain 376 for opening in self.openings: 377 opening(domain) 378 379 380 381 382 383 384 class Culvert_flow_energy: 385 """Culvert flow - transfer water from one hole to another 386 387 Using Momentum as Calculated by Culvert Flow !! 388 Could be Several Methods Investigated to do This !!! 389 390 2008_May_08 391 To Ole: 392 OK so here we need to get the Polygon Creating code to create 393 polygons for the culvert Based on 394 the two input Points (X0,Y0) and (X1,Y1) - need to be passed 395 to create polygon 396 397 The two centers are now passed on to create_polygon. 398 399 400 Input: Two points, pipe_size (either diameter or width, height), 401 mannings_rougness, 402 inlet/outlet energy_loss_coefficients, internal_bend_coefficent, 403 top-down_blockage_factor and bottom_up_blockage_factor 404 405 406 And the Delta H enquiry should be change from Openings in line 412 407 to the enquiry Polygons infront of the culvert 408 At the moment this script uses only Depth, later we can change it to 409 Energy... 410 411 Once we have Delta H can calculate a Flow Rate and from Flow Rate 412 an Outlet Velocity 413 The Outlet Velocity x Outlet Depth = Momentum to be applied at the Outlet... 414 415 Invert levels are optional. If left out they will default to the 416 elevation at the opening. 417 418 """ 419 420 def __init__(self, 421 domain, 24 422 label=None, 25 423 description=None, 26 424 end_point0=None, 27 425 end_point1=None, 28 29 30 426 width=None, 427 height=None, 428 diameter=None, 31 429 manning=None, # Mannings Roughness for Culvert 32 430 invert_level0=None, # Invert level at opening 0 … … 428 826 429 827 430 431 432 433 434 class Culvert_flow_energy: 435 """Culvert flow - transfer water from one hole to another 436 437 Using Momentum as Calculated by Culvert Flow !! 438 Could be Several Methods Investigated to do This !!! 439 440 2008_May_08 441 To Ole: 442 OK so here we need to get the Polygon Creating code to create 443 polygons for the culvert Based on 444 the two input Points (X0,Y0) and (X1,Y1) - need to be passed 445 to create polygon 446 447 The two centers are now passed on to create_polygon. 448 449 450 Input: Two points, pipe_size (either diameter or width, height), 451 mannings_rougness, 452 inlet/outlet energy_loss_coefficients, internal_bend_coefficent, 453 top-down_blockage_factor and bottom_up_blockage_factor 454 455 456 And the Delta H enquiry should be change from Openings in line 412 457 to the enquiry Polygons infront of the culvert 458 At the moment this script uses only Depth, later we can change it to 459 Energy... 460 461 Once we have Delta H can calculate a Flow Rate and from Flow Rate 462 an Outlet Velocity 463 The Outlet Velocity x Outlet Depth = Momentum to be applied at the Outlet... 464 465 Invert levels are optional. If left out they will default to the 466 elevation at the opening. 467 468 """ 469 470 def __init__(self, 471 domain, 472 label=None, 473 description=None, 474 end_point0=None, 475 end_point1=None, 476 width=None, 477 height=None, 478 diameter=None, 479 manning=None, # Mannings Roughness for Culvert 480 invert_level0=None, # Invert level at opening 0 481 invert_level1=None, # Invert level at opening 1 482 loss_exit=None, 483 loss_entry=None, 484 loss_bend=None, 485 loss_special=None, 486 blockage_topdwn=None, 487 blockage_bottup=None, 488 culvert_routine=None, 489 number_of_barrels=1, 490 update_interval=None, 491 verbose=False): 492 493 from Numeric import sqrt, sum 494 495 # Input check 496 if diameter is not None: 497 self.culvert_type = 'circle' 498 self.diameter = diameter 499 if height is not None or width is not None: 500 msg = 'Either diameter or width&height must be specified, ' 501 msg += 'but not both.' 502 raise Exception, msg 503 else: 504 if height is not None: 505 if width is None: 506 self.culvert_type = 'square' 507 width = height 508 else: 509 self.culvert_type = 'rectangle' 510 elif width is not None: 511 if height is None: 512 self.culvert_type = 'square' 513 height = width 514 else: 515 msg = 'Either diameter or width&height must be specified.' 516 raise Exception, msg 517 518 if height == width: 519 self.culvert_type = 'square' 520 521 self.height = height 522 self.width = width 523 524 525 assert self.culvert_type in ['circle', 'square', 'rectangle'] 526 527 assert number_of_barrels >= 1 528 self.number_of_barrels = number_of_barrels 529 530 531 # Set defaults 532 if manning is None: manning = 0.012 # Default roughness for pipe 533 if loss_exit is None: loss_exit = 1.00 534 if loss_entry is None: loss_entry = 0.50 535 if loss_bend is None: loss_bend=0.00 536 if loss_special is None: loss_special=0.00 537 if blockage_topdwn is None: blockage_topdwn=0.00 538 if blockage_bottup is None: blockage_bottup=0.00 539 if culvert_routine is None: 540 culvert_routine=boyd_generalised_culvert_model 541 542 if label is None: label = 'culvert_flow' 543 label += '_' + str(id(self)) 544 self.label = label 545 546 # File for storing culvert quantities 547 self.timeseries_filename = label + '_timeseries.csv' 548 fid = open(self.timeseries_filename, 'w') 549 fid.write('time, E0, E1, Velocity, Discharge\n') 550 fid.close() 551 552 # Log file for storing general textual output 553 self.log_filename = label + '.log' 554 log_to_file(self.log_filename, self.label) 555 log_to_file(self.log_filename, description) 556 log_to_file(self.log_filename, self.culvert_type) 557 558 559 # Create the fundamental culvert polygons from POLYGON 560 if self.culvert_type == 'circle': 561 # Redefine width and height for use with create_culvert_polygons 562 width = height = diameter 563 564 P = create_culvert_polygons(end_point0, 565 end_point1, 566 width=width, 567 height=height, 568 number_of_barrels=number_of_barrels) 569 570 if verbose is True: 571 pass 572 #plot_polygons([[end_point0, end_point1], 573 # P['exchange_polygon0'], 574 # P['exchange_polygon1'], 575 # P['enquiry_polygon0'], 576 # P['enquiry_polygon1']], 577 # figname='culvert_polygon_output') 578 #import sys; sys.exit() 579 580 581 # Check that all polygons lie within the mesh. 582 bounding_polygon = domain.get_boundary_polygon() 583 for key in P.keys(): 584 if key in ['exchange_polygon0', 585 'exchange_polygon1', 586 'enquiry_polygon0', 587 'enquiry_polygon1']: 588 for point in P[key]: 589 590 msg = 'Point %s in polygon %s for culvert %s did not'\ 591 %(str(point), key, self.label) 592 msg += 'fall within the domain boundary.' 593 assert is_inside_polygon(point, bounding_polygon), msg 594 595 596 # Create inflow object at each end of the culvert. 597 self.openings = [] 598 self.openings.append(Inflow(domain, 599 polygon=P['exchange_polygon0'])) 600 601 self.openings.append(Inflow(domain, 602 polygon=P['exchange_polygon1'])) 603 604 605 # Assume two openings for now: Referred to as 0 and 1 606 assert len(self.openings) == 2 607 608 # Store basic geometry 609 self.end_points = [end_point0, end_point1] 610 self.invert_levels = [invert_level0, invert_level1] 611 #self.enquiry_polygons = [P['enquiry_polygon0'], P['enquiry_polygon1']] 612 self.enquiry_polylines = [P['enquiry_polygon0'][:2], 613 P['enquiry_polygon1'][:2]] 614 self.vector = P['vector'] 615 self.length = P['length']; assert self.length > 0.0 616 self.verbose = verbose 617 self.last_time = 0.0 618 self.last_update = 0.0 # For use with update_interval 619 self.update_interval = update_interval 620 621 622 # Store hydraulic parameters 623 self.manning = manning 624 self.loss_exit = loss_exit 625 self.loss_entry = loss_entry 626 self.loss_bend = loss_bend 627 self.loss_special = loss_special 628 self.sum_loss = loss_exit + loss_entry + loss_bend + loss_special 629 self.blockage_topdwn = blockage_topdwn 630 self.blockage_bottup = blockage_bottup 631 632 # Store culvert routine 633 self.culvert_routine = culvert_routine 634 635 636 # Create objects to update momentum (a bit crude at this stage) 637 638 639 xmom0 = General_forcing(domain, 'xmomentum', 640 polygon=P['exchange_polygon0']) 641 642 xmom1 = General_forcing(domain, 'xmomentum', 643 polygon=P['exchange_polygon1']) 644 645 ymom0 = General_forcing(domain, 'ymomentum', 646 polygon=P['exchange_polygon0']) 647 648 ymom1 = General_forcing(domain, 'ymomentum', 649 polygon=P['exchange_polygon1']) 650 651 self.opening_momentum = [ [xmom0, ymom0], [xmom1, ymom1] ] 652 653 654 # Print something 655 s = 'Culvert Effective Length = %.2f m' %(self.length) 656 log_to_file(self.log_filename, s) 657 658 s = 'Culvert Direction is %s\n' %str(self.vector) 659 log_to_file(self.log_filename, s) 660 661 662 def __call__(self, domain): 663 from anuga.utilities.numerical_tools import mean 664 665 from anuga.config import g, epsilon 666 from Numeric import take, sqrt 667 from anuga.config import velocity_protection 668 669 670 log_filename = self.log_filename 671 672 # Time stuff 673 time = domain.get_time() 674 675 676 update = False 677 if self.update_interval is None: 678 update = True 679 else: 680 if time - self.last_update > self.update_interval or time == 0.0: 681 update = True 682 683 #print 'call', time, time - self.last_update 684 685 686 if update is True: 687 #print 'Updating', time, time - self.last_update 688 self.last_update = time 689 690 delta_t = time-self.last_time 691 s = '\nTime = %.2f, delta_t = %f' %(time, delta_t) 692 log_to_file(log_filename, s) 693 694 msg = 'Time did not advance' 695 if time > 0.0: assert delta_t > 0.0, msg 696 697 698 # Get average water depths at each opening 699 openings = self.openings # There are two Opening [0] and [1] 700 for i, opening in enumerate(openings): 701 dq = domain.quantities 702 703 # Compute mean values of selected quantitites in the 704 # exchange area in front of the culvert 705 # Stage and velocity comes from enquiry area 706 # and elevation from exchange area 707 708 stage = dq['stage'].get_values(location='centroids', 709 indices=opening.exchange_indices) 710 w = mean(stage) # Average stage 711 712 # Use invert level instead of elevation if specified 713 invert_level = self.invert_levels[i] 714 if invert_level is not None: 715 z = invert_level 716 else: 717 elevation = dq['elevation'].get_values(location='centroids', 718 indices=opening.exchange_indices) 719 z = mean(elevation) # Average elevation 720 721 # Estimated depth above the culvert inlet 722 d = w - z # Used for calculations involving depth 723 if d < 0.0: 724 # This is possible since w and z are taken at different locations 725 #msg = 'D < 0.0: %f' %d 726 #raise msg 727 d = 0.0 728 729 730 # Ratio of depth to culvert height. 731 # If ratio > 1 then culvert is running full 732 if self.culvert_type == 'circle': 733 ratio = d/self.diameter 734 else: 735 ratio = d/self.height 736 opening.ratio = ratio 737 738 739 # Average measures of energy in front of this opening 740 polyline = self.enquiry_polylines[i] 741 #print 't = %.4f, opening=%d,' %(domain.time, i), 742 opening.total_energy = domain.get_energy_through_cross_section(polyline, 743 kind='total') 744 #print 'Et = %.3f m' %opening.total_energy 745 746 # Store current average stage and depth with each opening object 747 opening.depth = d 748 opening.depth_trigger = d # Use this for now 749 opening.stage = w 750 opening.elevation = z 751 752 753 ################# End of the FOR loop ################################################ 754 755 # We now need to deal with each opening individually 756 757 # Determine flow direction based on total energy difference 758 delta_Et = openings[0].total_energy - openings[1].total_energy 759 760 if delta_Et > 0: 761 #print 'Flow U/S ---> D/S' 762 inlet = openings[0] 763 outlet = openings[1] 764 765 inlet.momentum = self.opening_momentum[0] 766 outlet.momentum = self.opening_momentum[1] 767 768 else: 769 #print 'Flow D/S ---> U/S' 770 inlet = openings[1] 771 outlet = openings[0] 772 773 inlet.momentum = self.opening_momentum[1] 774 outlet.momentum = self.opening_momentum[0] 775 776 delta_Et = -delta_Et 777 778 self.inlet = inlet 779 self.outlet = outlet 780 781 msg = 'Total energy difference is negative' 782 assert delta_Et > 0.0, msg 783 784 delta_h = inlet.stage - outlet.stage 785 delta_z = inlet.elevation - outlet.elevation 786 culvert_slope = (delta_z/self.length) 787 788 if culvert_slope < 0.0: 789 # Adverse gradient - flow is running uphill 790 # Flow will be purely controlled by uphill outlet face 791 if self.verbose is True: 792 print 'WARNING: Flow is running uphill. Watch Out!', inlet.elevation, outlet.elevation 793 794 795 s = 'Delta total energy = %.3f' %(delta_Et) 796 log_to_file(log_filename, s) 797 798 799 # Calculate discharge for one barrel and set inlet.rate and outlet.rate 800 Q, barrel_velocity, culvert_outlet_depth = self.culvert_routine(self, inlet, outlet, delta_Et, g) 801 802 # Adjust discharge for multiple barrels 803 Q *= self.number_of_barrels 804 805 # Compute barrel momentum 806 barrel_momentum = barrel_velocity*culvert_outlet_depth 807 808 s = 'Barrel velocity = %f' %barrel_velocity 809 log_to_file(log_filename, s) 810 811 # Compute momentum vector at outlet 812 outlet_mom_x, outlet_mom_y = self.vector * barrel_momentum 813 814 s = 'Directional momentum = (%f, %f)' %(outlet_mom_x, outlet_mom_y) 815 log_to_file(log_filename, s) 816 817 # Log timeseries to file 818 fid = open(self.timeseries_filename, 'a') 819 fid.write('%f, %f, %f, %f, %f\n'\ 820 %(time, 821 openings[0].total_energy, 822 openings[1].total_energy, 823 barrel_velocity, 824 Q)) 825 fid.close() 826 827 # Update momentum 828 delta_t = time - self.last_time 829 if delta_t > 0.0: 830 xmomentum_rate = outlet_mom_x - outlet.momentum[0].value 831 xmomentum_rate /= delta_t 832 833 ymomentum_rate = outlet_mom_y - outlet.momentum[1].value 834 ymomentum_rate /= delta_t 835 836 s = 'X Y MOM_RATE = (%f, %f) ' %(xmomentum_rate, ymomentum_rate) 837 log_to_file(log_filename, s) 838 else: 839 xmomentum_rate = ymomentum_rate = 0.0 840 841 842 # Set momentum rates for outlet jet 843 outlet.momentum[0].rate = xmomentum_rate 844 outlet.momentum[1].rate = ymomentum_rate 845 846 # Remember this value for next step (IMPORTANT) 847 outlet.momentum[0].value = outlet_mom_x 848 outlet.momentum[1].value = outlet_mom_y 849 850 if int(domain.time*100) % 100 == 0: 851 s = 'T=%.5f, Culvert Discharge = %.3f f'\ 852 %(time, Q) 853 s += ' Depth= %0.3f Momentum = (%0.3f, %0.3f)'\ 854 %(culvert_outlet_depth, outlet_mom_x,outlet_mom_y) 855 s += ' Momentum rate: (%.4f, %.4f)'\ 856 %(xmomentum_rate, ymomentum_rate) 857 s+='Outlet Vel= %.3f'\ 858 %(barrel_velocity) 859 log_to_file(log_filename, s) 860 861 862 863 864 # Execute flow term for each opening 865 # This is where Inflow objects are evaluated and update the domain 866 for opening in self.openings: 867 opening(domain) 868 869 # Execute momentum terms 870 # This is where Inflow objects are evaluated and update the domain 871 self.outlet.momentum[0](domain) 872 self.outlet.momentum[1](domain) 873 874 # Store value of time #FIXME(Ole): Maybe only every time we update 875 self.last_time = time 876 877 878 828 -
anuga_core/source/anuga/culvert_flows/example_rating_curve.csv
r6054 r6055 1 Name, type, width or diameter, height, length, losses, description2 Test Culvert, box, 3, 1.8, 5, 1, 50% blockedtest case1 Name, type, width or diameter, height, length, number of barrels, losses, description 2 Test Culvert, box, 3, 1.8, 4, 1, 1, 50% blocked single barrel test case 3 3 4 4 5 Rating Curve 6 Delta W (m), Q (m3/s), Velocity (m/s) 5 Rating Curve: Delta W (m), Q (m3/s), Velocity (m/s) 7 6 0,0,0 8 7 0.1,0.720243203,0.800270226 -
anuga_core/source/anuga/shallow_water/shallow_water_domain.py
r6045 r6055 1982 1982 #------------------------------------------------------------------------ 1983 1983 # This is the new element implemented by Ole and Rudy to allow direct 1984 # input of Inflowin mm/s1984 # input of Rainfall in mm/s 1985 1985 1986 1986 catchmentrainfall = Rainfall(rain=file_function('Q100_2hr_Rain.tms'))
Note: See TracChangeset
for help on using the changeset viewer.