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

Last change on this file since 7810 was 7810, checked in by hudson, 14 years ago

Added aggressive psyco optimisation, fixed benchmark app to work with new API.

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