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

Last change on this file since 7938 was 7938, checked in by steve, 9 years ago

adding in option for extrafractional steps.

File size: 73.4 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                       processor=0,
65                       numproc=1,
66                       number_of_full_nodes=None,
67                       number_of_full_triangles=None):
68
69        """Instantiate generic computational Domain.
70
71        Input:
72          source:    Either a mesh filename or coordinates of mesh vertices.
73                     If it is a filename values specified for triangles will
74                     be overridden.
75          triangles: Mesh connectivity (see mesh.py for more information)
76          boundary:  See mesh.py for more information
77
78          conserved_quantities: List of quantity names entering the
79                                conservation equations
80          evolved_quantities:   List of all quantities that evolve
81          other_quantities:     List of other quantity names
82
83          tagged_elements:
84          ...
85        """
86
87        number_of_full_nodes=None
88        number_of_full_triangles=None
89       
90        # Determine whether source is a mesh filename or coordinates
91        if 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
105        self.mesh = Mesh(coordinates, triangles,
106                         boundary=boundary,
107                         tagged_elements=tagged_elements,
108                         geo_reference=geo_reference,
109                         use_inscribed_circle=use_inscribed_circle,
110                         number_of_full_nodes=number_of_full_nodes,
111                         number_of_full_triangles=number_of_full_triangles,
112                         verbose=verbose)
113
114        # Expose Mesh attributes (FIXME: Maybe turn into methods)
115        self.triangles = self.mesh.triangles       
116        self.centroid_coordinates = self.mesh.centroid_coordinates
117        self.vertex_coordinates = self.mesh.vertex_coordinates
118        self.boundary = self.mesh.boundary
119        self.neighbours = self.mesh.neighbours
120        self.surrogate_neighbours = self.mesh.surrogate_neighbours
121        self.neighbour_edges = self.mesh.neighbour_edges
122        self.normals = self.mesh.normals
123        self.edgelengths = self.mesh.edgelengths
124        self.radii = self.mesh.radii
125        self.areas = self.mesh.areas
126
127        self.number_of_boundaries = self.mesh.number_of_boundaries
128        self.number_of_full_nodes = self.mesh.number_of_full_nodes
129        self.number_of_full_triangles = self.mesh.number_of_full_triangles
130        self.number_of_triangles_per_node = \
131                                    self.mesh.number_of_triangles_per_node
132
133        self.vertex_value_indices = self.mesh.vertex_value_indices
134        self.number_of_triangles = self.mesh.number_of_triangles
135
136        self.geo_reference = self.mesh.geo_reference
137
138        if verbose: log.critical('Initialising Domain')
139
140        # List of quantity names entering the conservation equations
141        if conserved_quantities is None:
142            self.conserved_quantities = []
143        else:
144            self.conserved_quantities = conserved_quantities
145
146        if evolved_quantities is None:
147            self.evolved_quantities = self.conserved_quantities
148        else:
149            self.evolved_quantities = evolved_quantities
150           
151        # List of other quantity names
152        if other_quantities is None:
153            self.other_quantities = []
154        else:
155            self.other_quantities = other_quantities
156
157        # Test that conserved_quantities are stored in the first entries of
158        # evolved_quantities
159        for i, quantity in enumerate(self.conserved_quantities):
160            msg = 'The conserved quantities must be the first entries of '
161            msg += 'evolved_quantities'
162            assert quantity == self.evolved_quantities[i], msg
163           
164
165        # Build dictionary of Quantity instances keyed by quantity names
166        self.quantities = {}
167
168        for name in self.evolved_quantities:
169            self.quantities[name] = Quantity(self)
170        for name in self.other_quantities:
171            self.quantities[name] = Quantity(self)
172
173        # Create an empty list for explicit forcing terms
174        self.forcing_terms = []
175
176        # Setup the ghost cell communication
177        if full_send_dict is None:
178            self.full_send_dict = {}
179        else:
180            self.full_send_dict = full_send_dict
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
192        if verbose: log.critical('Domain: Set up communication buffers '
193                                 '(parallel)')
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]
197            self.full_send_dict[key].append(num.zeros((buffer_shape, self.nsys),
198                                                      num.float))
199
200        for key in self.ghost_recv_dict:
201            buffer_shape = self.ghost_recv_dict[key][0].shape[0]
202            self.ghost_recv_dict[key].append( \
203                                            num.zeros((buffer_shape, self.nsys),
204                                             num.float))
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
211        self.tri_full_flag = num.ones(N, num.int)
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.
218        if not num.allclose(self.tri_full_flag[:self.number_of_full_nodes], 1):
219            if self.numproc>1:
220                log.critical('WARNING: Not all full triangles are store before '
221                             'ghost triangles')
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
228        from anuga.config import default_order
229        from anuga.config import max_timestep, min_timestep
230
231        self.beta_w = beta_w
232        self.epsilon = epsilon
233        self.protect_against_isolated_degenerate_timesteps = \
234                        protect_against_isolated_degenerate_timesteps
235
236
237        self.centroid_transmissive_bc = False
238        self.set_default_order(default_order)
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)
247        self.set_evolve_max_timestep(max_timestep)
248        self.set_evolve_min_timestep(min_timestep)
249        self.boundary_map = None  # Will be populated by set_boundary
250
251        # Model time
252        self.time = 0.0
253        self.finaltime = None
254        self.recorded_min_timestep = self.recorded_max_timestep = 0.0
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()
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 + self.starttime
499
500    ##
501    # @brief Set the default beta for limiting.
502    # @param beta The new beta value.
503    def set_beta(self, beta):
504        """Set default beta for limiting."""
505
506        self.beta = beta
507        for name in self.quantities:
508            Q = self.quantities[name]
509            Q.set_beta(beta)
510
511    ##
512    # @brief Get the beta value used for limiting.
513    # @return The beta value used for limiting.
514    def get_beta(self):
515        """Get default beta for limiting."""
516
517        return self.beta
518
519
520    ##
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    ##
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    ##
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.
587    def set_default_order(self, n):
588        """Set default (spatial) order to either 1 or 2."""
589
590        msg = 'Default order must be either 1 or 2. I got %s' % n
591        assert n in [1,2], msg
592
593        self.default_order = n
594        self._order_ = self.default_order
595
596    ##
597    # @brief Set values of named quantities.
598    # @param quantity_dict Dictionary containing name/value pairs.
599    def set_quantity_vertices_dict(self, quantity_dict):
600        """Set values for named quantities.
601        Supplied dictionary contains name/value pairs:
602
603        name:  Name of quantity
604        value: Compatible list, numeric array, const or function (see below)
605
606        The values will be stored in elements following their internal ordering.
607        """
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
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):
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)
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):
652        """Add values to a named quantity
653
654        E.g add_quantity('elevation', X)
655
656        Option are the same as in set_quantity.
657        """
658
659        # Do the expression stuff
660        if kwargs.has_key('expression'):
661            expression = kwargs['expression']
662            Q2 = self.create_quantity_from_expression(expression)
663        else:
664            # Create new temporary quantity
665            Q2 = Quantity(self)
666
667            # Assign specified values to temporary quantity
668            Q2.set_values(*args, **kwargs)
669
670        # Add temporary quantity to named quantity
671        Q1 = self.get_quantity(name)
672        self.set_quantity(name, Q1 + Q2)
673
674    ##
675    # @brief Get list of quantity names for the Domain.
676    # @return List of quantity names.
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        """
681
682        return self.quantities.keys()
683
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):
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
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.
710    def create_quantity_from_expression(self, expression):
711        """Create new quantity from other quantities using arbitrary expression.
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')
722            exp = '(xmomentum*xmomentum + ymomentum*ymomentum)**0.5'
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
728
729        return apply_expression_to_dictionary(expression, self.quantities)
730
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.
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.
765        FIXME: This would lead to implementation of a default boundary condition
766
767        Note: If a segment is listed in the boundary dictionary and if it is
768        not None, it *will* become a boundary - even if there is a neighbouring
769        triangle.  This would be the case for internal boundaries.
770
771        Boundary objects that are None will be skipped.
772
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.
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.
786            self.boundary_map = boundary_map
787        else:
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]
792
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
799        self.boundary_objects = []
800        for k, (vol_id, edge_id) in enumerate(x):
801            tag = self.boundary[(vol_id, edge_id)]
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:
807                    self.boundary_objects.append(((vol_id, edge_id), B))
808                    self.neighbours[vol_id, edge_id] = \
809                                        -len(self.boundary_objects)
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()
819                raise Exception, msg
820
821    ##
822    # @brief Set quantities based on a regional tag.
823    # @param args
824    # @param kwargs
825    def set_region(self, *args, **kwargs):
826        """Set quantities based on a regional tag.
827
828        It is most often called with the following parameters;
829        (self, tag, quantity, X, location='vertices')
830        tag:      the name of the regional tag used to specify the region
831        quantity: Name of quantity to change
832        X:        const or function - how the quantity is changed
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)
846
847    ##
848    # @brief ??
849    # @param functions A list or tuple of ??
850    def _set_region(self, functions):
851        # coerce to an iterable (list or tuple)
852        if type(functions) not in [types.ListType, types.TupleType]:
853            functions = [functions]
854
855        # The order of functions in the list is used.
856        tagged_elements = self.get_tagged_elements()
857        for function in functions:
858            for tag in tagged_elements.keys():
859                function(tag, tagged_elements[tag], self)
860
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.
867    def set_quantities_to_be_monitored(self, q,
868                                             polygon=None,
869                                             time_interval=None):
870        """Specify which quantities will be monitored for extrema.
871
872        q must be either:
873          - the name of a quantity or derived quantity such as 'stage-elevation'
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
879
880        If q is None, monitoring will be switched off altogether.
881
882        polygon (if specified) will only monitor triangles inside polygon.
883        If omitted all triangles will be included.
884
885        time_interval, if specified, will restrict monitoring to time steps in
886        that interval. If omitted all timesteps will be included.
887        """
888
889        from anuga.abstract_2d_finite_volumes.util import\
890             apply_expression_to_dictionary
891
892        if q is None:
893            self.quantities_to_be_monitored = None
894            self.monitor_polygon = None
895            self.monitor_time_interval = None
896            self.monitor_indices = None
897            return
898
899        # coerce 'q' to a list if it's a string
900        if isinstance(q, basestring):
901            q = [q]
902
903        # Check correctness and initialise
904        self.quantities_to_be_monitored = {}
905        for quantity_name in q:
906            msg = 'Quantity %s is not a valid conserved quantity' \
907                      % quantity_name
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)
918                          'min_time': None,     # Argmin (t)
919                          'max_time': None}     # Argmax (t)
920
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.
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
932                if polygon in self.quantities:
933                    raise Exception, msg
934
935                try:
936                    apply_expression_to_dictionary(polygon, self.quantities)
937                except:
938                    # At least polygon wasn't expression involving quantitites
939                    pass
940                else:
941                    raise Exception, msg
942
943                # In any case, we don't allow polygon to be a string
944                msg = ('argument "polygon" must not be a string: '
945                       'I got polygon="%s"') % polygon
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
956        self.monitor_time_interval = time_interval
957
958    ##
959    # @brief Check Domain integrity.
960    # @note Raises an exception if integrity breached.
961    def check_integrity(self):
962        self.mesh.check_integrity()
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
969        for i, quantity in enumerate(self.conserved_quantities):
970            msg = 'Conserved quantities must be the first entries '
971            msg += 'of evolved_quantities'
972            assert quantity == self.evolved_quantities[i], msg
973 
974
975    ##
976    # @brief Print timestep stats to stdout.
977    # @param track_speeds If True, print smallest track speed.
978    def write_time(self, track_speeds=False):
979        log.critical(self.timestepping_statistics(track_speeds))
980
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
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
995        triangle rather than the one with the largest speed.
996        """
997
998        from anuga.utilities.numerical_tools import histogram, create_bins
999
1000        # qwidth determines the the width of the text field used for quantities
1001        qwidth = self.qwidth = 12
1002
1003        msg = ''
1004
1005        model_time = self.get_time()
1006        if self.recorded_min_timestep == self.recorded_max_timestep:
1007            msg += 'Time = %.4f, delta t = %.8f, steps=%d' \
1008                       % (model_time, self.recorded_min_timestep, \
1009                                    self.number_of_steps)
1010        elif self.recorded_min_timestep > self.recorded_max_timestep:
1011            msg += 'Time = %.4f, steps=%d' \
1012                       % (model_time, self.number_of_steps)
1013        else:
1014            msg += 'Time = %.4f, delta t in [%.8f, %.8f], steps=%d' \
1015                       % (model_time, self.recorded_min_timestep,
1016                          self.recorded_max_timestep, self.number_of_steps)
1017
1018        msg += ' (%ds)' % (walltime() - self.last_walltime)
1019        self.last_walltime = walltime()
1020
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'
1029            msg += '  Speeds in [%f, %f]\n' % (num.min(self.max_speed),
1030                                               num.max(self.max_speed))
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):
1037                    # Open upper interval
1038                    hi = bins[i+1]
1039                    msg += '    [%f, %f[: %d\n' % (lo, hi, count)
1040                else:
1041                    # Closed upper interval
1042                    hi = num.max(self.max_speed)
1043                    msg += '    [%f, %f]: %d\n' % (lo, hi, count)
1044
1045            N = len(self.max_speed.flat)
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)
1053                for i, a in enumerate(speed):
1054                    if i % (N/10) == 0 and i != 0:
1055                        # For every 10% of the sorted speeds
1056                        msg += '    %d speeds in [%f, %f]\n' % (i-k, lower, a)
1057                        lower = a
1058                        k = i
1059
1060                msg += '    %d speeds in [%f, %f]\n'\
1061                           % (N-k, lower, max(speed))
1062
1063            # Find index of largest computed flux speed
1064            if triangle_id is None:
1065                k = self.k = num.argmax(self.max_speed)
1066            else:
1067                errmsg = 'Triangle_id %d does not exist in mesh: %s' \
1068                             % (triangle_id, str(self))
1069                assert 0 <= triangle_id < len(self), errmsg
1070                k = self.k = triangle_id
1071
1072            x, y = self.get_centroid_coordinates(absolute=True)[k]
1073            radius = self.get_radii()[k]
1074            area = self.get_areas()[k]
1075            max_speed = self.max_speed[k]
1076
1077            msg += '  Triangle #%d with centroid (%.4f, %.4f), ' % (k, x, y)
1078            msg += 'area = %.4f and radius = %.4f ' % (area, radius)
1079            if triangle_id is None:
1080                msg += 'had the largest computed speed: %.6f m/s ' % (max_speed)
1081            else:
1082                msg += 'had computed speed: %.6f m/s ' % (max_speed)
1083
1084            if max_speed > 0.0:
1085                msg += '(timestep=%.6f)\n' % (radius/max_speed)
1086            else:
1087                msg += '(timestep=%.6f)\n' % (0)
1088
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]
1097                C = q.get_values(location='centroids', indices=[k])
1098
1099                s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n' \
1100                         % (name.ljust(qwidth), V[0], V[1], V[2])
1101
1102                s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n' \
1103                         % (name.ljust(qwidth), E[0], E[1], E[2])
1104
1105                s += '    %s: centroid_value = %.4f\n' \
1106                         % (name.ljust(qwidth), C[0])
1107
1108                msg += s
1109
1110        return msg
1111
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):
1117        log.critical(self.boundary_statistics(quantities, tags))
1118
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):
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
1147        # Input checks
1148        if quantities is None:
1149            quantities = self.evolved_quantities
1150        elif type(quantities) == types.StringType:
1151            quantities = [quantities] #Turn it into a list
1152
1153        msg = ('Keyword argument quantities must be either None, '
1154               'string or list. I got %s') % str(quantities)
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
1162        msg = ('Keyword argument tags must be either None, '
1163               'string or list. I got %s') % str(tags)
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
1174        msg = 'Boundary values at time %.4f:\n' % self.get_time()
1175        for tag in tags:
1176            msg += '    %s:\n' % tag
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
1183                for i, ((vol_id,edge_id),B) in enumerate(self.boundary_objects):
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:
1190                    msg += ('        Sorry no information available about'
1191                            ' tag %s and quantity %s\n') % (tag, name)
1192                else:
1193                    msg += '        %s in [%12.8f, %12.8f]\n' \
1194                               % (string.ljust(name, maxwidth), minval, maxval)
1195
1196        return msg
1197
1198    ##
1199    # @brief Update extrema if requested by set_quantities_to_be_monitored.
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
1208        from anuga.config import single_precision as epsilon
1209
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)
1234            if info_block['max'] is None or \
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)
1243            if info_block['min'] is None or \
1244                   minval < info_block['min'] - epsilon:
1245                info_block['min'] = minval
1246                minloc = Q.get_minimum_location()
1247                info_block['min_location'] = minloc
1248                info_block['min_time'] = self.time
1249
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'):
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
1266        msg = 'Monitored quantities at time %.4f:\n' % self.get_time()
1267        if self.monitor_polygon is not None:
1268            p_str = str(self.monitor_polygon)
1269            msg += '- Restricted by polygon: %s' % p_str[:maxlen]
1270            if len(p_str) >= maxlen:
1271                msg += '...\n'
1272            else:
1273                msg += '\n'
1274
1275        if self.monitor_time_interval is not None:
1276            msg += '- Restricted by time interval: %s\n' \
1277                       % str(self.monitor_time_interval)
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():
1283            msg += '    %s:\n' % quantity_name
1284
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))
1289
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))
1293
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))
1297
1298        return msg
1299
1300    ##
1301    # @brief Get the timestep method.
1302    # @return The timestep method. One of 'euler', 'rk2' or 'rk3' or 1, 2, 3.
1303    def get_timestepping_method(self):
1304        return self.timestepping_method
1305
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):
1311        methods = ['euler', 'rk2', 'rk3']   
1312        if timestepping_method in methods:
1313            self.timestepping_method = timestepping_method
1314            return
1315        if timestepping_method in [1,2,3]:
1316            self.timetepping_method = methods[timestepping_method-1]
1317            return
1318
1319        msg = '%s is an incorrect timestepping type' % timestepping_method
1320        raise Exception, msg
1321
1322    ##
1323    # @brief Get the Domain simulation name.
1324    # @return The simulation name string.
1325    def get_name(self):
1326        return self.simulation_name
1327
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.
1332    def set_name(self, name):
1333        """Assign a name to this simulation.
1334        This will be used to identify the output sww file.
1335        """
1336
1337        # remove any '.sww' end
1338        if name.endswith('.sww'):
1339            name = name[:-4]
1340
1341        self.simulation_name = name
1342
1343    ##
1344    # @brief Get data directory path.
1345    # @return The data directory path string.
1346    def get_datadir(self):
1347        return self.datadir
1348
1349    ##
1350    # @brief Set data directory path.
1351    # @param name The data directory path string.
1352    def set_datadir(self, name):
1353        self.datadir = name
1354
1355    ##
1356    # @brief Get the start time value.
1357    # @return The start time value (float).
1358    def get_starttime(self):
1359        return self.starttime
1360
1361    ##
1362    # @brief Set the start time value.
1363    # @param time The start time value.
1364    def set_starttime(self, time):
1365        self.starttime = float(time)
1366
1367################################################################################
1368# Main components of evolve
1369################################################################################
1370
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):
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
1408        from anuga.config import epsilon
1409
1410        # FIXME: Maybe lump into a larger check prior to evolving
1411        msg = ('Boundary tags must be bound to boundary objects before '
1412               'evolving system, '
1413               'e.g. using the method set_boundary.\n'
1414               'This system has the boundary tags %s '
1415                   % self.get_boundary_tags())
1416        assert hasattr(self, 'boundary_objects'), msg
1417
1418        if yieldstep is None:
1419            yieldstep = self.evolve_max_timestep
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'
1427            raise Exception, msg
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
1434        N = len(self)                             # Number of triangles
1435        self.yieldtime = self.time + yieldstep    # set next yield time
1436
1437        # Initialise interval of timestep sizes (for reporting only)
1438        # Note that we set recorded_min_timestep to be large so that it comes
1439        # down through the evolution, similarly recorded_max_timestep
1440        self.recorded_min_timestep = self.evolve_max_timestep
1441        self.recorded_max_timestep = self.evolve_min_timestep
1442        self.number_of_steps = 0
1443        self.number_of_first_order_steps = 0
1444
1445        # Update ghosts
1446        self.update_ghosts()
1447
1448        # Initial update of vertex and edge values
1449        self.distribute_to_vertices_and_edges()
1450
1451        # Initial update boundary values
1452        self.update_boundary()
1453
1454        # Update extrema if necessary (for reporting)
1455        self.update_extrema()
1456
1457
1458
1459        # Or maybe restore from latest checkpoint
1460        if self.checkpoint is True:
1461            self.goto_latest_checkpoint()
1462
1463        if skip_initial_step is False:
1464            yield(self.time)      # Yield initial values
1465
1466        while True:
1467
1468            initial_time = self.time
1469
1470            #==========================================
1471            # Apply fluid flow fractional step
1472            #==========================================
1473            if self.get_timestepping_method() == 'euler':
1474                self.evolve_one_euler_step(yieldstep, finaltime)
1475
1476            elif self.get_timestepping_method() == 'rk2':
1477                self.evolve_one_rk2_step(yieldstep, finaltime)
1478
1479            elif self.get_timestepping_method() == 'rk3':
1480                self.evolve_one_rk3_step(yieldstep, finaltime)
1481
1482
1483            #==========================================
1484            # Apply other fractional steps
1485            #==========================================
1486            #self.apply_fractional_steps()
1487
1488            #==========================================
1489            # Centroid Values of variables should be ok,
1490            # so now setup quantites etc for output
1491            #==========================================
1492
1493            # Update time
1494            self.time = initial_time + self.timestep
1495
1496            # Update vertex and edge values
1497            self.distribute_to_vertices_and_edges()
1498
1499            # Update boundary values
1500            self.update_boundary()
1501
1502            # Update extrema if necessary (for reporting)
1503            self.update_extrema()           
1504
1505            self.number_of_steps += 1
1506            if self._order_ == 1:
1507                self.number_of_first_order_steps += 1
1508
1509            # Yield results
1510            if finaltime is not None and self.time >= finaltime-epsilon:
1511                if self.time > finaltime:
1512                    # FIXME (Ole, 30 April 2006): Do we need this check?
1513                    # Probably not (Ole, 18 September 2008).
1514                    # Now changed to Exception.
1515                    msg = ('WARNING (domain.py): time overshot finaltime. '
1516                           'Contact Ole.Nielsen@ga.gov.au')
1517                    raise Exception, msg
1518
1519                # Yield final time and stop
1520                self.time = finaltime
1521                yield(self.time)
1522                break
1523
1524            # if we are at the next yield point
1525            if self.time >= self.yieldtime:
1526                # Yield (intermediate) time and allow inspection of domain
1527                if self.checkpoint is True:
1528                    self.store_checkpoint()
1529                    self.delete_old_checkpoints()
1530
1531                # Pass control on to outer loop for more specific actions
1532                yield(self.time)
1533
1534                # Reinitialise
1535                self.yieldtime += yieldstep                 # move to next yield
1536                self.recorded_min_timestep = self.evolve_max_timestep
1537                self.recorded_max_timestep = self.evolve_min_timestep
1538                self.number_of_steps = 0
1539                self.number_of_first_order_steps = 0
1540                self.max_speed = num.zeros(N, num.float)
1541
1542    ##
1543    # @brief 'Euler' time step method.
1544    # @param yieldstep The reporting time step.
1545    # @param finaltime The simulation final time.
1546    def evolve_one_euler_step(self, yieldstep, finaltime):
1547        """One Euler Time Step
1548        Q^{n+1} = E(h) Q^n
1549
1550        Assumes that centroid values have been extrapolated to vertices and edges
1551        """
1552
1553        # Compute fluxes across each element edge
1554        self.compute_fluxes()
1555
1556        # Compute forcing terms
1557        self.compute_forcing_terms()
1558
1559        # Update timestep to fit yieldstep and finaltime
1560        self.update_timestep(yieldstep, finaltime)
1561
1562        # Update conserved quantities
1563        self.update_conserved_quantities()
1564
1565        # Update ghosts
1566        self.update_ghosts()
1567
1568
1569
1570
1571    ##
1572    # @brief 'rk2' time step method.
1573    # @param yieldstep The reporting time step.
1574    # @param finaltime The simulation final time.
1575    def evolve_one_rk2_step(self, yieldstep, finaltime):
1576        """One 2nd order RK timestep
1577        Q^{n+1} = 0.5 Q^n + 0.5 E(h)^2 Q^n
1578        """
1579
1580        # Save initial initial conserved quantities values
1581        self.backup_conserved_quantities()
1582
1583        ######
1584        # First euler step
1585        ######
1586
1587        # Compute fluxes across each element edge
1588        self.compute_fluxes()
1589
1590        # Compute forcing terms
1591        self.compute_forcing_terms()
1592
1593        # Update timestep to fit yieldstep and finaltime
1594        self.update_timestep(yieldstep, finaltime)
1595
1596        # Update conserved quantities
1597        self.update_conserved_quantities()
1598
1599        # Update ghosts
1600        self.update_ghosts()
1601
1602        # Update time
1603        self.time += self.timestep
1604
1605        # Update vertex and edge values
1606        self.distribute_to_vertices_and_edges()
1607
1608        # Update boundary values
1609        self.update_boundary()
1610
1611        ######
1612        # Second Euler step using the same timestep
1613        # calculated in the first step. Might lead to
1614        # stability problems but we have not seen any
1615        # example.
1616        ######
1617
1618        # Compute fluxes across each element edge
1619        self.compute_fluxes()
1620
1621        # Compute forcing terms
1622        self.compute_forcing_terms()
1623
1624        # Update conserved quantities
1625        self.update_conserved_quantities()
1626
1627        ######
1628        # Combine initial and final values
1629        # of conserved quantities and cleanup
1630        ######
1631
1632        # Combine steps
1633        self.saxpy_conserved_quantities(0.5, 0.5)
1634
1635        # Update ghosts
1636        self.update_ghosts()
1637
1638
1639    ##
1640    # @brief 'rk3' time step method.
1641    # @param yieldstep The reporting time step.
1642    # @param finaltime The simulation final time.
1643    def evolve_one_rk3_step(self, yieldstep, finaltime):
1644        """One 3rd order RK timestep
1645        Q^(1) = 3/4 Q^n + 1/4 E(h)^2 Q^n  (at time t^n + h/2)
1646        Q^{n+1} = 1/3 Q^n + 2/3 E(h) Q^(1) (at time t^{n+1})
1647        """
1648
1649        # Save initial initial conserved quantities values
1650        self.backup_conserved_quantities()
1651
1652        initial_time = self.time
1653
1654        ######
1655        # First euler step
1656        ######
1657
1658        # Compute fluxes across each element edge
1659        self.compute_fluxes()
1660
1661        # Compute forcing terms
1662        self.compute_forcing_terms()
1663
1664        # Update timestep to fit yieldstep and finaltime
1665        self.update_timestep(yieldstep, finaltime)
1666
1667        # Update conserved quantities
1668        self.update_conserved_quantities()
1669
1670        # Update ghosts
1671        self.update_ghosts()
1672
1673        # Update time
1674        self.time += self.timestep
1675
1676        # Update vertex and edge values
1677        self.distribute_to_vertices_and_edges()
1678
1679        # Update boundary values
1680        self.update_boundary()
1681
1682        ######
1683        # Second Euler step using the same timestep
1684        # calculated in the first step. Might lead to
1685        # stability problems but we have not seen any
1686        # example.
1687        ######
1688
1689        # Compute fluxes across each element edge
1690        self.compute_fluxes()
1691
1692        # Compute forcing terms
1693        self.compute_forcing_terms()
1694
1695        # Update conserved quantities
1696        self.update_conserved_quantities()
1697
1698        ######
1699        # Combine steps to obtain intermediate
1700        # solution at time t^n + 0.5 h
1701        ######
1702
1703        # Combine steps
1704        self.saxpy_conserved_quantities(0.25, 0.75)
1705
1706        # Update ghosts
1707        self.update_ghosts()
1708
1709        # Set substep time
1710        self.time = initial_time + self.timestep*0.5
1711
1712        # Update vertex and edge values
1713        self.distribute_to_vertices_and_edges()
1714
1715        # Update boundary values
1716        self.update_boundary()
1717
1718        ######
1719        # Third Euler step
1720        ######
1721
1722        # Compute fluxes across each element edge
1723        self.compute_fluxes()
1724
1725        # Compute forcing terms
1726        self.compute_forcing_terms()
1727
1728        # Update conserved quantities
1729        self.update_conserved_quantities()
1730
1731        ######
1732        # Combine final and initial values
1733        # and cleanup
1734        ######
1735
1736        # Combine steps
1737        self.saxpy_conserved_quantities(2.0/3.0, 1.0/3.0)
1738
1739        # Update ghosts
1740        self.update_ghosts()
1741
1742        # Set new time
1743        self.time = initial_time + self.timestep
1744
1745
1746    ##
1747    # @brief Evolve simulation to a final time.
1748    # @param finaltime Sinulation final time.
1749    def evolve_to_end(self, finaltime=1.0):
1750        """Iterate evolve all the way to the end."""
1751
1752        for _ in self.evolve(yieldstep=None, finaltime=finaltime):
1753            pass
1754
1755    ##
1756    # @brief Backup conserved quantities
1757    def backup_conserved_quantities(self):
1758
1759        # Backup conserved_quantities centroid values
1760        for name in self.conserved_quantities:
1761            Q = self.quantities[name]
1762            Q.backup_centroid_values()
1763
1764    ##
1765    # @brief Combines current C and saved centroid values S as C = aC + bS
1766    # @param a factor in combination
1767    # @param b factor in combination
1768    def saxpy_conserved_quantities(self, a, b):
1769
1770        # Backup conserved_quantities centroid values
1771        for name in self.conserved_quantities:
1772            Q = self.quantities[name]
1773            Q.saxpy_centroid_values(a, b)
1774
1775           
1776
1777
1778    ##
1779    # @brief Mapping between conserved quantites and evolved quantities
1780    # @param Input: q_cons array of conserved quantity values
1781    # @param Input: q_evol array of current evolved quantity values
1782    # @note  Output: Updated q_evol array
1783    def  conserved_values_to_evolved_values(self, q_cons, q_evol):
1784        """Needs to be overridden by Domain subclass
1785        """
1786
1787        if len(q_cons) == len(q_evol):
1788            q_evol[:] = q_cons
1789        else:
1790            msg = 'Method conserved_values_to_evolved_values must be overridden'
1791            msg += ' by Domain subclass'
1792            raise Exception, msg
1793
1794        return q_evol
1795   
1796    ##
1797    # @brief Update boundary values for all conserved quantities.
1798    def update_boundary(self):
1799        """Go through list of boundary objects and update boundary values
1800        for all conserved quantities on boundary.
1801        It is assumed that the ordering of conserved quantities is
1802        consistent between the domain and the boundary object, i.e.
1803        the jth element of vector q must correspond to the jth conserved
1804        quantity in domain.
1805        """
1806
1807        # FIXME: Update only those that change (if that can be worked out)
1808        # FIXME: Boundary objects should not include ghost nodes.
1809        for i, ((vol_id, edge_id), B) in enumerate(self.boundary_objects):
1810            if B is None:
1811                log.critical('WARNING: Ignored boundary segment (None)')
1812            else:
1813                q_bdry = B.evaluate(vol_id, edge_id)
1814
1815                if len(q_bdry) == len(self.evolved_quantities):
1816                    # conserved and evolved quantities are the same
1817                    q_evol = q_bdry
1818                elif len(q_bdry) == len(self.conserved_quantities):
1819                    # boundary just returns conserved quantities
1820                    # Need to calculate all the evolved quantities
1821                    # Use default conversion
1822
1823                    q_evol = self.get_evolved_quantities(vol_id, edge = edge_id)
1824
1825                    q_evol = self.conserved_values_to_evolved_values \
1826                                                            (q_bdry, q_evol)
1827                else:
1828                    msg = 'Boundary must return array of either conserved'
1829                    msg += ' or evolved quantities'
1830                    raise Exception, msg
1831               
1832                for j, name in enumerate(self.evolved_quantities):
1833                    Q = self.quantities[name]
1834                    Q.boundary_values[i] = q_evol[j]
1835
1836    ##
1837    # @brief Compute fluxes.
1838    # @note MUST BE OVERRIDEN IN SUBCLASS!
1839    def compute_fluxes(self):
1840        msg = 'Method compute_fluxes must be overridden by Domain subclass'
1841        raise Exception, msg
1842
1843    ##
1844    # @brief
1845    # @param yieldstep
1846    # @param finaltime
1847    def update_timestep(self, yieldstep, finaltime):
1848
1849        # Protect against degenerate timesteps arising from isolated
1850        # triangles
1851        self.apply_protection_against_isolated_degenerate_timesteps()
1852               
1853        # self.timestep is calculated from speed of characteristics
1854        # Apply CFL condition here
1855        timestep = min(self.CFL*self.flux_timestep, self.evolve_max_timestep)
1856
1857        # Record maximal and minimal values of timestep for reporting
1858        self.recorded_max_timestep = max(timestep, self.recorded_max_timestep)
1859        self.recorded_min_timestep = min(timestep, self.recorded_min_timestep)
1860
1861        # Protect against degenerate time steps
1862        if timestep < self.evolve_min_timestep:
1863            # Number of consecutive small steps taken b4 taking action
1864            self.smallsteps += 1
1865
1866            if self.smallsteps > self.max_smallsteps:
1867                self.smallsteps = 0 # Reset
1868
1869                if self._order_ == 1:
1870                    msg = 'WARNING: Too small timestep %.16f reached ' \
1871                              % timestep
1872                    msg += 'even after %d steps of 1 order scheme' \
1873                               % self.max_smallsteps
1874                    log.critical(msg)
1875                    timestep = self.evolve_min_timestep  # Try enforce min_step
1876
1877                    stats = self.timestepping_statistics(track_speeds=True)
1878                    log.critical(stats)
1879
1880                    raise Exception, msg
1881                else:
1882                    # Try to overcome situation by switching to 1 order
1883                    self._order_ = 1
1884        else:
1885            self.smallsteps = 0
1886            if self._order_ == 1 and self.default_order == 2:
1887                self._order_ = 2
1888
1889        # Ensure that final time is not exceeded
1890        if finaltime is not None and self.time + timestep > finaltime :
1891            timestep = finaltime-self.time
1892
1893        # Ensure that model time is aligned with yieldsteps
1894        if self.time + timestep > self.yieldtime:
1895            timestep = self.yieldtime - self.time
1896
1897        self.timestep = timestep
1898
1899    ##
1900    # @brief Compute forcing terms, if any.
1901    def compute_forcing_terms(self):
1902        """If there are any forcing functions driving the system
1903        they should be defined in Domain subclass and appended to
1904        the list self.forcing_terms
1905        """
1906
1907        # The parameter self.flux_timestep should be updated
1908        # by the forcing_terms to ensure stability
1909
1910        for f in self.forcing_terms:
1911            f(self)
1912
1913
1914    ##
1915    # @brief Update vectors of conserved quantities.
1916    def update_conserved_quantities(self):
1917        """Update vectors of conserved quantities using previously
1918        computed fluxes and specified forcing functions.
1919        """
1920
1921        N = len(self) # Number_of_triangles
1922        d = len(self.conserved_quantities)
1923
1924        timestep = self.timestep
1925
1926
1927        # Update conserved_quantities
1928        for name in self.conserved_quantities:
1929            Q = self.quantities[name]
1930            Q.update(timestep)
1931
1932            # Note that Q.explicit_update is reset by compute_fluxes
1933            # Where is Q.semi_implicit_update reset?
1934            # It is reset in quantity_ext.c
1935
1936    ##
1937    # @brief Sequential update of ghost cells
1938    def update_ghosts(self):
1939        # We must send the information from the full cells and
1940        # receive the information for the ghost cells
1941        # We have a list with ghosts expecting updates
1942
1943        #Update of ghost cells
1944        iproc = self.processor
1945        if self.full_send_dict.has_key(iproc):
1946
1947            # now store full as local id, global id, value
1948            Idf  = self.full_send_dict[iproc][0]
1949
1950            # now store ghost as local id, global id, value
1951            Idg = self.ghost_recv_dict[iproc][0]
1952
1953            for i, q in enumerate(self.conserved_quantities):
1954                Q_cv =  self.quantities[q].centroid_values
1955                num.put(Q_cv, Idg, num.take(Q_cv, Idf, axis=0))
1956
1957 
1958    ##
1959    # @brief Extrapolate conserved quantities from centroid to vertices
1960    #        and edge-midpoints for each volume.
1961    def distribute_to_vertices_and_edges(self):
1962        """Extrapolate conserved quantities from centroid to
1963        vertices and edge-midpoints for each volume
1964
1965        Default implementation is straight first order,
1966        i.e. constant values throughout each element and
1967        no reference to non-conserved quantities.
1968        """
1969
1970        for name in self.conserved_quantities:
1971            Q = self.quantities[name]
1972            if self._order_ == 1:
1973                Q.extrapolate_first_order()
1974            elif self._order_ == 2:
1975                Q.extrapolate_second_order()
1976            else:
1977                raise Exception, 'Unknown order'
1978
1979    ##
1980    # @brief Calculate the norm of the centroid values of a specific quantity,
1981    #        using normfunc.
1982    # @param quantity
1983    # @param normfunc
1984    def centroid_norm(self, quantity, normfunc):
1985        """Calculate the norm of the centroid values of a specific quantity,
1986        using normfunc.
1987
1988        normfunc should take a list to a float.
1989
1990        common normfuncs are provided in the module utilities.norms
1991        """
1992
1993        return normfunc(self.quantities[quantity].centroid_values)
1994
1995
1996
1997    def apply_protection_against_isolated_degenerate_timesteps(self):
1998
1999        # FIXME (Steve): This should be in shallow_water as it assumes x and y
2000        # momentum
2001        if self.protect_against_isolated_degenerate_timesteps is False:
2002            return
2003       
2004        # FIXME (Ole): Make this configurable
2005        if num.max(self.max_speed) < 10.0: 
2006            return
2007
2008        # Setup 10 bins for speed histogram
2009        from anuga.utilities.numerical_tools import histogram, create_bins
2010
2011        bins = create_bins(self.max_speed, 10)
2012        hist = histogram(self.max_speed, bins)
2013
2014        # Look for characteristic signature
2015        if len(hist) > 1 and hist[-1] > 0 and \
2016            hist[4] == hist[5] == hist[6] == hist[7] == hist[8] == 0:
2017            # Danger of isolated degenerate triangles
2018
2019            # Find triangles in last bin
2020            # FIXME - speed up using numeric package
2021            d = 0
2022            for i in range(self.number_of_full_triangles):
2023                if self.max_speed[i] > bins[-1]:
2024                    msg = 'Time=%f: Ignoring isolated high ' % self.time
2025                    msg += 'speed triangle '
2026                    msg += '#%d of %d with max speed=%f' \
2027                        % (i, self.number_of_full_triangles, self.max_speed[i])
2028
2029                    self.get_quantity('xmomentum').\
2030                        set_values(0.0, indices=[i])
2031                    self.get_quantity('ymomentum').\
2032                        set_values(0.0, indices=[i])
2033                    self.max_speed[i]=0.0
2034                    d += 1
2035
2036
2037######
2038# Initialise module
2039######
2040
2041# Optimisation with psyco
2042#from anuga.config import use_psyco
2043
2044#if use_psyco:
2045    #try:
2046        #import psyco
2047    #except:
2048        #import os
2049        #if os.name == 'posix' and os.uname()[4] in ['x86_64', 'ia64']:
2050            #pass
2051            ## Psyco isn't supported on 64 bit systems, but it doesn't matter
2052        #else:
2053            #log.critical('WARNING: psyco (speedup) could not be imported, '
2054                         #'you may want to consider installing it')
2055    #else:
2056        #psyco.bind(Generic_Domain.update_boundary)
2057        ##psyco.bind(Domain.update_timestep) # Not worth it
2058        #psyco.bind(Generic_Domain.update_conserved_quantities)
2059        #psyco.bind(Generic_Domain.distribute_to_vertices_and_edges)
2060
2061
2062if __name__ == "__main__":
2063    pass
Note: See TracBrowser for help on using the repository browser.