source: branches/numpy/anuga/abstract_2d_finite_volumes/domain.py @ 6304

Last change on this file since 6304 was 6304, checked in by rwilson, 15 years ago

Initial commit of numpy changes. Still a long way to go.

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