Changeset 6181
- Timestamp:
- Jan 16, 2009, 4:06:55 PM (16 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py
r6145 r6181 38 38 39 39 40 ## 41 # @brief Basic Domain class 42 # @note Subclasses Mesh. 40 43 class Domain(Mesh): 41 44 42 43 def __init__(self, 44 source=None, 45 triangles=None, 46 boundary=None, 47 conserved_quantities=None, 48 other_quantities=None, 49 tagged_elements=None, 50 geo_reference=None, 51 use_inscribed_circle=False, 52 mesh_filename=None, 53 use_cache=False, 54 verbose=False, 55 full_send_dict=None, 56 ghost_recv_dict=None, 57 processor=0, 58 numproc=1, 59 number_of_full_nodes=None, 60 number_of_full_triangles=None): 61 45 ## 46 # @brief Generic computational Domain constructor. 47 # @param source Name of mesh file or coords of mesh vertices. 48 # @param triangles Mesh connectivity (see mesh.py for more information). 49 # @param boundary (see mesh.py for more information) 50 # @param conserved_quantities List of names of quantities to be conserved. 51 # @param other_quantities List of names of other quantities. 52 # @param tagged_elements 53 # @param geo_reference 54 # @param use_inscribed_circle 55 # @param mesh_filename 56 # @param use_cache 57 # @param verbose 58 # @param full_send_dict 59 # @param ghost_recv_dict 60 # @param processor 61 # @param numproc 62 # @param number_of_full_nodes 63 # @param number_of_full_triangles 64 def __init__(self, source=None, 65 triangles=None, 66 boundary=None, 67 conserved_quantities=None, 68 other_quantities=None, 69 tagged_elements=None, 70 geo_reference=None, 71 use_inscribed_circle=False, 72 mesh_filename=None, 73 use_cache=False, 74 verbose=False, 75 full_send_dict=None, 76 ghost_recv_dict=None, 77 processor=0, 78 numproc=1, 79 number_of_full_nodes=None, 80 number_of_full_triangles=None): 62 81 63 82 """Instantiate generic computational Domain. … … 76 95 tagged_elements: 77 96 ... 78 79 80 97 """ 81 98 … … 85 102 else: 86 103 coordinates = source 87 88 104 89 105 # In case a filename has been specified, extract content … … 94 110 use_cache=use_cache, 95 111 verbose=verbose) 96 97 112 98 113 # Initialise underlying mesh structure … … 108 123 if verbose: print 'Initialising Domain' 109 124 110 # List of quantity names entering 111 # the conservation equations 125 # List of quantity names entering the conservation equations 112 126 if conserved_quantities is None: 113 127 self.conserved_quantities = [] … … 121 135 self.other_quantities = other_quantities 122 136 123 124 137 # Build dictionary of Quantity instances keyed by quantity names 125 138 self.quantities = {} … … 138 151 self.full_send_dict = {} 139 152 else: 140 self.full_send_dict 153 self.full_send_dict = full_send_dict 141 154 142 155 # List of other quantity names … … 149 162 self.numproc = numproc 150 163 151 152 164 # Setup Communication Buffers 153 165 if verbose: print 'Domain: Set up communication buffers (parallel)' … … 155 167 for key in self.full_send_dict: 156 168 buffer_shape = self.full_send_dict[key][0].shape[0] 157 self.full_send_dict[key].append(num.zeros((buffer_shape, self.nsys), num.Float)) 169 self.full_send_dict[key].append(num.zeros((buffer_shape, self.nsys), 170 num.Float)) 158 171 159 172 for key in self.ghost_recv_dict: 160 173 buffer_shape = self.ghost_recv_dict[key][0].shape[0] 161 self.ghost_recv_dict[key].append(num.zeros((buffer_shape, self.nsys), num.Float)) 174 self.ghost_recv_dict[key].append(num.zeros((buffer_shape, self.nsys), 175 num.Float)) 162 176 163 177 # Setup cell full flag … … 173 187 # Test the assumption that all full triangles are store before 174 188 # the ghost triangles. 175 if not num.allclose(self.tri_full_flag[:self.number_of_full_nodes], 1):176 print 'WARNING: Not all full triangles are store before ghost triangles'177 189 if not num.allclose(self.tri_full_flag[:self.number_of_full_nodes], 1): 190 print ('WARNING: ' 191 'Not all full triangles are store before ghost triangles') 178 192 179 193 # Defaults … … 182 196 from anuga.config import timestepping_method 183 197 from anuga.config import protect_against_isolated_degenerate_timesteps 184 from anuga.config import default_order 198 from anuga.config import default_order 199 185 200 self.beta_w = beta_w 186 201 self.epsilon = epsilon 187 self.protect_against_isolated_degenerate_timesteps = protect_against_isolated_degenerate_timesteps 188 202 self.protect_against_isolated_degenerate_timesteps = \ 203 protect_against_isolated_degenerate_timesteps 204 189 205 self.set_default_order(default_order) 190 206 … … 196 212 self.set_timestepping_method(timestepping_method) 197 213 self.set_beta(beta_w) 198 self.boundary_map = None # Will be populated by set_boundary 199 214 self.boundary_map = None # Will be populated by set_boundary 200 215 201 216 # Model time … … 209 224 210 225 self.last_walltime = walltime() 211 226 212 227 # Monitoring 213 228 self.quantities_to_be_monitored = None 214 229 self.monitor_polygon = None 215 self.monitor_time_interval = None 230 self.monitor_time_interval = None 216 231 self.monitor_indices = None 217 232 218 233 # Checkpointing and storage 219 234 from anuga.config import default_datadir 235 220 236 self.datadir = default_datadir 221 237 self.simulation_name = 'domain' 222 238 self.checkpoint = False 223 239 224 # To avoid calculating the flux across each edge twice, 225 # keep an integer (boolean) array, to be used during the flux 226 # calculation 240 # To avoid calculating the flux across each edge twice, keep an integer 241 # (boolean) array, to be used during the flux calculation. 227 242 N = len(self) # Number_of_triangles 228 243 self.already_computed_flux = num.zeros((N, 3), num.Int) 229 244 230 245 # Storage for maximal speeds computed for each triangle by 231 # compute_fluxes 246 # compute_fluxes. 232 247 # This is used for diagnostics only (reset at every yieldstep) 233 248 self.max_speed = num.zeros(N, num.Float) 234 249 235 250 if mesh_filename is not None: 236 # If the mesh file passed any quantity values 237 # ,initialise with these values.251 # If the mesh file passed any quantity values, 252 # initialise with these values. 238 253 if verbose: print 'Domain: Initialising quantity values' 239 254 self.set_quantity_vertices_dict(vertex_quantity_dict) 240 255 241 242 256 if verbose: print 'Domain: Done' 243 257 244 245 246 247 248 249 #--------------------------- 250 # Public interface to Domain 251 #--------------------------- 252 def get_conserved_quantities(self, vol_id, vertex=None, edge=None): 258 ## 259 # @brief Get conserved quantities for a volume. 260 # @param vol_id ID of the volume we want the conserved quantities for. 261 # @param vertex If specified, use as index for edge values. 262 # @param edge If specified, use as index for edge values. 263 # @return Vector of conserved quantities. 264 # @note If neither 'vertex' or 'edge' specified, use centroid values. 265 # @note If both 'vertex' and 'edge' specified, raise exception. 266 def get_conserved_quantities(self, vol_id, 267 vertex=None, 268 edge=None): 253 269 """Get conserved quantities at volume vol_id 254 270 … … 259 275 260 276 Return value: Vector of length == number_of_conserved quantities 261 262 277 """ 263 278 … … 265 280 msg = 'Values for both vertex and edge was specified.' 266 281 msg += 'Only one (or none) is allowed.' 267 raise msg282 raise Exception, msg 268 283 269 284 q = num.zeros(len(self.conserved_quantities), num.Float) … … 280 295 return q 281 296 297 ## 298 # @brief Set the relative model time. 299 # @param time The new model time (seconds). 282 300 def set_time(self, time=0.0): 283 301 """Set the model time (seconds)""" … … 287 305 self.time = time 288 306 307 ## 308 # @brief Get the model time. 309 # @return The absolute model time (seconds). 289 310 def get_time(self): 290 311 """Get the absolute model time (seconds)""" … … 292 313 return self.time + self.starttime 293 314 294 def set_beta(self,beta): 295 """Set default beta for limiting 296 """ 315 ## 316 # @brief Set the default beta for limiting. 317 # @param beta The new beta value. 318 def set_beta(self, beta): 319 """Set default beta for limiting""" 297 320 298 321 self.beta = beta … … 302 325 Q.set_beta(beta) 303 326 327 ## 328 # @brief Get the beta value used for limiting. 329 # @return The beta value used for limiting. 304 330 def get_beta(self): 305 """Get default beta for limiting 306 """ 331 """Get default beta for limiting""" 307 332 308 333 return self.beta 309 334 335 ## 336 # @brief Set default (spatial) order. 337 # @param n The new spatial order value. 338 # @note If 'n' is not 1 or 2, raise exception. 310 339 def set_default_order(self, n): 311 """Set default (spatial) order to either 1 or 2 312 """ 313 314 msg = 'Default order must be either 1 or 2. I got %s' %n 340 """Set default (spatial) order to either 1 or 2""" 341 342 msg = 'Default order must be either 1 or 2. I got %s' % n 315 343 assert n in [1,2], msg 316 344 … … 327 355 #self.set_timestepping_method('rk2') 328 356 #self.set_all_limiters(beta_rk2) 329 330 357 358 ## 359 # @brief Set values of named quantities. 360 # @param quantity_dict Dictionary containing name/value pairs. 331 361 def set_quantity_vertices_dict(self, quantity_dict): 332 362 """Set values for named quantities. 333 The index is the quantity 334 335 name: Name of quantity 336 X: Compatible list, Numeric array, const or function (see below) 337 338 The values will be stored in elements following their 339 internal ordering. 340 341 """ 342 363 Supplied dictionary contains name/value pairs: 364 365 name: Name of quantity 366 value: Compatible list, Numeric array, const or function (see below) 367 368 The values will be stored in elements following their internal ordering. 369 """ 370 343 371 # FIXME: Could we name this a bit more intuitively 344 372 # E.g. set_quantities_from_dictionary … … 346 374 self.set_quantity(key, quantity_dict[key], location='vertices') 347 375 348 349 def set_quantity(self, name, *args, **kwargs): 376 ## 377 # @brief Set value(s) for a named quantity. 378 # @param name Name of quantity to be updated. 379 # @param args Positional args. 380 # @param kwargs Keyword args. 381 # @note If 'kwargs' dict has 'expression' key, evaluate expression. 382 def set_quantity(self, name, 383 *args, **kwargs): 350 384 """Set values for named quantity 351 352 385 353 386 One keyword argument is documented here: … … 370 403 # Assign values 371 404 self.quantities[name].set_values(*args, **kwargs) 372 373 def add_quantity(self, name, *args, **kwargs): 405 406 ## 407 # @brief Add to a named quantity value. 408 # @param name Name of quantity to be added to. 409 # @param args Positional args. 410 # @param kwargs Keyword args. 411 # @note If 'kwargs' dict has 'expression' key, evaluate expression. 412 def add_quantity(self, name, 413 *args, **kwargs): 374 414 """Add values to a named quantity 375 415 376 416 E.g add_quantity('elevation', X) 377 417 378 418 Option are the same as in set_quantity. 379 419 """ 380 420 381 421 # Do the expression stuff 382 422 if kwargs.has_key('expression'): 383 423 expression = kwargs['expression'] 384 424 Q2 = self.create_quantity_from_expression(expression) 385 else: 425 else: 386 426 # Create new temporary quantity 387 427 Q2 = Quantity(self) 388 428 389 429 # Assign specified values to temporary quantity 390 430 Q2.set_values(*args, **kwargs) 391 431 392 432 # Add temporary quantity to named quantity 393 433 Q1 = self.get_quantity(name) 394 434 self.set_quantity(name, Q1 + Q2) 395 396 435 436 ## 437 # @brief Get list of quantity names for the Domain. 438 # @return List of quantity names. 397 439 def get_quantity_names(self): 398 440 """Get a list of all the quantity names that this domain is aware of. 399 441 Any value in the result should be a valid input to get_quantity. 400 442 """ 443 401 444 return self.quantities.keys() 402 445 403 def get_quantity(self, name, location='vertices', indices = None): 446 ## 447 # @brief Get a quantity object. 448 # @param name Name of the quantity value. 449 # @param location ?? 450 # @param indices ?? 451 # @return The quantity value object. 452 # @note 'location' and 'indices' are unused. 453 def get_quantity(self, name, 454 location='vertices', 455 indices = None): 404 456 """Get pointer to quantity object. 405 457 … … 413 465 return self.quantities[name] #.get_values( location, indices = indices) 414 466 415 416 467 ## 468 # @brief Create a quantity value from an expression. 469 # @param expression The expression (string) to be evaluated. 470 # @return The expression value, evaluated from this Domain's quantities. 471 # @note Valid expression operators are as defined in class Quantity. 417 472 def create_quantity_from_expression(self, expression): 418 """Create new quantity from other quantities using arbitrary expression 473 """Create new quantity from other quantities using arbitrary expression. 419 474 420 475 Combine existing quantities in domain using expression and return … … 426 481 427 482 Examples creating derived quantities: 428 429 483 Depth = domain.create_quantity_from_expression('stage-elevation') 430 431 exp = '(xmomentum*xmomentum + ymomentum*ymomentum)**0.5' 484 exp = '(xmomentum*xmomentum + ymomentum*ymomentum)**0.5' 432 485 Absolute_momentum = domain.create_quantity_from_expression(exp) 433 434 486 """ 435 487 436 488 from anuga.abstract_2d_finite_volumes.util import\ 437 489 apply_expression_to_dictionary 438 490 439 491 return apply_expression_to_dictionary(expression, self.quantities) 440 492 441 442 493 ## 494 # @brief Associate boundary objects with tagged boundary segments. 495 # @param boundary_map A dict of boundary objects keyed by symbolic tags to 496 # matched against tags in the internal dictionary 497 # self.boundary. 443 498 def set_boundary(self, boundary_map): 444 499 """Associate boundary objects with tagged boundary segments. … … 461 516 self.boundary_objects: ((vol_id, edge_id), boundary_object) 462 517 463 464 518 Pre-condition: 465 519 self.boundary has been built. … … 471 525 exception is raised. 472 526 However, if a tag is not used to the domain, no error is thrown. 473 FIXME: This would lead to implementation of a 474 default boundary condition 527 FIXME: This would lead to implementation of a default boundary condition 475 528 476 529 Note: If a segment is listed in the boundary dictionary and if it is 477 not None, it *will* become a boundary - 478 even if there is a neighbouring triangle. 479 This would be the case for internal boundaries 530 not None, it *will* become a boundary - even if there is a neighbouring 531 triangle. This would be the case for internal boundaries. 480 532 481 533 Boundary objects that are None will be skipped. 482 534 483 If a boundary_map has already been set 484 (i.e. set_boundary has been called before), the old boundary map485 will be updated with new values. The new map need not define all486 boundary tags, and can thus change onlythose that are needed.535 If a boundary_map has already been set (i.e. set_boundary has been 536 called before), the old boundary map will be updated with new values. 537 The new map need not define all boundary tags, and can thus change only 538 those that are needed. 487 539 488 540 FIXME: If set_boundary is called multiple times and if Boundary 489 541 object is changed into None, the neighbour structure will not be 490 542 restored!!! 491 492 493 543 """ 494 544 … … 496 546 # This the first call to set_boundary. Store 497 547 # map for later updates and for use with boundary_stats. 498 self.boundary_map = boundary_map 499 else: 548 self.boundary_map = boundary_map 549 else: 500 550 # This is a modification of an already existing map 501 551 # Update map an proceed normally 502 503 552 for key in boundary_map.keys(): 504 553 self.boundary_map[key] = boundary_map[key] 505 506 554 507 555 # FIXME (Ole): Try to remove the sorting and fix test_mesh.py 508 556 x = self.boundary.keys() … … 511 559 # Loop through edges that lie on the boundary and associate them with 512 560 # callable boundary objects depending on their tags 513 self.boundary_objects = [] 561 self.boundary_objects = [] 514 562 for k, (vol_id, edge_id) in enumerate(x): 515 tag = self.boundary[ (vol_id, edge_id)]563 tag = self.boundary[(vol_id, edge_id)] 516 564 517 565 if self.boundary_map.has_key(tag): … … 519 567 520 568 if B is not None: 521 self.boundary_objects.append( ((vol_id, edge_id), B))569 self.boundary_objects.append(((vol_id, edge_id), B)) 522 570 self.neighbours[vol_id, edge_id] = \ 523 571 -len(self.boundary_objects) 524 572 else: 525 573 pass 526 574 #FIXME: Check and perhaps fix neighbour structure 527 528 529 575 else: 530 576 msg = 'ERROR (domain.py): Tag "%s" has not been ' %tag … … 533 579 msg += 'in set_boundary.\n' 534 580 msg += 'The tags are: %s' %self.get_boundary_tags() 535 raise msg 536 537 581 raise Exception, msg 582 583 ## 584 # @brief Set quantities based on a regional tag. 585 # @param args 586 # @param kwargs 538 587 def set_region(self, *args, **kwargs): 539 """ 540 This method is used to set quantities based on a regional tag. 541 588 """Set quantities based on a regional tag. 589 542 590 It is most often called with the following parameters; 543 591 (self, tag, quantity, X, location='vertices') 544 tag: the name of the regional tag used to specify the region592 tag: the name of the regional tag used to specify the region 545 593 quantity: Name of quantity to change 546 X: const or function - how the quantity is changed594 X: const or function - how the quantity is changed 547 595 location: Where values are to be stored. 548 596 Permissible options are: vertices, centroid and unique vertices … … 558 606 func = region_set_region(*args, **kwargs) 559 607 self._set_region(func) 560 561 608 609 ## 610 # @brief ?? 611 # @param functions A list or tuple of ?? 562 612 def _set_region(self, functions): 613 # coerce to an iterable (list or tuple) 614 if type(functions) not in [types.ListType, types.TupleType]: 615 functions = [functions] 616 563 617 # The order of functions in the list is used. 564 if type(functions) not in [types.ListType,types.TupleType]:565 functions = [functions]566 618 for function in functions: 567 619 for tag in self.tagged_elements.keys(): 568 620 function(tag, self.tagged_elements[tag], self) 569 621 570 571 572 622 ## 623 # @brief Specify the quantities which will be monitored for extrema. 624 # @param q Single or list of quantity names to monitor. 625 # @param polygon If specified, monitor only triangles inside polygon. 626 # @param time_interval If specified, monitor only timesteps inside interval. 627 # @note If 'q' is None, do no monitoring. 573 628 def set_quantities_to_be_monitored(self, q, 574 polygon=None,575 time_interval=None):629 polygon=None, 630 time_interval=None): 576 631 """Specify which quantities will be monitored for extrema. 577 632 578 633 q must be either: 579 - the name of a quantity or aderived quantity such as 'stage-elevation'634 - the name of a quantity or derived quantity such as 'stage-elevation' 580 635 - a list of quantity names 581 636 - None … … 583 638 In the two first cases, the named quantities will be monitored at 584 639 each internal timestep 585 640 586 641 If q is None, monitoring will be switched off altogether. 587 642 588 polygon (if specified) will restrict monitoring totriangles inside polygon.643 polygon (if specified) will only monitor triangles inside polygon. 589 644 If omitted all triangles will be included. 590 645 591 time_interval (if specified)will restrict monitoring to time steps in646 time_interval, if specified, will restrict monitoring to time steps in 592 647 that interval. If omitted all timesteps will be included. 593 648 """ 594 649 595 650 from anuga.abstract_2d_finite_volumes.util import\ 596 apply_expression_to_dictionary 651 apply_expression_to_dictionary 597 652 598 653 if q is None: … … 600 655 self.monitor_polygon = None 601 656 self.monitor_time_interval = None 602 self.monitor_indices = None 657 self.monitor_indices = None 603 658 return 604 659 660 # coerce 'q' to a list if it's a string 605 661 if isinstance(q, basestring): 606 q = [q] # Turn argument into a list607 608 # Check correc ness and initialise662 q = [q] 663 664 # Check correctness and initialise 609 665 self.quantities_to_be_monitored = {} 610 666 for quantity_name in q: 611 msg = 'Quantity %s is not a valid conserved quantity'\ 612 %quantity_name 613 667 msg = 'Quantity %s is not a valid conserved quantity' \ 668 % quantity_name 614 669 615 670 if not quantity_name in self.quantities: … … 622 677 'min_location': None, # Argmin (x, y) 623 678 'max_location': None, # Argmax (x, y) 624 'min_time': None, # Argmin (t) 679 'min_time': None, # Argmin (t) 625 680 'max_time': None} # Argmax (t) 626 681 627 682 self.quantities_to_be_monitored[quantity_name] = info_block 628 629 630 683 631 684 if polygon is not None: 632 685 # Check input 633 686 if isinstance(polygon, basestring): 634 635 687 # Check if multiple quantities were accidentally 636 688 # given as separate argument rather than a list. 637 msg = 'Multiple quantities must be specified in a list. '638 msg +='Not as multiple arguments. '639 msg += 'I got "%s" as a second argument' %polygon640 689 msg = ('Multiple quantities must be specified in a list. ' 690 'Not as multiple arguments. ' 691 'I got "%s" as a second argument') % polygon 692 641 693 if polygon in self.quantities: 642 694 raise Exception, msg 643 695 644 696 try: 645 697 apply_expression_to_dictionary(polygon, self.quantities) 646 698 except: 647 # At least polygon wasn't anexpression involving quantitites699 # At least polygon wasn't expression involving quantitites 648 700 pass 649 701 else: … … 651 703 652 704 # In any case, we don't allow polygon to be a string 653 msg = 'argument "polygon" must not be a string: '654 msg += 'I got polygon=\'%s\' ' %polygon705 msg = ('argument "polygon" must not be a string: ' 706 'I got polygon="%s"') % polygon 655 707 raise Exception, msg 656 657 708 658 709 # Get indices for centroids that are inside polygon 659 710 points = self.get_centroid_coordinates(absolute=True) 660 711 self.monitor_indices = inside_polygon(points, polygon) 661 662 712 663 713 if time_interval is not None: 664 714 assert len(time_interval) == 2 665 715 666 667 716 self.monitor_polygon = polygon 668 self.monitor_time_interval = time_interval 669 670 671 #-------------------------- 672 # Miscellaneous diagnostics 673 #-------------------------- 717 self.monitor_time_interval = time_interval 718 719 ## 720 # @brief Check Domain integrity. 721 # @note Raises an exception if integrity breached. 674 722 def check_integrity(self): 675 723 Mesh.check_integrity(self) … … 680 728 681 729 ##assert hasattr(self, 'boundary_objects') 682 683 730 731 ## 732 # @brief Print timestep stats to stdout. 733 # @param track_speeds If True, print smallest track speed. 684 734 def write_time(self, track_speeds=False): 685 735 print self.timestepping_statistics(track_speeds) 686 736 687 688 def timestepping_statistics(self, 689 track_speeds=False, 690 triangle_id=None): 691 """Return string with time stepping statistics for printing or logging 737 ## 738 # @brief Get timestepping stats string. 739 # @param track_speeds If True, report location of smallest timestep. 740 # @param triangle_id If specified, use specific triangle. 741 # @return A string containing timestep stats. 742 def timestepping_statistics(self, track_speeds=False, 743 triangle_id=None): 744 """Return string with time stepping statistics 692 745 693 746 Optional boolean keyword track_speeds decides whether to report … … 696 749 697 750 Optional keyword triangle_id can be used to specify a particular 698 triangle rather than the one with the largest speed. 751 triangle rather than the one with the largest speed. 699 752 """ 700 753 701 754 from anuga.utilities.numerical_tools import histogram, create_bins 702 703 704 # qwidth determines the text field used for quantities 755 756 # qwidth determines the the width of the text field used for quantities 705 757 qwidth = self.qwidth = 12 706 707 758 708 759 msg = '' … … 723 774 model_time = self.get_time() 724 775 if self.min_timestep == self.max_timestep: 725 msg += 'Time = %.4f, delta t = %.8f, steps=%d' \726 %(model_time, self.min_timestep, self.number_of_steps)776 msg += 'Time = %.4f, delta t = %.8f, steps=%d' \ 777 % (model_time, self.min_timestep, self.number_of_steps) 727 778 elif self.min_timestep > self.max_timestep: 728 msg += 'Time = %.4f, steps=%d' \729 %(model_time, self.number_of_steps)779 msg += 'Time = %.4f, steps=%d' \ 780 % (model_time, self.number_of_steps) 730 781 else: 731 msg += 'Time = %.4f, delta t in [%.8f, %.8f], steps=%d' \732 %(model_time, self.min_timestep,733 self.max_timestep, self.number_of_steps)734 735 msg += ' (%ds)' % (walltime() - self.last_walltime)736 self.last_walltime = walltime() 737 782 msg += 'Time = %.4f, delta t in [%.8f, %.8f], steps=%d' \ 783 % (model_time, self.min_timestep, 784 self.max_timestep, self.number_of_steps) 785 786 msg += ' (%ds)' % (walltime() - self.last_walltime) 787 self.last_walltime = walltime() 788 738 789 if track_speeds is True: 739 790 msg += '\n' 740 741 791 742 792 # Setup 10 bins for speed histogram … … 745 795 746 796 msg += '------------------------------------------------\n' 747 msg += ' Speeds in [%f, %f]\n' % (min(self.max_speed),748 max(self.max_speed))797 msg += ' Speeds in [%f, %f]\n' % (min(self.max_speed), 798 max(self.max_speed)) 749 799 msg += ' Histogram:\n' 750 800 … … 753 803 lo = hi 754 804 if i+1 < len(bins): 755 # Open upper interval 805 # Open upper interval 756 806 hi = bins[i+1] 757 msg += ' [%f, %f[: %d\n' % (lo, hi, count)807 msg += ' [%f, %f[: %d\n' % (lo, hi, count) 758 808 else: 759 809 # Closed upper interval 760 810 hi = max(self.max_speed) 761 msg += ' [%f, %f]: %d\n' %(lo, hi, count) 762 811 msg += ' [%f, %f]: %d\n' % (lo, hi, count) 763 812 764 813 N = len(self.max_speed) … … 770 819 k = 0 771 820 lower = min(speed) 772 for i, a in enumerate(speed): 821 for i, a in enumerate(speed): 773 822 if i % (N/10) == 0 and i != 0: 774 823 # For every 10% of the sorted speeds 775 msg += ' %d speeds in [%f, %f]\n' % (i-k, lower, a)824 msg += ' %d speeds in [%f, %f]\n' % (i-k, lower, a) 776 825 lower = a 777 826 k = i 778 827 779 828 msg += ' %d speeds in [%f, %f]\n'\ 780 %(N-k, lower, max(speed)) 781 782 783 784 785 829 % (N-k, lower, max(speed)) 830 786 831 # Find index of largest computed flux speed 787 832 if triangle_id is None: 788 833 k = self.k = num.argmax(self.max_speed) 789 834 else: 790 errmsg = 'Triangle_id %d does not exist in mesh: %s' %(triangle_id,791 835 errmsg = 'Triangle_id %d does not exist in mesh: %s' \ 836 % (triangle_id, str(self)) 792 837 assert 0 <= triangle_id < len(self), errmsg 793 838 k = self.k = triangle_id 794 795 839 796 840 x, y = self.get_centroid_coordinates()[k] 797 841 radius = self.get_radii()[k] 798 area = self.get_areas()[k] 799 max_speed = self.max_speed[k] 800 801 msg += ' Triangle #%d with centroid (%.4f, %.4f), ' % (k, x, y)802 msg += 'area = %.4f and radius = %.4f ' % (area, radius)842 area = self.get_areas()[k] 843 max_speed = self.max_speed[k] 844 845 msg += ' Triangle #%d with centroid (%.4f, %.4f), ' % (k, x, y) 846 msg += 'area = %.4f and radius = %.4f ' % (area, radius) 803 847 if triangle_id is None: 804 msg += 'had the largest computed speed: %.6f m/s ' % (max_speed)848 msg += 'had the largest computed speed: %.6f m/s ' % (max_speed) 805 849 else: 806 msg += 'had computed speed: %.6f m/s ' % (max_speed)807 850 msg += 'had computed speed: %.6f m/s ' % (max_speed) 851 808 852 if max_speed > 0.0: 809 msg += '(timestep=%.6f)\n' % (radius/max_speed)853 msg += '(timestep=%.6f)\n' % (radius/max_speed) 810 854 else: 811 msg += '(timestep=%.6f)\n' % (0)812 855 msg += '(timestep=%.6f)\n' % (0) 856 813 857 # Report all quantity values at vertices, edges and centroid 814 858 msg += ' Quantity' … … 819 863 V = q.get_values(location='vertices', indices=[k])[0] 820 864 E = q.get_values(location='edges', indices=[k])[0] 821 C = q.get_values(location='centroids', indices=[k]) 822 823 s = ' %s: vertex_values = %.4f,\t %.4f,\t %.4f\n' \824 %(name.ljust(qwidth), V[0], V[1], V[2])825 826 s += ' %s: edge_values = %.4f,\t %.4f,\t %.4f\n' \827 %(name.ljust(qwidth), E[0], E[1], E[2])828 829 s += ' %s: centroid_value = %.4f\n' \830 %(name.ljust(qwidth), C[0])865 C = q.get_values(location='centroids', indices=[k]) 866 867 s = ' %s: vertex_values = %.4f,\t %.4f,\t %.4f\n' \ 868 % (name.ljust(qwidth), V[0], V[1], V[2]) 869 870 s += ' %s: edge_values = %.4f,\t %.4f,\t %.4f\n' \ 871 % (name.ljust(qwidth), E[0], E[1], E[2]) 872 873 s += ' %s: centroid_value = %.4f\n' \ 874 % (name.ljust(qwidth), C[0]) 831 875 832 876 msg += s … … 834 878 return msg 835 879 836 837 def write_boundary_statistics(self, quantities = None, tags = None): 880 ## 881 # @brief Print boundary forcing stats at each timestep to stdout. 882 # @param quantities A name or list of names of quantities to report. 883 # @param tags A name or list of names of tags to report. 884 def write_boundary_statistics(self, quantities=None, tags=None): 838 885 print self.boundary_statistics(quantities, tags) 839 886 840 def boundary_statistics(self, quantities = None, tags = None): 887 # @brief Get a string containing boundary forcing stats at each timestep. 888 # @param quantities A name or list of names of quantities to report. 889 # @param tags A name or list of names of tags to report. 890 # @note If 'quantities' is None, report all. Same for 'tags'. 891 def boundary_statistics(self, quantities=None, 892 tags=None): 841 893 """Output statistics about boundary forcing at each timestep 842 843 894 844 895 Input: … … 847 898 tags: either None, a string or a list of strings naming the 848 899 tags to be reported 849 850 900 851 901 Example output: … … 856 906 Tag 'ocean' 857 907 858 859 908 If quantities are specified only report on those. Otherwise take all 860 909 conserved quantities. 861 910 If tags are specified only report on those, otherwise take all tags. 862 863 """ 911 """ 912 913 import types, string 864 914 865 915 # Input checks 866 import types, string867 868 916 if quantities is None: 869 917 quantities = self.conserved_quantities … … 871 919 quantities = [quantities] #Turn it into a list 872 920 873 msg = 'Keyword argument quantities must be either None, '874 msg += 'string or list. I got %s' %str(quantities)921 msg = ('Keyword argument quantities must be either None, ' 922 'string or list. I got %s') % str(quantities) 875 923 assert type(quantities) == types.ListType, msg 876 877 924 878 925 if tags is None: … … 881 928 tags = [tags] #Turn it into a list 882 929 883 msg = 'Keyword argument tags must be either None, '884 msg += 'string or list. I got %s' %str(tags)930 msg = ('Keyword argument tags must be either None, ' 931 'string or list. I got %s') % str(tags) 885 932 assert type(tags) == types.ListType, msg 886 933 … … 893 940 894 941 # Output statistics 895 msg = 'Boundary values at time %.4f:\n' % self.get_time()942 msg = 'Boundary values at time %.4f:\n' % self.get_time() 896 943 for tag in tags: 897 msg += ' %s:\n' % tag944 msg += ' %s:\n' % tag 898 945 899 946 for name in quantities: … … 902 949 # Find range of boundary values for tag and q 903 950 maxval = minval = None 904 for i, ((vol_id, edge_id), B) in\ 905 enumerate(self.boundary_objects): 951 for i, ((vol_id,edge_id),B) in enumerate(self.boundary_objects): 906 952 if self.boundary[(vol_id, edge_id)] == tag: 907 953 v = q.boundary_values[i] … … 910 956 911 957 if minval is None or maxval is None: 912 msg += ' Sorry no information available about' +\913 ' tag %s and quantity %s\n' %(tag, name)958 msg += (' Sorry no information available about' 959 ' tag %s and quantity %s\n') % (tag, name) 914 960 else: 915 msg += ' %s in [%12.8f, %12.8f]\n'\ 916 %(string.ljust(name, maxwidth), minval, maxval) 917 961 msg += ' %s in [%12.8f, %12.8f]\n' \ 962 % (string.ljust(name, maxwidth), minval, maxval) 918 963 919 964 return msg 920 965 921 966 ## 967 # @brief Update extrema if requested by set_quantities_to_be_monitored. 922 968 def update_extrema(self): 923 969 """Update extrema if requested by set_quantities_to_be_monitored. … … 929 975 # Define a tolerance for extremum computations 930 976 from anuga.config import single_precision as epsilon 931 977 932 978 if self.quantities_to_be_monitored is None: 933 979 return … … 954 1000 # of the epsilon) 955 1001 maxval = Q.get_maximum_value(self.monitor_indices) 956 if info_block['max'] is None or \1002 if info_block['max'] is None or \ 957 1003 maxval > info_block['max'] + epsilon: 958 1004 info_block['max'] = maxval … … 961 1007 info_block['max_time'] = self.time 962 1008 963 964 1009 # Update minimum 965 1010 minval = Q.get_minimum_value(self.monitor_indices) 966 if info_block['min'] is None or \1011 if info_block['min'] is None or \ 967 1012 minval < info_block['min'] - epsilon: 968 info_block['min'] = minval 1013 info_block['min'] = minval 969 1014 minloc = Q.get_minimum_location() 970 1015 info_block['min_location'] = minloc 971 info_block['min_time'] = self.time 972 973 974 975 def quantity_statistics(self, precision = '%.4f'): 1016 info_block['min_time'] = self.time 1017 1018 ## 1019 # @brief Return string with statistics about quantities 1020 # @param precision A format string to use for float values. 1021 # @return The stats string. 1022 def quantity_statistics(self, precision='%.4f'): 976 1023 """Return string with statistics about quantities for 977 1024 printing or logging … … 980 1027 981 1028 set_quantities_to_be_monitored 982 983 1029 """ 984 1030 … … 986 1032 987 1033 # Output statistics 988 msg = 'Monitored quantities at time %.4f:\n' % self.get_time()1034 msg = 'Monitored quantities at time %.4f:\n' % self.get_time() 989 1035 if self.monitor_polygon is not None: 990 1036 p_str = str(self.monitor_polygon) 991 msg += '- Restricted by polygon: %s' % p_str[:maxlen]1037 msg += '- Restricted by polygon: %s' % p_str[:maxlen] 992 1038 if len(p_str) >= maxlen: 993 1039 msg += '...\n' … … 995 1041 msg += '\n' 996 1042 997 998 1043 if self.monitor_time_interval is not None: 999 msg += '- Restricted by time interval: %s\n' \1000 %str(self.monitor_time_interval)1044 msg += '- Restricted by time interval: %s\n' \ 1045 % str(self.monitor_time_interval) 1001 1046 time_interval_start = self.monitor_time_interval[0] 1002 1047 else: 1003 1048 time_interval_start = 0.0 1004 1049 1005 1006 1050 for quantity_name, info in self.quantities_to_be_monitored.items(): 1007 msg += ' %s:\n' %quantity_name 1008 1009 msg += ' values since time = %.2f in [%s, %s]\n'\ 1010 %(time_interval_start, 1011 get_textual_float(info['min'], precision), 1012 get_textual_float(info['max'], precision)) 1013 1014 msg += ' minimum attained at time = %s, location = %s\n'\ 1015 %(get_textual_float(info['min_time'], precision), 1016 get_textual_float(info['min_location'], precision)) 1017 1018 1019 msg += ' maximum attained at time = %s, location = %s\n'\ 1020 %(get_textual_float(info['max_time'], precision), 1021 get_textual_float(info['max_location'], precision)) 1022 1051 msg += ' %s:\n' % quantity_name 1052 1053 msg += ' values since time = %.2f in [%s, %s]\n' \ 1054 % (time_interval_start, 1055 get_textual_float(info['min'], precision), 1056 get_textual_float(info['max'], precision)) 1057 1058 msg += ' minimum attained at time = %s, location = %s\n' \ 1059 % (get_textual_float(info['min_time'], precision), 1060 get_textual_float(info['min_location'], precision)) 1061 1062 msg += ' maximum attained at time = %s, location = %s\n' \ 1063 % (get_textual_float(info['max_time'], precision), 1064 get_textual_float(info['max_location'], precision)) 1023 1065 1024 1066 return msg 1025 1067 1068 ## 1069 # @brief Get the timestep method. 1070 # @return The timestep method. One of 'euler', 'rk2' or 'rk3'. 1026 1071 def get_timestepping_method(self): 1027 1072 return self.timestepping_method 1028 1073 1029 def set_timestepping_method(self,timestepping_method): 1030 1074 ## 1075 # @brief Set the tmestep method to be used. 1076 # @param timestepping_method One of 'euler', 'rk2' or 'rk3'. 1077 # @note Raises exception of method not known. 1078 def set_timestepping_method(self, timestepping_method): 1031 1079 if timestepping_method in ['euler', 'rk2', 'rk3']: 1032 1080 self.timestepping_method = timestepping_method 1033 1081 return 1034 1082 1035 msg = '%s is an incorrect timestepping type' % timestepping_method1083 msg = '%s is an incorrect timestepping type' % timestepping_method 1036 1084 raise Exception, msg 1037 1085 1086 ## 1087 # @brief Get the Domain simulation name. 1088 # @return The simulation name string. 1038 1089 def get_name(self): 1039 1090 return self.simulation_name 1040 1091 1092 ## 1093 # @brief Set the simulation name. 1094 # @param name The name of the simulation. 1095 # @note The simulation name is also used for the output .sww file. 1041 1096 def set_name(self, name): 1042 1097 """Assign a name to this simulation. 1043 1098 This will be used to identify the output sww file. 1044 1045 """ 1099 """ 1100 1101 # remove any '.sww' end 1046 1102 if name.endswith('.sww'): 1047 1103 name = name[:-4] 1048 1104 1049 1105 self.simulation_name = name 1050 1106 1107 ## 1108 # @brief Get data directory path. 1109 # @return The data directory path string. 1051 1110 def get_datadir(self): 1052 1111 return self.datadir 1053 1112 1113 ## 1114 # @brief Set data directory path. 1115 # @param name The data directory path string. 1054 1116 def set_datadir(self, name): 1055 1117 self.datadir = name 1056 1118 1119 ## 1120 # @brief Get the start time value. 1121 # @return The start time value (float). 1057 1122 def get_starttime(self): 1058 1123 return self.starttime 1059 1124 1125 ## 1126 # @brief Set the start time value. 1127 # @param time The start time value. 1060 1128 def set_starttime(self, time): 1061 self.starttime = float(time) 1062 1063 1064 1065 #-------------------------- 1066 # Main components of evolve 1067 #-------------------------- 1068 1069 def evolve(self, 1070 yieldstep = None, 1071 finaltime = None, 1072 duration = None, 1073 skip_initial_step = False): 1129 self.starttime = float(time) 1130 1131 #-------------------------- 1132 # Main components of evolve 1133 #-------------------------- 1134 1135 ## 1136 # @brief Evolve the model through time. 1137 # @param yieldstep Interval between yields where results are stored, etc. 1138 # @param finaltime Time where simulation should end. 1139 # @param duration Duration of simulation. 1140 # @param skip_initial_step If True, skip the first yield step. 1141 def evolve(self, yieldstep=None, 1142 finaltime=None, 1143 duration=None, 1144 skip_initial_step=False): 1074 1145 """Evolve model through time starting from self.starttime. 1075 1076 1146 1077 1147 yieldstep: Interval between yields where results are stored, … … 1088 1158 If both duration and finaltime are given an exception is thrown. 1089 1159 1090 1091 1160 skip_initial_step: Boolean flag that decides whether the first 1092 1161 yield step is skipped or not. This is useful for example to avoid 1093 1162 duplicate steps when multiple evolve processes are dove tailed. 1094 1163 1095 1096 1164 Evolve is implemented as a generator and is to be called as such, e.g. 1097 1165 … … 1099 1167 <Do something with domain and t> 1100 1168 1101 1102 1169 All times are given in seconds 1103 1104 1170 """ 1105 1171 … … 1107 1173 1108 1174 # FIXME: Maybe lump into a larger check prior to evolving 1109 msg = 'Boundary tags must be bound to boundary objects before '1110 msg +='evolving system, '1111 msg +='e.g. using the method set_boundary.\n'1112 msg += 'This system has the boundary tags %s '\1113 %self.get_boundary_tags()1175 msg = ('Boundary tags must be bound to boundary objects before ' 1176 'evolving system, ' 1177 'e.g. using the method set_boundary.\n' 1178 'This system has the boundary tags %s ') \ 1179 % self.get_boundary_tags() 1114 1180 assert hasattr(self, 'boundary_objects'), msg 1115 1116 1181 1117 1182 if yieldstep is None: … … 1122 1187 self._order_ = self.default_order 1123 1188 1124 1125 1189 if finaltime is not None and duration is not None: 1126 1190 # print 'F', finaltime, duration 1127 1191 msg = 'Only one of finaltime and duration may be specified' 1128 raise msg1192 raise Exception, msg 1129 1193 else: 1130 1194 if finaltime is not None: … … 1133 1197 self.finaltime = self.starttime + float(duration) 1134 1198 1135 1136 1137 1199 N = len(self) # Number of triangles 1138 1200 self.yieldtime = 0.0 # Track time between 'yields' … … 1144 1206 self.number_of_first_order_steps = 0 1145 1207 1146 1147 1208 # Update ghosts 1148 1209 self.update_ghosts() … … 1153 1214 # Update extrema if necessary (for reporting) 1154 1215 self.update_extrema() 1155 1216 1156 1217 # Initial update boundary values 1157 1218 self.update_boundary() … … 1165 1226 1166 1227 while True: 1167 1168 1228 # Evolve One Step, using appropriate timestepping method 1169 1229 if self.get_timestepping_method() == 'euler': 1170 self.evolve_one_euler_step(yieldstep, finaltime)1171 1230 self.evolve_one_euler_step(yieldstep, finaltime) 1231 1172 1232 elif self.get_timestepping_method() == 'rk2': 1173 self.evolve_one_rk2_step(yieldstep, finaltime)1233 self.evolve_one_rk2_step(yieldstep, finaltime) 1174 1234 1175 1235 elif self.get_timestepping_method() == 'rk3': 1176 self.evolve_one_rk3_step(yieldstep, finaltime)1177 1236 self.evolve_one_rk3_step(yieldstep, finaltime) 1237 1178 1238 # Update extrema if necessary (for reporting) 1179 self.update_extrema() 1180 1239 self.update_extrema() 1181 1240 1182 1241 self.yieldtime += self.timestep … … 1187 1246 # Yield results 1188 1247 if finaltime is not None and self.time >= finaltime-epsilon: 1189 1190 1248 if self.time > finaltime: 1191 1249 # FIXME (Ole, 30 April 2006): Do we need this check? 1192 1250 # Probably not (Ole, 18 September 2008). Now changed to 1193 1251 # Exception 1194 msg = 'WARNING (domain.py): time overshot finaltime. '1195 msg += 'Contact Ole.Nielsen@ga.gov.au'1252 msg = ('WARNING (domain.py): time overshot finaltime. ' 1253 'Contact Ole.Nielsen@ga.gov.au') 1196 1254 raise Exception, msg 1197 1198 1255 1199 1256 # Yield final time and stop … … 1202 1259 break 1203 1260 1204 1205 1261 if self.yieldtime >= yieldstep: 1206 1262 # Yield (intermediate) time and allow inspection of domain 1207 1208 1263 if self.checkpoint is True: 1209 1264 self.store_checkpoint() … … 1221 1276 self.max_speed = num.zeros(N, num.Float) 1222 1277 1278 ## 1279 # @brief 'Euler' time step method. 1280 # @param yieldstep The reporting time step. 1281 # @param finaltime The simulation final time. 1223 1282 def evolve_one_euler_step(self, yieldstep, finaltime): 1224 """ 1225 One Euler Time Step 1283 """One Euler Time Step 1226 1284 Q^{n+1} = E(h) Q^n 1227 1285 """ … … 1248 1306 self.update_boundary() 1249 1307 1250 1251 1252 1253 1308 ## 1309 # @brief 'rk2' time step method. 1310 # @param yieldstep The reporting time step. 1311 # @param finaltime The simulation final time. 1254 1312 def evolve_one_rk2_step(self, yieldstep, finaltime): 1255 """ 1256 One 2nd order RK timestep 1313 """One 2nd order RK timestep 1257 1314 Q^{n+1} = 0.5 Q^n + 0.5 E(h)^2 Q^n 1258 1315 """ 1259 1316 1260 1317 # Save initial initial conserved quantities values 1261 self.backup_conserved_quantities() 1318 self.backup_conserved_quantities() 1262 1319 1263 1320 #-------------------------------------- … … 1289 1346 # Second Euler step 1290 1347 #------------------------------------ 1291 1348 1292 1349 # Compute fluxes across each element edge 1293 1350 self.compute_fluxes() … … 1300 1357 # of conserved quantities and cleanup 1301 1358 #------------------------------------ 1302 1359 1303 1360 # Combine steps 1304 1361 self.saxpy_conserved_quantities(0.5, 0.5) 1305 1362 1306 1363 # Update ghosts 1307 1364 self.update_ghosts() … … 1313 1370 self.update_boundary() 1314 1371 1315 1316 1372 ## 1373 # @brief 'rk3' time step method. 1374 # @param yieldstep The reporting time step. 1375 # @param finaltime The simulation final time. 1317 1376 def evolve_one_rk3_step(self, yieldstep, finaltime): 1318 """ 1319 One 3rd order RK timestep 1377 """One 3rd order RK timestep 1320 1378 Q^(1) = 3/4 Q^n + 1/4 E(h)^2 Q^n (at time t^n + h/2) 1321 1379 Q^{n+1} = 1/3 Q^n + 2/3 E(h) Q^(1) (at time t^{n+1}) … … 1323 1381 1324 1382 # Save initial initial conserved quantities values 1325 self.backup_conserved_quantities() 1383 self.backup_conserved_quantities() 1326 1384 1327 1385 initial_time = self.time 1328 1386 1329 1387 #-------------------------------------- 1330 1388 # First euler step … … 1352 1410 self.update_boundary() 1353 1411 1354 1355 1412 #------------------------------------ 1356 1413 # Second Euler step 1357 1414 #------------------------------------ 1358 1415 1359 1416 # Compute fluxes across each element edge 1360 1417 self.compute_fluxes() … … 1370 1427 # Combine steps 1371 1428 self.saxpy_conserved_quantities(0.25, 0.75) 1372 1429 1373 1430 # Update ghosts 1374 1431 self.update_ghosts() 1375 1432 1376 1377 1433 # Set substep time 1378 1434 self.time = initial_time + self.timestep*0.5 … … 1383 1439 # Update boundary values 1384 1440 self.update_boundary() 1385 1386 1441 1387 1442 #------------------------------------ 1388 1443 # Third Euler step 1389 1444 #------------------------------------ 1390 1445 1391 1446 # Compute fluxes across each element edge 1392 1447 self.compute_fluxes() … … 1399 1454 # and cleanup 1400 1455 #------------------------------------ 1401 1456 1402 1457 # Combine steps 1403 1458 self.saxpy_conserved_quantities(2.0/3.0, 1.0/3.0) 1404 1459 1405 1460 # Update ghosts 1406 1461 self.update_ghosts() 1407 1462 1408 1463 # Set new time 1409 self.time = initial_time + self.timestep 1464 self.time = initial_time + self.timestep 1410 1465 1411 1466 # Update vertex and edge values … … 1415 1470 self.update_boundary() 1416 1471 1417 1418 def evolve_to_end(self, finaltime = 1.0): 1419 """Iterate evolve all the way to the end """ 1472 ## 1473 # @brief Evolve simulation to a final time. 1474 # @param finaltime Sinulation final time. 1475 def evolve_to_end(self, finaltime=1.0): 1476 """Iterate evolve all the way to the end.""" 1420 1477 1421 1478 for _ in self.evolve(yieldstep=None, finaltime=finaltime): 1422 1479 pass 1423 1480 1424 1481 ## 1482 # @brief Backup conserved quantities. 1425 1483 def backup_conserved_quantities(self): 1426 1484 N = len(self) # Number_of_triangles … … 1429 1487 for name in self.conserved_quantities: 1430 1488 Q = self.quantities[name] 1431 Q.backup_centroid_values() 1432 1433 def saxpy_conserved_quantities(self,a,b): 1489 Q.backup_centroid_values() 1490 1491 ## 1492 # @brief ?? 1493 # @param a ?? 1494 # @param b ?? 1495 def saxpy_conserved_quantities(self, a, b): 1434 1496 N = len(self) #number_of_triangles 1435 1497 … … 1437 1499 for name in self.conserved_quantities: 1438 1500 Q = self.quantities[name] 1439 Q.saxpy_centroid_values(a,b) 1440 1441 1501 Q.saxpy_centroid_values(a, b) 1502 1503 ## 1504 # @brief Update boundary values for all conserved quantities. 1442 1505 def update_boundary(self): 1443 1506 """Go through list of boundary objects and update boundary values 1444 1507 for all conserved quantities on boundary. 1445 It is assumed that the ordering of conserved quantities is 1446 consistent between the domain and the boundary object, i.e. 1508 It is assumed that the ordering of conserved quantities is 1509 consistent between the domain and the boundary object, i.e. 1447 1510 the jth element of vector q must correspond to the jth conserved 1448 1511 quantity in domain. … … 1453 1516 for i, ((vol_id, edge_id), B) in enumerate(self.boundary_objects): 1454 1517 if B is None: 1455 print 'WARNING: Ignored boundary segment %d(None)'1518 print 'WARNING: Ignored boundary segment (None)' 1456 1519 else: 1520 print 'XXXX' 1457 1521 q = B.evaluate(vol_id, edge_id) 1458 1522 … … 1461 1525 Q.boundary_values[i] = q[j] 1462 1526 1463 1527 ## 1528 # @brief Compute fluxes. 1529 # @note MUST BE OVERRIDEN IN SUBCLASS! 1464 1530 def compute_fluxes(self): 1465 1531 msg = 'Method compute_fluxes must be overridden by Domain subclass' 1466 raise msg 1467 1468 1532 raise Exception, msg 1533 1534 ## 1535 # @brief 1536 # @param yieldstep 1537 # @param finaltime 1469 1538 def update_timestep(self, yieldstep, finaltime): 1470 1471 1539 from anuga.config import min_timestep, max_timestep 1472 1540 1473 1474 1475 1541 # Protect against degenerate timesteps arising from isolated 1476 1542 # triangles 1477 1543 # FIXME (Steve): This should be in shallow_water as it assumes x and y 1478 1544 # momentum 1479 if self.protect_against_isolated_degenerate_timesteps is True and \1545 if self.protect_against_isolated_degenerate_timesteps is True and \ 1480 1546 self.max_speed > 10.0: # FIXME (Ole): Make this configurable 1481 1547 1482 1548 # Setup 10 bins for speed histogram 1483 1549 from anuga.utilities.numerical_tools import histogram, create_bins 1484 1550 1485 1551 bins = create_bins(self.max_speed, 10) 1486 1552 hist = histogram(self.max_speed, bins) 1487 1553 1488 1554 # Look for characteristic signature 1489 if len(hist) > 1 and\ 1490 hist[-1] > 0 and\ 1555 if len(hist) > 1 and hist[-1] > 0 and \ 1491 1556 hist[4] == hist[5] == hist[6] == hist[7] == hist[8] == 0: 1492 1557 # Danger of isolated degenerate triangles 1493 # print self.timestepping_statistics(track_speeds=True) 1494 1558 # print self.timestepping_statistics(track_speeds=True) 1559 1495 1560 # Find triangles in last bin 1496 1561 # FIXME - speed up using Numeric … … 1498 1563 for i in range(self.number_of_full_triangles): 1499 1564 if self.max_speed[i] > bins[-1]: 1500 msg = 'Time=%f: Ignoring isolated high ' % self.time1501 msg += 'speed triangle ' 1502 msg += '#%d of %d with max speed=%f' \1503 %(i, self.number_of_full_triangles,1504 self.max_speed[i])1505 1565 msg = 'Time=%f: Ignoring isolated high ' % self.time 1566 msg += 'speed triangle ' 1567 msg += '#%d of %d with max speed=%f' \ 1568 % (i, self.number_of_full_triangles, 1569 self.max_speed[i]) 1570 1506 1571 # print 'Found offending triangle', i, 1507 # self.max_speed[i] 1508 self.get_quantity('xmomentum').set_values(0.0, indices=[i]) 1509 self.get_quantity('ymomentum').set_values(0.0, indices=[i]) 1572 # self.max_speed[i] 1573 self.get_quantity('xmomentum').\ 1574 set_values(0.0, indices=[i]) 1575 self.get_quantity('ymomentum').\ 1576 set_values(0.0, indices=[i]) 1510 1577 self.max_speed[i]=0.0 1511 1578 d += 1 1512 1513 #print 'Adjusted %d triangles' %d 1514 #print self.timestepping_statistics(track_speeds=True) 1515 1516 1517 1579 1518 1580 # self.timestep is calculated from speed of characteristics 1519 1581 # Apply CFL condition here … … 1524 1586 self.min_timestep = min(timestep, self.min_timestep) 1525 1587 1526 1527 1528 1588 # Protect against degenerate time steps 1529 1589 if timestep < min_timestep: 1530 1531 1590 # Number of consecutive small steps taken b4 taking action 1532 1591 self.smallsteps += 1 … … 1536 1595 1537 1596 if self._order_ == 1: 1538 msg = 'WARNING: Too small timestep %.16f reached ' \1539 %timestep1540 msg += 'even after %d steps of 1 order scheme' \1541 %self.max_smallsteps1597 msg = 'WARNING: Too small timestep %.16f reached ' \ 1598 % timestep 1599 msg += 'even after %d steps of 1 order scheme' \ 1600 % self.max_smallsteps 1542 1601 print msg 1543 1602 timestep = min_timestep # Try enforcing min_step … … 1549 1608 # Try to overcome situation by switching to 1 order 1550 1609 self._order_ = 1 1551 1552 1610 else: 1553 1611 self.smallsteps = 0 … … 1555 1613 self._order_ = 2 1556 1614 1557 1558 1615 # Ensure that final time is not exceeded 1559 1616 if finaltime is not None and self.time + timestep > finaltime : … … 1566 1623 self.timestep = timestep 1567 1624 1568 1569 1625 ## 1626 # @brief Compute forcing terms, if any. 1570 1627 def compute_forcing_terms(self): 1571 1628 """If there are any forcing functions driving the system … … 1577 1634 f(self) 1578 1635 1579 1636 ## 1637 # @brief Update vectors of conserved quantities. 1580 1638 def update_conserved_quantities(self): 1581 1639 """Update vectors of conserved quantities using previously … … 1597 1655 1598 1656 # Note that Q.explicit_update is reset by compute_fluxes 1599 # Where is Q.semi_implicit_update reset? 1657 # Where is Q.semi_implicit_update reset? 1600 1658 # It is reset in quantity_ext.c 1601 1602 1659 1660 ## 1661 # @brief ?? 1603 1662 def update_ghosts(self): 1604 1663 pass 1605 1664 1665 ## 1666 # @brief Extrapolate conserved quantities from centroid to vertices 1667 # and edge-midpoints for each volume. 1606 1668 def distribute_to_vertices_and_edges(self): 1607 1669 """Extrapolate conserved quantities from centroid to … … 1621 1683 #Q.limit() 1622 1684 else: 1623 raise 'Unknown order'1685 raise Exception, 'Unknown order' 1624 1686 #Q.interpolate_from_vertices_to_edges() 1625 1687 1626 1688 ## 1689 # @brief Calculate the norm of the centroid values of a specific quantity, 1690 # using normfunc. 1691 # @param quantity 1692 # @param normfunc 1627 1693 def centroid_norm(self, quantity, normfunc): 1628 1694 """Calculate the norm of the centroid values … … 1633 1699 common normfuncs are provided in the module utilities.norms 1634 1700 """ 1701 1635 1702 return normfunc(self.quantities[quantity].centroid_values) 1636 1637 1703 1638 1704 … … 1652 1718 # Psyco isn't supported on 64 bit systems, but it doesn't matter 1653 1719 else: 1654 msg = 'WARNING: psyco (speedup) could not import'+\1655 ', you may want to consider installing it'1720 msg = ('WARNING: psyco (speedup) could not be imported, ' 1721 'you may want to consider installing it') 1656 1722 print msg 1657 1723 else:
Note: See TracChangeset
for help on using the changeset viewer.
