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

Last change on this file since 8068 was 8068, checked in by habili, 12 years ago

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