source: trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/generic_domain.py @ 8125

Last change on this file since 8125 was 8125, checked in by wilsonr, 13 years ago

Changes to address ticket 360.

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