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

Last change on this file since 5306 was 5306, checked in by steve, 16 years ago

Working version of wave.py for Ole to work on.

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