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

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

Changed import of boundary conditions

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