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

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

revert back to previous until tests pass

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