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

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

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