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

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

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

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