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

Last change on this file since 5957 was 5957, checked in by ole, 15 years ago

Moved default_order into config.py

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