source: anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py @ 7562

Last change on this file since 7562 was 7562, checked in by steve, 14 years ago

Updating the balanced and parallel code

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