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

Last change on this file since 4217 was 4217, checked in by ole, 18 years ago

More diagnostics

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