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

Last change on this file since 4701 was 4701, checked in by ole, 17 years ago

Work towards ticket:192

File size: 41.0 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
11from anuga.config import epsilon
12
13from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
14from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
15     import Boundary
16from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
17     import File_boundary
18from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
19     import Dirichlet_boundary
20from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
21     import Time_boundary
22from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
23     import Transmissive_boundary
24
25from anuga.abstract_2d_finite_volumes.pmesh2domain import pmesh_to_domain
26from anuga.abstract_2d_finite_volumes.region\
27     import Set_region as region_set_region
28
29import types
30
31class Domain(Mesh):
32
33
34    def __init__(self,
35                 source=None,
36                 triangles=None,
37                 boundary=None,
38                 conserved_quantities=None,
39                 other_quantities=None,
40                 tagged_elements=None,
41                 geo_reference=None,
42                 use_inscribed_circle=False,
43                 mesh_filename=None,
44                 use_cache=False,
45                 verbose=False,
46                 full_send_dict=None,
47                 ghost_recv_dict=None,
48                 processor=0,
49                 numproc=1,
50                 number_of_full_nodes=None,
51                 number_of_full_triangles=None):
52
53
54        """Instantiate generic computational Domain.
55
56        Input:
57          source:    Either a mesh filename or coordinates of mesh vertices.
58                     If it is a filename values specified for triangles will
59                     be overridden.
60          triangles: Mesh connectivity (see mesh.py for more information)
61          boundary:  See mesh.py for more information
62
63          conserved_quantities: List of quantity names entering the
64                                conservation equations
65          other_quantities:     List of other quantity names
66
67          tagged_elements:
68          ...
69
70
71        """
72
73        # Determine whether source is a mesh filename or coordinates
74        if type(source) == types.StringType:
75            mesh_filename = source
76        else:
77            coordinates = source
78
79
80        # In case a filename has been specified, extract content
81        if mesh_filename is not None:
82            coordinates, triangles, boundary, vertex_quantity_dict, \
83                         tagged_elements, geo_reference = \
84                         pmesh_to_domain(file_name=mesh_filename,
85                                         use_cache=use_cache,
86                                         verbose=verbose)
87
88
89        # Initialise underlying mesh structure
90        Mesh.__init__(self, coordinates, triangles,
91                      boundary=boundary,
92                      tagged_elements=tagged_elements,
93                      geo_reference=geo_reference,
94                      use_inscribed_circle=use_inscribed_circle,
95                      number_of_full_nodes=number_of_full_nodes,
96                      number_of_full_triangles=number_of_full_triangles,
97                      verbose=verbose)
98
99        if verbose: print 'Initialising Domain'
100        from Numeric import zeros, Float, Int, ones
101        from quantity import Quantity, Conserved_quantity
102
103        # List of quantity names entering
104        # the conservation equations
105        if conserved_quantities is None:
106            self.conserved_quantities = []
107        else:
108            self.conserved_quantities = conserved_quantities
109
110        # List of other quantity names
111        if other_quantities is None:
112            self.other_quantities = []
113        else:
114            self.other_quantities = other_quantities
115
116
117        #Build dictionary of Quantity instances keyed by quantity names
118        self.quantities = {}
119
120        #FIXME: remove later - maybe OK, though....
121        for name in self.conserved_quantities:
122            self.quantities[name] = Conserved_quantity(self)
123        for name in self.other_quantities:
124            self.quantities[name] = Quantity(self)
125
126        #Create an empty list for explicit forcing terms
127        self.forcing_terms = []
128
129        #Setup the ghost cell communication
130        if full_send_dict is None:
131            self.full_send_dict = {}
132        else:
133            self.full_send_dict  = full_send_dict
134
135        # List of other quantity names
136        if ghost_recv_dict is None:
137            self.ghost_recv_dict = {}
138        else:
139            self.ghost_recv_dict = ghost_recv_dict
140
141        self.processor = processor
142        self.numproc = numproc
143
144
145        # Setup Communication Buffers
146        if verbose: print 'Domain: Set up communication buffers (parallel)'
147        self.nsys = len(self.conserved_quantities)
148        for key in self.full_send_dict:
149            buffer_shape = self.full_send_dict[key][0].shape[0]
150            self.full_send_dict[key].append(zeros( (buffer_shape,self.nsys) ,Float))
151
152
153        for key in self.ghost_recv_dict:
154            buffer_shape = self.ghost_recv_dict[key][0].shape[0]
155            self.ghost_recv_dict[key].append(zeros( (buffer_shape,self.nsys) ,Float))
156
157
158        # Setup cell full flag
159        # =1 for full
160        # =0 for ghost
161        N = len(self) #number_of_elements
162        self.tri_full_flag = ones(N, Int)
163        for i in self.ghost_recv_dict.keys():
164            for id in self.ghost_recv_dict[i][0]:
165                self.tri_full_flag[id] = 0
166
167        # Test the assumption that all full triangles are store before
168        # the ghost triangles.
169        assert allclose(self.tri_full_flag[:self.number_of_full_nodes],1)
170                       
171
172        # Defaults
173        from anuga.config import max_smallsteps, beta_w, beta_h, epsilon
174        from anuga.config import CFL
175        from anuga.config import protect_against_isolated_degenerate_timesteps
176        self.beta_w = beta_w
177        self.beta_h = beta_h
178        self.epsilon = epsilon
179        self.protect_against_isolated_degenerate_timesteps = protect_against_isolated_degenerate_timesteps
180       
181
182        # FIXME: Maybe have separate orders for h-limiter and w-limiter?
183        # Or maybe get rid of order altogether and use beta_w and beta_h
184        self.set_default_order(1)
185        #self.default_order = 1
186        #self._order_ = self.default_order
187
188        self.smallsteps = 0
189        self.max_smallsteps = max_smallsteps
190        self.number_of_steps = 0
191        self.number_of_first_order_steps = 0
192        self.CFL = CFL
193
194        self.boundary_map = None  # Will be populated by set_boundary       
195       
196
197        # Model time
198        self.time = 0.0
199        self.finaltime = None
200        self.min_timestep = self.max_timestep = 0.0
201        self.starttime = 0 #Physical starttime if any (0 is 1 Jan 1970 00:00:00)
202
203        # Monitoring
204        self.quantities_to_be_monitored = None
205        self.monitor_polygon = None
206        self.monitor_time_interval = None               
207
208
209        # Checkpointing and storage
210        from anuga.config import default_datadir
211        self.datadir = default_datadir
212        self.simulation_name = 'domain'
213        self.checkpoint = False
214
215        # MH310505 To avoid calculating the flux across each edge twice, keep an integer (boolean) array,
216        # to be used during the flux calculation
217        N = len(self) #number_of_triangles
218        self.already_computed_flux = zeros((N, 3), Int)
219
220        # Storage for maximal speeds computed for each triangle by compute_fluxes
221        # This is used for diagnostics only
222        self.max_speed = zeros(N, Float)
223
224        if mesh_filename is not None:
225            # If the mesh file passed any quantity values
226            # , initialise with these values.
227            if verbose: print 'Domain: Initialising quantity values'
228            self.set_quantity_vertices_dict(vertex_quantity_dict)
229
230
231        if verbose: print 'Domain: Done'
232
233
234
235
236    def set_default_order(self, n):
237        """Set default (spatial) order to either 1 or 2
238        """
239
240        msg = 'Default order must be either 1 or 2. I got %s' %n
241        assert n in [1,2], msg
242
243        self.default_order = n
244        self._order_ = self.default_order
245
246
247    #Public interface to Domain
248    def get_conserved_quantities(self, vol_id, vertex=None, edge=None):
249        """Get conserved quantities at volume vol_id
250
251        If vertex is specified use it as index for vertex values
252        If edge is specified use it as index for edge values
253        If neither are specified use centroid values
254        If both are specified an exeception is raised
255
256        Return value: Vector of length == number_of_conserved quantities
257
258        """
259
260        from Numeric import zeros, Float
261
262        if not (vertex is None or edge is None):
263            msg = 'Values for both vertex and edge was specified.'
264            msg += 'Only one (or none) is allowed.'
265            raise msg
266
267        q = zeros( len(self.conserved_quantities), Float)
268
269        for i, name in enumerate(self.conserved_quantities):
270            Q = self.quantities[name]
271            if vertex is not None:
272                q[i] = Q.vertex_values[vol_id, vertex]
273            elif edge is not None:
274                q[i] = Q.edge_values[vol_id, edge]
275            else:
276                q[i] = Q.centroid_values[vol_id]
277
278        return q
279
280    def set_time(self, time=0.0):
281        """Set the model time (seconds)"""
282
283        self.time = time
284
285    def get_time(self):
286        """Get the model time (seconds)"""
287
288        return self.time
289
290    def set_quantity_vertices_dict(self, quantity_dict):
291        """Set values for named quantities.
292        The index is the quantity
293
294        name: Name of quantity
295        X: Compatible list, Numeric array, const or function (see below)
296
297        The values will be stored in elements following their
298        internal ordering.
299
300        """
301       
302        # FIXME: Could we name this a bit more intuitively
303        # E.g. set_quantities_from_dictionary
304        for key in quantity_dict.keys():
305            self.set_quantity(key, quantity_dict[key], location='vertices')
306
307
308    def set_quantity(self, name, *args, **kwargs):
309        """Set values for named quantity
310
311
312        One keyword argument is documented here:
313        expression = None, # Arbitrary expression
314
315        expression:
316          Arbitrary expression involving quantity names
317
318        See Quantity.set_values for further documentation.
319        """
320
321        #FIXME (Ole): Allow new quantities here
322        #from quantity import Quantity, Conserved_quantity
323        #Create appropriate quantity object
324        ##if name in self.conserved_quantities:
325        ##    self.quantities[name] = Conserved_quantity(self)
326        ##else:
327        ##    self.quantities[name] = Quantity(self)
328
329
330        #Do the expression stuff
331        if kwargs.has_key('expression'):
332            expression = kwargs['expression']
333            del kwargs['expression']
334
335            Q = self.create_quantity_from_expression(expression)
336            kwargs['quantity'] = Q
337
338        #Assign values
339        self.quantities[name].set_values(*args, **kwargs)
340
341
342    def get_quantity_names(self):
343        """Get a list of all the quantity names that this domain is aware of.
344        Any value in the result should be a valid input to get_quantity.
345        """
346        return self.quantities.keys()
347
348    def get_quantity(self, name, location='vertices', indices = None):
349        """Get quantity object.
350
351        name: Name of quantity
352
353        See methods inside the quantity object for more options
354
355        FIXME: clean input args
356        """
357
358        return self.quantities[name] #.get_values( location, indices = indices)
359
360
361    def get_quantity_object(self, name):
362        """Get object for named quantity
363
364        name: Name of quantity
365
366        FIXME: Obsolete
367        """
368
369        print 'get_quantity_object has been deprecated. Please use get_quantity'
370        return self.quantities[name]
371
372
373    def create_quantity_from_expression(self, expression):
374        """Create new quantity from other quantities using arbitrary expression
375
376        Combine existing quantities in domain using expression and return
377        result as a new quantity.
378
379        Note, the new quantity could e.g. be used in set_quantity
380
381        Valid expressions are limited to operators defined in class Quantity
382
383        Examples creating derived quantities:
384
385            Depth = domain.create_quantity_from_expression('stage-elevation')
386
387            exp = '(xmomentum*xmomentum + ymomentum*ymomentum)**0.5')       
388            Absolute_momentum = domain.create_quantity_from_expression(exp)
389
390        """
391
392        from anuga.abstract_2d_finite_volumes.util import\
393             apply_expression_to_dictionary
394       
395        return apply_expression_to_dictionary(expression, self.quantities)
396
397
398
399    #def modify_boundary(self, boundary_map):
400    #    """Modify existing boundary by elements in boundary map#
401    #
402    #    Input:#
403    #
404    #    boundary_map: Dictionary mapping tags to boundary objects
405    #
406    #    See set_boundary for more details on how this works
407    #
408    #      OBSOLETE
409    #    """
410    #
411    #    for key in boundary_map.keys():
412    #        self.boundary_map[key] = boundary_map[key]
413    #
414    #    self.set_boundary(self.boundary_map)
415       
416       
417
418    def set_boundary(self, boundary_map):
419        """Associate boundary objects with tagged boundary segments.
420
421        Input boundary_map is a dictionary of boundary objects keyed
422        by symbolic tags to matched against tags in the internal dictionary
423        self.boundary.
424
425        As result one pointer to a boundary object is stored for each vertex
426        in the list self.boundary_objects.
427        More entries may point to the same boundary object
428
429        Schematically the mapping is from two dictionaries to one list
430        where the index is used as pointer to the boundary_values arrays
431        within each quantity.
432
433        self.boundary:          (vol_id, edge_id): tag
434        boundary_map (input):   tag: boundary_object
435        ----------------------------------------------
436        self.boundary_objects:  ((vol_id, edge_id), boundary_object)
437
438
439        Pre-condition:
440          self.boundary has been built.
441
442        Post-condition:
443          self.boundary_objects is built
444
445        If a tag from the domain doesn't appear in the input dictionary an
446        exception is raised.
447        However, if a tag is not used to the domain, no error is thrown.
448        FIXME: This would lead to implementation of a
449        default boundary condition
450
451        Note: If a segment is listed in the boundary dictionary and if it is
452        not None, it *will* become a boundary -
453        even if there is a neighbouring triangle.
454        This would be the case for internal boundaries
455
456        Boundary objects that are None will be skipped.
457
458        If a boundary_map has already been set
459        (i.e. set_boundary has been called before), the old boundary map
460        will be updated with new values. The new map need not define all
461        boundary tags, and can thus change only those that are needed.
462
463        FIXME: If set_boundary is called multiple times and if Boundary
464        object is changed into None, the neighbour structure will not be
465        restored!!!
466
467
468        """
469
470        if self.boundary_map is None:
471            # This the first call to set_boundary. Store
472            # map for later updates and for use with boundary_stats.
473            self.boundary_map = boundary_map
474        else:   
475            # This is a modification of an already existing map
476            # Update map an proceed normally
477
478            for key in boundary_map.keys():
479                self.boundary_map[key] = boundary_map[key]
480               
481       
482        #FIXME: Try to remove the sorting and fix test_mesh.py
483        x = self.boundary.keys()
484        x.sort()
485
486        #Loop through edges that lie on the boundary and associate them with
487        #callable boundary objects depending on their tags
488        self.boundary_objects = []       
489        for k, (vol_id, edge_id) in enumerate(x):
490            tag = self.boundary[ (vol_id, edge_id) ]
491
492            if self.boundary_map.has_key(tag):
493                B = self.boundary_map[tag]  #Get callable boundary object
494
495                if B is not None:
496                    self.boundary_objects.append( ((vol_id, edge_id), B) )
497                    self.neighbours[vol_id, edge_id] = \
498                                            -len(self.boundary_objects)
499                else:
500                    pass
501                    #FIXME: Check and perhaps fix neighbour structure
502
503
504            else:
505                msg = 'ERROR (domain.py): Tag "%s" has not been ' %tag
506                msg += 'bound to a boundary object.\n'
507                msg += 'All boundary tags defined in domain must appear '
508                msg += 'in set_boundary.\n'
509                msg += 'The tags are: %s' %self.get_boundary_tags()
510                raise msg
511
512
513    def set_region(self, *args, **kwargs):
514        """
515        This method is used to set quantities based on a regional tag.
516       
517        It is most often called with the following parameters;
518        (self, tag, quantity, X, location='vertices')
519        tag: the name of the regional tag used to specify the region
520        quantity: Name of quantity to change
521        X: const or function - how the quantity is changed
522        location: Where values are to be stored.
523            Permissible options are: vertices, centroid and unique vertices
524
525        A callable region class or a list of callable region classes
526        can also be passed into this function.
527        """
528        #print "*args", args
529        #print "**kwargs", kwargs
530        if len(args) == 1:
531            self._set_region(*args, **kwargs)
532        else:
533            #Assume it is arguments for the region.set_region function
534            func = region_set_region(*args, **kwargs)
535            self._set_region(func)
536           
537       
538    def _set_region(self, functions):
539        # The order of functions in the list is used.
540        if type(functions) not in [types.ListType,types.TupleType]:
541            functions = [functions]
542        for function in functions:
543            for tag in self.tagged_elements.keys():
544                function(tag, self.tagged_elements[tag], self)
545
546
547
548
549    def set_quantities_to_be_monitored(self, q,
550                                       polygon=None,
551                                       time_interval=None):
552        """Specify which quantities will be monitored for extrema.
553
554        q must be either:
555          - the name of a quantity
556          - a list of quantity names
557          - None
558
559        In the two first cases, the named quantities will be monitored at
560        each internal timestep
561       
562        If q is None, monitoring will be switched off altogether.
563
564        polygon (if specified) will restrict monitoring to triangles inside polygon.
565        If omitted all triangles will be included.
566
567        time_interval (if specified) will restrict monitoring to time steps in
568        that interval. If omitted all timesteps will be included.
569
570        FIXME: Derived quantities such as 'stage-elevation' will appear later
571        """
572
573        # FIXME (Ole): This is under construction. See ticket:192
574
575        if q is None:
576            self.quantities_to_be_monitored = None
577            self.monitor_polygon = None
578            self.monitor_time_interval = None                   
579            return
580
581        if isinstance(q, basestring):
582            q = [q] # Turn argument into a list
583
584        # Check correcness
585        for quantity_name in q:
586            msg = 'Quantity %s is not a valid conserved quantity'\
587                  %quantity_name
588           
589            assert quantity_name in self.conserved_quantities, msg
590
591        if polygon is not None:
592            # FIXME Check input
593            pass
594
595        if time_interval is not None:
596            # FIXME Check input
597            pass       
598
599
600        self.quantities_to_be_monitored = q
601        self.monitor_polygon = polygon
602        self.monitor_time_interval = time_interval       
603       
604
605
606    #MISC
607    def check_integrity(self):
608        Mesh.check_integrity(self)
609
610        for quantity in self.conserved_quantities:
611            msg = 'Conserved quantities must be a subset of all quantities'
612            assert quantity in self.quantities, msg
613
614        ##assert hasattr(self, 'boundary_objects')
615
616    def write_time(self, track_speeds=False):
617        print self.timestepping_statistics(track_speeds)
618
619
620    def timestepping_statistics(self, track_speeds=False):
621        """Return string with time stepping statistics for printing or logging
622
623        Optional boolean keyword track_speeds decides whether to report location of
624        smallest timestep as well as a histogram and percentile report.
625        """
626
627        from anuga.utilities.numerical_tools import histogram, create_bins
628       
629
630        msg = ''
631        if self.min_timestep == self.max_timestep:
632            msg += 'Time = %.4f, delta t = %.8f, steps=%d (%d)'\
633                   %(self.time, self.min_timestep, self.number_of_steps,
634                     self.number_of_first_order_steps)
635        elif self.min_timestep > self.max_timestep:
636            msg += 'Time = %.4f, steps=%d (%d)'\
637                   %(self.time, self.number_of_steps,
638                     self.number_of_first_order_steps)
639        else:
640            msg += 'Time = %.4f, delta t in [%.8f, %.8f], steps=%d (%d)'\
641                   %(self.time, self.min_timestep,
642                     self.max_timestep, self.number_of_steps,
643                     self.number_of_first_order_steps)
644
645        if track_speeds is True:
646            msg += '\n'
647
648
649            #Setup 10 bins for speed histogram
650            bins = create_bins(self.max_speed, 10)
651            hist = histogram(self.max_speed, bins)
652
653            msg += '------------------------------------------------\n'
654            msg += '  Speeds in [%f, %f]\n' %(min(self.max_speed), max(self.max_speed))
655            msg += '  Histogram:\n'
656
657            hi = bins[0]
658            for i, count in enumerate(hist):
659                lo = hi
660                if i+1 < len(bins):
661                    #Open upper interval               
662                    hi = bins[i+1]
663                    msg += '    [%f, %f[: %d\n' %(lo, hi, count)               
664                else:
665                    #Closed upper interval
666                    hi = max(self.max_speed)
667                    msg += '    [%f, %f]: %d\n' %(lo, hi, count)
668
669
670            N = len(self.max_speed)
671            if N > 10:
672                msg += '  Percentiles (10%):\n'
673                speed = self.max_speed.tolist()
674                speed.sort()
675
676                k = 0
677                lower = min(speed)
678                for i, a in enumerate(speed):       
679                    if i % (N/10) == 0 and i != 0: #For every 10% of the sorted speeds
680                        msg += '    %d speeds in [%f, %f]\n' %(i-k, lower, a)
681                        lower = a
682                        k = i
683                       
684                msg += '    %d speeds in [%f, %f]\n'\
685                       %(N-k, lower, max(speed))                   
686               
687                     
688           
689
690               
691           
692            # Find index of largest computed flux speed
693            k = argmax(self.max_speed)
694
695            x, y = self.get_centroid_coordinates()[k]
696            radius = self.get_radii()[k]
697            area = self.get_areas()[k]     
698            max_speed = self.max_speed[k]           
699
700            msg += '  Triangle #%d with centroid (%.4f, %.4f), ' %(k, x, y)
701            msg += 'area = %.4f and radius = %.4f ' %(area, radius)
702            msg += 'had the largest computed speed: %.6f m/s ' %(max_speed)
703            if max_speed > 0.0:
704                msg += '(timestep=%.6f)\n' %(radius/max_speed)
705            else:
706                msg += '(timestep=%.6f)\n' %(0)     
707           
708            # Report all quantity values at vertices
709            msg += '    Quantity \t vertex values\t\t\t\t\t centroid values\n'
710            for name in self.quantities:
711                q = self.quantities[name]
712
713                V = q.get_values(location='vertices', indices=[k])[0]
714                C = q.get_values(location='centroids', indices=[k])                 
715               
716                s = '    %s:\t %.4f,\t %.4f,\t %.4f,\t %.4f\n'\
717                    %(name, V[0], V[1], V[2], C[0])                 
718
719                msg += s
720
721        return msg
722
723
724    def write_boundary_statistics(self, quantities = None, tags = None):
725        print self.boundary_statistics(quantities, tags)
726
727    def boundary_statistics(self, quantities = None, tags = None):
728        """Output statistics about boundary forcing at each timestep
729
730
731        Input:
732          quantities: either None, a string or a list of strings naming the quantities to be reported
733          tags:       either None, a string or a list of strings naming the tags to be reported
734
735
736        Example output:
737        Tag 'wall':
738            stage in [2, 5.5]
739            xmomentum in []
740            ymomentum in []
741        Tag 'ocean'
742
743
744        If quantities are specified only report on those. Otherwise take all conserved quantities.
745        If tags are specified only report on those, otherwise take all tags.
746
747        """
748
749        #Input checks
750        import types, string
751
752        if quantities is None:
753            quantities = self.conserved_quantities
754        elif type(quantities) == types.StringType:
755            quantities = [quantities] #Turn it into a list
756
757        msg = 'Keyword argument quantities must be either None, '
758        msg += 'string or list. I got %s' %str(quantities)
759        assert type(quantities) == types.ListType, msg
760
761
762        if tags is None:
763            tags = self.get_boundary_tags()
764        elif type(tags) == types.StringType:
765            tags = [tags] #Turn it into a list
766
767        msg = 'Keyword argument tags must be either None, '
768        msg += 'string or list. I got %s' %str(tags)
769        assert type(tags) == types.ListType, msg
770
771        #Determine width of longest quantity name (for cosmetic purposes)
772        maxwidth = 0
773        for name in quantities:
774            w = len(name)
775            if w > maxwidth:
776                maxwidth = w
777
778        #Output stats
779        msg = 'Boundary values at time %.4f:\n' %self.time
780        for tag in tags:
781            msg += '    %s:\n' %tag
782
783            for name in quantities:
784                q = self.quantities[name]
785
786                #Find range of boundary values for tag and q
787                maxval = minval = None
788                for i, ((vol_id, edge_id), B) in\
789                        enumerate(self.boundary_objects):
790                    if self.boundary[(vol_id, edge_id)] == tag:
791                        v = q.boundary_values[i]
792                        if minval is None or v < minval: minval = v
793                        if maxval is None or v > maxval: maxval = v
794
795                if minval is None or maxval is None:
796                    msg += '        Sorry no information available about' +\
797                           ' tag %s and quantity %s\n' %(tag, name)
798                else:
799                    msg += '        %s in [%12.8f, %12.8f]\n'\
800                           %(string.ljust(name, maxwidth), minval, maxval)
801
802
803        return msg
804
805
806
807    def quantity_statistics(self):
808        """Return string with statistics about quantities for printing or logging
809
810        Quantities reported are specified through method
811
812           set_quantities_to_be_monitored
813           
814        """
815
816        pass
817
818
819
820
821    def get_name(self):
822        return self.simulation_name
823
824    def set_name(self, name):
825        """Assign a name to this simulation.
826        This will be used to identify the output sww file.
827
828        """
829        if name.endswith('.sww'):
830            name = name[:-4]
831           
832        self.simulation_name = name
833
834    def get_datadir(self):
835        return self.datadir
836
837    def set_datadir(self, name):
838        self.datadir = name
839
840
841
842    #def set_defaults(self):
843    #    """Set default values for uninitialised quantities.
844    #    Should be overridden or specialised by specific modules
845    #    """#
846    #
847    #    for name in self.conserved_quantities + self.other_quantities:
848    #        self.set_quantity(name, 0.0)
849
850
851    ###########################
852    #Main components of evolve
853
854    def evolve(self,
855               yieldstep = None,
856               finaltime = None,
857               duration = None,
858               skip_initial_step = False):
859        """Evolve model through time starting from self.starttime.
860
861
862        yieldstep: Interval between yields where results are stored,
863                   statistics written and domain inspected or
864                   possibly modified. If omitted the internal predefined
865                   max timestep is used.
866                   Internally, smaller timesteps may be taken.
867
868        duration: Duration of simulation
869
870        finaltime: Time where simulation should end. This is currently
871        relative time.  So it's the same as duration.
872
873        If both duration and finaltime are given an exception is thrown.
874
875
876        skip_initial_step: Boolean flag that decides whether the first
877        yield step is skipped or not. This is useful for example to avoid
878        duplicate steps when multiple evolve processes are dove tailed.
879
880
881        Evolve is implemented as a generator and is to be called as such, e.g.
882
883        for t in domain.evolve(yieldstep, finaltime):
884            <Do something with domain and t>
885
886
887        All times are given in seconds
888
889        """
890
891        from anuga.config import min_timestep, max_timestep, epsilon
892
893        #FIXME: Maybe lump into a larger check prior to evolving
894        msg = 'Boundary tags must be bound to boundary objects before '
895        msg += 'evolving system, '
896        msg += 'e.g. using the method set_boundary.\n'
897        msg += 'This system has the boundary tags %s '\
898               %self.get_boundary_tags()
899        assert hasattr(self, 'boundary_objects'), msg
900
901
902        if yieldstep is None:
903            yieldstep = max_timestep
904        else:
905            yieldstep = float(yieldstep)
906
907        self._order_ = self.default_order
908
909
910        if finaltime is not None and duration is not None:
911            #print 'F', finaltime, duration
912            msg = 'Only one of finaltime and duration may be specified'
913            raise msg
914        else:
915            if finaltime is not None:
916                self.finaltime = float(finaltime)
917            if duration is not None:
918                self.finaltime = self.starttime + float(duration)
919
920
921
922
923        self.yieldtime = 0.0 #Time between 'yields'
924
925        #Initialise interval of timestep sizes (for reporting only)
926        self.min_timestep = max_timestep
927        self.max_timestep = min_timestep
928        self.number_of_steps = 0
929        self.number_of_first_order_steps = 0
930
931        #update ghosts
932        self.update_ghosts()
933
934        #Initial update of vertex and edge values
935        self.distribute_to_vertices_and_edges()
936
937        #Initial update boundary values
938        self.update_boundary()
939
940        #Or maybe restore from latest checkpoint
941        if self.checkpoint is True:
942            self.goto_latest_checkpoint()
943
944        if skip_initial_step is False:
945            yield(self.time)  #Yield initial values
946
947        while True:
948            #Compute fluxes across each element edge
949            self.compute_fluxes()
950
951            #Update timestep to fit yieldstep and finaltime
952            self.update_timestep(yieldstep, finaltime)
953
954            #Update conserved quantities
955            self.update_conserved_quantities()
956
957            #update ghosts
958            self.update_ghosts()
959
960            #Update vertex and edge values
961            self.distribute_to_vertices_and_edges()
962
963            #Update boundary values
964            self.update_boundary()
965
966            #Update time
967            self.time += self.timestep
968            self.yieldtime += self.timestep
969            self.number_of_steps += 1
970            if self._order_ == 1:
971                self.number_of_first_order_steps += 1
972
973            #Yield results
974            if finaltime is not None and self.time >= finaltime-epsilon:
975
976                if self.time > finaltime:
977                    #FIXME (Ole, 30 April 2006): Do we need this check?
978                    print 'WARNING (domain.py): time overshot finaltime. Contact Ole.Nielsen@ga.gov.au'
979                    self.time = finaltime
980
981                # Yield final time and stop
982                self.time = finaltime
983                yield(self.time)
984                break
985
986
987            if self.yieldtime >= yieldstep:
988                # Yield (intermediate) time and allow inspection of domain
989
990                if self.checkpoint is True:
991                    self.store_checkpoint()
992                    self.delete_old_checkpoints()
993
994                # Pass control on to outer loop for more specific actions
995                yield(self.time)
996
997                # Reinitialise
998                self.yieldtime = 0.0
999                self.min_timestep = max_timestep
1000                self.max_timestep = min_timestep
1001                self.number_of_steps = 0
1002                self.number_of_first_order_steps = 0
1003
1004
1005    def evolve_to_end(self, finaltime = 1.0):
1006        """Iterate evolve all the way to the end
1007        """
1008
1009        for _ in self.evolve(yieldstep=None, finaltime=finaltime):
1010            pass
1011
1012
1013
1014    def update_boundary(self):
1015        """Go through list of boundary objects and update boundary values
1016        for all conserved quantities on boundary.
1017        It is assumed that the ordering of conserved quantities is
1018        consistent between the domain and the boundary object, i.e.
1019        the jth element of vector q must correspond to the jth conserved
1020        quantity in domain.
1021        """
1022
1023        #FIXME: Update only those that change (if that can be worked out)
1024        #FIXME: Boundary objects should not include ghost nodes.
1025        for i, ((vol_id, edge_id), B) in enumerate(self.boundary_objects):
1026            if B is None:
1027                print 'WARNING: Ignored boundary segment %d (None)'
1028            else:
1029                q = B.evaluate(vol_id, edge_id)
1030
1031                for j, name in enumerate(self.conserved_quantities):
1032                    Q = self.quantities[name]
1033                    Q.boundary_values[i] = q[j]
1034
1035
1036    def compute_fluxes(self):
1037        msg = 'Method compute_fluxes must be overridden by Domain subclass'
1038        raise msg
1039
1040
1041    def update_timestep(self, yieldstep, finaltime):
1042
1043        from anuga.config import min_timestep, max_timestep
1044
1045       
1046       
1047        # Protect against degenerate timesteps arising from isolated
1048        # triangles
1049        if self.protect_against_isolated_degenerate_timesteps is True and\
1050               self.max_speed > 10.0:
1051
1052            # Setup 10 bins for speed histogram
1053            from anuga.utilities.numerical_tools import histogram, create_bins
1054       
1055            bins = create_bins(self.max_speed, 10)
1056            hist = histogram(self.max_speed, bins)
1057
1058            # Look for characteristic signature
1059            if len(hist) > 1 and\
1060                hist[-1] > 0 and\
1061                hist[4] == hist[5] == hist[6] == hist[7] == hist[8] == 0:
1062                    # Danger of isolated degenerate triangles
1063                    # print self.timestepping_statistics(track_speeds=True) 
1064                   
1065                    # Find triangles in last bin
1066                    # FIXME - speed up using Numeric
1067                    d = 0
1068                    for i in range(self.number_of_full_triangles):
1069                        if self.max_speed[i] > bins[-1]:
1070                            msg = 'Time=%f: Ignoring isolated high speed triangle ' %self.time
1071                            msg += '#%d of %d with max speed=%f'\
1072                                  %(i, self.number_of_full_triangles, self.max_speed[i])
1073                            #print msg
1074                   
1075                            # print 'Found offending triangle', i, self.max_speed[i]
1076                            self.get_quantity('xmomentum').set_values(0.0, indices=[i])
1077                            self.get_quantity('ymomentum').set_values(0.0, indices=[i])
1078                            self.max_speed[i]=0.0
1079                            d += 1
1080                   
1081                    #print 'Adjusted %d triangles' %d       
1082                    #print self.timestepping_statistics(track_speeds=True)     
1083               
1084
1085                       
1086        # self.timestep is calculated from speed of characteristics
1087        # Apply CFL condition here
1088        timestep = min(self.CFL*self.timestep, max_timestep)
1089
1090        #Record maximal and minimal values of timestep for reporting
1091        self.max_timestep = max(timestep, self.max_timestep)
1092        self.min_timestep = min(timestep, self.min_timestep)
1093
1094           
1095 
1096        #Protect against degenerate time steps
1097        if timestep < min_timestep:
1098
1099            #Number of consecutive small steps taken b4 taking action
1100            self.smallsteps += 1
1101
1102            if self.smallsteps > self.max_smallsteps:
1103                self.smallsteps = 0 #Reset
1104
1105                if self._order_ == 1:
1106                    msg = 'WARNING: Too small timestep %.16f reached '\
1107                          %timestep
1108                    msg += 'even after %d steps of 1 order scheme'\
1109                           %self.max_smallsteps
1110                    print msg
1111                    timestep = min_timestep  #Try enforcing min_step
1112
1113
1114                    #print self.timestepping_statistics(track_speeds=True)
1115
1116
1117                    raise Exception, msg
1118                else:
1119                    #Try to overcome situation by switching to 1 order
1120                    self._order_ = 1
1121
1122        else:
1123            self.smallsteps = 0
1124            if self._order_ == 1 and self.default_order == 2:
1125                self._order_ = 2
1126
1127
1128        #Ensure that final time is not exceeded
1129        if finaltime is not None and self.time + timestep > finaltime :
1130            timestep = finaltime-self.time
1131
1132        #Ensure that model time is aligned with yieldsteps
1133        if self.yieldtime + timestep > yieldstep:
1134            timestep = yieldstep-self.yieldtime
1135
1136        self.timestep = timestep
1137
1138
1139
1140    def compute_forcing_terms(self):
1141        """If there are any forcing functions driving the system
1142        they should be defined in Domain subclass and appended to
1143        the list self.forcing_terms
1144        """
1145
1146        for f in self.forcing_terms:
1147            f(self)
1148
1149
1150
1151    def update_conserved_quantities(self):
1152        """Update vectors of conserved quantities using previously
1153        computed fluxes specified forcing functions.
1154        """
1155
1156        from Numeric import ones, sum, equal, Float
1157
1158        N = len(self) #number_of_triangles
1159        d = len(self.conserved_quantities)
1160
1161        timestep = self.timestep
1162
1163        #Compute forcing terms
1164        self.compute_forcing_terms()
1165
1166        #Update conserved_quantities
1167        for name in self.conserved_quantities:
1168            Q = self.quantities[name]
1169            Q.update(timestep)
1170
1171            #Clean up
1172            #Note that Q.explicit_update is reset by compute_fluxes
1173
1174            #MH090605 commented out the following since semi_implicit_update is now re-initialized
1175            #at the end of the _update function in quantity_ext.c (This is called by the
1176            #preceeding Q.update(timestep) statement above).
1177            #For run_profile.py with N=128, the time of update_conserved_quantities is cut from 14.00 secs
1178            #to 8.35 secs
1179
1180            #Q.semi_implicit_update[:] = 0.0
1181
1182    def update_ghosts(self):
1183        pass
1184
1185    def distribute_to_vertices_and_edges(self):
1186        """Extrapolate conserved quantities from centroid to
1187        vertices and edge-midpoints for each volume
1188
1189        Default implementation is straight first order,
1190        i.e. constant values throughout each element and
1191        no reference to non-conserved quantities.
1192        """
1193
1194        for name in self.conserved_quantities:
1195            Q = self.quantities[name]
1196            if self._order_ == 1:
1197                Q.extrapolate_first_order()
1198            elif self._order_ == 2:
1199                Q.extrapolate_second_order()
1200                Q.limit()
1201            else:
1202                raise 'Unknown order'
1203            Q.interpolate_from_vertices_to_edges()
1204
1205
1206    def centroid_norm(self, quantity, normfunc):
1207        """Calculate the norm of the centroid values
1208        of a specific quantity, using normfunc.
1209
1210        normfunc should take a list to a float.
1211
1212        common normfuncs are provided in the module utilities.norms
1213        """
1214        return normfunc(self.quantities[quantity].centroid_values)
1215
1216
1217
1218##############################################
1219#Initialise module
1220
1221#Optimisation with psyco
1222from anuga.config import use_psyco
1223if use_psyco:
1224    try:
1225        import psyco
1226    except:
1227        import os
1228        if os.name == 'posix' and os.uname()[4] == 'x86_64':
1229            pass
1230            #Psyco isn't supported on 64 bit systems, but it doesn't matter
1231        else:
1232            msg = 'WARNING: psyco (speedup) could not import'+\
1233                  ', you may want to consider installing it'
1234            print msg
1235    else:
1236        psyco.bind(Domain.update_boundary)
1237        #psyco.bind(Domain.update_timestep)     #Not worth it
1238        psyco.bind(Domain.update_conserved_quantities)
1239        psyco.bind(Domain.distribute_to_vertices_and_edges)
1240
1241
1242if __name__ == "__main__":
1243    pass
Note: See TracBrowser for help on using the repository browser.