- Timestamp:
- Jun 30, 2009, 2:07:41 PM (15 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/shallow_water/shallow_water_domain.py
r7090 r7276 1 """Finite-volume computations of the shallow water wave equation. 1 """ 2 Finite-volume computations of the shallow water wave equation. 2 3 3 4 Title: ANGUA shallow_water_domain - 2D triangular domains for finite-volume … … 53 54 Journal of Hydraulic Engineering, vol. 127, No. 7 July 1999 54 55 55 Hydrodynamic modelling of coastal inundation. 56 Hydrodynamic modelling of coastal inundation. 56 57 Nielsen, O., S. Roberts, D. Gray, A. McPherson and A. Hitchman 57 58 In Zerger, A. and Argent, R.M. (eds) MODSIM 2005 International Congress on … … 71 72 $Author$ 72 73 $Date$ 73 74 74 """ 75 75 … … 80 80 # $LastChangedBy$ 81 81 82 import Numeric as num 82 83 import numpy as num 83 84 84 85 from anuga.abstract_2d_finite_volumes.neighbour_mesh import segment_midpoints … … 96 97 from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\ 97 98 import AWI_boundary 98 from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\99 import AWI_boundary100 99 101 100 from anuga.pmesh.mesh_interface import create_mesh_from_regions … … 114 113 from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a 115 114 116 from anuga.fit_interpolate.interpolate import Modeltime_too_late, Modeltime_too_early 117 118 from anuga.utilities.polygon import inside_polygon, polygon_area, is_inside_polygon 115 from anuga.fit_interpolate.interpolate import Modeltime_too_late, \ 116 Modeltime_too_early 117 118 from anuga.utilities.polygon import inside_polygon, polygon_area, \ 119 is_inside_polygon 119 120 120 121 from types import IntType, FloatType 121 122 from warnings import warn 122 123 123 import math 124 125 # 124 125 ################################################################################ 126 126 # Shallow water domain 127 #--------------------- 127 ################################################################################ 128 129 ## 130 # @brief Class for a shallow water domain. 128 131 class Domain(Generic_Domain): 129 132 130 133 conserved_quantities = ['stage', 'xmomentum', 'ymomentum'] 131 134 other_quantities = ['elevation', 'friction'] 132 135 136 ## 137 # @brief Instantiate a shallow water domain. 138 # @param coordinates 139 # @param vertices 140 # @param boundary 141 # @param tagged_elements 142 # @param geo_reference 143 # @param use_inscribed_circle 144 # @param mesh_filename 145 # @param use_cache 146 # @param verbose 147 # @param full_send_dict 148 # @param ghost_recv_dict 149 # @param processor 150 # @param numproc 151 # @param number_of_full_nodes 152 # @param number_of_full_triangles 133 153 def __init__(self, 134 154 coordinates=None, … … 147 167 number_of_full_nodes=None, 148 168 number_of_full_triangles=None): 149 150 169 151 170 other_quantities = ['elevation', 'friction'] … … 167 186 numproc, 168 187 number_of_full_nodes=number_of_full_nodes, 169 number_of_full_triangles=number_of_full_triangles) 170 188 number_of_full_triangles=number_of_full_triangles) 171 189 172 190 self.set_minimum_allowed_height(minimum_allowed_height) 173 191 174 192 self.maximum_allowed_speed = maximum_allowed_speed 175 193 self.g = g 176 self.beta_w 177 self.beta_w_dry 178 self.beta_uh 194 self.beta_w = beta_w 195 self.beta_w_dry = beta_w_dry 196 self.beta_uh = beta_uh 179 197 self.beta_uh_dry = beta_uh_dry 180 self.beta_vh 198 self.beta_vh = beta_vh 181 199 self.beta_vh_dry = beta_vh_dry 182 200 self.alpha_balance = alpha_balance … … 193 211 self.set_store_vertices_uniquely(False) 194 212 self.minimum_storable_height = minimum_storable_height 195 self.quantities_to_be_stored = ['stage', 'xmomentum','ymomentum']213 self.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum'] 196 214 197 215 # Limiters … … 200 218 self.use_centroid_velocities = use_centroid_velocities 201 219 202 220 ## 221 # @brief 222 # @param beta 203 223 def set_all_limiters(self, beta): 204 """Shorthand to assign one constant value [0,1 [to all limiters.224 """Shorthand to assign one constant value [0,1] to all limiters. 205 225 0 Corresponds to first order, where as larger values make use of 206 the second order scheme. 207 """ 208 209 self.beta_w 210 self.beta_w_dry 226 the second order scheme. 227 """ 228 229 self.beta_w = beta 230 self.beta_w_dry = beta 211 231 self.quantities['stage'].beta = beta 212 213 self.beta_uh 232 233 self.beta_uh = beta 214 234 self.beta_uh_dry = beta 215 235 self.quantities['xmomentum'].beta = beta 216 217 self.beta_vh 236 237 self.beta_vh = beta 218 238 self.beta_vh_dry = beta 219 239 self.quantities['ymomentum'].beta = beta 220 221 222 240 241 ## 242 # @brief 243 # @param flag 244 # @param reduction 223 245 def set_store_vertices_uniquely(self, flag, reduction=None): 224 246 """Decide whether vertex values should be stored uniquely as … … 235 257 #self.reduction = min #Looks better near steep slopes 236 258 237 259 ## 260 # @brief Set the minimum depth that will be written to an SWW file. 261 # @param minimum_storable_height The minimum stored height (in m). 238 262 def set_minimum_storable_height(self, minimum_storable_height): 239 """ 240 Set the minimum depth that will be recognised when writing 263 """Set the minimum depth that will be recognised when writing 241 264 to an sww file. This is useful for removing thin water layers 242 265 that seems to be caused by friction creep. … … 244 267 The minimum allowed sww depth is in meters. 245 268 """ 269 246 270 self.minimum_storable_height = minimum_storable_height 247 271 248 272 ## 273 # @brief 274 # @param minimum_allowed_height 249 275 def set_minimum_allowed_height(self, minimum_allowed_height): 250 """ 251 Set the minimum depth that will be recognised in the numerical 252 scheme 276 """Set minimum depth that will be recognised in the numerical scheme. 253 277 254 278 The minimum allowed depth is in meters. … … 263 287 #and deal with them adaptively similarly to how we used to use 1 order 264 288 #steps to recover. 289 265 290 self.minimum_allowed_height = minimum_allowed_height 266 self.H0 = minimum_allowed_height 267 268 291 self.H0 = minimum_allowed_height 292 293 ## 294 # @brief 295 # @param maximum_allowed_speed 269 296 def set_maximum_allowed_speed(self, maximum_allowed_speed): 270 """ 271 Set the maximum particle speed that is allowed in water 297 """Set the maximum particle speed that is allowed in water 272 298 shallower than minimum_allowed_height. This is useful for 273 299 controlling speeds in very thin layers of water and at the same time 274 300 allow some movement avoiding pooling of water. 275 276 """ 301 """ 302 277 303 self.maximum_allowed_speed = maximum_allowed_speed 278 304 279 280 def set_points_file_block_line_size(self,points_file_block_line_size): 281 """ 282 Set the minimum depth that will be recognised when writing 305 ## 306 # @brief 307 # @param points_file_block_line_size 308 def set_points_file_block_line_size(self, points_file_block_line_size): 309 """Set the minimum depth that will be recognised when writing 283 310 to an sww file. This is useful for removing thin water layers 284 311 that seems to be caused by friction creep. … … 287 314 """ 288 315 self.points_file_block_line_size = points_file_block_line_size 289 290 316 317 ## 318 # @brief Set the quantities that will be written to an SWW file. 319 # @param q The quantities to be written. 320 # @note Param 'q' may be None, single quantity or list of quantity strings. 321 # @note If 'q' is None, no quantities will be stored in the SWW file. 291 322 def set_quantities_to_be_stored(self, q): 292 323 """Specify which quantities will be stored in the sww file. … … 300 331 each yieldstep (This is in addition to the quantities elevation 301 332 and friction) 302 333 303 334 If q is None, storage will be switched off altogether. 304 335 """ 305 306 336 307 337 if q is None: … … 315 345 # Check correcness 316 346 for quantity_name in q: 317 msg = 'Quantity %s is not a valid conserved quantity'\ 318 %quantity_name 319 347 msg = ('Quantity %s is not a valid conserved quantity' 348 % quantity_name) 320 349 assert quantity_name in self.conserved_quantities, msg 321 350 322 351 self.quantities_to_be_stored = q 323 352 324 325 353 ## 354 # @brief 355 # @param indices 326 356 def get_wet_elements(self, indices=None): 327 357 """Return indices for elements where h > minimum_allowed_height … … 333 363 indices = get_wet_elements() 334 364 335 Note, centroid values are used for this operation 365 Note, centroid values are used for this operation 336 366 """ 337 367 … … 340 370 from anuga.config import minimum_allowed_height 341 371 342 343 372 elevation = self.get_quantity('elevation').\ 373 get_values(location='centroids', indices=indices) 374 stage = self.get_quantity('stage').\ 344 375 get_values(location='centroids', indices=indices) 345 stage = self.get_quantity('stage').\346 get_values(location='centroids', indices=indices)347 376 depth = stage - elevation 348 377 … … 350 379 wet_indices = num.compress(depth > minimum_allowed_height, 351 380 num.arange(len(depth))) 352 return wet_indices 353 354 381 return wet_indices 382 383 ## 384 # @brief 385 # @param indices 355 386 def get_maximum_inundation_elevation(self, indices=None): 356 387 """Return highest elevation where h > 0 … … 362 393 q = get_maximum_inundation_elevation() 363 394 364 Note, centroid values are used for this operation 395 Note, centroid values are used for this operation 365 396 """ 366 397 367 398 wet_elements = self.get_wet_elements(indices) 368 399 return self.get_quantity('elevation').\ 369 get_maximum_value(indices=wet_elements) 370 371 400 get_maximum_value(indices=wet_elements) 401 402 ## 403 # @brief 404 # @param indices 372 405 def get_maximum_inundation_location(self, indices=None): 373 406 """Return location of highest elevation where h > 0 … … 379 412 q = get_maximum_inundation_location() 380 413 381 Note, centroid values are used for this operation 414 Note, centroid values are used for this operation 382 415 """ 383 416 384 417 wet_elements = self.get_wet_elements(indices) 385 418 return self.get_quantity('elevation').\ 386 get_maximum_location(indices=wet_elements) 387 388 389 390 391 def get_flow_through_cross_section(self, polyline, 392 verbose=False): 393 """Get the total flow through an arbitrary poly line. 394 395 This is a run-time equivalent of the function with same name 419 get_maximum_location(indices=wet_elements) 420 421 ## 422 # @brief Get the total flow through an arbitrary poly line. 423 # @param polyline Representation of desired cross section. 424 # @param verbose True if this method is to be verbose. 425 # @note 'polyline' may contain multiple sections allowing complex shapes. 426 # @note Assume absolute UTM coordinates. 427 def get_flow_through_cross_section(self, polyline, verbose=False): 428 """Get the total flow through an arbitrary poly line. 429 430 This is a run-time equivalent of the function with same name 396 431 in data_manager.py 397 432 398 433 Input: 399 polyline: Representation of desired cross section - it may contain 400 multiple sections allowing for complex shapes. Assume 434 polyline: Representation of desired cross section - it may contain 435 multiple sections allowing for complex shapes. Assume 401 436 absolute UTM coordinates. 402 Format [[x0, y0], [x1, y1], ...] 403 404 Output: 437 Format [[x0, y0], [x1, y1], ...] 438 439 Output: 405 440 Q: Total flow [m^3/s] across given segments. 406 407 408 """ 409 441 """ 442 410 443 # Find all intersections and associated triangles. 411 segments = self.get_intersecting_segments(polyline, 412 use_cache=True, 444 segments = self.get_intersecting_segments(polyline, use_cache=True, 413 445 verbose=verbose) 414 446 415 447 # Get midpoints 416 midpoints = segment_midpoints(segments) 417 448 midpoints = segment_midpoints(segments) 449 418 450 # Make midpoints Geospatial instances 419 midpoints = ensure_geospatial(midpoints, self.geo_reference) 420 421 # Compute flow 451 midpoints = ensure_geospatial(midpoints, self.geo_reference) 452 453 # Compute flow 422 454 if verbose: print 'Computing flow through specified cross section' 423 455 424 456 # Get interpolated values 425 457 xmomentum = self.get_quantity('xmomentum') 426 ymomentum = self.get_quantity('ymomentum') 427 428 uh = xmomentum.get_values(interpolation_points=midpoints, use_cache=True) 429 vh = ymomentum.get_values(interpolation_points=midpoints, use_cache=True) 430 458 ymomentum = self.get_quantity('ymomentum') 459 460 uh = xmomentum.get_values(interpolation_points=midpoints, 461 use_cache=True) 462 vh = ymomentum.get_values(interpolation_points=midpoints, 463 use_cache=True) 464 431 465 # Compute and sum flows across each segment 432 total_flow =0466 total_flow = 0 433 467 for i in range(len(uh)): 434 435 # Inner product of momentum vector with segment normal [m^2/s] 468 # Inner product of momentum vector with segment normal [m^2/s] 436 469 normal = segments[i].normal 437 normal_momentum = uh[i]*normal[0] + vh[i]*normal[1] 438 470 normal_momentum = uh[i]*normal[0] + vh[i]*normal[1] 471 439 472 # Flow across this segment [m^3/s] 440 473 segment_flow = normal_momentum*segments[i].length … … 442 475 # Accumulate 443 476 total_flow += segment_flow 444 477 445 478 return total_flow 446 479 447 448 480 ## 481 # @brief 482 # @param polyline Representation of desired cross section. 483 # @param kind Select energy type to compute ('specific' or 'total'). 484 # @param verbose True if this method is to be verbose. 485 # @note 'polyline' may contain multiple sections allowing complex shapes. 486 # @note Assume absolute UTM coordinates. 449 487 def get_energy_through_cross_section(self, polyline, 450 488 kind='total', 451 verbose=False): 489 verbose=False): 452 490 """Obtain average energy head [m] across specified cross section. 453 491 454 492 Inputs: 455 polyline: Representation of desired cross section - it may contain 456 multiple sections allowing for complex shapes. Assume 493 polyline: Representation of desired cross section - it may contain 494 multiple sections allowing for complex shapes. Assume 457 495 absolute UTM coordinates. 458 496 Format [[x0, y0], [x1, y1], ...] 459 kind: Select which energy to compute. 497 kind: Select which energy to compute. 460 498 Options are 'specific' and 'total' (default) 461 499 … … 463 501 E: Average energy [m] across given segments for all stored times. 464 502 465 The average velocity is computed for each triangle intersected by 466 the polyline and averaged weighted by segment lengths. 467 468 The typical usage of this function would be to get average energy of 469 flow in a channel, and the polyline would then be a cross section 503 The average velocity is computed for each triangle intersected by 504 the polyline and averaged weighted by segment lengths. 505 506 The typical usage of this function would be to get average energy of 507 flow in a channel, and the polyline would then be a cross section 470 508 perpendicular to the flow. 471 509 472 #FIXME (Ole) - need name for this energy reflecting that its dimension 510 #FIXME (Ole) - need name for this energy reflecting that its dimension 473 511 is [m]. 474 512 """ 475 513 476 from anuga.config import g, epsilon, velocity_protection as h0 477 478 514 from anuga.config import g, epsilon, velocity_protection as h0 515 479 516 # Find all intersections and associated triangles. 480 segments = self.get_intersecting_segments(polyline, 481 use_cache=True, 517 segments = self.get_intersecting_segments(polyline, use_cache=True, 482 518 verbose=verbose) 483 519 484 520 # Get midpoints 485 midpoints = segment_midpoints(segments) 486 521 midpoints = segment_midpoints(segments) 522 487 523 # Make midpoints Geospatial instances 488 midpoints = ensure_geospatial(midpoints, self.geo_reference) 489 490 # Compute energy 491 if verbose: print 'Computing %s energy' %kind 492 524 midpoints = ensure_geospatial(midpoints, self.geo_reference) 525 526 # Compute energy 527 if verbose: print 'Computing %s energy' %kind 528 493 529 # Get interpolated values 494 stage = self.get_quantity('stage') 495 elevation = self.get_quantity('elevation') 530 stage = self.get_quantity('stage') 531 elevation = self.get_quantity('elevation') 496 532 xmomentum = self.get_quantity('xmomentum') 497 ymomentum = self.get_quantity('ymomentum') 533 ymomentum = self.get_quantity('ymomentum') 498 534 499 535 w = stage.get_values(interpolation_points=midpoints, use_cache=True) 500 z = elevation.get_values(interpolation_points=midpoints, use_cache=True) 501 uh = xmomentum.get_values(interpolation_points=midpoints, use_cache=True) 502 vh = ymomentum.get_values(interpolation_points=midpoints, use_cache=True) 503 h = w-z # Depth 504 536 z = elevation.get_values(interpolation_points=midpoints, use_cache=True) 537 uh = xmomentum.get_values(interpolation_points=midpoints, 538 use_cache=True) 539 vh = ymomentum.get_values(interpolation_points=midpoints, 540 use_cache=True) 541 h = w-z # Depth 542 505 543 # Compute total length of polyline for use with weighted averages 506 544 total_line_length = 0.0 507 545 for segment in segments: 508 546 total_line_length += segment.length 509 547 510 548 # Compute and sum flows across each segment 511 average_energy =0.0549 average_energy = 0.0 512 550 for i in range(len(w)): 513 514 551 # Average velocity across this segment 515 552 if h[i] > epsilon: … … 519 556 else: 520 557 u = v = 0.0 521 522 speed_squared = u*u + v*v 558 559 speed_squared = u*u + v*v 523 560 kinetic_energy = 0.5*speed_squared/g 524 561 525 562 if kind == 'specific': 526 563 segment_energy = h[i] + kinetic_energy 527 564 elif kind == 'total': 528 segment_energy = w[i] + kinetic_energy 565 segment_energy = w[i] + kinetic_energy 529 566 else: 530 567 msg = 'Energy kind must be either "specific" or "total".' 531 568 msg += ' I got %s' %kind 532 533 569 534 570 # Add to weighted average 535 571 weigth = segments[i].length/total_line_length 536 572 average_energy += segment_energy*weigth 537 538 573 539 574 return average_energy 540 541 542 543 575 576 ## 577 # @brief Run integrity checks on shallow water domain. 544 578 def check_integrity(self): 545 579 Generic_Domain.check_integrity(self) 546 580 547 581 #Check that we are solving the shallow water wave equation 548 549 582 msg = 'First conserved quantity must be "stage"' 550 583 assert self.conserved_quantities[0] == 'stage', msg … … 554 587 assert self.conserved_quantities[2] == 'ymomentum', msg 555 588 589 ## 590 # @brief 556 591 def extrapolate_second_order_sw(self): 557 #Call correct module function 558 #(either from this module or C-extension) 592 #Call correct module function (either from this module or C-extension) 559 593 extrapolate_second_order_sw(self) 560 594 595 ## 596 # @brief 561 597 def compute_fluxes(self): 562 #Call correct module function 563 #(either from this module or C-extension) 598 #Call correct module function (either from this module or C-extension) 564 599 compute_fluxes(self) 565 600 601 ## 602 # @brief 566 603 def distribute_to_vertices_and_edges(self): 567 604 # Call correct module function 568 605 if self.use_edge_limiter: 569 distribute_using_edge_limiter(self) 606 distribute_using_edge_limiter(self) 570 607 else: 571 608 distribute_using_vertex_limiter(self) 572 609 573 574 575 610 ## 611 # @brief Evolve the model by one step. 612 # @param yieldstep 613 # @param finaltime 614 # @param duration 615 # @param skip_initial_step 576 616 def evolve(self, 577 yieldstep = None, 578 finaltime = None, 579 duration = None, 580 skip_initial_step = False): 581 """Specialisation of basic evolve method from parent class 582 """ 617 yieldstep=None, 618 finaltime=None, 619 duration=None, 620 skip_initial_step=False): 621 """Specialisation of basic evolve method from parent class""" 583 622 584 623 # Call check integrity here rather than from user scripts 585 624 # self.check_integrity() 586 625 587 msg = ' Parameter beta_w must be in the interval [0, 2['626 msg = 'Attribute self.beta_w must be in the interval [0, 2]' 588 627 assert 0 <= self.beta_w <= 2.0, msg 589 628 590 591 629 # Initial update of vertex and edge values before any STORAGE 592 # and or visualisation 630 # and or visualisation. 593 631 # This is done again in the initialisation of the Generic_Domain 594 # evolve loop but we do it here to ensure the values are ok for storage 632 # evolve loop but we do it here to ensure the values are ok for storage. 595 633 self.distribute_to_vertices_and_edges() 596 634 597 635 if self.store is True and self.time == 0.0: 598 636 self.initialise_storage() 599 # print 'Storing results in ' + self.writer.filename600 637 else: 601 638 pass … … 606 643 607 644 # Call basic machinery from parent class 608 for t in Generic_Domain.evolve(self, 609 yieldstep=yieldstep, 610 finaltime=finaltime, 611 duration=duration, 645 for t in Generic_Domain.evolve(self, yieldstep=yieldstep, 646 finaltime=finaltime, duration=duration, 612 647 skip_initial_step=skip_initial_step): 613 614 648 # Store model data, e.g. for subsequent visualisation 615 649 if self.store is True: … … 622 656 yield(t) 623 657 624 658 ## 659 # @brief 625 660 def initialise_storage(self): 626 661 """Create and initialise self.writer object for storing data. … … 636 671 self.writer.store_connectivity() 637 672 638 673 ## 674 # @brief 675 # @param name 639 676 def store_timestep(self, name): 640 677 """Store named quantity and time. … … 643 680 self.write has been initialised 644 681 """ 682 645 683 self.writer.store_timestep(name) 646 684 647 685 ## 686 # @brief Get time stepping statistics string for printing. 687 # @param track_speeds 688 # @param triangle_id 648 689 def timestepping_statistics(self, 649 690 track_speeds=False, 650 triangle_id=None): 691 triangle_id=None): 651 692 """Return string with time stepping statistics for printing or logging 652 693 … … 656 697 """ 657 698 658 from anuga.config import epsilon, g 659 699 from anuga.config import epsilon, g 660 700 661 701 # Call basic machinery from parent class 662 msg = Generic_Domain.timestepping_statistics(self, 663 track_speeds, 702 msg = Generic_Domain.timestepping_statistics(self, track_speeds, 664 703 triangle_id) 665 704 666 705 if track_speeds is True: 667 668 706 # qwidth determines the text field used for quantities 669 707 qwidth = self.qwidth 670 708 671 709 # Selected triangle 672 710 k = self.k … … 674 712 # Report some derived quantities at vertices, edges and centroid 675 713 # specific to the shallow water wave equation 676 677 714 z = self.quantities['elevation'] 678 w = self.quantities['stage'] 715 w = self.quantities['stage'] 679 716 680 717 Vw = w.get_values(location='vertices', indices=[k])[0] … … 684 721 Vz = z.get_values(location='vertices', indices=[k])[0] 685 722 Ez = z.get_values(location='edges', indices=[k])[0] 686 Cz = z.get_values(location='centroids', indices=[k]) 687 723 Cz = z.get_values(location='centroids', indices=[k]) 688 724 689 725 name = 'depth' … … 691 727 Eh = Ew-Ez 692 728 Ch = Cw-Cz 693 729 694 730 s = ' %s: vertex_values = %.4f,\t %.4f,\t %.4f\n'\ 695 731 %(name.ljust(qwidth), Vh[0], Vh[1], Vh[2]) 696 732 697 733 s += ' %s: edge_values = %.4f,\t %.4f,\t %.4f\n'\ 698 734 %(name.ljust(qwidth), Eh[0], Eh[1], Eh[2]) 699 735 700 736 s += ' %s: centroid_value = %.4f\n'\ 701 %(name.ljust(qwidth), Ch[0]) 702 737 %(name.ljust(qwidth), Ch[0]) 738 703 739 msg += s 704 740 … … 709 745 Euh = uh.get_values(location='edges', indices=[k])[0] 710 746 Cuh = uh.get_values(location='centroids', indices=[k]) 711 747 712 748 Vvh = vh.get_values(location='vertices', indices=[k])[0] 713 749 Evh = vh.get_values(location='edges', indices=[k])[0] … … 717 753 Vu = Vuh/(Vh + epsilon) 718 754 Eu = Euh/(Eh + epsilon) 719 Cu = Cuh/(Ch + epsilon) 755 Cu = Cuh/(Ch + epsilon) 720 756 name = 'U' 721 757 s = ' %s: vertex_values = %.4f,\t %.4f,\t %.4f\n'\ 722 758 %(name.ljust(qwidth), Vu[0], Vu[1], Vu[2]) 723 759 724 760 s += ' %s: edge_values = %.4f,\t %.4f,\t %.4f\n'\ 725 761 %(name.ljust(qwidth), Eu[0], Eu[1], Eu[2]) 726 762 727 763 s += ' %s: centroid_value = %.4f\n'\ 728 %(name.ljust(qwidth), Cu[0]) 729 764 %(name.ljust(qwidth), Cu[0]) 765 730 766 msg += s 731 767 732 768 Vv = Vvh/(Vh + epsilon) 733 769 Ev = Evh/(Eh + epsilon) 734 Cv = Cvh/(Ch + epsilon) 770 Cv = Cvh/(Ch + epsilon) 735 771 name = 'V' 736 772 s = ' %s: vertex_values = %.4f,\t %.4f,\t %.4f\n'\ 737 773 %(name.ljust(qwidth), Vv[0], Vv[1], Vv[2]) 738 774 739 775 s += ' %s: edge_values = %.4f,\t %.4f,\t %.4f\n'\ 740 776 %(name.ljust(qwidth), Ev[0], Ev[1], Ev[2]) 741 777 742 778 s += ' %s: centroid_value = %.4f\n'\ 743 %(name.ljust(qwidth), Cv[0]) 744 779 %(name.ljust(qwidth), Cv[0]) 780 745 781 msg += s 746 747 782 748 783 # Froude number in each direction … … 751 786 Efx = Eu/(num.sqrt(g*Eh) + epsilon) 752 787 Cfx = Cu/(num.sqrt(g*Ch) + epsilon) 753 788 754 789 s = ' %s: vertex_values = %.4f,\t %.4f,\t %.4f\n'\ 755 790 %(name.ljust(qwidth), Vfx[0], Vfx[1], Vfx[2]) 756 791 757 792 s += ' %s: edge_values = %.4f,\t %.4f,\t %.4f\n'\ 758 793 %(name.ljust(qwidth), Efx[0], Efx[1], Efx[2]) 759 794 760 795 s += ' %s: centroid_value = %.4f\n'\ 761 %(name.ljust(qwidth), Cfx[0]) 762 796 %(name.ljust(qwidth), Cfx[0]) 797 763 798 msg += s 764 765 799 766 800 name = 'Froude (y)' … … 768 802 Efy = Ev/(num.sqrt(g*Eh) + epsilon) 769 803 Cfy = Cv/(num.sqrt(g*Ch) + epsilon) 770 804 771 805 s = ' %s: vertex_values = %.4f,\t %.4f,\t %.4f\n'\ 772 806 %(name.ljust(qwidth), Vfy[0], Vfy[1], Vfy[2]) 773 807 774 808 s += ' %s: edge_values = %.4f,\t %.4f,\t %.4f\n'\ 775 809 %(name.ljust(qwidth), Efy[0], Efy[1], Efy[2]) 776 810 777 811 s += ' %s: centroid_value = %.4f\n'\ 778 %(name.ljust(qwidth), Cfy[0]) 779 780 msg += s 781 782 812 %(name.ljust(qwidth), Cfy[0]) 813 814 msg += s 783 815 784 816 return msg 785 817 786 818 787 819 … … 880 912 """Create volumetric balance report suitable for printing or logging. 881 913 """ 882 883 X = self.compute_boundary_flows()884 boundary_flows, total_boundary_inflow, total_boundary_outflow = X914 915 (boundary_flows, total_boundary_inflow, 916 total_boundary_outflow) = self.compute_boundary_flows() 885 917 886 918 s = '---------------------------\n' 887 919 s += 'Volumetric balance report:\n' 888 920 s += '--------------------------\n' 889 s += 'Total boundary inflow [m^3/s]: %.2 e\n' % total_boundary_inflow890 s += 'Total boundary outflow [m^3/s]: %.2 e\n' % total_boundary_outflow921 s += 'Total boundary inflow [m^3/s]: %.2f\n' % total_boundary_inflow 922 s += 'Total boundary outflow [m^3/s]: %.2f\n' % total_boundary_outflow 891 923 s += 'Net boundary flow by tags [m^3/s]\n' 892 924 for tag in boundary_flows: 893 s += ' %s [m^3/s]: %.2 e\n' % (tag, boundary_flows[tag])894 895 s += 'Total net boundary flow [m^3/s]: %.2 e\n' % (total_boundary_inflow + total_boundary_outflow)896 s += 'Total volume in domain [m^3]: %.2 e\n' % self.compute_total_volume()925 s += ' %s [m^3/s]: %.2f\n' % (tag, boundary_flows[tag]) 926 927 s += 'Total net boundary flow [m^3/s]: %.2f\n' % (total_boundary_inflow + total_boundary_outflow) 928 s += 'Total volume in domain [m^3]: %.2f\n' % self.compute_total_volume() 897 929 898 930 # The go through explicit forcing update and record the rate of change for stage and … … 904 936 return s 905 937 906 907 908 909 910 #=============== End of class Shallow Water Domain =============================== 911 938 ################################################################################ 939 # End of class Shallow Water Domain 940 ################################################################################ 912 941 913 942 #----------------- … … 915 944 #----------------- 916 945 946 ## @brief Compute fluxes and timestep suitable for all volumes in domain. 947 # @param domain The domain to calculate fluxes for. 917 948 def compute_fluxes(domain): 918 """Compute all fluxes and the timestep suitable for all volumes 919 in domain. 949 """Compute fluxes and timestep suitable for all volumes in domain. 920 950 921 951 Compute total flux for each conserved quantity using "flux_function" … … 933 963 domain.explicit_update is reset to computed flux values 934 964 domain.timestep is set to the largest step satisfying all volumes. 935 936 965 937 966 This wrapper calls the underlying C version of compute fluxes … … 939 968 940 969 import sys 941 942 N = len(domain) # number_of_triangles 970 from shallow_water_ext import compute_fluxes_ext_central \ 971 as compute_fluxes_ext 972 973 N = len(domain) # number_of_triangles 943 974 944 975 # Shortcuts … … 949 980 950 981 timestep = float(sys.maxint) 951 from shallow_water_ext import\952 compute_fluxes_ext_central as compute_fluxes_ext953 954 982 955 983 flux_timestep = compute_fluxes_ext(timestep, … … 980 1008 domain.flux_timestep = flux_timestep 981 1009 982 983 984 #--------------------------------------- 1010 ################################################################################ 985 1011 # Module functions for gradient limiting 986 #--------------------------------------- 987 988 989 # MH090605 The following method belongs to the shallow_water domain class 990 # see comments in the corresponding method in shallow_water_ext.c 1012 ################################################################################ 1013 1014 ## 1015 # @brief Wrapper for C version of extrapolate_second_order_sw. 1016 # @param domain The domain to operate on. 1017 # @note MH090605 The following method belongs to the shallow_water domain class 1018 # see comments in the corresponding method in shallow_water_ext.c 991 1019 def extrapolate_second_order_sw(domain): 992 """Wrapper calling C version of extrapolate_second_order_sw 993 """ 1020 """Wrapper calling C version of extrapolate_second_order_sw""" 1021 994 1022 import sys 1023 from shallow_water_ext import extrapolate_second_order_sw as extrapol2 995 1024 996 1025 N = len(domain) # number_of_triangles … … 1002 1031 Elevation = domain.quantities['elevation'] 1003 1032 1004 from shallow_water_ext import extrapolate_second_order_sw as extrapol21005 1033 extrapol2(domain, 1006 1034 domain.surrogate_neighbours, … … 1018 1046 int(domain.optimise_dry_cells)) 1019 1047 1020 1048 ## 1049 # @brief Distribution from centroids to vertices specific to the SWW eqn. 1050 # @param domain The domain to operate on. 1021 1051 def distribute_using_vertex_limiter(domain): 1022 """Distribution from centroids to vertices specific to the 1023 shallow water wave 1024 equation. 1052 """Distribution from centroids to vertices specific to the SWW equation. 1025 1053 1026 1054 It will ensure that h (w-z) is always non-negative even in the … … 1039 1067 Postcondition 1040 1068 Conserved quantities defined at vertices 1041 1042 1069 """ 1043 1044 1045 1070 1046 1071 # Remove very thin layers of water … … 1073 1098 raise 'Unknown order' 1074 1099 1075 1076 1100 # Take bed elevation into account when water heights are small 1077 1101 balance_deep_and_shallow(domain) … … 1082 1106 Q.interpolate_from_vertices_to_edges() 1083 1107 1084 1085 1108 ## 1109 # @brief Distribution from centroids to edges specific to the SWW eqn. 1110 # @param domain The domain to operate on. 1086 1111 def distribute_using_edge_limiter(domain): 1087 """Distribution from centroids to edges specific to the 1088 shallow water wave 1089 equation. 1112 """Distribution from centroids to edges specific to the SWW eqn. 1090 1113 1091 1114 It will ensure that h (w-z) is always non-negative even in the … … 1103 1126 Postcondition 1104 1127 Conserved quantities defined at vertices 1105 1106 1128 """ 1107 1129 1108 1130 # Remove very thin layers of water 1109 1131 protect_against_infinitesimal_and_negative_heights(domain) 1110 1111 1132 1112 1133 for name in domain.conserved_quantities: … … 1126 1147 Q.interpolate_from_vertices_to_edges() 1127 1148 1128 1149 ## 1150 # @brief Protect against infinitesimal heights and associated high velocities. 1151 # @param domain The domain to operate on. 1129 1152 def protect_against_infinitesimal_and_negative_heights(domain): 1130 """Protect against infinitesimal heights and associated high velocities 1131 """ 1153 """Protect against infinitesimal heights and associated high velocities""" 1154 1155 from shallow_water_ext import protect 1132 1156 1133 1157 # Shortcuts … … 1137 1161 ymomc = domain.quantities['ymomentum'].centroid_values 1138 1162 1139 from shallow_water_ext import protect1140 1141 1163 protect(domain.minimum_allowed_height, domain.maximum_allowed_speed, 1142 1164 domain.epsilon, wc, zc, xmomc, ymomc) 1143 1165 1144 1145 1166 ## 1167 # @brief Wrapper for C function balance_deep_and_shallow_c(). 1168 # @param domain The domain to operate on. 1146 1169 def balance_deep_and_shallow(domain): 1147 1170 """Compute linear combination between stage as computed by … … 1158 1181 """ 1159 1182 1160 from shallow_water_ext import balance_deep_and_shallow as balance_deep_and_shallow_c1161 1183 from shallow_water_ext import balance_deep_and_shallow \ 1184 as balance_deep_and_shallow_c 1162 1185 1163 1186 # Shortcuts 1164 1187 wc = domain.quantities['stage'].centroid_values 1165 1188 zc = domain.quantities['elevation'].centroid_values 1166 1167 1189 wv = domain.quantities['stage'].vertex_values 1168 1190 zv = domain.quantities['elevation'].vertex_values … … 1177 1199 1178 1200 balance_deep_and_shallow_c(domain, 1179 wc, zc, wv, zv, wc, 1201 wc, zc, wv, zv, wc, 1180 1202 xmomc, ymomc, xmomv, ymomv) 1181 1203 1182 1204 1183 1184 1185 #------------------------------------------------------------------ 1205 ################################################################################ 1186 1206 # Boundary conditions - specific to the shallow water wave equation 1187 #------------------------------------------------------------------ 1207 ################################################################################ 1208 1209 ## 1210 # @brief Class for a reflective boundary. 1211 # @note Inherits from Boundary. 1188 1212 class Reflective_boundary(Boundary): 1189 1213 """Reflective boundary returns same conserved quantities as … … 1195 1219 """ 1196 1220 1197 def __init__(self, domain = None): 1221 ## 1222 # @brief Instantiate a Reflective_boundary. 1223 # @param domain 1224 def __init__(self, domain=None): 1198 1225 Boundary.__init__(self) 1199 1226 1200 1227 if domain is None: 1201 1228 msg = 'Domain must be specified for reflective boundary' 1202 raise msg1229 raise Exception, msg 1203 1230 1204 1231 # Handy shorthands 1205 self.stage 1206 self.xmom 1207 self.ymom 1232 self.stage = domain.quantities['stage'].edge_values 1233 self.xmom = domain.quantities['xmomentum'].edge_values 1234 self.ymom = domain.quantities['ymomentum'].edge_values 1208 1235 self.normals = domain.normals 1209 1236 1210 self.conserved_quantities = num.zeros(3, num.Float) 1211 1237 self.conserved_quantities = num.zeros(3, num.float) 1238 1239 ## 1240 # @brief Return a representation of this instance. 1212 1241 def __repr__(self): 1213 1242 return 'Reflective_boundary' 1214 1243 1215 1244 ## 1245 # @brief Calculate reflections (reverse outward momentum). 1246 # @param vol_id 1247 # @param edge_id 1216 1248 def evaluate(self, vol_id, edge_id): 1217 1249 """Reflective boundaries reverses the outward momentum … … 1226 1258 normal = self.normals[vol_id, 2*edge_id:2*edge_id+2] 1227 1259 1228 1229 1260 r = rotate(q, normal, direction = 1) 1230 1261 r[1] = -r[1] … … 1234 1265 1235 1266 1236 1237 1267 ## 1268 # @brief Class for a transmissive boundary. 1269 # @note Inherits from Boundary. 1238 1270 class Transmissive_momentum_set_stage_boundary(Boundary): 1239 1271 """Returns same momentum conserved quantities as … … 1244 1276 Example: 1245 1277 1246 def waveform(t): 1278 def waveform(t): 1247 1279 return sea_level + normalized_amplitude/cosh(t-25)**2 1248 1280 1249 1281 Bts = Transmissive_momentum_set_stage_boundary(domain, waveform) 1250 1251 1282 1252 1283 Underlying domain must be specified when boundary is instantiated 1253 1284 """ 1254 1285 1255 def __init__(self, domain = None, function=None): 1286 ## 1287 # @brief Instantiate a Reflective_boundary. 1288 # @param domain 1289 # @param function 1290 def __init__(self, domain=None, function=None): 1256 1291 Boundary.__init__(self) 1257 1292 1258 1293 if domain is None: 1259 1294 msg = 'Domain must be specified for this type boundary' 1260 raise msg1295 raise Exception, msg 1261 1296 1262 1297 if function is None: 1263 1298 msg = 'Function must be specified for this type boundary' 1264 raise msg1265 1266 self.domain 1299 raise Exception, msg 1300 1301 self.domain = domain 1267 1302 self.function = function 1268 1303 1304 ## 1305 # @brief Return a representation of this instance. 1269 1306 def __repr__(self): 1270 1307 return 'Transmissive_momentum_set_stage_boundary(%s)' %self.domain 1271 1308 1309 ## 1310 # @brief Calculate transmissive results. 1311 # @param vol_id 1312 # @param edge_id 1272 1313 def evaluate(self, vol_id, edge_id): 1273 1314 """Transmissive momentum set stage boundaries return the edge momentum … … 1277 1318 q = self.domain.get_conserved_quantities(vol_id, edge = edge_id) 1278 1319 1279 1280 1320 t = self.domain.get_time() 1281 1321 1282 1322 if hasattr(self.function, 'time'): 1283 # Roll boundary over if time exceeds 1323 # Roll boundary over if time exceeds 1284 1324 while t > self.function.time[-1]: 1285 msg = 'WARNING: domain time %.2f has exceeded' % t1325 msg = 'WARNING: domain time %.2f has exceeded' % t 1286 1326 msg += 'time provided in ' 1287 1327 msg += 'transmissive_momentum_set_stage boundary object.\n' … … 1290 1330 t -= self.function.time[-1] 1291 1331 1292 1293 1332 value = self.function(t) 1294 1333 try: … … 1296 1335 except: 1297 1336 x = float(value[0]) 1298 1337 1299 1338 q[0] = x 1300 1339 return q 1301 1302 1340 1303 1341 # FIXME: Consider this (taken from File_boundary) to allow … … 1310 1348 1311 1349 1312 # Backward compatibility1313 # FIXME(Ole): Deprecate1350 ## 1351 # @brief Deprecated boundary class. 1314 1352 class Transmissive_Momentum_Set_Stage_boundary(Transmissive_momentum_set_stage_boundary): 1315 1353 pass 1316 1354 1317 1318 1355 1356 ## 1357 # @brief A transmissive boundary, momentum set to zero. 1358 # @note Inherits from Bouondary. 1319 1359 class Transmissive_stage_zero_momentum_boundary(Boundary): 1320 """Return same stage as those present in its neighbour volume. Set momentum to zero. 1360 """Return same stage as those present in its neighbour volume. 1361 Set momentum to zero. 1321 1362 1322 1363 Underlying domain must be specified when boundary is instantiated 1323 1364 """ 1324 1365 1366 ## 1367 # @brief Instantiate a Transmissive (zero momentum) boundary. 1368 # @param domain 1325 1369 def __init__(self, domain=None): 1326 1370 Boundary.__init__(self) 1327 1371 1328 1372 if domain is None: 1329 msg = 'Domain must be specified for '1330 msg += 'Transmissive_stage_zero_momentum boundary'1373 msg = ('Domain must be specified for ' 1374 'Transmissive_stage_zero_momentum boundary') 1331 1375 raise Exception, msg 1332 1376 1333 1377 self.domain = domain 1334 1378 1379 ## 1380 # @brief Return a representation of this instance. 1335 1381 def __repr__(self): 1336 return 'Transmissive_stage_zero_momentum_boundary(%s)' %self.domain 1337 1382 return 'Transmissive_stage_zero_momentum_boundary(%s)' % self.domain 1383 1384 ## 1385 # @brief Calculate transmissive (zero momentum) results. 1386 # @param vol_id 1387 # @param edge_id 1338 1388 def evaluate(self, vol_id, edge_id): 1339 1389 """Transmissive boundaries return the edge values … … 1342 1392 1343 1393 q = self.domain.get_conserved_quantities(vol_id, edge=edge_id) 1344 1394 1345 1395 q[1] = q[2] = 0.0 1346 1396 return q 1347 1397 1348 1398 1349 1399 ## 1400 # @brief Class for a Dirichlet discharge boundary. 1401 # @note Inherits from Boundary. 1350 1402 class Dirichlet_discharge_boundary(Boundary): 1351 1403 """ … … 1356 1408 """ 1357 1409 1410 ## 1411 # @brief Instantiate a Dirichlet discharge boundary. 1412 # @param domain 1413 # @param stage0 1414 # @param wh0 1358 1415 def __init__(self, domain=None, stage0=None, wh0=None): 1359 1416 Boundary.__init__(self) 1360 1417 1361 1418 if domain is None: 1362 msg = 'Domain must be specified for this type boundary'1363 raise msg1419 msg = 'Domain must be specified for this type of boundary' 1420 raise Exception, msg 1364 1421 1365 1422 if stage0 is None: 1366 raise 'set stage'1423 raise Exception, 'Stage must be specified for this type of boundary' 1367 1424 1368 1425 if wh0 is None: 1369 1426 wh0 = 0.0 1370 1427 1371 self.domain 1372 self.stage0 1428 self.domain = domain 1429 self.stage0 = stage0 1373 1430 self.wh0 = wh0 1374 1431 1432 ## 1433 # @brief Return a representation of this instance. 1375 1434 def __repr__(self): 1376 return 'Dirichlet_Discharge_boundary(%s)' %self.domain 1377 1435 return 'Dirichlet_Discharge_boundary(%s)' % self.domain 1436 1437 ## 1438 # @brief Calculate Dirichlet discharge boundary results. 1439 # @param vol_id 1440 # @param edge_id 1378 1441 def evaluate(self, vol_id, edge_id): 1379 """Set discharge in the (inward) normal direction 1380 """ 1442 """Set discharge in the (inward) normal direction""" 1381 1443 1382 1444 normal = self.domain.get_normal(vol_id,edge_id) 1383 1445 q = [self.stage0, -self.wh0*normal[0], -self.wh0*normal[1]] 1384 1446 return q 1385 1386 1447 1387 1448 # FIXME: Consider this (taken from File_boundary) to allow … … 1394 1455 1395 1456 1396 1397 # Backward compatibility 1457 # Backward compatibility 1398 1458 # FIXME(Ole): Deprecate 1459 ## 1460 # @brief Deprecated 1399 1461 class Dirichlet_Discharge_boundary(Dirichlet_discharge_boundary): 1400 1462 pass 1401 1402 1403 1404 1405 1463 1464 1406 1465 class Inflow_boundary(Boundary): 1407 1466 """Apply given flow in m^3/s to boundary segment. … … 1491 1550 1492 1551 class Field_boundary(Boundary): 1493 """Set boundary from given field represented in an sww file containing values 1494 for stage, xmomentum and ymomentum. 1495 Optionally, the user can specify mean_stage to offset the stage provided in the 1496 sww file. 1497 1498 This function is a thin wrapper around the generic File_boundary. The 1552 """Set boundary from given field represented in an sww file containing 1553 values for stage, xmomentum and ymomentum. 1554 1555 Optionally, the user can specify mean_stage to offset the stage provided 1556 in the sww file. 1557 1558 This function is a thin wrapper around the generic File_boundary. The 1499 1559 difference between the file_boundary and field_boundary is only that the 1500 1560 field_boundary will allow you to change the level of the stage height when 1501 you read in the boundary condition. This is very useful when running 1502 different tide heights in the same area as you need only to convert one 1503 boundary condition to a SWW file, ideally for tide height of 0 m 1561 you read in the boundary condition. This is very useful when running 1562 different tide heights in the same area as you need only to convert one 1563 boundary condition to a SWW file, ideally for tide height of 0 m 1504 1564 (saving disk space). Then you can use field_boundary to read this SWW file 1505 1565 and change the stage height (tide) on the fly depending on the scenario. 1506 1507 1566 """ 1508 1567 1509 1510 def __init__(self, filename, domain, 1568 ## 1569 # @brief Construct an instance of a 'field' boundary. 1570 # @param filename Name of SWW file containing stage, x and ymomentum. 1571 # @param domain Shallow water domain for which the boundary applies. 1572 # @param mean_stage Mean water level added to stage derived from SWW file. 1573 # @param time_thinning Time step thinning factor. 1574 # @param time_limit 1575 # @param boundary_polygon 1576 # @param default_boundary None or an instance of Boundary. 1577 # @param use_cache True if caching is to be used. 1578 # @param verbose True if this method is to be verbose. 1579 def __init__(self, 1580 filename, 1581 domain, 1511 1582 mean_stage=0.0, 1512 1583 time_thinning=1, 1513 1584 time_limit=None, 1514 boundary_polygon=None, 1515 default_boundary=None, 1585 boundary_polygon=None, 1586 default_boundary=None, 1516 1587 use_cache=False, 1517 1588 verbose=False): … … 1523 1594 from the boundary condition 1524 1595 time_thinning: Will set how many time steps from the sww file read in 1525 will be interpolated to the boundary. For example if 1596 will be interpolated to the boundary. For example if 1526 1597 the sww file has 1 second time steps and is 24 hours 1527 in length it has 86400 time steps. If you set 1528 time_thinning to 1 it will read all these steps. 1598 in length it has 86400 time steps. If you set 1599 time_thinning to 1 it will read all these steps. 1529 1600 If you set it to 100 it will read every 100th step eg 1530 1601 only 864 step. This parameter is very useful to increase 1531 the speed of a model run that you are setting up 1602 the speed of a model run that you are setting up 1532 1603 and testing. 1533 1534 default_boundary: Must be either None or an instance of a 1604 1605 default_boundary: Must be either None or an instance of a 1535 1606 class descending from class Boundary. 1536 This will be used in case model time exceeds 1607 This will be used in case model time exceeds 1537 1608 that available in the underlying data. 1609 1538 1610 Note that mean_stage will also be added to this. 1539 1611 1540 1612 use_cache: 1541 1613 verbose: 1542 1543 1614 """ 1544 1615 1545 1616 # Create generic file_boundary object 1546 self.file_boundary = File_boundary(filename, domain, 1617 self.file_boundary = File_boundary(filename, 1618 domain, 1547 1619 time_thinning=time_thinning, 1548 1620 time_limit=time_limit, … … 1552 1624 verbose=verbose) 1553 1625 1554 1555 1626 # Record information from File_boundary 1556 1627 self.F = self.file_boundary.F 1557 1628 self.domain = self.file_boundary.domain 1558 1629 1559 1630 # Record mean stage 1560 1631 self.mean_stage = mean_stage 1561 1632 1562 1633 ## 1634 # @note Generate a string representation of this instance. 1563 1635 def __repr__(self): 1564 1636 return 'Field boundary' 1565 1637 1566 1638 ## 1639 # @brief Calculate 'field' boundary results. 1640 # @param vol_id 1641 # @param edge_id 1567 1642 def evaluate(self, vol_id=None, edge_id=None): 1568 1643 """Return linearly interpolated values based on domain.time … … 1570 1645 vol_id and edge_id are ignored 1571 1646 """ 1572 1647 1573 1648 # Evaluate file boundary 1574 1649 q = self.file_boundary.evaluate(vol_id, edge_id) … … 1580 1655 return q 1581 1656 1582 1583 1584 #----------------------- 1657 1658 ################################################################################ 1585 1659 # Standard forcing terms 1586 #----------------------- 1587 1660 ################################################################################ 1661 1662 ## 1663 # @brief Apply gravitational pull in the presence of bed slope. 1664 # @param domain The domain to apply gravity to. 1665 # @note Wrapper for C function gravity_c(). 1588 1666 def gravity(domain): 1589 1667 """Apply gravitational pull in the presence of bed slope … … 1591 1669 """ 1592 1670 1671 from shallow_water_ext import gravity as gravity_c 1672 1593 1673 xmom = domain.quantities['xmomentum'].explicit_update 1594 1674 ymom = domain.quantities['ymomentum'].explicit_update … … 1602 1682 x = domain.get_vertex_coordinates() 1603 1683 g = domain.g 1604 1605 1606 from shallow_water_ext import gravity as gravity_c 1607 gravity_c(g, h, z, x, xmom, ymom) #, 1.0e-6) 1608 1609 1610 1684 1685 gravity_c(g, h, z, x, xmom, ymom) #, 1.0e-6) 1686 1687 ## 1688 # @brief Apply friction to a surface (implicit). 1689 # @param domain The domain to apply Manning friction to. 1690 # @note Wrapper for C function manning_friction_c(). 1611 1691 def manning_friction_implicit(domain): 1612 """Apply (Manning) friction to water momentum 1692 """Apply (Manning) friction to water momentum 1613 1693 Wrapper for c version 1614 1694 """ 1615 1695 1616 1617 #print 'Implicit friction' 1696 from shallow_water_ext import manning_friction as manning_friction_c 1618 1697 1619 1698 xmom = domain.quantities['xmomentum'] … … 1634 1713 g = domain.g 1635 1714 1636 from shallow_water_ext import manning_friction as manning_friction_c1637 1715 manning_friction_c(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update) 1638 1716 1639 1717 ## 1718 # @brief Apply friction to a surface (explicit). 1719 # @param domain The domain to apply Manning friction to. 1720 # @note Wrapper for C function manning_friction_c(). 1640 1721 def manning_friction_explicit(domain): 1641 """Apply (Manning) friction to water momentum 1722 """Apply (Manning) friction to water momentum 1642 1723 Wrapper for c version 1643 1724 """ 1644 1725 1645 # print 'Explicit friction'1726 from shallow_water_ext import manning_friction as manning_friction_c 1646 1727 1647 1728 xmom = domain.quantities['xmomentum'] … … 1662 1743 g = domain.g 1663 1744 1664 from shallow_water_ext import manning_friction as manning_friction_c1665 1745 manning_friction_c(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update) 1666 1746 1667 1747 1668 1748 # FIXME (Ole): This was implemented for use with one of the analytical solutions (Sampson?) 1669 # Is it still needed (30 Oct 2007)? 1749 ## 1750 # @brief Apply linear friction to a surface. 1751 # @param domain The domain to apply Manning friction to. 1752 # @note Is this still used (30 Oct 2007)? 1670 1753 def linear_friction(domain): 1671 1754 """Apply linear friction to water momentum … … 1699 1782 xmom_update[k] += S*uh[k] 1700 1783 ymom_update[k] += S*vh[k] 1701 1702 1703 1704 1784 1705 1785 def depth_dependent_friction(domain, default_friction, … … 1722 1802 """ 1723 1803 1724 import Numericas num1804 import numpy as num 1725 1805 1726 1806 # Create a temp array to store updated depth dependent friction for wet elements 1727 1807 # EHR this is outwardly inneficient but not obvious how to avoid recreating each call?????? 1728 1808 N=len(domain) 1729 wet_friction = num.zeros(N, num. Float)1809 wet_friction = num.zeros(N, num.float) 1730 1810 wet_friction[:] = default_n0 # Initially assign default_n0 to all array so sure have no zeros values 1731 1811 … … 1774 1854 1775 1855 1776 1777 1778 1779 1780 #--------------------------------- 1856 ################################################################################ 1781 1857 # Experimental auxiliary functions 1782 #--------------------------------- 1858 ################################################################################ 1859 1860 ## 1861 # @brief Check forcefield parameter. 1862 # @param f Object to check. 1863 # @note 'f' may be a callable object or a scalar value. 1783 1864 def check_forcefield(f): 1784 """Check that f is either 1865 """Check that force object is as expected. 1866 1867 Check that f is either: 1785 1868 1: a callable object f(t,x,y), where x and y are vectors 1786 1869 and that it returns an array or a list of same length … … 1791 1874 if callable(f): 1792 1875 N = 3 1793 x = num.ones(3, num. Float)1794 y = num.ones(3, num. Float)1876 x = num.ones(3, num.float) 1877 y = num.ones(3, num.float) 1795 1878 try: 1796 1879 q = f(1.0, x=x, y=y) … … 1798 1881 msg = 'Function %s could not be executed:\n%s' %(f, e) 1799 1882 # FIXME: Reconsider this semantics 1800 raise msg1883 raise Exception, msg 1801 1884 1802 1885 try: 1803 q = num.array(q, num. Float)1886 q = num.array(q, num.float) 1804 1887 except: 1805 msg = 'Return value from vector function %s could ' %f1806 msg += 'not be converted into a Numeric array of floats.\n'1807 msg += 'Specified function should return either list or array.'1808 raise msg1888 msg = ('Return value from vector function %s could not ' 1889 'be converted into a numeric array of floats.\nSpecified ' 1890 'function should return either list or array.' % f) 1891 raise Exception, msg 1809 1892 1810 1893 # Is this really what we want? 1811 msg = 'Return vector from function %s ' %f 1812 msg += 'must have same lenght as input vectors' 1813 assert len(q) == N, msg 1814 1894 # info is "(func name, filename, defining line)" 1895 func_info = (f.func_name, f.func_code.co_filename, 1896 f.func_code.co_firstlineno) 1897 func_msg = 'Function %s (defined in %s, line %d)' % func_info 1898 try: 1899 result_len = len(q) 1900 except: 1901 msg = '%s must return vector' % func_msg 1902 self.fail(msg) 1903 msg = '%s must return vector of length %d' % (func_msg, N) 1904 assert result_len == N, msg 1815 1905 else: 1816 1906 try: 1817 1907 f = float(f) 1818 1908 except: 1819 msg = 'Force field %s must be either a scalar' %f 1820 msg += ' or a vector function' 1821 raise Exception(msg) 1909 msg = ('Force field %s must be a scalar value coercible to float.' 1910 % str(f)) 1911 raise Exception, msg 1912 1822 1913 return f 1823 1914 1824 1915 1916 ## 1917 # Class to apply a wind stress to a domain. 1825 1918 class Wind_stress: 1826 1919 """Apply wind stress to water momentum in terms of … … 1828 1921 """ 1829 1922 1923 ## 1924 # @brief Create an instance of Wind_stress. 1925 # @param *args 1926 # @param **kwargs 1830 1927 def __init__(self, *args, **kwargs): 1831 1928 """Initialise windfield from wind speed s [m/s] … … 1880 1977 else: 1881 1978 # Assume info is in 2 keyword arguments 1882 1883 1979 if len(kwargs) == 2: 1884 1980 s = kwargs['s'] 1885 1981 phi = kwargs['phi'] 1886 1982 else: 1887 raise 'Assumes two keyword arguments: s=..., phi=....'1983 raise Exception, 'Assumes two keyword arguments: s=..., phi=....' 1888 1984 1889 1985 self.speed = check_forcefield(s) … … 1892 1988 self.const = eta_w*rho_a/rho_w 1893 1989 1894 1990 ## 1991 # @brief 'execute' this class instance. 1992 # @param domain 1895 1993 def __call__(self, domain): 1896 """Evaluate windfield based on values found in domain 1897 """ 1994 """Evaluate windfield based on values found in domain""" 1898 1995 1899 1996 from math import pi, cos, sin, sqrt … … 1902 1999 ymom_update = domain.quantities['ymomentum'].explicit_update 1903 2000 1904 N = len(domain) # number_of_triangles2001 N = len(domain) # number_of_triangles 1905 2002 t = domain.time 1906 2003 … … 1910 2007 else: 1911 2008 # Assume s is a scalar 1912 1913 2009 try: 1914 s_vec = self.speed * num.ones(N, num. Float)2010 s_vec = self.speed * num.ones(N, num.float) 1915 2011 except: 1916 2012 msg = 'Speed must be either callable or a scalar: %s' %self.s 1917 2013 raise msg 1918 1919 2014 1920 2015 if callable(self.phi): … … 1925 2020 1926 2021 try: 1927 phi_vec = self.phi * num.ones(N, num. Float)2022 phi_vec = self.phi * num.ones(N, num.float) 1928 2023 except: 1929 2024 msg = 'Angle must be either callable or a scalar: %s' %self.phi … … 1934 2029 1935 2030 2031 ## 2032 # @brief Assign wind field values 2033 # @param xmom_update 2034 # @param ymom_update 2035 # @param s_vec 2036 # @param phi_vec 2037 # @param const 1936 2038 def assign_windfield_values(xmom_update, ymom_update, 1937 2039 s_vec, phi_vec, const): 1938 2040 """Python version of assigning wind field to update vectors. 1939 A cversion also exists (for speed)2041 A C version also exists (for speed) 1940 2042 """ 2043 1941 2044 from math import pi, cos, sin, sqrt 1942 2045 … … 1959 2062 1960 2063 1961 1962 1963 2064 ## 2065 # @brief A class for a general explicit forcing term. 1964 2066 class General_forcing: 1965 2067 """General explicit forcing term for update of quantity 1966 2068 1967 2069 This is used by Inflow and Rainfall for instance 1968 2070 1969 2071 1970 2072 General_forcing(quantity_name, rate, center, radius, polygon) 1971 2073 1972 2074 domain: ANUGA computational domain 1973 quantity_name: Name of quantity to update. 2075 quantity_name: Name of quantity to update. 1974 2076 It must be a known conserved quantity. 1975 1976 rate [?/s]: Total rate of change over the specified area. 2077 2078 rate [?/s]: Total rate of change over the specified area. 1977 2079 This parameter can be either a constant or a 1978 function of time. Positive values indicate increases, 2080 function of time. Positive values indicate increases, 1979 2081 negative values indicate decreases. 1980 2082 Rate can be None at initialisation but must be specified … … 1990 2092 Either center, radius or polygon can be specified but not both. 1991 2093 If neither are specified the entire domain gets updated. 1992 All coordinates to be specified in absolute UTM coordinates (x, y) assuming the zone of domain. 1993 2094 All coordinates to be specified in absolute UTM coordinates (x, y) assuming the zone of domain. 2095 1994 2096 Inflow or Rainfall for examples of use 1995 2097 """ … … 1998 2100 # FIXME (AnyOne) : Add various methods to allow spatial variations 1999 2101 2102 ## 2103 # @brief Create an instance of this forcing term. 2104 # @param domain 2105 # @param quantity_name 2106 # @param rate 2107 # @param center 2108 # @param radius 2109 # @param polygon 2110 # @param default_rate 2111 # @param verbose 2000 2112 def __init__(self, 2001 2113 domain, 2002 2114 quantity_name, 2003 2115 rate=0.0, 2004 center=None, radius=None, 2116 center=None, 2117 radius=None, 2005 2118 polygon=None, 2006 2119 default_rate=None, 2007 2120 verbose=False): 2008 2121 2122 from math import pi, cos, sin 2123 2009 2124 if center is None: 2010 msg = 'I got radius but no center.' 2125 msg = 'I got radius but no center.' 2011 2126 assert radius is None, msg 2012 2127 2013 2128 if radius is None: 2014 msg += 'I got center but no radius.' 2129 msg += 'I got center but no radius.' 2015 2130 assert center is None, msg 2016 2017 2018 2019 from math import pi, cos, sin2020 2131 2021 2132 self.domain = domain … … 2024 2135 self.center = ensure_numeric(center) 2025 2136 self.radius = radius 2026 self.polygon = polygon 2137 self.polygon = polygon 2027 2138 self.verbose = verbose 2028 self.value = 0.0 # Can be used to remember value at2029 # previous timestep in order to obtain rate2139 self.value = 0.0 # Can be used to remember value at 2140 # previous timestep in order to obtain rate 2030 2141 2031 2142 # Get boundary (in absolute coordinates) 2032 2143 bounding_polygon = domain.get_boundary_polygon() 2033 2034 2144 2035 2145 # Update area if applicable … … 2039 2149 assert polygon is None, msg 2040 2150 2041 2042 2151 # Check that circle center lies within the mesh. 2043 msg = 'Center %s specified for forcing term did not' % (str(center))2152 msg = 'Center %s specified for forcing term did not' % str(center) 2044 2153 msg += 'fall within the domain boundary.' 2045 2154 assert is_inside_polygon(center, bounding_polygon), msg … … 2049 2158 periphery_points = [] 2050 2159 for i in range(N): 2051 2052 2160 theta = 2*pi*i/100 2053 2161 2054 2162 x = center[0] + radius*cos(theta) 2055 2163 y = center[1] + radius*sin(theta) … … 2057 2165 periphery_points.append([x,y]) 2058 2166 2059 2060 2167 for point in periphery_points: 2061 msg = 'Point %s on periphery for forcing term' % (str(point))2168 msg = 'Point %s on periphery for forcing term' % str(point) 2062 2169 msg += ' did not fall within the domain boundary.' 2063 2170 assert is_inside_polygon(point, bounding_polygon), msg 2064 2171 2065 2066 2172 if polygon is not None: 2067 2068 2173 # Check that polygon lies within the mesh. 2069 2174 for point in self.polygon: 2070 msg = 'Point %s in polygon for forcing term' % (point)2175 msg = 'Point %s in polygon for forcing term' % str(point) 2071 2176 msg += ' did not fall within the domain boundary.' 2072 2177 assert is_inside_polygon(point, bounding_polygon), msg 2073 2074 2075 2178 2076 2179 # Pointer to update vector … … 2078 2181 2079 2182 # Determine indices in flow area 2080 N = len(domain) 2183 N = len(domain) 2081 2184 points = domain.get_centroid_coordinates(absolute=True) 2082 2185 … … 2085 2188 if self.center is not None and self.radius is not None: 2086 2189 # Inlet is circular 2087 2088 inlet_region = 'center=%s, radius=%s' %(self.center, self.radius) 2089 2190 inlet_region = 'center=%s, radius=%s' % (self.center, self.radius) 2191 2090 2192 self.exchange_indices = [] 2091 2193 for k in range(N): 2092 x, y = points[k,:] # Centroid2093 2194 x, y = points[k,:] # Centroid 2195 2094 2196 c = self.center 2095 2197 if ((x-c[0])**2+(y-c[1])**2) < self.radius**2: 2096 2198 self.exchange_indices.append(k) 2097 2098 if self.polygon is not None: 2199 2200 if self.polygon is not None: 2099 2201 # Inlet is polygon 2100 2101 2202 inlet_region = 'polygon=%s' % (self.polygon) 2102 2203 self.exchange_indices = inside_polygon(points, self.polygon) 2103 2104 2105 2204 2106 2205 if self.exchange_indices is None: 2107 2206 self.exchange_area = polygon_area(bounding_polygon) … … 2109 2208 if len(self.exchange_indices) == 0: 2110 2209 msg = 'No triangles have been identified in ' 2111 msg += 'specified region: %s' % inlet_region2210 msg += 'specified region: %s' % inlet_region 2112 2211 raise Exception, msg 2113 2212 … … 2127 2226 2128 2227 # Check and store default_rate 2129 msg = 'Keyword argument default_rate must be either None '2130 msg += 'or a function of time.\n I got %s' %(str(default_rate))2131 assert default_rate is None or \2132 type(default_rate) in [IntType, FloatType] or \2133 callable(default_rate), msg2134 2228 msg = ('Keyword argument default_rate must be either None ' 2229 'or a function of time.\nI got %s.' % str(default_rate)) 2230 assert (default_rate is None or 2231 type(default_rate) in [IntType, FloatType] or 2232 callable(default_rate)), msg 2233 2135 2234 if default_rate is not None: 2136 2137 2235 # If it is a constant, make it a function 2138 2236 if not callable(default_rate): … … 2140 2238 default_rate = lambda t: tmp 2141 2239 2142 2143 2240 # Check that default_rate is a function of one argument 2144 2241 try: … … 2148 2245 2149 2246 self.default_rate = default_rate 2150 self.default_rate_invoked = False # Flag 2151 2152 2247 self.default_rate_invoked = False # Flag 2248 2249 ## 2250 # @brief Execute this instance. 2251 # @param domain 2153 2252 def __call__(self, domain): 2154 """Apply inflow function at time specified in domain and update stage 2155 """ 2253 """Apply inflow function at time specified in domain, update stage""" 2156 2254 2157 2255 # Call virtual method allowing local modifications 2158 2159 2256 t = domain.get_time() 2160 2257 try: … … 2164 2261 except Modeltime_too_late, e: 2165 2262 if self.default_rate is None: 2166 raise Exception, e # Reraise exception2263 raise Exception, e # Reraise exception 2167 2264 else: 2168 2265 # Pass control to default rate function 2169 2266 rate = self.default_rate(t) 2170 2267 2171 2268 if self.default_rate_invoked is False: 2172 2269 # Issue warning the first time 2173 msg = '%s' %str(e)2174 msg += 'Instead I will use the default rate: %s\n'\2175 %str(self.default_rate)2176 msg += 'Note: Further warnings will be supressed'2270 msg = ('%s\n' 2271 'Instead I will use the default rate: %s\n' 2272 'Note: Further warnings will be supressed' 2273 % (str(e), str(self.default_rate))) 2177 2274 warn(msg) 2178 2275 2179 2276 # FIXME (Ole): Replace this crude flag with 2180 2277 # Python's ability to print warnings only once. 2181 2278 # See http://docs.python.org/lib/warning-filter.html 2182 2279 self.default_rate_invoked = True 2183 2184 2185 2186 2187 2280 2188 2281 if rate is None: 2189 msg = 'Attribute rate must be specified in General_forcing'2190 msg += ' or its descendants before attempting to call it'2282 msg = ('Attribute rate must be specified in General_forcing ' 2283 'or its descendants before attempting to call it') 2191 2284 raise Exception, msg 2192 2193 2285 2194 2286 # Now rate is a number 2195 2287 if self.verbose is True: 2196 print 'Rate of %s at time = %.2f = %f' %(self.quantity_name, 2197 domain.get_time(), 2198 rate) 2199 2288 print 'Rate of %s at time = %.2f = %f' % (self.quantity_name, 2289 domain.get_time(), 2290 rate) 2200 2291 2201 2292 if self.exchange_indices is None: … … 2206 2297 self.update[k] += rate 2207 2298 2208 2299 ## 2300 # @brief Update the internal rate. 2301 # @param t A callable or scalar used to set the rate. 2302 # @return The new rate. 2209 2303 def update_rate(self, t): 2210 2304 """Virtual method allowing local modifications by writing an 2211 2305 overriding version in descendant 2212 2213 """ 2214 2215 2216 2217 2306 """ 2307 2308 if callable(self.rate): 2309 rate = self.rate(t) 2310 else: 2311 rate = self.rate 2218 2312 2219 2313 return rate 2220 2314 2221 2315 ## 2316 # @brief Get values for the specified quantity. 2317 # @param quantity_name Name of the quantity of interest. 2318 # @return The value(s) of the quantity. 2319 # @note If 'quantity_name' is None, use self.quantity_name. 2222 2320 def get_quantity_values(self, quantity_name=None): 2223 """Return values for specified quantity restricted to opening 2224 2321 """Return values for specified quantity restricted to opening 2322 2225 2323 Optionally a quantity name can be specified if values from another 2226 2324 quantity is sought 2227 2325 """ 2228 2326 2229 2327 if quantity_name is None: 2230 2328 quantity_name = self.quantity_name 2231 2329 2232 2330 q = self.domain.quantities[quantity_name] 2233 return q.get_values(location='centroids', 2331 return q.get_values(location='centroids', 2234 2332 indices=self.exchange_indices) 2235 2236 2333 2334 ## 2335 # @brief Set value for the specified quantity. 2336 # @param val The value object used to set value. 2337 # @param quantity_name Name of the quantity of interest. 2338 # @note If 'quantity_name' is None, use self.quantity_name. 2237 2339 def set_quantity_values(self, val, quantity_name=None): 2238 """Set values for specified quantity restricted to opening 2239 2340 """Set values for specified quantity restricted to opening 2341 2240 2342 Optionally a quantity name can be specified if values from another 2241 quantity is sought 2343 quantity is sought 2242 2344 """ 2243 2345 2244 2346 if quantity_name is None: 2245 2347 quantity_name = self.quantity_name 2246 2247 q = self.domain.quantities[self.quantity_name] 2248 q.set_values(val, 2249 location='centroids', 2250 indices=self.exchange_indices) 2251 2252 2253 2348 2349 q = self.domain.quantities[self.quantity_name] 2350 q.set_values(val, 2351 location='centroids', 2352 indices=self.exchange_indices) 2353 2354 2355 ## 2356 # @brief A class for rainfall forcing function. 2357 # @note Inherits from General_forcing. 2254 2358 class Rainfall(General_forcing): 2255 2359 """Class Rainfall - general 'rain over entire domain' forcing term. 2256 2360 2257 2361 Used for implementing Rainfall over the entire domain. 2258 2259 2260 2261 Need to add Spatial Varying Capability 2262 2263 2362 2363 Current Limited to only One Gauge.. 2364 2365 Need to add Spatial Varying Capability 2366 (This module came from copying and amending the Inflow Code) 2367 2264 2368 Rainfall(rain) 2265 2369 2266 domain 2267 rain [mm/s]: Total rain rate over the specified domain. 2268 2269 2370 domain 2371 rain [mm/s]: Total rain rate over the specified domain. 2372 NOTE: Raingauge Data needs to reflect the time step. 2373 IE: if Gauge is mm read at a time step, then the input 2270 2374 here is as mm/(timeStep) so 10mm in 5minutes becomes 2271 2375 10/(5x60) = 0.0333mm/s. 2272 2273 2376 2274 2377 This parameter can be either a constant or a 2275 function of time. Positive values indicate inflow, 2378 function of time. Positive values indicate inflow, 2276 2379 negative values indicate outflow. 2277 2380 (and be used for Infiltration - Write Seperate Module) … … 2279 2382 the inflow region and then applied to update the 2280 2383 stage quantity. 2281 2282 Note also that any reference to the internal attribute 'rate'2283 later will use the one converted to m/s.2284 2384 2285 2385 polygon: Specifies a polygon to restrict the rainfall. 2286 2386 2287 2387 Examples 2288 2388 How to put them in a run File... 2289 2389 2290 2390 #------------------------------------------------------------------------ 2291 2391 # Setup specialised forcing terms … … 2294 2394 # input of Rainfall in mm/s 2295 2395 2296 catchmentrainfall = Rainfall(rain=file_function('Q100_2hr_Rain.tms')) 2396 catchmentrainfall = Rainfall(rain=file_function('Q100_2hr_Rain.tms')) 2297 2397 # Note need path to File in String. 2298 2398 # Else assumed in same directory … … 2301 2401 """ 2302 2402 2303 2403 ## 2404 # @brief Create an instance of the class. 2405 # @param domain Domain of interest. 2406 # @param rate Total rain rate over the specified domain (mm/s). 2407 # @param center 2408 # @param radius 2409 # @param polygon Polygon to restrict rainfall. 2410 # @param default_rate 2411 # @param verbose True if this instance is to be verbose. 2304 2412 def __init__(self, 2305 2413 domain, 2306 rate=0.0, 2307 center=None, radius=None, 2414 rate=0.0, 2415 center=None, 2416 radius=None, 2308 2417 polygon=None, 2309 default_rate=None, 2418 default_rate=None, 2310 2419 verbose=False): 2311 2420 … … 2316 2425 rain = rate/1000.0 2317 2426 2318 if default_rate is not None: 2427 if default_rate is not None: 2319 2428 if callable(default_rate): 2320 2429 default_rain = lambda t: default_rate(t)/1000.0 … … 2324 2433 default_rain = None 2325 2434 2326 2327 2435 2328 2436 … … 2331 2439 'stage', 2332 2440 rate=rain, 2333 center=center, radius=radius, 2441 center=center, 2442 radius=radius, 2334 2443 polygon=polygon, 2335 2444 default_rate=default_rain, 2336 2445 verbose=verbose) 2337 2446 2338 2339 2340 2341 2342 2447 2448 ## 2449 # @brief A class for inflow (rain and drain) forcing function. 2450 # @note Inherits from General_forcing. 2343 2451 class Inflow(General_forcing): 2344 2452 """Class Inflow - general 'rain and drain' forcing term. 2345 2453 2346 2454 Useful for implementing flows in and out of the domain. 2347 2455 2348 2456 Inflow(flow, center, radius, polygon) 2349 2457 2350 2458 domain 2351 rate [m^3/s]: Total flow rate over the specified area. 2459 rate [m^3/s]: Total flow rate over the specified area. 2352 2460 This parameter can be either a constant or a 2353 function of time. Positive values indicate inflow, 2461 function of time. Positive values indicate inflow, 2354 2462 negative values indicate outflow. 2355 2463 The specified flow will be divided by the area of 2356 the inflow region and then applied to update stage. 2464 the inflow region and then applied to update stage. 2357 2465 center [m]: Coordinates at center of flow point 2358 2466 radius [m]: Size of circular area … … 2360 2468 2361 2469 Either center, radius or polygon must be specified 2362 2470 2363 2471 Examples 2364 2472 2365 2473 # Constant drain at 0.003 m^3/s. 2366 2474 # The outflow area is 0.07**2*pi=0.0154 m^2 2367 # This corresponds to a rate of change of 0.003/0.0154 = 0.2 m/s 2368 # 2475 # This corresponds to a rate of change of 0.003/0.0154 = 0.2 m/s 2476 # 2369 2477 Inflow((0.7, 0.4), 0.07, -0.003) 2370 2478 … … 2372 2480 # Tap turning up to a maximum inflow of 0.0142 m^3/s. 2373 2481 # The inflow area is 0.03**2*pi = 0.00283 m^2 2374 # This corresponds to a rate of change of 0.0142/0.00283 = 5 m/s 2482 # This corresponds to a rate of change of 0.0142/0.00283 = 5 m/s 2375 2483 # over the specified area 2376 2484 Inflow((0.5, 0.5), 0.03, lambda t: min(0.01*t, 0.0142)) … … 2387 2495 2388 2496 domain.forcing_terms.append(hydrograph) 2389 2390 2497 """ 2391 2498 2392 2499 ## 2500 # @brief Create an instance of the class. 2501 # @param domain Domain of interest. 2502 # @param rate Total rain rate over the specified domain (mm/s). 2503 # @param center 2504 # @param radius 2505 # @param polygon Polygon to restrict rainfall. 2506 # @param default_rate 2507 # @param verbose True if this instance is to be verbose. 2393 2508 def __init__(self, 2394 2509 domain, 2395 rate=0.0, 2396 center=None, radius=None, 2510 rate=0.0, 2511 center=None, 2512 radius=None, 2397 2513 polygon=None, 2398 2514 default_rate=None, 2399 verbose=False): 2400 2401 2515 verbose=False): 2402 2516 # Create object first to make area is available 2403 2517 General_forcing.__init__(self, … … 2405 2519 'stage', 2406 2520 rate=rate, 2407 center=center, radius=radius, 2521 center=center, 2522 radius=radius, 2408 2523 polygon=polygon, 2409 2524 default_rate=default_rate, 2410 2525 verbose=verbose) 2411 2526 2527 ## 2528 # @brief Update the instance rate. 2529 # @param t New rate object. 2412 2530 def update_rate(self, t): 2413 2531 """Virtual method allowing local modifications by writing an 2414 2532 overriding version in descendant 2415 2533 2416 This one converts m^3/s to m/s which can be added directly 2534 This one converts m^3/s to m/s which can be added directly 2417 2535 to 'stage' in ANUGA 2418 2536 """ 2419 2537 2420 2421 2422 2423 2538 if callable(self.rate): 2539 _rate = self.rate(t)/self.exchange_area 2540 else: 2541 _rate = self.rate/self.exchange_area 2424 2542 2425 2543 return _rate 2426 2544 2427 2545 2428 2429 2430 #------------------ 2546 ################################################################################ 2431 2547 # Initialise module 2432 #------------------ 2433 2548 ################################################################################ 2434 2549 2435 2550 from anuga.utilities import compile 2436 2551 if compile.can_use_C_extension('shallow_water_ext.c'): 2437 # Underlying C implementations can be accessed 2438 2552 # Underlying C implementations can be accessed 2439 2553 from shallow_water_ext import rotate, assign_windfield_values 2440 2554 else: 2441 msg = 'C implementations could not be accessed by %s.\n ' % __file__2555 msg = 'C implementations could not be accessed by %s.\n ' % __file__ 2442 2556 msg += 'Make sure compile_all.py has been run as described in ' 2443 2557 msg += 'the ANUGA installation guide.' 2444 2558 raise Exception, msg 2445 2446 2559 2447 2560 # Optimisation with psyco … … 2456 2569 #Psyco isn't supported on 64 bit systems, but it doesn't matter 2457 2570 else: 2458 msg = 'WARNING: psyco (speedup) could not import'+\2459 ', you may want to consider installing it'2571 msg = ('WARNING: psyco (speedup) could not be imported, ' 2572 'you may want to consider installing it') 2460 2573 print msg 2461 2574 else: … … 2463 2576 psyco.bind(Domain.compute_fluxes) 2464 2577 2578 2465 2579 if __name__ == "__main__": 2466 2580 pass 2467 2468
Note: See TracChangeset
for help on using the changeset viewer.
