Changeset 6181


Ignore:
Timestamp:
Jan 16, 2009, 4:06:55 PM (10 years ago)
Author:
rwilson
Message:

PEP8'd file, plus @brief.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py

    r6145 r6181  
    3838
    3939
     40##
     41# @brief Basic Domain class
     42# @note Subclasses Mesh.
    4043class Domain(Mesh):
    4144
    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):
    6281
    6382        """Instantiate generic computational Domain.
     
    7695          tagged_elements:
    7796          ...
    78 
    79 
    8097        """
    8198
     
    85102        else:
    86103            coordinates = source
    87 
    88104
    89105        # In case a filename has been specified, extract content
     
    94110                                         use_cache=use_cache,
    95111                                         verbose=verbose)
    96 
    97112
    98113        # Initialise underlying mesh structure
     
    108123        if verbose: print 'Initialising Domain'
    109124
    110         # List of quantity names entering
    111         # the conservation equations
     125        # List of quantity names entering the conservation equations
    112126        if conserved_quantities is None:
    113127            self.conserved_quantities = []
     
    121135            self.other_quantities = other_quantities
    122136
    123 
    124137        # Build dictionary of Quantity instances keyed by quantity names
    125138        self.quantities = {}
     
    138151            self.full_send_dict = {}
    139152        else:
    140             self.full_send_dict  = full_send_dict
     153            self.full_send_dict = full_send_dict
    141154
    142155        # List of other quantity names
     
    149162        self.numproc = numproc
    150163
    151 
    152164        # Setup Communication Buffers
    153165        if verbose: print 'Domain: Set up communication buffers (parallel)'
     
    155167        for key in self.full_send_dict:
    156168            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))
    158171
    159172        for key in self.ghost_recv_dict:
    160173            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))
    162176
    163177        # Setup cell full flag
     
    173187        # Test the assumption that all full triangles are store before
    174188        # 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')
    178192
    179193        # Defaults
     
    182196        from anuga.config import timestepping_method
    183197        from anuga.config import protect_against_isolated_degenerate_timesteps
    184         from anuga.config import default_order       
     198        from anuga.config import default_order
     199
    185200        self.beta_w = beta_w
    186201        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
    189205        self.set_default_order(default_order)
    190206
     
    196212        self.set_timestepping_method(timestepping_method)
    197213        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
    200215
    201216        # Model time
     
    209224
    210225        self.last_walltime = walltime()
    211        
     226
    212227        # Monitoring
    213228        self.quantities_to_be_monitored = None
    214229        self.monitor_polygon = None
    215         self.monitor_time_interval = None               
     230        self.monitor_time_interval = None
    216231        self.monitor_indices = None
    217232
    218233        # Checkpointing and storage
    219234        from anuga.config import default_datadir
     235
    220236        self.datadir = default_datadir
    221237        self.simulation_name = 'domain'
    222238        self.checkpoint = False
    223239
    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.
    227242        N = len(self) # Number_of_triangles
    228243        self.already_computed_flux = num.zeros((N, 3), num.Int)
    229244
    230245        # Storage for maximal speeds computed for each triangle by
    231         # compute_fluxes
     246        # compute_fluxes.
    232247        # This is used for diagnostics only (reset at every yieldstep)
    233248        self.max_speed = num.zeros(N, num.Float)
    234249
    235250        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.
    238253            if verbose: print 'Domain: Initialising quantity values'
    239254            self.set_quantity_vertices_dict(vertex_quantity_dict)
    240255
    241 
    242256        if verbose: print 'Domain: Done'
    243257
    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):
    253269        """Get conserved quantities at volume vol_id
    254270
     
    259275
    260276        Return value: Vector of length == number_of_conserved quantities
    261 
    262277        """
    263278
     
    265280            msg = 'Values for both vertex and edge was specified.'
    266281            msg += 'Only one (or none) is allowed.'
    267             raise msg
     282            raise Exception, msg
    268283
    269284        q = num.zeros(len(self.conserved_quantities), num.Float)
     
    280295        return q
    281296
     297    ##
     298    # @brief Set the relative model time.
     299    # @param time The new model time (seconds).
    282300    def set_time(self, time=0.0):
    283301        """Set the model time (seconds)"""
     
    287305        self.time = time
    288306
     307    ##
     308    # @brief Get the model time.
     309    # @return The absolute model time (seconds).
    289310    def get_time(self):
    290311        """Get the absolute model time (seconds)"""
     
    292313        return self.time + self.starttime
    293314
    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"""
    297320
    298321        self.beta = beta
     
    302325            Q.set_beta(beta)
    303326
     327    ##
     328    # @brief Get the beta value used for limiting.
     329    # @return The beta value used for limiting.
    304330    def get_beta(self):
    305         """Get default beta for limiting
    306         """
     331        """Get default beta for limiting"""
    307332
    308333        return self.beta
    309334
     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.
    310339    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
    315343        assert n in [1,2], msg
    316344
     
    327355            #self.set_timestepping_method('rk2')
    328356            #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.
    331361    def set_quantity_vertices_dict(self, quantity_dict):
    332362        """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
    343371        # FIXME: Could we name this a bit more intuitively
    344372        # E.g. set_quantities_from_dictionary
     
    346374            self.set_quantity(key, quantity_dict[key], location='vertices')
    347375
    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):
    350384        """Set values for named quantity
    351 
    352385
    353386        One keyword argument is documented here:
     
    370403        # Assign values
    371404        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):
    374414        """Add values to a named quantity
    375        
     415
    376416        E.g add_quantity('elevation', X)
    377        
     417
    378418        Option are the same as in set_quantity.
    379419        """
    380            
     420
    381421        # Do the expression stuff
    382422        if kwargs.has_key('expression'):
    383423            expression = kwargs['expression']
    384424            Q2 = self.create_quantity_from_expression(expression)
    385         else:   
     425        else:
    386426            # Create new temporary quantity
    387427            Q2 = Quantity(self)
    388        
     428
    389429            # Assign specified values to temporary quantity
    390430            Q2.set_values(*args, **kwargs)
    391            
     431
    392432        # Add temporary quantity to named quantity
    393433        Q1 = self.get_quantity(name)
    394434        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.
    397439    def get_quantity_names(self):
    398440        """Get a list of all the quantity names that this domain is aware of.
    399441        Any value in the result should be a valid input to get_quantity.
    400442        """
     443
    401444        return self.quantities.keys()
    402445
    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):
    404456        """Get pointer to quantity object.
    405457
     
    413465        return self.quantities[name] #.get_values( location, indices = indices)
    414466
    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.
    417472    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.
    419474
    420475        Combine existing quantities in domain using expression and return
     
    426481
    427482        Examples creating derived quantities:
    428 
    429483            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'
    432485            Absolute_momentum = domain.create_quantity_from_expression(exp)
    433 
    434486        """
    435487
    436488        from anuga.abstract_2d_finite_volumes.util import\
    437489             apply_expression_to_dictionary
    438        
     490
    439491        return apply_expression_to_dictionary(expression, self.quantities)
    440492
    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.
    443498    def set_boundary(self, boundary_map):
    444499        """Associate boundary objects with tagged boundary segments.
     
    461516        self.boundary_objects:  ((vol_id, edge_id), boundary_object)
    462517
    463 
    464518        Pre-condition:
    465519          self.boundary has been built.
     
    471525        exception is raised.
    472526        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
    475528
    476529        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.
    480532
    481533        Boundary objects that are None will be skipped.
    482534
    483         If a boundary_map has already been set
    484         (i.e. set_boundary has been called before), the old boundary map
    485         will be updated with new values. The new map need not define all
    486         boundary tags, and can thus change only those 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.
    487539
    488540        FIXME: If set_boundary is called multiple times and if Boundary
    489541        object is changed into None, the neighbour structure will not be
    490542        restored!!!
    491 
    492 
    493543        """
    494544
     
    496546            # This the first call to set_boundary. Store
    497547            # 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:
    500550            # This is a modification of an already existing map
    501551            # Update map an proceed normally
    502 
    503552            for key in boundary_map.keys():
    504553                self.boundary_map[key] = boundary_map[key]
    505                
    506        
     554
    507555        # FIXME (Ole): Try to remove the sorting and fix test_mesh.py
    508556        x = self.boundary.keys()
     
    511559        # Loop through edges that lie on the boundary and associate them with
    512560        # callable boundary objects depending on their tags
    513         self.boundary_objects = []       
     561        self.boundary_objects = []
    514562        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)]
    516564
    517565            if self.boundary_map.has_key(tag):
     
    519567
    520568                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))
    522570                    self.neighbours[vol_id, edge_id] = \
    523                                             -len(self.boundary_objects)
     571                                        -len(self.boundary_objects)
    524572                else:
    525573                    pass
    526574                    #FIXME: Check and perhaps fix neighbour structure
    527 
    528 
    529575            else:
    530576                msg = 'ERROR (domain.py): Tag "%s" has not been ' %tag
     
    533579                msg += 'in set_boundary.\n'
    534580                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
    538587    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
    542590        It is most often called with the following parameters;
    543591        (self, tag, quantity, X, location='vertices')
    544         tag: the name of the regional tag used to specify the region
     592        tag:      the name of the regional tag used to specify the region
    545593        quantity: Name of quantity to change
    546         X: const or function - how the quantity is changed
     594        X:        const or function - how the quantity is changed
    547595        location: Where values are to be stored.
    548596            Permissible options are: vertices, centroid and unique vertices
     
    558606            func = region_set_region(*args, **kwargs)
    559607            self._set_region(func)
    560            
    561        
     608
     609    ##
     610    # @brief ??
     611    # @param functions A list or tuple of ??
    562612    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
    563617        # The order of functions in the list is used.
    564         if type(functions) not in [types.ListType,types.TupleType]:
    565             functions = [functions]
    566618        for function in functions:
    567619            for tag in self.tagged_elements.keys():
    568620                function(tag, self.tagged_elements[tag], self)
    569621
    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.
    573628    def set_quantities_to_be_monitored(self, q,
    574                                        polygon=None,
    575                                        time_interval=None):
     629                                             polygon=None,
     630                                             time_interval=None):
    576631        """Specify which quantities will be monitored for extrema.
    577632
    578633        q must be either:
    579           - the name of a quantity or a derived quantity such as 'stage-elevation'
     634          - the name of a quantity or derived quantity such as 'stage-elevation'
    580635          - a list of quantity names
    581636          - None
     
    583638        In the two first cases, the named quantities will be monitored at
    584639        each internal timestep
    585        
     640
    586641        If q is None, monitoring will be switched off altogether.
    587642
    588         polygon (if specified) will restrict monitoring to triangles inside polygon.
     643        polygon (if specified) will only monitor triangles inside polygon.
    589644        If omitted all triangles will be included.
    590645
    591         time_interval (if specified) will restrict monitoring to time steps in
     646        time_interval, if specified, will restrict monitoring to time steps in
    592647        that interval. If omitted all timesteps will be included.
    593648        """
    594649
    595650        from anuga.abstract_2d_finite_volumes.util import\
    596              apply_expression_to_dictionary       
     651             apply_expression_to_dictionary
    597652
    598653        if q is None:
     
    600655            self.monitor_polygon = None
    601656            self.monitor_time_interval = None
    602             self.monitor_indices = None           
     657            self.monitor_indices = None
    603658            return
    604659
     660        # coerce 'q' to a list if it's a string
    605661        if isinstance(q, basestring):
    606             q = [q] # Turn argument into a list
    607 
    608         # Check correcness and initialise
     662            q = [q]
     663
     664        # Check correctness and initialise
    609665        self.quantities_to_be_monitored = {}
    610666        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
    614669
    615670            if not quantity_name in self.quantities:
     
    622677                          'min_location': None, # Argmin (x, y)
    623678                          'max_location': None, # Argmax (x, y)
    624                           'min_time': None,     # Argmin (t) 
     679                          'min_time': None,     # Argmin (t)
    625680                          'max_time': None}     # Argmax (t)
    626            
     681
    627682            self.quantities_to_be_monitored[quantity_name] = info_block
    628 
    629            
    630683
    631684        if polygon is not None:
    632685            # Check input
    633686            if isinstance(polygon, basestring):
    634 
    635687                # Check if multiple quantities were accidentally
    636688                # 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' %polygon
    640                
     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
    641693                if polygon in self.quantities:
    642694                    raise Exception, msg
    643                
     695
    644696                try:
    645697                    apply_expression_to_dictionary(polygon, self.quantities)
    646698                except:
    647                     # At least polygon wasn't an expression involving quantitites
     699                    # At least polygon wasn't expression involving quantitites
    648700                    pass
    649701                else:
     
    651703
    652704                # 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\' ' %polygon
     705                msg = ('argument "polygon" must not be a string: '
     706                       'I got polygon="%s"') % polygon
    655707                raise Exception, msg
    656 
    657708
    658709            # Get indices for centroids that are inside polygon
    659710            points = self.get_centroid_coordinates(absolute=True)
    660711            self.monitor_indices = inside_polygon(points, polygon)
    661            
    662712
    663713        if time_interval is not None:
    664714            assert len(time_interval) == 2
    665715
    666        
    667716        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.
    674722    def check_integrity(self):
    675723        Mesh.check_integrity(self)
     
    680728
    681729        ##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.
    684734    def write_time(self, track_speeds=False):
    685735        print self.timestepping_statistics(track_speeds)
    686736
    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
    692745
    693746        Optional boolean keyword track_speeds decides whether to report
     
    696749
    697750        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.
    699752        """
    700753
    701754        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
    705757        qwidth = self.qwidth = 12
    706 
    707758
    708759        msg = ''
     
    723774        model_time = self.get_time()
    724775        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)
    727778        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)
    730781        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
    738789        if track_speeds is True:
    739790            msg += '\n'
    740 
    741791
    742792            # Setup 10 bins for speed histogram
     
    745795
    746796            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))
    749799            msg += '  Histogram:\n'
    750800
     
    753803                lo = hi
    754804                if i+1 < len(bins):
    755                     # Open upper interval               
     805                    # Open upper interval
    756806                    hi = bins[i+1]
    757                     msg += '    [%f, %f[: %d\n' %(lo, hi, count)               
     807                    msg += '    [%f, %f[: %d\n' % (lo, hi, count)
    758808                else:
    759809                    # Closed upper interval
    760810                    hi = max(self.max_speed)
    761                     msg += '    [%f, %f]: %d\n' %(lo, hi, count)
    762 
     811                    msg += '    [%f, %f]: %d\n' % (lo, hi, count)
    763812
    764813            N = len(self.max_speed)
     
    770819                k = 0
    771820                lower = min(speed)
    772                 for i, a in enumerate(speed):       
     821                for i, a in enumerate(speed):
    773822                    if i % (N/10) == 0 and i != 0:
    774823                        # 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)
    776825                        lower = a
    777826                        k = i
    778                        
     827
    779828                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
    786831            # Find index of largest computed flux speed
    787832            if triangle_id is None:
    788833                k = self.k = num.argmax(self.max_speed)
    789834            else:
    790                 errmsg = 'Triangle_id %d does not exist in mesh: %s' %(triangle_id,
    791                                                                    str(self))
     835                errmsg = 'Triangle_id %d does not exist in mesh: %s' \
     836                             % (triangle_id, str(self))
    792837                assert 0 <= triangle_id < len(self), errmsg
    793838                k = self.k = triangle_id
    794            
    795839
    796840            x, y = self.get_centroid_coordinates()[k]
    797841            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)
    803847            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)
    805849            else:
    806                 msg += 'had computed speed: %.6f m/s ' %(max_speed)
    807                
     850                msg += 'had computed speed: %.6f m/s ' % (max_speed)
     851
    808852            if max_speed > 0.0:
    809                 msg += '(timestep=%.6f)\n' %(radius/max_speed)
     853                msg += '(timestep=%.6f)\n' % (radius/max_speed)
    810854            else:
    811                 msg += '(timestep=%.6f)\n' %(0)     
    812            
     855                msg += '(timestep=%.6f)\n' % (0)
     856
    813857            # Report all quantity values at vertices, edges and centroid
    814858            msg += '    Quantity'
     
    819863                V = q.get_values(location='vertices', indices=[k])[0]
    820864                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])
    831875
    832876                msg += s
     
    834878        return msg
    835879
    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):
    838885        print self.boundary_statistics(quantities, tags)
    839886
    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):
    841893        """Output statistics about boundary forcing at each timestep
    842 
    843894
    844895        Input:
     
    847898          tags:       either None, a string or a list of strings naming the
    848899                      tags to be reported
    849 
    850900
    851901        Example output:
     
    856906        Tag 'ocean'
    857907
    858 
    859908        If quantities are specified only report on those. Otherwise take all
    860909        conserved quantities.
    861910        If tags are specified only report on those, otherwise take all tags.
    862 
    863         """
     911        """
     912
     913        import types, string
    864914
    865915        # Input checks
    866         import types, string
    867 
    868916        if quantities is None:
    869917            quantities = self.conserved_quantities
     
    871919            quantities = [quantities] #Turn it into a list
    872920
    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)
    875923        assert type(quantities) == types.ListType, msg
    876 
    877924
    878925        if tags is None:
     
    881928            tags = [tags] #Turn it into a list
    882929
    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)
    885932        assert type(tags) == types.ListType, msg
    886933
     
    893940
    894941        # Output statistics
    895         msg = 'Boundary values at time %.4f:\n' %self.get_time()
     942        msg = 'Boundary values at time %.4f:\n' % self.get_time()
    896943        for tag in tags:
    897             msg += '    %s:\n' %tag
     944            msg += '    %s:\n' % tag
    898945
    899946            for name in quantities:
     
    902949                # Find range of boundary values for tag and q
    903950                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):
    906952                    if self.boundary[(vol_id, edge_id)] == tag:
    907953                        v = q.boundary_values[i]
     
    910956
    911957                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)
    914960                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)
    918963
    919964        return msg
    920965
    921 
     966    ##
     967    # @brief Update extrema if requested by set_quantities_to_be_monitored.
    922968    def update_extrema(self):
    923969        """Update extrema if requested by set_quantities_to_be_monitored.
     
    929975        # Define a tolerance for extremum computations
    930976        from anuga.config import single_precision as epsilon
    931        
     977
    932978        if self.quantities_to_be_monitored is None:
    933979            return
     
    9541000            # of the epsilon)
    9551001            maxval = Q.get_maximum_value(self.monitor_indices)
    956             if info_block['max'] is None or\
     1002            if info_block['max'] is None or \
    9571003                   maxval > info_block['max'] + epsilon:
    9581004                info_block['max'] = maxval
     
    9611007                info_block['max_time'] = self.time
    9621008
    963 
    9641009            # Update minimum
    9651010            minval = Q.get_minimum_value(self.monitor_indices)
    966             if info_block['min'] is None or\
     1011            if info_block['min'] is None or \
    9671012                   minval < info_block['min'] - epsilon:
    968                 info_block['min'] = minval               
     1013                info_block['min'] = minval
    9691014                minloc = Q.get_minimum_location()
    9701015                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'):
    9761023        """Return string with statistics about quantities for
    9771024        printing or logging
     
    9801027
    9811028           set_quantities_to_be_monitored
    982            
    9831029        """
    9841030
     
    9861032
    9871033        # Output statistics
    988         msg = 'Monitored quantities at time %.4f:\n' %self.get_time()
     1034        msg = 'Monitored quantities at time %.4f:\n' % self.get_time()
    9891035        if self.monitor_polygon is not None:
    9901036            p_str = str(self.monitor_polygon)
    991             msg += '- Restricted by polygon: %s' %p_str[:maxlen]
     1037            msg += '- Restricted by polygon: %s' % p_str[:maxlen]
    9921038            if len(p_str) >= maxlen:
    9931039                msg += '...\n'
     
    9951041                msg += '\n'
    9961042
    997 
    9981043        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)
    10011046            time_interval_start = self.monitor_time_interval[0]
    10021047        else:
    10031048            time_interval_start = 0.0
    10041049
    1005            
    10061050        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))
    10231065
    10241066        return msg
    10251067
     1068    ##
     1069    # @brief Get the timestep method.
     1070    # @return The timestep method. One of 'euler', 'rk2' or 'rk3'.
    10261071    def get_timestepping_method(self):
    10271072        return self.timestepping_method
    10281073
    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):
    10311079        if timestepping_method in ['euler', 'rk2', 'rk3']:
    10321080            self.timestepping_method = timestepping_method
    10331081            return
    10341082
    1035         msg = '%s is an incorrect timestepping type'% timestepping_method
     1083        msg = '%s is an incorrect timestepping type' % timestepping_method
    10361084        raise Exception, msg
    10371085
     1086    ##
     1087    # @brief Get the Domain simulation name.
     1088    # @return The simulation name string.
    10381089    def get_name(self):
    10391090        return self.simulation_name
    10401091
     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.
    10411096    def set_name(self, name):
    10421097        """Assign a name to this simulation.
    10431098        This will be used to identify the output sww file.
    1044 
    1045         """
     1099        """
     1100
     1101        # remove any '.sww' end
    10461102        if name.endswith('.sww'):
    10471103            name = name[:-4]
    1048            
     1104
    10491105        self.simulation_name = name
    10501106
     1107    ##
     1108    # @brief Get data directory path.
     1109    # @return The data directory path string.
    10511110    def get_datadir(self):
    10521111        return self.datadir
    10531112
     1113    ##
     1114    # @brief Set data directory path.
     1115    # @param name The data directory path string.
    10541116    def set_datadir(self, name):
    10551117        self.datadir = name
    10561118
     1119    ##
     1120    # @brief Get the start time value.
     1121    # @return The start time value (float).
    10571122    def get_starttime(self):
    10581123        return self.starttime
    10591124
     1125    ##
     1126    # @brief Set the start time value.
     1127    # @param time The start time value.
    10601128    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):
    10741145        """Evolve model through time starting from self.starttime.
    1075 
    10761146
    10771147        yieldstep: Interval between yields where results are stored,
     
    10881158        If both duration and finaltime are given an exception is thrown.
    10891159
    1090 
    10911160        skip_initial_step: Boolean flag that decides whether the first
    10921161        yield step is skipped or not. This is useful for example to avoid
    10931162        duplicate steps when multiple evolve processes are dove tailed.
    10941163
    1095 
    10961164        Evolve is implemented as a generator and is to be called as such, e.g.
    10971165
     
    10991167            <Do something with domain and t>
    11001168
    1101 
    11021169        All times are given in seconds
    1103 
    11041170        """
    11051171
     
    11071173
    11081174        # 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()
    11141180        assert hasattr(self, 'boundary_objects'), msg
    1115 
    11161181
    11171182        if yieldstep is None:
     
    11221187        self._order_ = self.default_order
    11231188
    1124 
    11251189        if finaltime is not None and duration is not None:
    11261190            # print 'F', finaltime, duration
    11271191            msg = 'Only one of finaltime and duration may be specified'
    1128             raise msg
     1192            raise Exception, msg
    11291193        else:
    11301194            if finaltime is not None:
     
    11331197                self.finaltime = self.starttime + float(duration)
    11341198
    1135 
    1136 
    11371199        N = len(self) # Number of triangles
    11381200        self.yieldtime = 0.0 # Track time between 'yields'
     
    11441206        self.number_of_first_order_steps = 0
    11451207
    1146 
    11471208        # Update ghosts
    11481209        self.update_ghosts()
     
    11531214        # Update extrema if necessary (for reporting)
    11541215        self.update_extrema()
    1155        
     1216
    11561217        # Initial update boundary values
    11571218        self.update_boundary()
     
    11651226
    11661227        while True:
    1167 
    11681228            # Evolve One Step, using appropriate timestepping method
    11691229            if self.get_timestepping_method() == 'euler':
    1170                 self.evolve_one_euler_step(yieldstep,finaltime)
    1171                
     1230                self.evolve_one_euler_step(yieldstep, finaltime)
     1231
    11721232            elif self.get_timestepping_method() == 'rk2':
    1173                 self.evolve_one_rk2_step(yieldstep,finaltime)
     1233                self.evolve_one_rk2_step(yieldstep, finaltime)
    11741234
    11751235            elif self.get_timestepping_method() == 'rk3':
    1176                 self.evolve_one_rk3_step(yieldstep,finaltime)               
    1177            
     1236                self.evolve_one_rk3_step(yieldstep, finaltime)
     1237
    11781238            # Update extrema if necessary (for reporting)
    1179             self.update_extrema()           
    1180 
     1239            self.update_extrema()
    11811240
    11821241            self.yieldtime += self.timestep
     
    11871246            # Yield results
    11881247            if finaltime is not None and self.time >= finaltime-epsilon:
    1189 
    11901248                if self.time > finaltime:
    11911249                    # FIXME (Ole, 30 April 2006): Do we need this check?
    11921250                    # Probably not (Ole, 18 September 2008). Now changed to
    11931251                    # 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')
    11961254                    raise Exception, msg
    1197                    
    11981255
    11991256                # Yield final time and stop
     
    12021259                break
    12031260
    1204 
    12051261            if self.yieldtime >= yieldstep:
    12061262                # Yield (intermediate) time and allow inspection of domain
    1207 
    12081263                if self.checkpoint is True:
    12091264                    self.store_checkpoint()
     
    12211276                self.max_speed = num.zeros(N, num.Float)
    12221277
     1278    ##
     1279    # @brief 'Euler' time step method.
     1280    # @param yieldstep The reporting time step.
     1281    # @param finaltime The simulation final time.
    12231282    def evolve_one_euler_step(self, yieldstep, finaltime):
    1224         """
    1225         One Euler Time Step
     1283        """One Euler Time Step
    12261284        Q^{n+1} = E(h) Q^n
    12271285        """
     
    12481306        self.update_boundary()
    12491307
    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.
    12541312    def evolve_one_rk2_step(self, yieldstep, finaltime):
    1255         """
    1256         One 2nd order RK timestep
     1313        """One 2nd order RK timestep
    12571314        Q^{n+1} = 0.5 Q^n + 0.5 E(h)^2 Q^n
    12581315        """
    12591316
    12601317        # Save initial initial conserved quantities values
    1261         self.backup_conserved_quantities()           
     1318        self.backup_conserved_quantities()
    12621319
    12631320        #--------------------------------------
     
    12891346        # Second Euler step
    12901347        #------------------------------------
    1291            
     1348
    12921349        # Compute fluxes across each element edge
    12931350        self.compute_fluxes()
     
    13001357        # of conserved quantities and cleanup
    13011358        #------------------------------------
    1302        
     1359
    13031360        # Combine steps
    13041361        self.saxpy_conserved_quantities(0.5, 0.5)
    1305  
     1362
    13061363        # Update ghosts
    13071364        self.update_ghosts()
     
    13131370        self.update_boundary()
    13141371
    1315 
    1316 
     1372    ##
     1373    # @brief 'rk3' time step method.
     1374    # @param yieldstep The reporting time step.
     1375    # @param finaltime The simulation final time.
    13171376    def evolve_one_rk3_step(self, yieldstep, finaltime):
    1318         """
    1319         One 3rd order RK timestep
     1377        """One 3rd order RK timestep
    13201378        Q^(1) = 3/4 Q^n + 1/4 E(h)^2 Q^n  (at time t^n + h/2)
    13211379        Q^{n+1} = 1/3 Q^n + 2/3 E(h) Q^(1) (at time t^{n+1})
     
    13231381
    13241382        # Save initial initial conserved quantities values
    1325         self.backup_conserved_quantities()           
     1383        self.backup_conserved_quantities()
    13261384
    13271385        initial_time = self.time
    1328        
     1386
    13291387        #--------------------------------------
    13301388        # First euler step
     
    13521410        self.update_boundary()
    13531411
    1354 
    13551412        #------------------------------------
    13561413        # Second Euler step
    13571414        #------------------------------------
    1358            
     1415
    13591416        # Compute fluxes across each element edge
    13601417        self.compute_fluxes()
     
    13701427        # Combine steps
    13711428        self.saxpy_conserved_quantities(0.25, 0.75)
    1372  
     1429
    13731430        # Update ghosts
    13741431        self.update_ghosts()
    13751432
    1376 
    13771433        # Set substep time
    13781434        self.time = initial_time + self.timestep*0.5
     
    13831439        # Update boundary values
    13841440        self.update_boundary()
    1385 
    13861441
    13871442        #------------------------------------
    13881443        # Third Euler step
    13891444        #------------------------------------
    1390            
     1445
    13911446        # Compute fluxes across each element edge
    13921447        self.compute_fluxes()
     
    13991454        # and cleanup
    14001455        #------------------------------------
    1401        
     1456
    14021457        # Combine steps
    14031458        self.saxpy_conserved_quantities(2.0/3.0, 1.0/3.0)
    1404  
     1459
    14051460        # Update ghosts
    14061461        self.update_ghosts()
    14071462
    14081463        # Set new time
    1409         self.time = initial_time + self.timestep       
     1464        self.time = initial_time + self.timestep
    14101465
    14111466        # Update vertex and edge values
     
    14151470        self.update_boundary()
    14161471
    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."""
    14201477
    14211478        for _ in self.evolve(yieldstep=None, finaltime=finaltime):
    14221479            pass
    14231480
    1424 
     1481    ##
     1482    # @brief Backup conserved quantities.
    14251483    def backup_conserved_quantities(self):
    14261484        N = len(self) # Number_of_triangles
     
    14291487        for name in self.conserved_quantities:
    14301488            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):
    14341496        N = len(self) #number_of_triangles
    14351497
     
    14371499        for name in self.conserved_quantities:
    14381500            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.
    14421505    def update_boundary(self):
    14431506        """Go through list of boundary objects and update boundary values
    14441507        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.
    14471510        the jth element of vector q must correspond to the jth conserved
    14481511        quantity in domain.
     
    14531516        for i, ((vol_id, edge_id), B) in enumerate(self.boundary_objects):
    14541517            if B is None:
    1455                 print 'WARNING: Ignored boundary segment %d (None)'
     1518                print 'WARNING: Ignored boundary segment (None)'
    14561519            else:
     1520                print 'XXXX'
    14571521                q = B.evaluate(vol_id, edge_id)
    14581522
     
    14611525                    Q.boundary_values[i] = q[j]
    14621526
    1463 
     1527    ##
     1528    # @brief Compute fluxes.
     1529    # @note MUST BE OVERRIDEN IN SUBCLASS!
    14641530    def compute_fluxes(self):
    14651531        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
    14691538    def update_timestep(self, yieldstep, finaltime):
    1470 
    14711539        from anuga.config import min_timestep, max_timestep
    14721540
    1473        
    1474        
    14751541        # Protect against degenerate timesteps arising from isolated
    14761542        # triangles
    14771543        # FIXME (Steve): This should be in shallow_water as it assumes x and y
    14781544        # momentum
    1479         if self.protect_against_isolated_degenerate_timesteps is True and\
     1545        if self.protect_against_isolated_degenerate_timesteps is True and \
    14801546               self.max_speed > 10.0: # FIXME (Ole): Make this configurable
    14811547
    14821548            # Setup 10 bins for speed histogram
    14831549            from anuga.utilities.numerical_tools import histogram, create_bins
    1484        
     1550
    14851551            bins = create_bins(self.max_speed, 10)
    14861552            hist = histogram(self.max_speed, bins)
    14871553
    14881554            # 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 \
    14911556                hist[4] == hist[5] == hist[6] == hist[7] == hist[8] == 0:
    14921557                    # Danger of isolated degenerate triangles
    1493                     # print self.timestepping_statistics(track_speeds=True) 
    1494                    
     1558                    # print self.timestepping_statistics(track_speeds=True)
     1559
    14951560                    # Find triangles in last bin
    14961561                    # FIXME - speed up using Numeric
     
    14981563                    for i in range(self.number_of_full_triangles):
    14991564                        if self.max_speed[i] > bins[-1]:
    1500                             msg = 'Time=%f: Ignoring isolated high ' %self.time
    1501                             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
    15061571                            # 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])
    15101577                            self.max_speed[i]=0.0
    15111578                            d += 1
    1512                    
    1513                     #print 'Adjusted %d triangles' %d       
    1514                     #print self.timestepping_statistics(track_speeds=True)     
    1515                
    1516 
    1517                        
     1579
    15181580        # self.timestep is calculated from speed of characteristics
    15191581        # Apply CFL condition here
     
    15241586        self.min_timestep = min(timestep, self.min_timestep)
    15251587
    1526            
    1527  
    15281588        # Protect against degenerate time steps
    15291589        if timestep < min_timestep:
    1530 
    15311590            # Number of consecutive small steps taken b4 taking action
    15321591            self.smallsteps += 1
     
    15361595
    15371596                if self._order_ == 1:
    1538                     msg = 'WARNING: Too small timestep %.16f reached '\
    1539                           %timestep
    1540                     msg += 'even after %d steps of 1 order scheme'\
    1541                            %self.max_smallsteps
     1597                    msg = 'WARNING: Too small timestep %.16f reached ' \
     1598                              % timestep
     1599                    msg += 'even after %d steps of 1 order scheme' \
     1600                               % self.max_smallsteps
    15421601                    print msg
    15431602                    timestep = min_timestep  # Try enforcing min_step
     
    15491608                    # Try to overcome situation by switching to 1 order
    15501609                    self._order_ = 1
    1551 
    15521610        else:
    15531611            self.smallsteps = 0
     
    15551613                self._order_ = 2
    15561614
    1557 
    15581615        # Ensure that final time is not exceeded
    15591616        if finaltime is not None and self.time + timestep > finaltime :
     
    15661623        self.timestep = timestep
    15671624
    1568 
    1569 
     1625    ##
     1626    # @brief Compute forcing terms, if any.
    15701627    def compute_forcing_terms(self):
    15711628        """If there are any forcing functions driving the system
     
    15771634            f(self)
    15781635
    1579 
     1636    ##
     1637    # @brief Update vectors of conserved quantities.
    15801638    def update_conserved_quantities(self):
    15811639        """Update vectors of conserved quantities using previously
     
    15971655
    15981656            # 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?
    16001658            # It is reset in quantity_ext.c
    1601            
    1602 
     1659
     1660    ##
     1661    # @brief ??
    16031662    def update_ghosts(self):
    16041663        pass
    16051664
     1665    ##
     1666    # @brief Extrapolate conserved quantities from centroid to vertices
     1667    #        and edge-midpoints for each volume.
    16061668    def distribute_to_vertices_and_edges(self):
    16071669        """Extrapolate conserved quantities from centroid to
     
    16211683                #Q.limit()
    16221684            else:
    1623                 raise 'Unknown order'
     1685                raise Exception, 'Unknown order'
    16241686            #Q.interpolate_from_vertices_to_edges()
    16251687
    1626 
     1688    ##
     1689    # @brief Calculate the norm of the centroid values of a specific quantity,
     1690    #        using normfunc.
     1691    # @param quantity
     1692    # @param normfunc
    16271693    def centroid_norm(self, quantity, normfunc):
    16281694        """Calculate the norm of the centroid values
     
    16331699        common normfuncs are provided in the module utilities.norms
    16341700        """
     1701
    16351702        return normfunc(self.quantities[quantity].centroid_values)
    1636 
    16371703
    16381704
     
    16521718            # Psyco isn't supported on 64 bit systems, but it doesn't matter
    16531719        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')
    16561722            print msg
    16571723    else:
Note: See TracChangeset for help on using the changeset viewer.