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

Last change on this file since 7780 was 7780, checked in by hudson, 13 years ago

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