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

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

Refactored get_flow_through_cross_section and added a runtime version to
shallow_water_domain. It still needs some work and testing.

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