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

Last change on this file since 5571 was 5442, checked in by ole, 16 years ago

Retired h-limiter and beta_h as per ticket:194.
All unit tests and validation tests pass.

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