source: anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py @ 7053

Last change on this file since 7053 was 7053, checked in by ole, 15 years ago

Made timestepping_statistics report absolute coordinates for centroid of
triangle with fastest speeds (track_speeds=True).

File size: 66.0 KB
Line 
1"""Class Domain - 2D triangular domains for finite-volume computations of
2   conservation laws.
3
4
5   Copyright 2004
6   Ole Nielsen, Stephen Roberts, Duncan Gray
7   Geoscience Australia
8"""
9
10import types
11from time import time as walltime
12
13from anuga.config import epsilon
14from anuga.config import beta_euler, beta_rk2
15from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
16from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
17     import Boundary
18from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
19     import File_boundary
20from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
21     import AWI_boundary
22from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
23     import Dirichlet_boundary
24from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
25     import Time_boundary
26from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
27     import Transmissive_boundary
28
29from anuga.abstract_2d_finite_volumes.pmesh2domain import pmesh_to_domain
30from anuga.abstract_2d_finite_volumes.region\
31     import Set_region as region_set_region
32from anuga.utilities.polygon import inside_polygon
33from anuga.abstract_2d_finite_volumes.util import get_textual_float
34from quantity import Quantity
35
36import Numeric as num
37
38
39##
40# @brief Generic Domain class
41class Domain:
42
43    ##
44    # @brief Generic computational Domain constructor.
45    # @param source Name of mesh file or coords of mesh vertices.
46    # @param triangles Mesh connectivity (see mesh.py for more information).
47    # @param boundary (see mesh.py for more information)
48    # @param conserved_quantities List of names of quantities to be conserved.
49    # @param other_quantities List of names of other quantities.
50    # @param tagged_elements ??
51    # @param geo_reference ??
52    # @param use_inscribed_circle ??
53    # @param mesh_filename ??
54    # @param use_cache ??
55    # @param verbose True if this method is to be verbose.
56    # @param full_send_dict ??
57    # @param ghost_recv_dict ??
58    # @param processor ??
59    # @param numproc ??
60    # @param number_of_full_nodes ??
61    # @param number_of_full_triangles ??
62    def __init__(self, source=None,
63                       triangles=None,
64                       boundary=None,
65                       conserved_quantities=None,
66                       other_quantities=None,
67                       tagged_elements=None,
68                       geo_reference=None,
69                       use_inscribed_circle=False,
70                       mesh_filename=None,
71                       use_cache=False,
72                       verbose=False,
73                       full_send_dict=None,
74                       ghost_recv_dict=None,
75                       processor=0,
76                       numproc=1,
77                       number_of_full_nodes=None,
78                       number_of_full_triangles=None):
79
80        """Instantiate generic computational Domain.
81
82        Input:
83          source:    Either a mesh filename or coordinates of mesh vertices.
84                     If it is a filename values specified for triangles will
85                     be overridden.
86          triangles: Mesh connectivity (see mesh.py for more information)
87          boundary:  See mesh.py for more information
88
89          conserved_quantities: List of quantity names entering the
90                                conservation equations
91          other_quantities:     List of other quantity names
92
93          tagged_elements:
94          ...
95        """
96
97
98        number_of_full_nodes=None
99        number_of_full_triangles=None
100       
101        # Determine whether source is a mesh filename or coordinates
102        if type(source) == types.StringType:
103            mesh_filename = source
104        else:
105            coordinates = source
106
107        # In case a filename has been specified, extract content
108        if mesh_filename is not None:
109            coordinates, triangles, boundary, vertex_quantity_dict, \
110                         tagged_elements, geo_reference = \
111                         pmesh_to_domain(file_name=mesh_filename,
112                                         use_cache=use_cache,
113                                         verbose=verbose)
114
115        # Initialise underlying mesh structure
116        self.mesh = Mesh(coordinates, triangles,
117                         boundary=boundary,
118                         tagged_elements=tagged_elements,
119                         geo_reference=geo_reference,
120                         use_inscribed_circle=use_inscribed_circle,
121                         number_of_full_nodes=number_of_full_nodes,
122                         number_of_full_triangles=number_of_full_triangles,
123                         verbose=verbose)
124
125        # Expose Mesh attributes (FIXME: Maybe turn into methods)
126        self.triangles = self.mesh.triangles       
127        self.centroid_coordinates = self.mesh.centroid_coordinates
128        self.vertex_coordinates = self.mesh.vertex_coordinates
129        self.boundary = self.mesh.boundary
130        self.neighbours = self.mesh.neighbours
131        self.surrogate_neighbours = self.mesh.surrogate_neighbours
132        self.neighbour_edges = self.mesh.neighbour_edges
133        self.normals = self.mesh.normals
134        self.edgelengths = self.mesh.edgelengths
135        self.radii = self.mesh.radii
136        self.areas = self.mesh.areas
137
138        self.number_of_boundaries = self.mesh.number_of_boundaries
139        self.number_of_full_nodes = self.mesh.number_of_full_nodes
140        self.number_of_full_triangles = self.mesh.number_of_full_triangles
141        self.number_of_triangles_per_node = self.mesh.number_of_triangles_per_node
142
143        self.vertex_value_indices = self.mesh.vertex_value_indices
144        self.number_of_triangles = self.mesh.number_of_triangles
145
146        self.geo_reference = self.mesh.geo_reference
147
148        if verbose: print 'Initialising Domain'
149
150        # List of quantity names entering the conservation equations
151        if conserved_quantities is None:
152            self.conserved_quantities = []
153        else:
154            self.conserved_quantities = conserved_quantities
155
156        # List of other quantity names
157        if other_quantities is None:
158            self.other_quantities = []
159        else:
160            self.other_quantities = other_quantities
161
162        # Build dictionary of Quantity instances keyed by quantity names
163        self.quantities = {}
164
165        # FIXME: remove later - maybe OK, though....
166        for name in self.conserved_quantities:
167            self.quantities[name] = Quantity(self)
168        for name in self.other_quantities:
169            self.quantities[name] = Quantity(self)
170
171        # Create an empty list for explicit forcing terms
172        self.forcing_terms = []
173
174        # Setup the ghost cell communication
175        if full_send_dict is None:
176            self.full_send_dict = {}
177        else:
178            self.full_send_dict = full_send_dict
179
180        # List of other quantity names
181        if ghost_recv_dict is None:
182            self.ghost_recv_dict = {}
183        else:
184            self.ghost_recv_dict = ghost_recv_dict
185
186        self.processor = processor
187        self.numproc = numproc
188
189        # Setup Communication Buffers
190        if verbose: print 'Domain: Set up communication buffers (parallel)'
191        self.nsys = len(self.conserved_quantities)
192        for key in self.full_send_dict:
193            buffer_shape = self.full_send_dict[key][0].shape[0]
194            self.full_send_dict[key].append(num.zeros((buffer_shape, self.nsys),
195                                                      num.Float))
196
197        for key in self.ghost_recv_dict:
198            buffer_shape = self.ghost_recv_dict[key][0].shape[0]
199            self.ghost_recv_dict[key].append(num.zeros((buffer_shape, self.nsys),
200                                             num.Float))
201
202        # Setup cell full flag
203        # =1 for full
204        # =0 for ghost
205        N = len(self) #number_of_elements
206        self.number_of_elements = N
207        self.tri_full_flag = num.ones(N, num.Int)
208        for i in self.ghost_recv_dict.keys():
209            for id in self.ghost_recv_dict[i][0]:
210                self.tri_full_flag[id] = 0
211
212        # Test the assumption that all full triangles are store before
213        # the ghost triangles.
214        if not num.allclose(self.tri_full_flag[:self.number_of_full_nodes], 1):
215            if self.numproc>1:
216                print ('WARNING:  '
217                       'Not all full triangles are store before ghost triangles')
218
219        # Defaults
220        from anuga.config import max_smallsteps, beta_w, epsilon
221        from anuga.config import CFL
222        from anuga.config import timestepping_method
223        from anuga.config import protect_against_isolated_degenerate_timesteps
224        from anuga.config import default_order
225
226        self.beta_w = beta_w
227        self.epsilon = epsilon
228        self.protect_against_isolated_degenerate_timesteps = \
229                        protect_against_isolated_degenerate_timesteps
230
231        self.set_default_order(default_order)
232
233        self.smallsteps = 0
234        self.max_smallsteps = max_smallsteps
235        self.number_of_steps = 0
236        self.number_of_first_order_steps = 0
237        self.CFL = CFL
238        self.set_timestepping_method(timestepping_method)
239        self.set_beta(beta_w)
240        self.boundary_map = None  # Will be populated by set_boundary
241
242        # Model time
243        self.time = 0.0
244        self.finaltime = None
245        self.min_timestep = self.max_timestep = 0.0
246        self.starttime = 0 # Physical starttime if any
247                           # (0 is 1 Jan 1970 00:00:00)
248        self.timestep = 0.0
249        self.flux_timestep = 0.0
250
251        self.last_walltime = walltime()
252
253        # Monitoring
254        self.quantities_to_be_monitored = None
255        self.monitor_polygon = None
256        self.monitor_time_interval = None
257        self.monitor_indices = None
258
259        # Checkpointing and storage
260        from anuga.config import default_datadir
261
262        self.datadir = default_datadir
263        self.simulation_name = 'domain'
264        self.checkpoint = False
265
266        # To avoid calculating the flux across each edge twice, keep an integer
267        # (boolean) array, to be used during the flux calculation.
268        N = len(self) # Number_of_triangles
269        self.already_computed_flux = num.zeros((N, 3), num.Int)
270
271        # Storage for maximal speeds computed for each triangle by
272        # compute_fluxes.
273        # This is used for diagnostics only (reset at every yieldstep)
274        self.max_speed = num.zeros(N, num.Float)
275
276        if mesh_filename is not None:
277            # If the mesh file passed any quantity values,
278            # initialise with these values.
279            if verbose: print 'Domain: Initialising quantity values'
280            self.set_quantity_vertices_dict(vertex_quantity_dict)
281
282        if verbose: print 'Domain: Done'
283
284    ######
285    # Expose underlying Mesh functionality
286    ######
287
288    def __len__(self):
289        return len(self.mesh)
290
291    def get_centroid_coordinates(self, *args, **kwargs):
292        return self.mesh.get_centroid_coordinates(*args, **kwargs)
293
294    def get_radii(self, *args, **kwargs):
295        return self.mesh.get_radii(*args, **kwargs)
296
297    def get_areas(self, *args, **kwargs):
298        return self.mesh.get_areas(*args, **kwargs)
299
300    def get_area(self, *args, **kwargs):
301        return self.mesh.get_area(*args, **kwargs)
302
303    def get_vertex_coordinates(self, *args, **kwargs):
304        return self.mesh.get_vertex_coordinates(*args, **kwargs)
305
306    def get_triangles(self, *args, **kwargs):
307        return self.mesh.get_triangles(*args, **kwargs)
308
309    def get_nodes(self, *args, **kwargs):
310        return self.mesh.get_nodes(*args, **kwargs)
311
312    def get_number_of_nodes(self, *args, **kwargs):
313        return self.mesh.get_number_of_nodes(*args, **kwargs)
314
315    def get_normal(self, *args, **kwargs):
316        return self.mesh.get_normal(*args, **kwargs)
317
318    def get_triangle_containing_point(self, *args, **kwargs):
319        return self.mesh.get_triangle_containing_point(*args, **kwargs)
320
321    def get_intersecting_segments(self, *args, **kwargs):
322        return self.mesh.get_intersecting_segments(*args, **kwargs)
323
324    def get_disconnected_triangles(self, *args, **kwargs):
325        return self.mesh.get_disconnected_triangles(*args, **kwargs)
326
327    def get_boundary_tags(self, *args, **kwargs):
328        return self.mesh.get_boundary_tags(*args, **kwargs)
329
330    def get_boundary_polygon(self, *args, **kwargs):
331        return self.mesh.get_boundary_polygon(*args, **kwargs)
332
333    # FIXME(Ole): This doesn't seem to be required
334    def get_number_of_triangles_per_node(self, *args, **kwargs):
335        return self.mesh.get_number_of_triangles_per_node(*args, **kwargs)
336
337    def get_triangles_and_vertices_per_node(self, *args, **kwargs):
338        return self.mesh.get_triangles_and_vertices_per_node(*args, **kwargs)
339
340    def get_interpolation_object(self, *args, **kwargs):
341        return self.mesh.get_interpolation_object(*args, **kwargs)
342
343    def get_tagged_elements(self, *args, **kwargs):
344        return self.mesh.get_tagged_elements(*args, **kwargs)
345
346    def get_lone_vertices(self, *args, **kwargs):
347        return self.mesh.get_lone_vertices(*args, **kwargs)
348
349    def get_unique_vertices(self, *args, **kwargs):
350        return self.mesh.get_unique_vertices(*args, **kwargs)
351
352    def get_georeference(self, *args, **kwargs):
353        return self.mesh.get_georeference(*args, **kwargs)
354
355    def set_georeference(self, *args, **kwargs):
356        self.mesh.set_georeference(*args, **kwargs)
357
358    def build_tagged_elements_dictionary(self, *args, **kwargs):
359        self.mesh.build_tagged_elements_dictionary(*args, **kwargs)
360
361    def statistics(self, *args, **kwargs):
362        return self.mesh.statistics(*args, **kwargs)
363       
364    def get_extent(self, *args, **kwargs):
365        return self.mesh.get_extent(*args, **kwargs)   
366
367    ##
368    # @brief Get conserved quantities for a volume.
369    # @param vol_id ID of the volume we want the conserved quantities for.
370    # @param vertex If specified, use as index for edge values.
371    # @param edge If specified, use as index for edge values.
372    # @return Vector of conserved quantities.
373    # @note If neither 'vertex' or 'edge' specified, use centroid values.
374    # @note If both 'vertex' and 'edge' specified, raise exception.
375    def get_conserved_quantities(self, vol_id,
376                                       vertex=None,
377                                       edge=None):
378        """Get conserved quantities at volume vol_id.
379
380        If vertex is specified use it as index for vertex values
381        If edge is specified use it as index for edge values
382        If neither are specified use centroid values
383        If both are specified an exeception is raised
384
385        Return value: Vector of length == number_of_conserved quantities
386        """
387
388        if not (vertex is None or edge is None):
389            msg = 'Values for both vertex and edge was specified.'
390            msg += 'Only one (or none) is allowed.'
391            raise Exception, msg
392
393        q = num.zeros(len(self.conserved_quantities), num.Float)
394
395        for i, name in enumerate(self.conserved_quantities):
396            Q = self.quantities[name]
397            if vertex is not None:
398                q[i] = Q.vertex_values[vol_id, vertex]
399            elif edge is not None:
400                q[i] = Q.edge_values[vol_id, edge]
401            else:
402                q[i] = Q.centroid_values[vol_id]
403
404        return q
405
406    ##
407    # @brief Set the relative model time.
408    # @param time The new model time (seconds).
409    def set_time(self, time=0.0):
410        """Set the model time (seconds)."""
411
412        # FIXME: this is setting the relative time
413        # Note that get_time and set_time are now not symmetric
414
415        self.time = time
416
417    ##
418    # @brief Get the model time.
419    # @return The absolute model time (seconds).
420    def get_time(self):
421        """Get the absolute model time (seconds)."""
422
423        return self.time + self.starttime
424
425    ##
426    # @brief Set the default beta for limiting.
427    # @param beta The new beta value.
428    def set_beta(self, beta):
429        """Set default beta for limiting."""
430
431        self.beta = beta
432        for name in self.quantities:
433            Q = self.quantities[name]
434            Q.set_beta(beta)
435
436    ##
437    # @brief Get the beta value used for limiting.
438    # @return The beta value used for limiting.
439    def get_beta(self):
440        """Get default beta for limiting."""
441
442        return self.beta
443
444    ##
445    # @brief Set default (spatial) order.
446    # @param n The new spatial order value.
447    # @note If 'n' is not 1 or 2, raise exception.
448    def set_default_order(self, n):
449        """Set default (spatial) order to either 1 or 2."""
450
451        msg = 'Default order must be either 1 or 2. I got %s' % n
452        assert n in [1,2], msg
453
454        self.default_order = n
455        self._order_ = self.default_order
456
457    ##
458    # @brief Set values of named quantities.
459    # @param quantity_dict Dictionary containing name/value pairs.
460    def set_quantity_vertices_dict(self, quantity_dict):
461        """Set values for named quantities.
462        Supplied dictionary contains name/value pairs:
463
464        name:  Name of quantity
465        value: Compatible list, Numeric array, const or function (see below)
466
467        The values will be stored in elements following their internal ordering.
468        """
469
470        # FIXME: Could we name this a bit more intuitively
471        # E.g. set_quantities_from_dictionary
472        for key in quantity_dict.keys():
473            self.set_quantity(key, quantity_dict[key], location='vertices')
474
475    ##
476    # @brief Set value(s) for a named quantity.
477    # @param name Name of quantity to be updated.
478    # @param args Positional args.
479    # @param kwargs Keyword args.
480    # @note If 'kwargs' dict has 'expression' key, evaluate expression.
481    def set_quantity(self, name,
482                           *args, **kwargs):
483        """Set values for named quantity
484
485        One keyword argument is documented here:
486        expression = None, # Arbitrary expression
487
488        expression:
489          Arbitrary expression involving quantity names
490
491        See Quantity.set_values for further documentation.
492        """
493
494        # Do the expression stuff
495        if kwargs.has_key('expression'):
496            expression = kwargs['expression']
497            del kwargs['expression']
498
499            Q = self.create_quantity_from_expression(expression)
500            kwargs['quantity'] = Q
501
502        # Assign values
503        self.quantities[name].set_values(*args, **kwargs)
504
505    ##
506    # @brief Add to a named quantity value.
507    # @param name Name of quantity to be added to.
508    # @param args Positional args.
509    # @param kwargs Keyword args.
510    # @note If 'kwargs' dict has 'expression' key, evaluate expression.
511    def add_quantity(self, name,
512                           *args, **kwargs):
513        """Add values to a named quantity
514
515        E.g add_quantity('elevation', X)
516
517        Option are the same as in set_quantity.
518        """
519
520        # Do the expression stuff
521        if kwargs.has_key('expression'):
522            expression = kwargs['expression']
523            Q2 = self.create_quantity_from_expression(expression)
524        else:
525            # Create new temporary quantity
526            Q2 = Quantity(self)
527
528            # Assign specified values to temporary quantity
529            Q2.set_values(*args, **kwargs)
530
531        # Add temporary quantity to named quantity
532        Q1 = self.get_quantity(name)
533        self.set_quantity(name, Q1 + Q2)
534
535    ##
536    # @brief Get list of quantity names for the Domain.
537    # @return List of quantity names.
538    def get_quantity_names(self):
539        """Get a list of all the quantity names that this domain is aware of.
540        Any value in the result should be a valid input to get_quantity.
541        """
542
543        return self.quantities.keys()
544
545    ##
546    # @brief Get a quantity object.
547    # @param name Name of the quantity value.
548    # @param location ??
549    # @param indices ??
550    # @return The quantity value object.
551    # @note 'location' and 'indices' are unused.
552    def get_quantity(self, name,
553                           location='vertices',
554                           indices = None):
555        """Get pointer to quantity object.
556
557        name: Name of quantity
558
559        See methods inside the quantity object for more options
560
561        FIXME: clean input args
562        """
563
564        return self.quantities[name] #.get_values( location, indices = indices)
565
566    ##
567    # @brief Create a quantity value from an expression.
568    # @param expression The expression (string) to be evaluated.
569    # @return The expression value, evaluated from this Domain's quantities.
570    # @note Valid expression operators are as defined in class Quantity.
571    def create_quantity_from_expression(self, expression):
572        """Create new quantity from other quantities using arbitrary expression.
573
574        Combine existing quantities in domain using expression and return
575        result as a new quantity.
576
577        Note, the new quantity could e.g. be used in set_quantity
578
579        Valid expressions are limited to operators defined in class Quantity
580
581        Examples creating derived quantities:
582            Depth = domain.create_quantity_from_expression('stage-elevation')
583            exp = '(xmomentum*xmomentum + ymomentum*ymomentum)**0.5'
584            Absolute_momentum = domain.create_quantity_from_expression(exp)
585        """
586
587        from anuga.abstract_2d_finite_volumes.util import\
588             apply_expression_to_dictionary
589
590        return apply_expression_to_dictionary(expression, self.quantities)
591
592    ##
593    # @brief Associate boundary objects with tagged boundary segments.
594    # @param boundary_map A dict of boundary objects keyed by symbolic tags to
595    #                     matched against tags in the internal dictionary
596    #                     self.boundary.
597    def set_boundary(self, boundary_map):
598        """Associate boundary objects with tagged boundary segments.
599
600        Input boundary_map is a dictionary of boundary objects keyed
601        by symbolic tags to matched against tags in the internal dictionary
602        self.boundary.
603
604        As result one pointer to a boundary object is stored for each vertex
605        in the list self.boundary_objects.
606        More entries may point to the same boundary object
607
608        Schematically the mapping is from two dictionaries to one list
609        where the index is used as pointer to the boundary_values arrays
610        within each quantity.
611
612        self.boundary:          (vol_id, edge_id): tag
613        boundary_map (input):   tag: boundary_object
614        ----------------------------------------------
615        self.boundary_objects:  ((vol_id, edge_id), boundary_object)
616
617        Pre-condition:
618          self.boundary has been built.
619
620        Post-condition:
621          self.boundary_objects is built
622
623        If a tag from the domain doesn't appear in the input dictionary an
624        exception is raised.
625        However, if a tag is not used to the domain, no error is thrown.
626        FIXME: This would lead to implementation of a default boundary condition
627
628        Note: If a segment is listed in the boundary dictionary and if it is
629        not None, it *will* become a boundary - even if there is a neighbouring
630        triangle.  This would be the case for internal boundaries.
631
632        Boundary objects that are None will be skipped.
633
634        If a boundary_map has already been set (i.e. set_boundary has been
635        called before), the old boundary map will be updated with new values.
636        The new map need not define all boundary tags, and can thus change only
637        those that are needed.
638
639        FIXME: If set_boundary is called multiple times and if Boundary
640        object is changed into None, the neighbour structure will not be
641        restored!!!
642        """
643
644        if self.boundary_map is None:
645            # This the first call to set_boundary. Store
646            # map for later updates and for use with boundary_stats.
647            self.boundary_map = boundary_map
648        else:
649            # This is a modification of an already existing map
650            # Update map an proceed normally
651            for key in boundary_map.keys():
652                self.boundary_map[key] = boundary_map[key]
653
654        # FIXME (Ole): Try to remove the sorting and fix test_mesh.py
655        x = self.boundary.keys()
656        x.sort()
657
658        # Loop through edges that lie on the boundary and associate them with
659        # callable boundary objects depending on their tags
660        self.boundary_objects = []
661        for k, (vol_id, edge_id) in enumerate(x):
662            tag = self.boundary[(vol_id, edge_id)]
663
664            if self.boundary_map.has_key(tag):
665                B = self.boundary_map[tag]  # Get callable boundary object
666
667                if B is not None:
668                    self.boundary_objects.append(((vol_id, edge_id), B))
669                    self.neighbours[vol_id, edge_id] = \
670                                        -len(self.boundary_objects)
671                else:
672                    pass
673                    #FIXME: Check and perhaps fix neighbour structure
674            else:
675                msg = 'ERROR (domain.py): Tag "%s" has not been ' %tag
676                msg += 'bound to a boundary object.\n'
677                msg += 'All boundary tags defined in domain must appear '
678                msg += 'in set_boundary.\n'
679                msg += 'The tags are: %s' %self.get_boundary_tags()
680                raise Exception, msg
681
682    ##
683    # @brief Set quantities based on a regional tag.
684    # @param args
685    # @param kwargs
686    def set_region(self, *args, **kwargs):
687        """Set quantities based on a regional tag.
688
689        It is most often called with the following parameters;
690        (self, tag, quantity, X, location='vertices')
691        tag:      the name of the regional tag used to specify the region
692        quantity: Name of quantity to change
693        X:        const or function - how the quantity is changed
694        location: Where values are to be stored.
695            Permissible options are: vertices, centroid and unique vertices
696
697        A callable region class or a list of callable region classes
698        can also be passed into this function.
699        """
700
701        if len(args) == 1:
702            self._set_region(*args, **kwargs)
703        else:
704            # Assume it is arguments for the region.set_region function
705            func = region_set_region(*args, **kwargs)
706            self._set_region(func)
707
708    ##
709    # @brief ??
710    # @param functions A list or tuple of ??
711    def _set_region(self, functions):
712        # coerce to an iterable (list or tuple)
713        if type(functions) not in [types.ListType, types.TupleType]:
714            functions = [functions]
715
716        # The order of functions in the list is used.
717        tagged_elements = self.get_tagged_elements()
718        for function in functions:
719            for tag in tagged_elements.keys():
720                function(tag, tagged_elements[tag], self)
721
722    ##
723    # @brief Specify the quantities which will be monitored for extrema.
724    # @param q Single or list of quantity names to monitor.
725    # @param polygon If specified, monitor only triangles inside polygon.
726    # @param time_interval If specified, monitor only timesteps inside interval.
727    # @note If 'q' is None, do no monitoring.
728    def set_quantities_to_be_monitored(self, q,
729                                             polygon=None,
730                                             time_interval=None):
731        """Specify which quantities will be monitored for extrema.
732
733        q must be either:
734          - the name of a quantity or derived quantity such as 'stage-elevation'
735          - a list of quantity names
736          - None
737
738        In the two first cases, the named quantities will be monitored at
739        each internal timestep
740
741        If q is None, monitoring will be switched off altogether.
742
743        polygon (if specified) will only monitor triangles inside polygon.
744        If omitted all triangles will be included.
745
746        time_interval, if specified, will restrict monitoring to time steps in
747        that interval. If omitted all timesteps will be included.
748        """
749
750        from anuga.abstract_2d_finite_volumes.util import\
751             apply_expression_to_dictionary
752
753        if q is None:
754            self.quantities_to_be_monitored = None
755            self.monitor_polygon = None
756            self.monitor_time_interval = None
757            self.monitor_indices = None
758            return
759
760        # coerce 'q' to a list if it's a string
761        if isinstance(q, basestring):
762            q = [q]
763
764        # Check correctness and initialise
765        self.quantities_to_be_monitored = {}
766        for quantity_name in q:
767            msg = 'Quantity %s is not a valid conserved quantity' \
768                      % quantity_name
769
770            if not quantity_name in self.quantities:
771                # See if this expression is valid
772                apply_expression_to_dictionary(quantity_name, self.quantities)
773
774            # Initialise extrema information
775            info_block = {'min': None,          # Min value
776                          'max': None,          # Max value
777                          'min_location': None, # Argmin (x, y)
778                          'max_location': None, # Argmax (x, y)
779                          'min_time': None,     # Argmin (t)
780                          'max_time': None}     # Argmax (t)
781
782            self.quantities_to_be_monitored[quantity_name] = info_block
783
784        if polygon is not None:
785            # Check input
786            if isinstance(polygon, basestring):
787                # Check if multiple quantities were accidentally
788                # given as separate argument rather than a list.
789                msg = ('Multiple quantities must be specified in a list. '
790                       'Not as multiple arguments. '
791                       'I got "%s" as a second argument') % polygon
792
793                if polygon in self.quantities:
794                    raise Exception, msg
795
796                try:
797                    apply_expression_to_dictionary(polygon, self.quantities)
798                except:
799                    # At least polygon wasn't expression involving quantitites
800                    pass
801                else:
802                    raise Exception, msg
803
804                # In any case, we don't allow polygon to be a string
805                msg = ('argument "polygon" must not be a string: '
806                       'I got polygon="%s"') % polygon
807                raise Exception, msg
808
809            # Get indices for centroids that are inside polygon
810            points = self.get_centroid_coordinates(absolute=True)
811            self.monitor_indices = inside_polygon(points, polygon)
812
813        if time_interval is not None:
814            assert len(time_interval) == 2
815
816        self.monitor_polygon = polygon
817        self.monitor_time_interval = time_interval
818
819    ##
820    # @brief Check Domain integrity.
821    # @note Raises an exception if integrity breached.
822    def check_integrity(self):
823        self.mesh.check_integrity()
824
825        for quantity in self.conserved_quantities:
826            msg = 'Conserved quantities must be a subset of all quantities'
827            assert quantity in self.quantities, msg
828
829        ##assert hasattr(self, 'boundary_objects')
830
831    ##
832    # @brief Print timestep stats to stdout.
833    # @param track_speeds If True, print smallest track speed.
834    def write_time(self, track_speeds=False):
835        print self.timestepping_statistics(track_speeds)
836
837    ##
838    # @brief Get timestepping stats string.
839    # @param track_speeds If True, report location of smallest timestep.
840    # @param triangle_id If specified, use specific triangle.
841    # @return A string containing timestep stats.
842    def timestepping_statistics(self, track_speeds=False,
843                                      triangle_id=None):
844        """Return string with time stepping statistics
845
846        Optional boolean keyword track_speeds decides whether to report
847        location of smallest timestep as well as a histogram and percentile
848        report.
849
850        Optional keyword triangle_id can be used to specify a particular
851        triangle rather than the one with the largest speed.
852        """
853
854        from anuga.utilities.numerical_tools import histogram, create_bins
855
856        # qwidth determines the the width of the text field used for quantities
857        qwidth = self.qwidth = 12
858
859        msg = ''
860        #if self.min_timestep == self.max_timestep:
861        #    msg += 'Time = %.4f, delta t = %.8f, steps=%d (%d)'\
862        #           %(self.time, self.min_timestep, self.number_of_steps,
863        #             self.number_of_first_order_steps)
864        #elif self.min_timestep > self.max_timestep:
865        #    msg += 'Time = %.4f, steps=%d (%d)'\
866        #           %(self.time, self.number_of_steps,
867        #             self.number_of_first_order_steps)
868        #else:
869        #    msg += 'Time = %.4f, delta t in [%.8f, %.8f], steps=%d (%d)'\
870        #           %(self.time, self.min_timestep,
871        #             self.max_timestep, self.number_of_steps,
872        #             self.number_of_first_order_steps)
873
874        model_time = self.get_time()
875        if self.min_timestep == self.max_timestep:
876            msg += 'Time = %.4f, delta t = %.8f, steps=%d' \
877                       % (model_time, self.min_timestep, self.number_of_steps)
878        elif self.min_timestep > self.max_timestep:
879            msg += 'Time = %.4f, steps=%d' \
880                       % (model_time, self.number_of_steps)
881        else:
882            msg += 'Time = %.4f, delta t in [%.8f, %.8f], steps=%d' \
883                       % (model_time, self.min_timestep,
884                          self.max_timestep, self.number_of_steps)
885
886        msg += ' (%ds)' % (walltime() - self.last_walltime)
887        self.last_walltime = walltime()
888
889        if track_speeds is True:
890            msg += '\n'
891
892            # Setup 10 bins for speed histogram
893            bins = create_bins(self.max_speed, 10)
894            hist = histogram(self.max_speed, bins)
895
896            msg += '------------------------------------------------\n'
897            msg += '  Speeds in [%f, %f]\n' % (min(self.max_speed),
898                                               max(self.max_speed))
899            msg += '  Histogram:\n'
900
901            hi = bins[0]
902            for i, count in enumerate(hist):
903                lo = hi
904                if i+1 < len(bins):
905                    # Open upper interval
906                    hi = bins[i+1]
907                    msg += '    [%f, %f[: %d\n' % (lo, hi, count)
908                else:
909                    # Closed upper interval
910                    hi = max(self.max_speed)
911                    msg += '    [%f, %f]: %d\n' % (lo, hi, count)
912
913            N = len(self.max_speed)
914            if N > 10:
915                msg += '  Percentiles (10%):\n'
916                speed = self.max_speed.tolist()
917                speed.sort()
918
919                k = 0
920                lower = min(speed)
921                for i, a in enumerate(speed):
922                    if i % (N/10) == 0 and i != 0:
923                        # For every 10% of the sorted speeds
924                        msg += '    %d speeds in [%f, %f]\n' % (i-k, lower, a)
925                        lower = a
926                        k = i
927
928                msg += '    %d speeds in [%f, %f]\n'\
929                           % (N-k, lower, max(speed))
930
931            # Find index of largest computed flux speed
932            if triangle_id is None:
933                k = self.k = num.argmax(self.max_speed)
934            else:
935                errmsg = 'Triangle_id %d does not exist in mesh: %s' \
936                             % (triangle_id, str(self))
937                assert 0 <= triangle_id < len(self), errmsg
938                k = self.k = triangle_id
939
940            x, y = self.get_centroid_coordinates(absolute=True)[k]
941            radius = self.get_radii()[k]
942            area = self.get_areas()[k]
943            max_speed = self.max_speed[k]
944
945            msg += '  Triangle #%d with centroid (%.4f, %.4f), ' % (k, x, y)
946            msg += 'area = %.4f and radius = %.4f ' % (area, radius)
947            if triangle_id is None:
948                msg += 'had the largest computed speed: %.6f m/s ' % (max_speed)
949            else:
950                msg += 'had computed speed: %.6f m/s ' % (max_speed)
951
952            if max_speed > 0.0:
953                msg += '(timestep=%.6f)\n' % (radius/max_speed)
954            else:
955                msg += '(timestep=%.6f)\n' % (0)
956
957            # Report all quantity values at vertices, edges and centroid
958            msg += '    Quantity'
959            msg += '------------\n'
960            for name in self.quantities:
961                q = self.quantities[name]
962
963                V = q.get_values(location='vertices', indices=[k])[0]
964                E = q.get_values(location='edges', indices=[k])[0]
965                C = q.get_values(location='centroids', indices=[k])
966
967                s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n' \
968                         % (name.ljust(qwidth), V[0], V[1], V[2])
969
970                s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n' \
971                         % (name.ljust(qwidth), E[0], E[1], E[2])
972
973                s += '    %s: centroid_value = %.4f\n' \
974                         % (name.ljust(qwidth), C[0])
975
976                msg += s
977
978        return msg
979
980    ##
981    # @brief Print boundary forcing stats at each timestep to stdout.
982    # @param quantities A name or list of names of quantities to report.
983    # @param tags A name or list of names of tags to report.
984    def write_boundary_statistics(self, quantities=None, tags=None):
985        print self.boundary_statistics(quantities, tags)
986
987    # @brief Get a string containing boundary forcing stats at each timestep.
988    # @param quantities A name or list of names of quantities to report.
989    # @param tags A name or list of names of tags to report.
990    # @note If 'quantities' is None, report all.  Same for 'tags'.
991    def boundary_statistics(self, quantities=None,
992                                  tags=None):
993        """Output statistics about boundary forcing at each timestep
994
995        Input:
996          quantities: either None, a string or a list of strings naming the
997                      quantities to be reported
998          tags:       either None, a string or a list of strings naming the
999                      tags to be reported
1000
1001        Example output:
1002        Tag 'wall':
1003            stage in [2, 5.5]
1004            xmomentum in []
1005            ymomentum in []
1006        Tag 'ocean'
1007
1008        If quantities are specified only report on those. Otherwise take all
1009        conserved quantities.
1010        If tags are specified only report on those, otherwise take all tags.
1011        """
1012
1013        import types, string
1014
1015        # Input checks
1016        if quantities is None:
1017            quantities = self.conserved_quantities
1018        elif type(quantities) == types.StringType:
1019            quantities = [quantities] #Turn it into a list
1020
1021        msg = ('Keyword argument quantities must be either None, '
1022               'string or list. I got %s') % str(quantities)
1023        assert type(quantities) == types.ListType, msg
1024
1025        if tags is None:
1026            tags = self.get_boundary_tags()
1027        elif type(tags) == types.StringType:
1028            tags = [tags] #Turn it into a list
1029
1030        msg = ('Keyword argument tags must be either None, '
1031               'string or list. I got %s') % str(tags)
1032        assert type(tags) == types.ListType, msg
1033
1034        # Determine width of longest quantity name (for cosmetic purposes)
1035        maxwidth = 0
1036        for name in quantities:
1037            w = len(name)
1038            if w > maxwidth:
1039                maxwidth = w
1040
1041        # Output statistics
1042        msg = 'Boundary values at time %.4f:\n' % self.get_time()
1043        for tag in tags:
1044            msg += '    %s:\n' % tag
1045
1046            for name in quantities:
1047                q = self.quantities[name]
1048
1049                # Find range of boundary values for tag and q
1050                maxval = minval = None
1051                for i, ((vol_id,edge_id),B) in enumerate(self.boundary_objects):
1052                    if self.boundary[(vol_id, edge_id)] == tag:
1053                        v = q.boundary_values[i]
1054                        if minval is None or v < minval: minval = v
1055                        if maxval is None or v > maxval: maxval = v
1056
1057                if minval is None or maxval is None:
1058                    msg += ('        Sorry no information available about'
1059                            ' tag %s and quantity %s\n') % (tag, name)
1060                else:
1061                    msg += '        %s in [%12.8f, %12.8f]\n' \
1062                               % (string.ljust(name, maxwidth), minval, maxval)
1063
1064        return msg
1065
1066    ##
1067    # @brief Update extrema if requested by set_quantities_to_be_monitored.
1068    def update_extrema(self):
1069        """Update extrema if requested by set_quantities_to_be_monitored.
1070        This data is used for reporting e.g. by running
1071        print domain.quantity_statistics()
1072        and may also stored in output files (see data_manager in shallow_water)
1073        """
1074
1075        # Define a tolerance for extremum computations
1076        from anuga.config import single_precision as epsilon
1077
1078        if self.quantities_to_be_monitored is None:
1079            return
1080
1081        # Observe time interval restriction if any
1082        if self.monitor_time_interval is not None and\
1083               (self.time < self.monitor_time_interval[0] or\
1084               self.time > self.monitor_time_interval[1]):
1085            return
1086
1087        # Update extrema for each specified quantity subject to
1088        # polygon restriction (via monitor_indices).
1089        for quantity_name in self.quantities_to_be_monitored:
1090
1091            if quantity_name in self.quantities:
1092                Q = self.get_quantity(quantity_name)
1093            else:
1094                Q = self.create_quantity_from_expression(quantity_name)
1095
1096            info_block = self.quantities_to_be_monitored[quantity_name]
1097
1098            # Update maximum
1099            # (n > None is always True, but we check explicitly because
1100            # of the epsilon)
1101            maxval = Q.get_maximum_value(self.monitor_indices)
1102            if info_block['max'] is None or \
1103                   maxval > info_block['max'] + epsilon:
1104                info_block['max'] = maxval
1105                maxloc = Q.get_maximum_location()
1106                info_block['max_location'] = maxloc
1107                info_block['max_time'] = self.time
1108
1109            # Update minimum
1110            minval = Q.get_minimum_value(self.monitor_indices)
1111            if info_block['min'] is None or \
1112                   minval < info_block['min'] - epsilon:
1113                info_block['min'] = minval
1114                minloc = Q.get_minimum_location()
1115                info_block['min_location'] = minloc
1116                info_block['min_time'] = self.time
1117
1118    ##
1119    # @brief Return string with statistics about quantities
1120    # @param precision A format string to use for float values.
1121    # @return The stats string.
1122    def quantity_statistics(self, precision='%.4f'):
1123        """Return string with statistics about quantities for
1124        printing or logging
1125
1126        Quantities reported are specified through method
1127
1128           set_quantities_to_be_monitored
1129        """
1130
1131        maxlen = 128 # Max length of polygon string representation
1132
1133        # Output statistics
1134        msg = 'Monitored quantities at time %.4f:\n' % self.get_time()
1135        if self.monitor_polygon is not None:
1136            p_str = str(self.monitor_polygon)
1137            msg += '- Restricted by polygon: %s' % p_str[:maxlen]
1138            if len(p_str) >= maxlen:
1139                msg += '...\n'
1140            else:
1141                msg += '\n'
1142
1143        if self.monitor_time_interval is not None:
1144            msg += '- Restricted by time interval: %s\n' \
1145                       % str(self.monitor_time_interval)
1146            time_interval_start = self.monitor_time_interval[0]
1147        else:
1148            time_interval_start = 0.0
1149
1150        for quantity_name, info in self.quantities_to_be_monitored.items():
1151            msg += '    %s:\n' % quantity_name
1152
1153            msg += '      values since time = %.2f in [%s, %s]\n' \
1154                       % (time_interval_start,
1155                          get_textual_float(info['min'], precision),
1156                          get_textual_float(info['max'], precision))
1157
1158            msg += '      minimum attained at time = %s, location = %s\n' \
1159                       % (get_textual_float(info['min_time'], precision),
1160                          get_textual_float(info['min_location'], precision))
1161
1162            msg += '      maximum attained at time = %s, location = %s\n' \
1163                       % (get_textual_float(info['max_time'], precision),
1164                          get_textual_float(info['max_location'], precision))
1165
1166        return msg
1167
1168    ##
1169    # @brief Get the timestep method.
1170    # @return The timestep method. One of 'euler', 'rk2' or 'rk3'.
1171    def get_timestepping_method(self):
1172        return self.timestepping_method
1173
1174    ##
1175    # @brief Set the tmestep method to be used.
1176    # @param timestepping_method One of 'euler', 'rk2' or 'rk3'.
1177    # @note Raises exception of method not known.
1178    def set_timestepping_method(self, timestepping_method):
1179        if timestepping_method in ['euler', 'rk2', 'rk3']:
1180            self.timestepping_method = timestepping_method
1181            return
1182
1183        msg = '%s is an incorrect timestepping type' % timestepping_method
1184        raise Exception, msg
1185
1186    ##
1187    # @brief Get the Domain simulation name.
1188    # @return The simulation name string.
1189    def get_name(self):
1190        return self.simulation_name
1191
1192    ##
1193    # @brief Set the simulation name.
1194    # @param name The name of the simulation.
1195    # @note The simulation name is also used for the output .sww file.
1196    def set_name(self, name):
1197        """Assign a name to this simulation.
1198        This will be used to identify the output sww file.
1199        """
1200
1201        # remove any '.sww' end
1202        if name.endswith('.sww'):
1203            name = name[:-4]
1204
1205        self.simulation_name = name
1206
1207    ##
1208    # @brief Get data directory path.
1209    # @return The data directory path string.
1210    def get_datadir(self):
1211        return self.datadir
1212
1213    ##
1214    # @brief Set data directory path.
1215    # @param name The data directory path string.
1216    def set_datadir(self, name):
1217        self.datadir = name
1218
1219    ##
1220    # @brief Get the start time value.
1221    # @return The start time value (float).
1222    def get_starttime(self):
1223        return self.starttime
1224
1225    ##
1226    # @brief Set the start time value.
1227    # @param time The start time value.
1228    def set_starttime(self, time):
1229        self.starttime = float(time)
1230
1231################################################################################
1232# Main components of evolve
1233################################################################################
1234
1235    ##
1236    # @brief Evolve the model through time.
1237    # @param yieldstep Interval between yields where results are stored, etc.
1238    # @param finaltime Time where simulation should end.
1239    # @param duration Duration of simulation.
1240    # @param skip_initial_step If True, skip the first yield step.
1241    def evolve(self, yieldstep=None,
1242                     finaltime=None,
1243                     duration=None,
1244                     skip_initial_step=False):
1245        """Evolve model through time starting from self.starttime.
1246
1247        yieldstep: Interval between yields where results are stored,
1248                   statistics written and domain inspected or
1249                   possibly modified. If omitted the internal predefined
1250                   max timestep is used.
1251                   Internally, smaller timesteps may be taken.
1252
1253        duration: Duration of simulation
1254
1255        finaltime: Time where simulation should end. This is currently
1256        relative time.  So it's the same as duration.
1257
1258        If both duration and finaltime are given an exception is thrown.
1259
1260        skip_initial_step: Boolean flag that decides whether the first
1261        yield step is skipped or not. This is useful for example to avoid
1262        duplicate steps when multiple evolve processes are dove tailed.
1263
1264        Evolve is implemented as a generator and is to be called as such, e.g.
1265
1266        for t in domain.evolve(yieldstep, finaltime):
1267            <Do something with domain and t>
1268
1269        All times are given in seconds
1270        """
1271
1272        from anuga.config import min_timestep, max_timestep, epsilon
1273
1274        # FIXME: Maybe lump into a larger check prior to evolving
1275        msg = ('Boundary tags must be bound to boundary objects before '
1276               'evolving system, '
1277               'e.g. using the method set_boundary.\n'
1278               'This system has the boundary tags %s '
1279                   % self.get_boundary_tags())
1280        assert hasattr(self, 'boundary_objects'), msg
1281
1282        if yieldstep is None:
1283            yieldstep = max_timestep
1284        else:
1285            yieldstep = float(yieldstep)
1286
1287        self._order_ = self.default_order
1288
1289        if finaltime is not None and duration is not None:
1290            msg = 'Only one of finaltime and duration may be specified'
1291            raise Exception, msg
1292        else:
1293            if finaltime is not None:
1294                self.finaltime = float(finaltime)
1295            if duration is not None:
1296                self.finaltime = self.starttime + float(duration)
1297
1298        N = len(self)                             # Number of triangles
1299        self.yieldtime = self.time + yieldstep    # set next yield time
1300
1301        # Initialise interval of timestep sizes (for reporting only)
1302        self.min_timestep = max_timestep
1303        self.max_timestep = min_timestep
1304        self.number_of_steps = 0
1305        self.number_of_first_order_steps = 0
1306
1307        # Update ghosts
1308        self.update_ghosts()
1309
1310        # Initial update of vertex and edge values
1311        self.distribute_to_vertices_and_edges()
1312
1313        # Update extrema if necessary (for reporting)
1314        self.update_extrema()
1315
1316        # Initial update boundary values
1317        self.update_boundary()
1318
1319        # Or maybe restore from latest checkpoint
1320        if self.checkpoint is True:
1321            self.goto_latest_checkpoint()
1322
1323        if skip_initial_step is False:
1324            yield(self.time)      # Yield initial values
1325
1326        while True:
1327            # Evolve One Step, using appropriate timestepping method
1328            if self.get_timestepping_method() == 'euler':
1329                self.evolve_one_euler_step(yieldstep, finaltime)
1330
1331            elif self.get_timestepping_method() == 'rk2':
1332                self.evolve_one_rk2_step(yieldstep, finaltime)
1333
1334            elif self.get_timestepping_method() == 'rk3':
1335                self.evolve_one_rk3_step(yieldstep, finaltime)
1336
1337            # Update extrema if necessary (for reporting)
1338            self.update_extrema()           
1339
1340            self.number_of_steps += 1
1341            if self._order_ == 1:
1342                self.number_of_first_order_steps += 1
1343
1344            # Yield results
1345            if finaltime is not None and self.time >= finaltime-epsilon:
1346                if self.time > finaltime:
1347                    # FIXME (Ole, 30 April 2006): Do we need this check?
1348                    # Probably not (Ole, 18 September 2008).
1349                    # Now changed to Exception.
1350                    msg = ('WARNING (domain.py): time overshot finaltime. '
1351                           'Contact Ole.Nielsen@ga.gov.au')
1352                    raise Exception, msg
1353
1354                # Yield final time and stop
1355                self.time = finaltime
1356                yield(self.time)
1357                break
1358
1359            # if we are at the next yield point
1360            if self.time >= self.yieldtime:
1361                # Yield (intermediate) time and allow inspection of domain
1362                if self.checkpoint is True:
1363                    self.store_checkpoint()
1364                    self.delete_old_checkpoints()
1365
1366                # Pass control on to outer loop for more specific actions
1367                yield(self.time)
1368
1369                # Reinitialise
1370                self.yieldtime += yieldstep                 # move to next yield
1371                self.min_timestep = max_timestep
1372                self.max_timestep = min_timestep
1373                self.number_of_steps = 0
1374                self.number_of_first_order_steps = 0
1375                self.max_speed = num.zeros(N, num.Float)
1376
1377    ##
1378    # @brief 'Euler' time step method.
1379    # @param yieldstep The reporting time step.
1380    # @param finaltime The simulation final time.
1381    def evolve_one_euler_step(self, yieldstep, finaltime):
1382        """One Euler Time Step
1383        Q^{n+1} = E(h) Q^n
1384        """
1385
1386        # Compute fluxes across each element edge
1387        self.compute_fluxes()
1388
1389        # Update timestep to fit yieldstep and finaltime
1390        self.update_timestep(yieldstep, finaltime)
1391
1392        # Update conserved quantities
1393        self.update_conserved_quantities()
1394
1395        # Update ghosts
1396        self.update_ghosts()
1397
1398        # Update time
1399        self.time += self.timestep
1400
1401        # Update vertex and edge values
1402        self.distribute_to_vertices_and_edges()
1403
1404        # Update boundary values
1405        self.update_boundary()
1406
1407    ##
1408    # @brief 'rk2' time step method.
1409    # @param yieldstep The reporting time step.
1410    # @param finaltime The simulation final time.
1411    def evolve_one_rk2_step(self, yieldstep, finaltime):
1412        """One 2nd order RK timestep
1413        Q^{n+1} = 0.5 Q^n + 0.5 E(h)^2 Q^n
1414        """
1415
1416        # Save initial initial conserved quantities values
1417        self.backup_conserved_quantities()
1418
1419        ######
1420        # First euler step
1421        ######
1422
1423        # Compute fluxes across each element edge
1424        self.compute_fluxes()
1425
1426        # Update timestep to fit yieldstep and finaltime
1427        self.update_timestep(yieldstep, finaltime)
1428
1429        # Update conserved quantities
1430        self.update_conserved_quantities()
1431
1432        # Update ghosts
1433        self.update_ghosts()
1434
1435        # Update time
1436        self.time += self.timestep
1437
1438        # Update vertex and edge values
1439        self.distribute_to_vertices_and_edges()
1440
1441        # Update boundary values
1442        self.update_boundary()
1443
1444        ######
1445        # Second Euler step
1446        ######
1447
1448        # Compute fluxes across each element edge
1449        self.compute_fluxes()
1450
1451        # Update conserved quantities
1452        self.update_conserved_quantities()
1453
1454        ######
1455        # Combine initial and final values
1456        # of conserved quantities and cleanup
1457        ######
1458
1459        # Combine steps
1460        self.saxpy_conserved_quantities(0.5, 0.5)
1461
1462        # Update ghosts
1463        self.update_ghosts()
1464
1465        # Update vertex and edge values
1466        self.distribute_to_vertices_and_edges()
1467
1468        # Update boundary values
1469        self.update_boundary()
1470
1471    ##
1472    # @brief 'rk3' time step method.
1473    # @param yieldstep The reporting time step.
1474    # @param finaltime The simulation final time.
1475    def evolve_one_rk3_step(self, yieldstep, finaltime):
1476        """One 3rd order RK timestep
1477        Q^(1) = 3/4 Q^n + 1/4 E(h)^2 Q^n  (at time t^n + h/2)
1478        Q^{n+1} = 1/3 Q^n + 2/3 E(h) Q^(1) (at time t^{n+1})
1479        """
1480
1481        # Save initial initial conserved quantities values
1482        self.backup_conserved_quantities()
1483
1484        initial_time = self.time
1485
1486        ######
1487        # First euler step
1488        ######
1489
1490        # Compute fluxes across each element edge
1491        self.compute_fluxes()
1492
1493        # Update timestep to fit yieldstep and finaltime
1494        self.update_timestep(yieldstep, finaltime)
1495
1496        # Update conserved quantities
1497        self.update_conserved_quantities()
1498
1499        # Update ghosts
1500        self.update_ghosts()
1501
1502        # Update time
1503        self.time += self.timestep
1504
1505        # Update vertex and edge values
1506        self.distribute_to_vertices_and_edges()
1507
1508        # Update boundary values
1509        self.update_boundary()
1510
1511        ######
1512        # Second Euler step
1513        ######
1514
1515        # Compute fluxes across each element edge
1516        self.compute_fluxes()
1517
1518        # Update conserved quantities
1519        self.update_conserved_quantities()
1520
1521        ######
1522        # Combine steps to obtain intermediate
1523        # solution at time t^n + 0.5 h
1524        ######
1525
1526        # Combine steps
1527        self.saxpy_conserved_quantities(0.25, 0.75)
1528
1529        # Update ghosts
1530        self.update_ghosts()
1531
1532        # Set substep time
1533        self.time = initial_time + self.timestep*0.5
1534
1535        # Update vertex and edge values
1536        self.distribute_to_vertices_and_edges()
1537
1538        # Update boundary values
1539        self.update_boundary()
1540
1541        ######
1542        # Third Euler step
1543        ######
1544
1545        # Compute fluxes across each element edge
1546        self.compute_fluxes()
1547
1548        # Update conserved quantities
1549        self.update_conserved_quantities()
1550
1551        ######
1552        # Combine final and initial values
1553        # and cleanup
1554        ######
1555
1556        # Combine steps
1557        self.saxpy_conserved_quantities(2.0/3.0, 1.0/3.0)
1558
1559        # Update ghosts
1560        self.update_ghosts()
1561
1562        # Set new time
1563        self.time = initial_time + self.timestep
1564
1565        # Update vertex and edge values
1566        self.distribute_to_vertices_and_edges()
1567
1568        # Update boundary values
1569        self.update_boundary()
1570
1571    ##
1572    # @brief Evolve simulation to a final time.
1573    # @param finaltime Sinulation final time.
1574    def evolve_to_end(self, finaltime=1.0):
1575        """Iterate evolve all the way to the end."""
1576
1577        for _ in self.evolve(yieldstep=None, finaltime=finaltime):
1578            pass
1579
1580    ##
1581    # @brief Backup conserved quantities.
1582    def backup_conserved_quantities(self):
1583        N = len(self) # Number_of_triangles
1584
1585        # Backup conserved_quantities centroid values
1586        for name in self.conserved_quantities:
1587            Q = self.quantities[name]
1588            Q.backup_centroid_values()
1589
1590    ##
1591    # @brief ??
1592    # @param a ??
1593    # @param b ??
1594    def saxpy_conserved_quantities(self, a, b):
1595        N = len(self) #number_of_triangles
1596
1597        # Backup conserved_quantities centroid values
1598        for name in self.conserved_quantities:
1599            Q = self.quantities[name]
1600            Q.saxpy_centroid_values(a, b)
1601
1602    ##
1603    # @brief Update boundary values for all conserved quantities.
1604    def update_boundary(self):
1605        """Go through list of boundary objects and update boundary values
1606        for all conserved quantities on boundary.
1607        It is assumed that the ordering of conserved quantities is
1608        consistent between the domain and the boundary object, i.e.
1609        the jth element of vector q must correspond to the jth conserved
1610        quantity in domain.
1611        """
1612
1613        # FIXME: Update only those that change (if that can be worked out)
1614        # FIXME: Boundary objects should not include ghost nodes.
1615        for i, ((vol_id, edge_id), B) in enumerate(self.boundary_objects):
1616            if B is None:
1617                print 'WARNING: Ignored boundary segment (None)'
1618            else:
1619                q = B.evaluate(vol_id, edge_id)
1620
1621                for j, name in enumerate(self.conserved_quantities):
1622                    Q = self.quantities[name]
1623                    Q.boundary_values[i] = q[j]
1624
1625    ##
1626    # @brief Compute fluxes.
1627    # @note MUST BE OVERRIDEN IN SUBCLASS!
1628    def compute_fluxes(self):
1629        msg = 'Method compute_fluxes must be overridden by Domain subclass'
1630        raise Exception, msg
1631
1632    ##
1633    # @brief
1634    # @param yieldstep
1635    # @param finaltime
1636    def update_timestep(self, yieldstep, finaltime):
1637        from anuga.config import min_timestep, max_timestep
1638
1639        # Protect against degenerate timesteps arising from isolated
1640        # triangles
1641        # FIXME (Steve): This should be in shallow_water as it assumes x and y
1642        # momentum
1643        if self.protect_against_isolated_degenerate_timesteps is True and \
1644               self.max_speed > 10.0: # FIXME (Ole): Make this configurable
1645
1646            # Setup 10 bins for speed histogram
1647            from anuga.utilities.numerical_tools import histogram, create_bins
1648
1649            bins = create_bins(self.max_speed, 10)
1650            hist = histogram(self.max_speed, bins)
1651
1652            # Look for characteristic signature
1653            if len(hist) > 1 and hist[-1] > 0 and \
1654                hist[4] == hist[5] == hist[6] == hist[7] == hist[8] == 0:
1655                    # Danger of isolated degenerate triangles
1656                    # print self.timestepping_statistics(track_speeds=True)
1657
1658                    # Find triangles in last bin
1659                    # FIXME - speed up using Numeric
1660                    d = 0
1661                    for i in range(self.number_of_full_triangles):
1662                        if self.max_speed[i] > bins[-1]:
1663                            msg = 'Time=%f: Ignoring isolated high ' % self.time
1664                            msg += 'speed triangle '
1665                            msg += '#%d of %d with max speed=%f' \
1666                                      % (i, self.number_of_full_triangles,
1667                                         self.max_speed[i])
1668
1669                            # print 'Found offending triangle', i,
1670                            # self.max_speed[i]
1671                            self.get_quantity('xmomentum').\
1672                                            set_values(0.0, indices=[i])
1673                            self.get_quantity('ymomentum').\
1674                                            set_values(0.0, indices=[i])
1675                            self.max_speed[i]=0.0
1676                            d += 1
1677
1678        # self.timestep is calculated from speed of characteristics
1679        # Apply CFL condition here
1680        timestep = min(self.CFL*self.flux_timestep, max_timestep)
1681
1682        # Record maximal and minimal values of timestep for reporting
1683        self.max_timestep = max(timestep, self.max_timestep)
1684        self.min_timestep = min(timestep, self.min_timestep)
1685
1686        # Protect against degenerate time steps
1687        if timestep < min_timestep:
1688            # Number of consecutive small steps taken b4 taking action
1689            self.smallsteps += 1
1690
1691            if self.smallsteps > self.max_smallsteps:
1692                self.smallsteps = 0 # Reset
1693
1694                if self._order_ == 1:
1695                    msg = 'WARNING: Too small timestep %.16f reached ' \
1696                              % timestep
1697                    msg += 'even after %d steps of 1 order scheme' \
1698                               % self.max_smallsteps
1699                    print msg
1700                    timestep = min_timestep  # Try enforcing min_step
1701
1702                    print self.timestepping_statistics(track_speeds=True)
1703
1704                    raise Exception, msg
1705                else:
1706                    # Try to overcome situation by switching to 1 order
1707                    self._order_ = 1
1708        else:
1709            self.smallsteps = 0
1710            if self._order_ == 1 and self.default_order == 2:
1711                self._order_ = 2
1712
1713        # Ensure that final time is not exceeded
1714        if finaltime is not None and self.time + timestep > finaltime :
1715            timestep = finaltime-self.time
1716
1717        # Ensure that model time is aligned with yieldsteps
1718        if self.time + timestep > self.yieldtime:
1719            timestep = self.yieldtime - self.time
1720
1721        self.timestep = timestep
1722
1723    ##
1724    # @brief Compute forcing terms, if any.
1725    def compute_forcing_terms(self):
1726        """If there are any forcing functions driving the system
1727        they should be defined in Domain subclass and appended to
1728        the list self.forcing_terms
1729        """
1730
1731        for f in self.forcing_terms:
1732            f(self)
1733
1734    ##
1735    # @brief Update vectors of conserved quantities.
1736    def update_conserved_quantities(self):
1737        """Update vectors of conserved quantities using previously
1738        computed fluxes specified forcing functions.
1739        """
1740
1741        N = len(self) # Number_of_triangles
1742        d = len(self.conserved_quantities)
1743
1744        timestep = self.timestep
1745
1746        # Compute forcing terms
1747        self.compute_forcing_terms()
1748
1749        # Update conserved_quantities
1750        for name in self.conserved_quantities:
1751            Q = self.quantities[name]
1752            Q.update(timestep)
1753
1754            # Note that Q.explicit_update is reset by compute_fluxes
1755            # Where is Q.semi_implicit_update reset?
1756            # It is reset in quantity_ext.c
1757
1758    ##
1759    # @brief Sequential update of ghost cells
1760    def update_ghosts(self):
1761        # We must send the information from the full cells and
1762        # receive the information for the ghost cells
1763        # We have a list with ghosts expecting updates
1764
1765
1766        from Numeric import take,put
1767
1768
1769        #Update of ghost cells
1770        iproc = self.processor
1771        if self.full_send_dict.has_key(iproc):
1772
1773            # now store full as local id, global id, value
1774            Idf  = self.full_send_dict[iproc][0]
1775
1776            # now store ghost as local id, global id, value
1777            Idg = self.ghost_recv_dict[iproc][0]
1778
1779            for i, q in enumerate(self.conserved_quantities):
1780                Q_cv =  self.quantities[q].centroid_values
1781                put(Q_cv,     Idg, take(Q_cv,     Idf))
1782
1783 
1784    ##
1785    # @brief Extrapolate conserved quantities from centroid to vertices
1786    #        and edge-midpoints for each volume.
1787    def distribute_to_vertices_and_edges(self):
1788        """Extrapolate conserved quantities from centroid to
1789        vertices and edge-midpoints for each volume
1790
1791        Default implementation is straight first order,
1792        i.e. constant values throughout each element and
1793        no reference to non-conserved quantities.
1794        """
1795
1796        for name in self.conserved_quantities:
1797            Q = self.quantities[name]
1798            if self._order_ == 1:
1799                Q.extrapolate_first_order()
1800            elif self._order_ == 2:
1801                Q.extrapolate_second_order()
1802            else:
1803                raise Exception, 'Unknown order'
1804
1805    ##
1806    # @brief Calculate the norm of the centroid values of a specific quantity,
1807    #        using normfunc.
1808    # @param quantity
1809    # @param normfunc
1810    def centroid_norm(self, quantity, normfunc):
1811        """Calculate the norm of the centroid values of a specific quantity,
1812        using normfunc.
1813
1814        normfunc should take a list to a float.
1815
1816        common normfuncs are provided in the module utilities.norms
1817        """
1818
1819        return normfunc(self.quantities[quantity].centroid_values)
1820
1821
1822######
1823# Initialise module
1824######
1825
1826# Optimisation with psyco
1827from anuga.config import use_psyco
1828
1829if use_psyco:
1830    try:
1831        import psyco
1832    except:
1833        import os
1834        if os.name == 'posix' and os.uname()[4] in ['x86_64', 'ia64']:
1835            pass
1836            # Psyco isn't supported on 64 bit systems, but it doesn't matter
1837        else:
1838            msg = ('WARNING: psyco (speedup) could not be imported, '
1839                   'you may want to consider installing it')
1840            print msg
1841    else:
1842        psyco.bind(Domain.update_boundary)
1843        #psyco.bind(Domain.update_timestep) # Not worth it
1844        psyco.bind(Domain.update_conserved_quantities)
1845        psyco.bind(Domain.distribute_to_vertices_and_edges)
1846
1847
1848if __name__ == "__main__":
1849    pass
Note: See TracBrowser for help on using the repository browser.