source: branches/anuga_1_2_1/anuga_core/source/anuga/abstract_2d_finite_volumes/generic_domain.py @ 8083

Last change on this file since 8083 was 8083, checked in by habili, 13 years ago

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