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

Last change on this file since 4711 was 4711, checked in by ole, 17 years ago

Added tolerance for extremum computations

File size: 47.0 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
10from Numeric import allclose, argmax
11from anuga.config import epsilon
12
13from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
14from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
15     import Boundary
16from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
17     import File_boundary
18from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
19     import Dirichlet_boundary
20from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
21     import Time_boundary
22from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
23     import Transmissive_boundary
24
25from anuga.abstract_2d_finite_volumes.pmesh2domain import pmesh_to_domain
26from anuga.abstract_2d_finite_volumes.region\
27     import Set_region as region_set_region
28
29from anuga.utilities.polygon import inside_polygon
30from anuga.abstract_2d_finite_volumes.util import get_textual_float
31
32import types
33
34class Domain(Mesh):
35
36
37    def __init__(self,
38                 source=None,
39                 triangles=None,
40                 boundary=None,
41                 conserved_quantities=None,
42                 other_quantities=None,
43                 tagged_elements=None,
44                 geo_reference=None,
45                 use_inscribed_circle=False,
46                 mesh_filename=None,
47                 use_cache=False,
48                 verbose=False,
49                 full_send_dict=None,
50                 ghost_recv_dict=None,
51                 processor=0,
52                 numproc=1,
53                 number_of_full_nodes=None,
54                 number_of_full_triangles=None):
55
56
57        """Instantiate generic computational Domain.
58
59        Input:
60          source:    Either a mesh filename or coordinates of mesh vertices.
61                     If it is a filename values specified for triangles will
62                     be overridden.
63          triangles: Mesh connectivity (see mesh.py for more information)
64          boundary:  See mesh.py for more information
65
66          conserved_quantities: List of quantity names entering the
67                                conservation equations
68          other_quantities:     List of other quantity names
69
70          tagged_elements:
71          ...
72
73
74        """
75
76        # Determine whether source is a mesh filename or coordinates
77        if type(source) == types.StringType:
78            mesh_filename = source
79        else:
80            coordinates = source
81
82
83        # In case a filename has been specified, extract content
84        if mesh_filename is not None:
85            coordinates, triangles, boundary, vertex_quantity_dict, \
86                         tagged_elements, geo_reference = \
87                         pmesh_to_domain(file_name=mesh_filename,
88                                         use_cache=use_cache,
89                                         verbose=verbose)
90
91
92        # Initialise underlying mesh structure
93        Mesh.__init__(self, coordinates, triangles,
94                      boundary=boundary,
95                      tagged_elements=tagged_elements,
96                      geo_reference=geo_reference,
97                      use_inscribed_circle=use_inscribed_circle,
98                      number_of_full_nodes=number_of_full_nodes,
99                      number_of_full_triangles=number_of_full_triangles,
100                      verbose=verbose)
101
102        if verbose: print 'Initialising Domain'
103        from Numeric import zeros, Float, Int, ones
104        from quantity import Quantity, Conserved_quantity
105
106        # List of quantity names entering
107        # the conservation equations
108        if conserved_quantities is None:
109            self.conserved_quantities = []
110        else:
111            self.conserved_quantities = conserved_quantities
112
113        # List of other quantity names
114        if other_quantities is None:
115            self.other_quantities = []
116        else:
117            self.other_quantities = other_quantities
118
119
120        #Build dictionary of Quantity instances keyed by quantity names
121        self.quantities = {}
122
123        #FIXME: remove later - maybe OK, though....
124        for name in self.conserved_quantities:
125            self.quantities[name] = Conserved_quantity(self)
126        for name in self.other_quantities:
127            self.quantities[name] = Quantity(self)
128
129        #Create an empty list for explicit forcing terms
130        self.forcing_terms = []
131
132        #Setup the ghost cell communication
133        if full_send_dict is None:
134            self.full_send_dict = {}
135        else:
136            self.full_send_dict  = full_send_dict
137
138        # List of other quantity names
139        if ghost_recv_dict is None:
140            self.ghost_recv_dict = {}
141        else:
142            self.ghost_recv_dict = ghost_recv_dict
143
144        self.processor = processor
145        self.numproc = numproc
146
147
148        # Setup Communication Buffers
149        if verbose: print 'Domain: Set up communication buffers (parallel)'
150        self.nsys = len(self.conserved_quantities)
151        for key in self.full_send_dict:
152            buffer_shape = self.full_send_dict[key][0].shape[0]
153            self.full_send_dict[key].append(zeros( (buffer_shape,self.nsys) ,Float))
154
155
156        for key in self.ghost_recv_dict:
157            buffer_shape = self.ghost_recv_dict[key][0].shape[0]
158            self.ghost_recv_dict[key].append(zeros( (buffer_shape,self.nsys) ,Float))
159
160
161        # Setup cell full flag
162        # =1 for full
163        # =0 for ghost
164        N = len(self) #number_of_elements
165        self.tri_full_flag = ones(N, Int)
166        for i in self.ghost_recv_dict.keys():
167            for id in self.ghost_recv_dict[i][0]:
168                self.tri_full_flag[id] = 0
169
170        # Test the assumption that all full triangles are store before
171        # the ghost triangles.
172        assert allclose(self.tri_full_flag[:self.number_of_full_nodes],1)
173                       
174
175        # Defaults
176        from anuga.config import max_smallsteps, beta_w, beta_h, epsilon
177        from anuga.config import CFL
178        from anuga.config import protect_against_isolated_degenerate_timesteps
179        self.beta_w = beta_w
180        self.beta_h = beta_h
181        self.epsilon = epsilon
182        self.protect_against_isolated_degenerate_timesteps = protect_against_isolated_degenerate_timesteps
183       
184
185        # FIXME: Maybe have separate orders for h-limiter and w-limiter?
186        # Or maybe get rid of order altogether and use beta_w and beta_h
187        self.set_default_order(1)
188        #self.default_order = 1
189        #self._order_ = self.default_order
190
191        self.smallsteps = 0
192        self.max_smallsteps = max_smallsteps
193        self.number_of_steps = 0
194        self.number_of_first_order_steps = 0
195        self.CFL = CFL
196
197        self.boundary_map = None  # Will be populated by set_boundary       
198       
199
200        # Model time
201        self.time = 0.0
202        self.finaltime = None
203        self.min_timestep = self.max_timestep = 0.0
204        self.starttime = 0 #Physical starttime if any (0 is 1 Jan 1970 00:00:00)
205
206        # Monitoring
207        self.quantities_to_be_monitored = None
208        self.monitor_polygon = None
209        self.monitor_time_interval = None               
210        self.monitor_indices = None
211
212        # Checkpointing and storage
213        from anuga.config import default_datadir
214        self.datadir = default_datadir
215        self.simulation_name = 'domain'
216        self.checkpoint = False
217
218        # MH310505 To avoid calculating the flux across each edge twice, keep an integer (boolean) array,
219        # to be used during the flux calculation
220        N = len(self) #number_of_triangles
221        self.already_computed_flux = zeros((N, 3), Int)
222
223        # Storage for maximal speeds computed for each triangle by compute_fluxes
224        # This is used for diagnostics only
225        self.max_speed = zeros(N, Float)
226
227        if mesh_filename is not None:
228            # If the mesh file passed any quantity values
229            # , initialise with these values.
230            if verbose: print 'Domain: Initialising quantity values'
231            self.set_quantity_vertices_dict(vertex_quantity_dict)
232
233
234        if verbose: print 'Domain: Done'
235
236
237
238
239    def set_default_order(self, n):
240        """Set default (spatial) order to either 1 or 2
241        """
242
243        msg = 'Default order must be either 1 or 2. I got %s' %n
244        assert n in [1,2], msg
245
246        self.default_order = n
247        self._order_ = self.default_order
248
249
250    #Public interface to Domain
251    def get_conserved_quantities(self, vol_id, vertex=None, edge=None):
252        """Get conserved quantities at volume vol_id
253
254        If vertex is specified use it as index for vertex values
255        If edge is specified use it as index for edge values
256        If neither are specified use centroid values
257        If both are specified an exeception is raised
258
259        Return value: Vector of length == number_of_conserved quantities
260
261        """
262
263        from Numeric import zeros, Float
264
265        if not (vertex is None or edge is None):
266            msg = 'Values for both vertex and edge was specified.'
267            msg += 'Only one (or none) is allowed.'
268            raise msg
269
270        q = zeros( len(self.conserved_quantities), Float)
271
272        for i, name in enumerate(self.conserved_quantities):
273            Q = self.quantities[name]
274            if vertex is not None:
275                q[i] = Q.vertex_values[vol_id, vertex]
276            elif edge is not None:
277                q[i] = Q.edge_values[vol_id, edge]
278            else:
279                q[i] = Q.centroid_values[vol_id]
280
281        return q
282
283    def set_time(self, time=0.0):
284        """Set the model time (seconds)"""
285
286        self.time = time
287
288    def get_time(self):
289        """Get the model time (seconds)"""
290
291        return self.time
292
293    def set_quantity_vertices_dict(self, quantity_dict):
294        """Set values for named quantities.
295        The index is the quantity
296
297        name: Name of quantity
298        X: Compatible list, Numeric array, const or function (see below)
299
300        The values will be stored in elements following their
301        internal ordering.
302
303        """
304       
305        # FIXME: Could we name this a bit more intuitively
306        # E.g. set_quantities_from_dictionary
307        for key in quantity_dict.keys():
308            self.set_quantity(key, quantity_dict[key], location='vertices')
309
310
311    def set_quantity(self, name, *args, **kwargs):
312        """Set values for named quantity
313
314
315        One keyword argument is documented here:
316        expression = None, # Arbitrary expression
317
318        expression:
319          Arbitrary expression involving quantity names
320
321        See Quantity.set_values for further documentation.
322        """
323
324        #FIXME (Ole): Allow new quantities here
325        #from quantity import Quantity, Conserved_quantity
326        #Create appropriate quantity object
327        ##if name in self.conserved_quantities:
328        ##    self.quantities[name] = Conserved_quantity(self)
329        ##else:
330        ##    self.quantities[name] = Quantity(self)
331
332
333        #Do the expression stuff
334        if kwargs.has_key('expression'):
335            expression = kwargs['expression']
336            del kwargs['expression']
337
338            Q = self.create_quantity_from_expression(expression)
339            kwargs['quantity'] = Q
340
341        #Assign values
342        self.quantities[name].set_values(*args, **kwargs)
343
344
345    def get_quantity_names(self):
346        """Get a list of all the quantity names that this domain is aware of.
347        Any value in the result should be a valid input to get_quantity.
348        """
349        return self.quantities.keys()
350
351    def get_quantity(self, name, location='vertices', indices = None):
352        """Get quantity object.
353
354        name: Name of quantity
355
356        See methods inside the quantity object for more options
357
358        FIXME: clean input args
359        """
360
361        return self.quantities[name] #.get_values( location, indices = indices)
362
363
364    def get_quantity_object(self, name):
365        """Get object for named quantity
366
367        name: Name of quantity
368
369        FIXME: Obsolete
370        """
371
372        print 'get_quantity_object has been deprecated. Please use get_quantity'
373        return self.quantities[name]
374
375
376    def create_quantity_from_expression(self, expression):
377        """Create new quantity from other quantities using arbitrary expression
378
379        Combine existing quantities in domain using expression and return
380        result as a new quantity.
381
382        Note, the new quantity could e.g. be used in set_quantity
383
384        Valid expressions are limited to operators defined in class Quantity
385
386        Examples creating derived quantities:
387
388            Depth = domain.create_quantity_from_expression('stage-elevation')
389
390            exp = '(xmomentum*xmomentum + ymomentum*ymomentum)**0.5')       
391            Absolute_momentum = domain.create_quantity_from_expression(exp)
392
393        """
394
395        from anuga.abstract_2d_finite_volumes.util import\
396             apply_expression_to_dictionary
397       
398        return apply_expression_to_dictionary(expression, self.quantities)
399
400
401
402    #def modify_boundary(self, boundary_map):
403    #    """Modify existing boundary by elements in boundary map#
404    #
405    #    Input:#
406    #
407    #    boundary_map: Dictionary mapping tags to boundary objects
408    #
409    #    See set_boundary for more details on how this works
410    #
411    #      OBSOLETE
412    #    """
413    #
414    #    for key in boundary_map.keys():
415    #        self.boundary_map[key] = boundary_map[key]
416    #
417    #    self.set_boundary(self.boundary_map)
418       
419       
420
421    def set_boundary(self, boundary_map):
422        """Associate boundary objects with tagged boundary segments.
423
424        Input boundary_map is a dictionary of boundary objects keyed
425        by symbolic tags to matched against tags in the internal dictionary
426        self.boundary.
427
428        As result one pointer to a boundary object is stored for each vertex
429        in the list self.boundary_objects.
430        More entries may point to the same boundary object
431
432        Schematically the mapping is from two dictionaries to one list
433        where the index is used as pointer to the boundary_values arrays
434        within each quantity.
435
436        self.boundary:          (vol_id, edge_id): tag
437        boundary_map (input):   tag: boundary_object
438        ----------------------------------------------
439        self.boundary_objects:  ((vol_id, edge_id), boundary_object)
440
441
442        Pre-condition:
443          self.boundary has been built.
444
445        Post-condition:
446          self.boundary_objects is built
447
448        If a tag from the domain doesn't appear in the input dictionary an
449        exception is raised.
450        However, if a tag is not used to the domain, no error is thrown.
451        FIXME: This would lead to implementation of a
452        default boundary condition
453
454        Note: If a segment is listed in the boundary dictionary and if it is
455        not None, it *will* become a boundary -
456        even if there is a neighbouring triangle.
457        This would be the case for internal boundaries
458
459        Boundary objects that are None will be skipped.
460
461        If a boundary_map has already been set
462        (i.e. set_boundary has been called before), the old boundary map
463        will be updated with new values. The new map need not define all
464        boundary tags, and can thus change only those that are needed.
465
466        FIXME: If set_boundary is called multiple times and if Boundary
467        object is changed into None, the neighbour structure will not be
468        restored!!!
469
470
471        """
472
473        if self.boundary_map is None:
474            # This the first call to set_boundary. Store
475            # map for later updates and for use with boundary_stats.
476            self.boundary_map = boundary_map
477        else:   
478            # This is a modification of an already existing map
479            # Update map an proceed normally
480
481            for key in boundary_map.keys():
482                self.boundary_map[key] = boundary_map[key]
483               
484       
485        #FIXME: Try to remove the sorting and fix test_mesh.py
486        x = self.boundary.keys()
487        x.sort()
488
489        #Loop through edges that lie on the boundary and associate them with
490        #callable boundary objects depending on their tags
491        self.boundary_objects = []       
492        for k, (vol_id, edge_id) in enumerate(x):
493            tag = self.boundary[ (vol_id, edge_id) ]
494
495            if self.boundary_map.has_key(tag):
496                B = self.boundary_map[tag]  #Get callable boundary object
497
498                if B is not None:
499                    self.boundary_objects.append( ((vol_id, edge_id), B) )
500                    self.neighbours[vol_id, edge_id] = \
501                                            -len(self.boundary_objects)
502                else:
503                    pass
504                    #FIXME: Check and perhaps fix neighbour structure
505
506
507            else:
508                msg = 'ERROR (domain.py): Tag "%s" has not been ' %tag
509                msg += 'bound to a boundary object.\n'
510                msg += 'All boundary tags defined in domain must appear '
511                msg += 'in set_boundary.\n'
512                msg += 'The tags are: %s' %self.get_boundary_tags()
513                raise msg
514
515
516    def set_region(self, *args, **kwargs):
517        """
518        This method is used to set quantities based on a regional tag.
519       
520        It is most often called with the following parameters;
521        (self, tag, quantity, X, location='vertices')
522        tag: the name of the regional tag used to specify the region
523        quantity: Name of quantity to change
524        X: const or function - how the quantity is changed
525        location: Where values are to be stored.
526            Permissible options are: vertices, centroid and unique vertices
527
528        A callable region class or a list of callable region classes
529        can also be passed into this function.
530        """
531        #print "*args", args
532        #print "**kwargs", kwargs
533        if len(args) == 1:
534            self._set_region(*args, **kwargs)
535        else:
536            #Assume it is arguments for the region.set_region function
537            func = region_set_region(*args, **kwargs)
538            self._set_region(func)
539           
540       
541    def _set_region(self, functions):
542        # The order of functions in the list is used.
543        if type(functions) not in [types.ListType,types.TupleType]:
544            functions = [functions]
545        for function in functions:
546            for tag in self.tagged_elements.keys():
547                function(tag, self.tagged_elements[tag], self)
548
549
550
551
552    def set_quantities_to_be_monitored(self, q,
553                                       polygon=None,
554                                       time_interval=None):
555        """Specify which quantities will be monitored for extrema.
556
557        q must be either:
558          - the name of a quantity
559          - a list of quantity names
560          - None
561
562        In the two first cases, the named quantities will be monitored at
563        each internal timestep
564       
565        If q is None, monitoring will be switched off altogether.
566
567        polygon (if specified) will restrict monitoring to triangles inside polygon.
568        If omitted all triangles will be included.
569
570        time_interval (if specified) will restrict monitoring to time steps in
571        that interval. If omitted all timesteps will be included.
572
573        FIXME: Derived quantities such as 'stage-elevation' will appear later
574        """
575
576        # FIXME (Ole): This is under construction. See ticket:192
577
578        from anuga.abstract_2d_finite_volumes.util import\
579             apply_expression_to_dictionary       
580
581        if q is None:
582            self.quantities_to_be_monitored = None
583            self.monitor_polygon = None
584            self.monitor_time_interval = None
585            self.monitor_indices = None           
586            return
587
588        if isinstance(q, basestring):
589            q = [q] # Turn argument into a list
590
591        # Check correcness and initialise
592        self.quantities_to_be_monitored = {}
593        for quantity_name in q:
594            msg = 'Quantity %s is not a valid conserved quantity'\
595                  %quantity_name
596           
597
598            if not quantity_name in self.quantities:
599                # See if this expression is valid
600                apply_expression_to_dictionary(quantity_name, self.quantities)
601
602            # Initialise extrema information
603            info_block = {'min': None,          # Min value
604                          'max': None,          # Max value
605                          'min_location': None, # Argmin (x, y)
606                          'max_location': None, # Argmax (x, y)
607                          'min_time': None,     # Argmin (t)
608                          'max_time': None}     # Argmax (t)
609           
610            self.quantities_to_be_monitored[quantity_name] = info_block
611
612           
613
614        if polygon is not None:
615            # Check input
616            if isinstance(polygon, basestring):
617
618                # Check if multiple quantities were accidentally
619                # given as separate argument rather than a list.
620                msg = 'Multiple quantities must be specified in a list. '
621                msg += 'Not as multiple arguments. '
622                msg += 'I got "%s" as a second argument' %polygon
623               
624                if polygon in self.quantities:
625                    raise Exception, msg
626               
627                try:
628                    apply_expression_to_dictionary(polygon, self.quantities)
629                except:
630                    # At least polygon wasn't an expression involving quantitites
631                    pass
632                else:
633                    raise Exception, msg
634
635                # In any case, we don't allow polygon to be a string
636                msg = 'argument "polygon" must not be a string: '
637                msg += 'I got polygon=\'%s\' ' %polygon
638                raise Exception, msg
639
640
641            # Get indices for centroids that are inside polygon
642            points = self.get_centroid_coordinates(absolute=True)
643            self.monitor_indices = inside_polygon(points, polygon)
644           
645
646        if time_interval is not None:
647            assert len(time_interval) == 2
648
649       
650        self.monitor_polygon = polygon
651        self.monitor_time_interval = time_interval       
652       
653
654
655    #MISC
656    def check_integrity(self):
657        Mesh.check_integrity(self)
658
659        for quantity in self.conserved_quantities:
660            msg = 'Conserved quantities must be a subset of all quantities'
661            assert quantity in self.quantities, msg
662
663        ##assert hasattr(self, 'boundary_objects')
664
665    def write_time(self, track_speeds=False):
666        print self.timestepping_statistics(track_speeds)
667
668
669    def timestepping_statistics(self, track_speeds=False):
670        """Return string with time stepping statistics for printing or logging
671
672        Optional boolean keyword track_speeds decides whether to report location of
673        smallest timestep as well as a histogram and percentile report.
674        """
675
676        from anuga.utilities.numerical_tools import histogram, create_bins
677       
678
679        msg = ''
680        if self.min_timestep == self.max_timestep:
681            msg += 'Time = %.4f, delta t = %.8f, steps=%d (%d)'\
682                   %(self.time, self.min_timestep, self.number_of_steps,
683                     self.number_of_first_order_steps)
684        elif self.min_timestep > self.max_timestep:
685            msg += 'Time = %.4f, steps=%d (%d)'\
686                   %(self.time, self.number_of_steps,
687                     self.number_of_first_order_steps)
688        else:
689            msg += 'Time = %.4f, delta t in [%.8f, %.8f], steps=%d (%d)'\
690                   %(self.time, self.min_timestep,
691                     self.max_timestep, self.number_of_steps,
692                     self.number_of_first_order_steps)
693
694        if track_speeds is True:
695            msg += '\n'
696
697
698            #Setup 10 bins for speed histogram
699            bins = create_bins(self.max_speed, 10)
700            hist = histogram(self.max_speed, bins)
701
702            msg += '------------------------------------------------\n'
703            msg += '  Speeds in [%f, %f]\n' %(min(self.max_speed), max(self.max_speed))
704            msg += '  Histogram:\n'
705
706            hi = bins[0]
707            for i, count in enumerate(hist):
708                lo = hi
709                if i+1 < len(bins):
710                    #Open upper interval               
711                    hi = bins[i+1]
712                    msg += '    [%f, %f[: %d\n' %(lo, hi, count)               
713                else:
714                    #Closed upper interval
715                    hi = max(self.max_speed)
716                    msg += '    [%f, %f]: %d\n' %(lo, hi, count)
717
718
719            N = len(self.max_speed)
720            if N > 10:
721                msg += '  Percentiles (10%):\n'
722                speed = self.max_speed.tolist()
723                speed.sort()
724
725                k = 0
726                lower = min(speed)
727                for i, a in enumerate(speed):       
728                    if i % (N/10) == 0 and i != 0: #For every 10% of the sorted speeds
729                        msg += '    %d speeds in [%f, %f]\n' %(i-k, lower, a)
730                        lower = a
731                        k = i
732                       
733                msg += '    %d speeds in [%f, %f]\n'\
734                       %(N-k, lower, max(speed))                   
735               
736                     
737           
738
739               
740           
741            # Find index of largest computed flux speed
742            k = argmax(self.max_speed)
743
744            x, y = self.get_centroid_coordinates()[k]
745            radius = self.get_radii()[k]
746            area = self.get_areas()[k]     
747            max_speed = self.max_speed[k]           
748
749            msg += '  Triangle #%d with centroid (%.4f, %.4f), ' %(k, x, y)
750            msg += 'area = %.4f and radius = %.4f ' %(area, radius)
751            msg += 'had the largest computed speed: %.6f m/s ' %(max_speed)
752            if max_speed > 0.0:
753                msg += '(timestep=%.6f)\n' %(radius/max_speed)
754            else:
755                msg += '(timestep=%.6f)\n' %(0)     
756           
757            # Report all quantity values at vertices
758            msg += '    Quantity \t vertex values\t\t\t\t\t centroid values\n'
759            for name in self.quantities:
760                q = self.quantities[name]
761
762                V = q.get_values(location='vertices', indices=[k])[0]
763                C = q.get_values(location='centroids', indices=[k])                 
764               
765                s = '    %s:\t %.4f,\t %.4f,\t %.4f,\t %.4f\n'\
766                    %(name, V[0], V[1], V[2], C[0])                 
767
768                msg += s
769
770        return msg
771
772
773    def write_boundary_statistics(self, quantities = None, tags = None):
774        print self.boundary_statistics(quantities, tags)
775
776    def boundary_statistics(self, quantities = None, tags = None):
777        """Output statistics about boundary forcing at each timestep
778
779
780        Input:
781          quantities: either None, a string or a list of strings naming the quantities to be reported
782          tags:       either None, a string or a list of strings naming the tags to be reported
783
784
785        Example output:
786        Tag 'wall':
787            stage in [2, 5.5]
788            xmomentum in []
789            ymomentum in []
790        Tag 'ocean'
791
792
793        If quantities are specified only report on those. Otherwise take all conserved quantities.
794        If tags are specified only report on those, otherwise take all tags.
795
796        """
797
798        # Input checks
799        import types, string
800
801        if quantities is None:
802            quantities = self.conserved_quantities
803        elif type(quantities) == types.StringType:
804            quantities = [quantities] #Turn it into a list
805
806        msg = 'Keyword argument quantities must be either None, '
807        msg += 'string or list. I got %s' %str(quantities)
808        assert type(quantities) == types.ListType, msg
809
810
811        if tags is None:
812            tags = self.get_boundary_tags()
813        elif type(tags) == types.StringType:
814            tags = [tags] #Turn it into a list
815
816        msg = 'Keyword argument tags must be either None, '
817        msg += 'string or list. I got %s' %str(tags)
818        assert type(tags) == types.ListType, msg
819
820        # Determine width of longest quantity name (for cosmetic purposes)
821        maxwidth = 0
822        for name in quantities:
823            w = len(name)
824            if w > maxwidth:
825                maxwidth = w
826
827        # Output stats
828        msg = 'Boundary values at time %.4f:\n' %self.time
829        for tag in tags:
830            msg += '    %s:\n' %tag
831
832            for name in quantities:
833                q = self.quantities[name]
834
835                # Find range of boundary values for tag and q
836                maxval = minval = None
837                for i, ((vol_id, edge_id), B) in\
838                        enumerate(self.boundary_objects):
839                    if self.boundary[(vol_id, edge_id)] == tag:
840                        v = q.boundary_values[i]
841                        if minval is None or v < minval: minval = v
842                        if maxval is None or v > maxval: maxval = v
843
844                if minval is None or maxval is None:
845                    msg += '        Sorry no information available about' +\
846                           ' tag %s and quantity %s\n' %(tag, name)
847                else:
848                    msg += '        %s in [%12.8f, %12.8f]\n'\
849                           %(string.ljust(name, maxwidth), minval, maxval)
850
851
852        return msg
853
854
855    def update_extrema(self):
856        """Update extrema if requested by set_quantities_to_be_monitored.
857        This data is used for reporting e.g. by running
858        print domain.quantity_statistics()
859        and may also stored in output files (see data_manager in shallow_water)
860        """
861
862        # Define a tolerance for extremum computations
863        epsilon = 1.0e-6 # Import 'single_precision' from config
864       
865        if self.quantities_to_be_monitored is None:
866            return
867
868        # Observe time interval restriction if any
869        if self.monitor_time_interval is not None and\
870               (self.time < self.monitor_time_interval[0] or\
871               self.time > self.monitor_time_interval[1]):
872            return
873
874        # Update extrema for each specified quantity subject to
875        # polygon restriction (via monitor_indices).
876        for quantity_name in self.quantities_to_be_monitored:
877
878            if quantity_name in self.quantities:
879                Q = self.get_quantity(quantity_name)
880            else:
881                Q = self.create_quantity_from_expression(quantity_name)
882
883            info_block = self.quantities_to_be_monitored[quantity_name]
884
885            # Update maximum
886            # (n > None is always True, but we check explicitly because of the epsilon)
887            maxval = Q.get_maximum_value(self.monitor_indices)
888            if info_block['max'] is None or\
889                   maxval > info_block['max'] + epsilon:
890                info_block['max'] = maxval
891                maxloc = Q.get_maximum_location()
892                info_block['max_location'] = maxloc
893                info_block['max_time'] = self.time
894
895
896            # Update minimum
897            minval = Q.get_minimum_value(self.monitor_indices)
898            if info_block['min'] is None or\
899                   minval < info_block['min'] - epsilon:
900                info_block['min'] = minval               
901                minloc = Q.get_minimum_location()
902                info_block['min_location'] = minloc
903                info_block['min_time'] = self.time               
904       
905
906
907    def quantity_statistics(self, precision = '%.4f'):
908        """Return string with statistics about quantities for printing or logging
909
910        Quantities reported are specified through method
911
912           set_quantities_to_be_monitored
913           
914        """
915
916        maxlen = 128 # Max length of polygon string representation
917
918        # Output stats
919        msg = 'Monitored quantities at time %.4f:\n' %self.time
920        if self.monitor_polygon is not None:
921            p_str = str(self.monitor_polygon)
922            msg += '- Restricted by polygon: %s' %p_str[:maxlen]
923            if len(p_str) >= maxlen:
924                msg += '...\n'
925            else:
926                msg += '\n'
927
928
929        if self.monitor_time_interval is not None:
930            msg += '- Restricted by time interval: %s\n' %str(self.monitor_time_interval)
931            time_interval_start = self.monitor_time_interval[0]
932        else:
933            time_interval_start = 0.0
934
935           
936        for quantity_name, info in self.quantities_to_be_monitored.items():
937            msg += '    %s:\n' %quantity_name
938
939            msg += '      values since time = %.2f in [%s, %s]\n'\
940                   %(time_interval_start,
941                     get_textual_float(info['min'], precision),
942                     get_textual_float(info['max'], precision))                     
943                     
944            msg += '      minimum attained at time = %s, location = %s\n'\
945                   %(get_textual_float(info['min_time'], precision),
946                     get_textual_float(info['min_location'], precision))                     
947                     
948
949            msg += '      maximum attained at time = %s, location = %s\n'\
950                   %(get_textual_float(info['max_time'], precision),
951                     get_textual_float(info['max_location'], precision))                   
952
953
954        return msg
955
956
957    def get_name(self):
958        return self.simulation_name
959
960    def set_name(self, name):
961        """Assign a name to this simulation.
962        This will be used to identify the output sww file.
963
964        """
965        if name.endswith('.sww'):
966            name = name[:-4]
967           
968        self.simulation_name = name
969
970    def get_datadir(self):
971        return self.datadir
972
973    def set_datadir(self, name):
974        self.datadir = name
975
976
977
978    #def set_defaults(self):
979    #    """Set default values for uninitialised quantities.
980    #    Should be overridden or specialised by specific modules
981    #    """#
982    #
983    #    for name in self.conserved_quantities + self.other_quantities:
984    #        self.set_quantity(name, 0.0)
985
986
987    ###########################
988    #Main components of evolve
989
990    def evolve(self,
991               yieldstep = None,
992               finaltime = None,
993               duration = None,
994               skip_initial_step = False):
995        """Evolve model through time starting from self.starttime.
996
997
998        yieldstep: Interval between yields where results are stored,
999                   statistics written and domain inspected or
1000                   possibly modified. If omitted the internal predefined
1001                   max timestep is used.
1002                   Internally, smaller timesteps may be taken.
1003
1004        duration: Duration of simulation
1005
1006        finaltime: Time where simulation should end. This is currently
1007        relative time.  So it's the same as duration.
1008
1009        If both duration and finaltime are given an exception is thrown.
1010
1011
1012        skip_initial_step: Boolean flag that decides whether the first
1013        yield step is skipped or not. This is useful for example to avoid
1014        duplicate steps when multiple evolve processes are dove tailed.
1015
1016
1017        Evolve is implemented as a generator and is to be called as such, e.g.
1018
1019        for t in domain.evolve(yieldstep, finaltime):
1020            <Do something with domain and t>
1021
1022
1023        All times are given in seconds
1024
1025        """
1026
1027        from anuga.config import min_timestep, max_timestep, epsilon
1028
1029        # FIXME: Maybe lump into a larger check prior to evolving
1030        msg = 'Boundary tags must be bound to boundary objects before '
1031        msg += 'evolving system, '
1032        msg += 'e.g. using the method set_boundary.\n'
1033        msg += 'This system has the boundary tags %s '\
1034               %self.get_boundary_tags()
1035        assert hasattr(self, 'boundary_objects'), msg
1036
1037
1038        if yieldstep is None:
1039            yieldstep = max_timestep
1040        else:
1041            yieldstep = float(yieldstep)
1042
1043        self._order_ = self.default_order
1044
1045
1046        if finaltime is not None and duration is not None:
1047            # print 'F', finaltime, duration
1048            msg = 'Only one of finaltime and duration may be specified'
1049            raise msg
1050        else:
1051            if finaltime is not None:
1052                self.finaltime = float(finaltime)
1053            if duration is not None:
1054                self.finaltime = self.starttime + float(duration)
1055
1056
1057
1058
1059        self.yieldtime = 0.0 # Track time between 'yields'
1060
1061        # Initialise interval of timestep sizes (for reporting only)
1062        self.min_timestep = max_timestep
1063        self.max_timestep = min_timestep
1064        self.number_of_steps = 0
1065        self.number_of_first_order_steps = 0
1066
1067        # Update ghosts
1068        self.update_ghosts()
1069
1070        # Initial update of vertex and edge values
1071        self.distribute_to_vertices_and_edges()
1072
1073        # Update extrema if necessary (for reporting)
1074        self.update_extrema()
1075       
1076        # Initial update boundary values
1077        self.update_boundary()
1078
1079        # Or maybe restore from latest checkpoint
1080        if self.checkpoint is True:
1081            self.goto_latest_checkpoint()
1082
1083        if skip_initial_step is False:
1084            yield(self.time)  # Yield initial values
1085
1086        while True:
1087            # Compute fluxes across each element edge
1088            self.compute_fluxes()
1089
1090            # Update timestep to fit yieldstep and finaltime
1091            self.update_timestep(yieldstep, finaltime)
1092
1093            # Update conserved quantities
1094            self.update_conserved_quantities()
1095
1096            # Update ghosts
1097            self.update_ghosts()
1098
1099            # Update vertex and edge values
1100            self.distribute_to_vertices_and_edges()
1101
1102            # Update extrema if necessary (for reporting)
1103            self.update_extrema()           
1104
1105            # Update boundary values
1106            self.update_boundary()
1107
1108            # Update time
1109            self.time += self.timestep
1110            self.yieldtime += self.timestep
1111            self.number_of_steps += 1
1112            if self._order_ == 1:
1113                self.number_of_first_order_steps += 1
1114
1115            # Yield results
1116            if finaltime is not None and self.time >= finaltime-epsilon:
1117
1118                if self.time > finaltime:
1119                    # FIXME (Ole, 30 April 2006): Do we need this check?
1120                    print 'WARNING (domain.py): time overshot finaltime. Contact Ole.Nielsen@ga.gov.au'
1121                    self.time = finaltime
1122
1123                # Yield final time and stop
1124                self.time = finaltime
1125                yield(self.time)
1126                break
1127
1128
1129            if self.yieldtime >= yieldstep:
1130                # Yield (intermediate) time and allow inspection of domain
1131
1132                if self.checkpoint is True:
1133                    self.store_checkpoint()
1134                    self.delete_old_checkpoints()
1135
1136                # Pass control on to outer loop for more specific actions
1137                yield(self.time)
1138
1139                # Reinitialise
1140                self.yieldtime = 0.0
1141                self.min_timestep = max_timestep
1142                self.max_timestep = min_timestep
1143                self.number_of_steps = 0
1144                self.number_of_first_order_steps = 0
1145
1146
1147    def evolve_to_end(self, finaltime = 1.0):
1148        """Iterate evolve all the way to the end
1149        """
1150
1151        for _ in self.evolve(yieldstep=None, finaltime=finaltime):
1152            pass
1153
1154
1155
1156    def update_boundary(self):
1157        """Go through list of boundary objects and update boundary values
1158        for all conserved quantities on boundary.
1159        It is assumed that the ordering of conserved quantities is
1160        consistent between the domain and the boundary object, i.e.
1161        the jth element of vector q must correspond to the jth conserved
1162        quantity in domain.
1163        """
1164
1165        #FIXME: Update only those that change (if that can be worked out)
1166        #FIXME: Boundary objects should not include ghost nodes.
1167        for i, ((vol_id, edge_id), B) in enumerate(self.boundary_objects):
1168            if B is None:
1169                print 'WARNING: Ignored boundary segment %d (None)'
1170            else:
1171                q = B.evaluate(vol_id, edge_id)
1172
1173                for j, name in enumerate(self.conserved_quantities):
1174                    Q = self.quantities[name]
1175                    Q.boundary_values[i] = q[j]
1176
1177
1178    def compute_fluxes(self):
1179        msg = 'Method compute_fluxes must be overridden by Domain subclass'
1180        raise msg
1181
1182
1183    def update_timestep(self, yieldstep, finaltime):
1184
1185        from anuga.config import min_timestep, max_timestep
1186
1187       
1188       
1189        # Protect against degenerate timesteps arising from isolated
1190        # triangles
1191        if self.protect_against_isolated_degenerate_timesteps is True and\
1192               self.max_speed > 10.0:
1193
1194            # Setup 10 bins for speed histogram
1195            from anuga.utilities.numerical_tools import histogram, create_bins
1196       
1197            bins = create_bins(self.max_speed, 10)
1198            hist = histogram(self.max_speed, bins)
1199
1200            # Look for characteristic signature
1201            if len(hist) > 1 and\
1202                hist[-1] > 0 and\
1203                hist[4] == hist[5] == hist[6] == hist[7] == hist[8] == 0:
1204                    # Danger of isolated degenerate triangles
1205                    # print self.timestepping_statistics(track_speeds=True) 
1206                   
1207                    # Find triangles in last bin
1208                    # FIXME - speed up using Numeric
1209                    d = 0
1210                    for i in range(self.number_of_full_triangles):
1211                        if self.max_speed[i] > bins[-1]:
1212                            msg = 'Time=%f: Ignoring isolated high speed triangle ' %self.time
1213                            msg += '#%d of %d with max speed=%f'\
1214                                  %(i, self.number_of_full_triangles, self.max_speed[i])
1215                            #print msg
1216                   
1217                            # print 'Found offending triangle', i, self.max_speed[i]
1218                            self.get_quantity('xmomentum').set_values(0.0, indices=[i])
1219                            self.get_quantity('ymomentum').set_values(0.0, indices=[i])
1220                            self.max_speed[i]=0.0
1221                            d += 1
1222                   
1223                    #print 'Adjusted %d triangles' %d       
1224                    #print self.timestepping_statistics(track_speeds=True)     
1225               
1226
1227                       
1228        # self.timestep is calculated from speed of characteristics
1229        # Apply CFL condition here
1230        timestep = min(self.CFL*self.timestep, max_timestep)
1231
1232        #Record maximal and minimal values of timestep for reporting
1233        self.max_timestep = max(timestep, self.max_timestep)
1234        self.min_timestep = min(timestep, self.min_timestep)
1235
1236           
1237 
1238        #Protect against degenerate time steps
1239        if timestep < min_timestep:
1240
1241            #Number of consecutive small steps taken b4 taking action
1242            self.smallsteps += 1
1243
1244            if self.smallsteps > self.max_smallsteps:
1245                self.smallsteps = 0 #Reset
1246
1247                if self._order_ == 1:
1248                    msg = 'WARNING: Too small timestep %.16f reached '\
1249                          %timestep
1250                    msg += 'even after %d steps of 1 order scheme'\
1251                           %self.max_smallsteps
1252                    print msg
1253                    timestep = min_timestep  #Try enforcing min_step
1254
1255
1256                    #print self.timestepping_statistics(track_speeds=True)
1257
1258
1259                    raise Exception, msg
1260                else:
1261                    #Try to overcome situation by switching to 1 order
1262                    self._order_ = 1
1263
1264        else:
1265            self.smallsteps = 0
1266            if self._order_ == 1 and self.default_order == 2:
1267                self._order_ = 2
1268
1269
1270        #Ensure that final time is not exceeded
1271        if finaltime is not None and self.time + timestep > finaltime :
1272            timestep = finaltime-self.time
1273
1274        #Ensure that model time is aligned with yieldsteps
1275        if self.yieldtime + timestep > yieldstep:
1276            timestep = yieldstep-self.yieldtime
1277
1278        self.timestep = timestep
1279
1280
1281
1282    def compute_forcing_terms(self):
1283        """If there are any forcing functions driving the system
1284        they should be defined in Domain subclass and appended to
1285        the list self.forcing_terms
1286        """
1287
1288        for f in self.forcing_terms:
1289            f(self)
1290
1291
1292
1293    def update_conserved_quantities(self):
1294        """Update vectors of conserved quantities using previously
1295        computed fluxes specified forcing functions.
1296        """
1297
1298        from Numeric import ones, sum, equal, Float
1299
1300        N = len(self) #number_of_triangles
1301        d = len(self.conserved_quantities)
1302
1303        timestep = self.timestep
1304
1305        #Compute forcing terms
1306        self.compute_forcing_terms()
1307
1308        #Update conserved_quantities
1309        for name in self.conserved_quantities:
1310            Q = self.quantities[name]
1311            Q.update(timestep)
1312
1313            #Clean up
1314            #Note that Q.explicit_update is reset by compute_fluxes
1315
1316            #MH090605 commented out the following since semi_implicit_update is now re-initialized
1317            #at the end of the _update function in quantity_ext.c (This is called by the
1318            #preceeding Q.update(timestep) statement above).
1319            #For run_profile.py with N=128, the time of update_conserved_quantities is cut from 14.00 secs
1320            #to 8.35 secs
1321
1322            #Q.semi_implicit_update[:] = 0.0
1323
1324    def update_ghosts(self):
1325        pass
1326
1327    def distribute_to_vertices_and_edges(self):
1328        """Extrapolate conserved quantities from centroid to
1329        vertices and edge-midpoints for each volume
1330
1331        Default implementation is straight first order,
1332        i.e. constant values throughout each element and
1333        no reference to non-conserved quantities.
1334        """
1335
1336        for name in self.conserved_quantities:
1337            Q = self.quantities[name]
1338            if self._order_ == 1:
1339                Q.extrapolate_first_order()
1340            elif self._order_ == 2:
1341                Q.extrapolate_second_order()
1342                Q.limit()
1343            else:
1344                raise 'Unknown order'
1345            Q.interpolate_from_vertices_to_edges()
1346
1347
1348    def centroid_norm(self, quantity, normfunc):
1349        """Calculate the norm of the centroid values
1350        of a specific quantity, using normfunc.
1351
1352        normfunc should take a list to a float.
1353
1354        common normfuncs are provided in the module utilities.norms
1355        """
1356        return normfunc(self.quantities[quantity].centroid_values)
1357
1358
1359
1360##############################################
1361#Initialise module
1362
1363#Optimisation with psyco
1364from anuga.config import use_psyco
1365if use_psyco:
1366    try:
1367        import psyco
1368    except:
1369        import os
1370        if os.name == 'posix' and os.uname()[4] == 'x86_64':
1371            pass
1372            #Psyco isn't supported on 64 bit systems, but it doesn't matter
1373        else:
1374            msg = 'WARNING: psyco (speedup) could not import'+\
1375                  ', you may want to consider installing it'
1376            print msg
1377    else:
1378        psyco.bind(Domain.update_boundary)
1379        #psyco.bind(Domain.update_timestep)     #Not worth it
1380        psyco.bind(Domain.update_conserved_quantities)
1381        psyco.bind(Domain.distribute_to_vertices_and_edges)
1382
1383
1384if __name__ == "__main__":
1385    pass
Note: See TracBrowser for help on using the repository browser.