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

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

Moved shallow water out from the old pyvolution directory.
All tests pass, most examples and okushiri works again.

File size: 32.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, Christopher Zoppou
7   Geoscience Australia
8"""
9
10from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
11from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
12     import Boundary
13from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
14     import File_boundary
15from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
16     import Dirichlet_boundary
17from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
18     import Time_boundary
19from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
20     import Transmissive_boundary
21
22from anuga.abstract_2d_finite_volumes.pmesh2domain import pmesh_to_domain
23from anuga.abstract_2d_finite_volumes.region\
24     import Set_region as region_set_region
25
26import types
27
28class Domain(Mesh):
29
30
31    def __init__(self,
32                 source=None,
33                 triangles=None,
34                 boundary=None,
35                 conserved_quantities=None,
36                 other_quantities=None,
37                 tagged_elements=None,
38                 geo_reference=None,
39                 use_inscribed_circle=False,
40                 mesh_filename=None,
41                 use_cache=False,
42                 verbose=False,
43                 full_send_dict=None,
44                 ghost_recv_dict=None,
45                 processor=0,
46                 numproc=1):
47
48
49        """Instantiate generic computational Domain.
50
51        Input:
52          source:    Either a mesh filename or coordinates of mesh vertices.
53                     If it is a filename values specified for triangles will
54                     be overridden.
55          triangles: Mesh connectivity (see mesh.py for more information)
56          boundary:  See mesh.py for more information
57
58          conserved_quantities: List of quantity names entering the
59                                conservation equations
60          other_quantities:     List of other quantity names
61
62          tagged_elements:
63          ...
64
65
66        """
67
68        # Determine whether source is a mesh filename or coordinates
69        if type(source) == types.StringType:
70            mesh_filename = source
71        else:
72            coordinates = source
73
74
75        # In case a filename has been specified, extract content
76        if mesh_filename is not None:
77            coordinates, triangles, boundary, vertex_quantity_dict, \
78                         tagged_elements, geo_reference = \
79                         pmesh_to_domain(file_name=mesh_filename,
80                                         use_cache=use_cache,
81                                         verbose=verbose)
82
83
84        # Initialise underlying mesh structure
85        Mesh.__init__(self, coordinates, triangles, boundary,
86                      tagged_elements, geo_reference, use_inscribed_circle,
87                      verbose=verbose)
88
89        if verbose: print 'Initialising Domain'
90        from Numeric import zeros, Float, Int, ones
91        from quantity import Quantity, Conserved_quantity
92
93        # List of quantity names entering
94        # the conservation equations
95        if conserved_quantities is None:
96            self.conserved_quantities = []
97        else:
98            self.conserved_quantities = conserved_quantities
99
100        # List of other quantity names
101        if other_quantities is None:
102            self.other_quantities = []
103        else:
104            self.other_quantities = other_quantities
105
106
107        #Build dictionary of Quantity instances keyed by quantity names
108        self.quantities = {}
109
110        #FIXME: remove later - maybe OK, though....
111        for name in self.conserved_quantities:
112            self.quantities[name] = Conserved_quantity(self)
113        for name in self.other_quantities:
114            self.quantities[name] = Quantity(self)
115
116        #Create an empty list for explicit forcing terms
117        self.forcing_terms = []
118
119        #Setup the ghost cell communication
120        if full_send_dict is None:
121            self.full_send_dict = {}
122        else:
123            self.full_send_dict  = full_send_dict
124
125        # List of other quantity names
126        if ghost_recv_dict  is None:
127            self.ghost_recv_dict  = {}
128        else:
129            self.ghost_recv_dict  = ghost_recv_dict
130
131        self.processor = processor
132        self.numproc   = numproc
133
134        # Setup Communication Buffers
135
136        if verbose: print 'Domain: Set up communication buffers (parallel)'
137        self.nsys = len(self.conserved_quantities)
138        for key in self.full_send_dict:
139            buffer_shape = self.full_send_dict[key][0].shape[0]
140            self.full_send_dict[key].append(zeros( (buffer_shape,self.nsys) ,Float))
141
142
143        for key in self.ghost_recv_dict:
144            buffer_shape = self.ghost_recv_dict[key][0].shape[0]
145            self.ghost_recv_dict[key].append(zeros( (buffer_shape,self.nsys) ,Float))
146
147
148        # Setup cell full flag
149        # =1 for full
150        # =0 for ghost
151        N=self.number_of_elements
152        self.tri_full_flag = ones(N, Int)
153        for i in self.ghost_recv_dict.keys():
154            for id in self.ghost_recv_dict[i][0]:
155                self.tri_full_flag[id] = 0
156
157
158        #Defaults
159        from anuga.config import max_smallsteps, beta_w, beta_h, epsilon, CFL
160        self.beta_w = beta_w
161        self.beta_h = beta_h
162        self.epsilon = epsilon
163
164        #FIXME: Maybe have separate orders for h-limiter and w-limiter?
165        #Or maybe get rid of order altogether and use beta_w and beta_h
166        self.set_default_order(1)
167        #self.default_order = 1
168        #self.order = self.default_order
169
170        self.smallsteps = 0
171        self.max_smallsteps = max_smallsteps
172        self.number_of_steps = 0
173        self.number_of_first_order_steps = 0
174        self.CFL = CFL
175
176        #Model time
177        self.time = 0.0
178        self.finaltime = None
179        self.min_timestep = self.max_timestep = 0.0
180        self.starttime = 0 #Physical starttime if any (0 is 1 Jan 1970 00:00:00)
181
182        ######OBSOLETE
183        #Origin in UTM coordinates
184        #FIXME: This should be set if read by a msh file
185        #self.zone = zone
186        #self.xllcorner = xllcorner
187        #self.yllcorner = yllcorner
188
189
190        #Checkpointing and storage
191        from anuga.config import default_datadir
192        self.datadir = default_datadir
193        self.filename = 'domain'
194        self.checkpoint = False
195
196        #MH310505 To avoid calculating the flux across each edge twice, keep an integer (boolean) array,
197        #to be used during the flux calculation
198        N=self.number_of_elements
199        self.already_computed_flux = zeros((N, 3), Int)
200
201        if mesh_filename is not None:
202            # If the mesh file passed any quantity values
203            # , initialise with these values.
204            if verbose: print 'Domain: Initialising quantity values'
205            self.set_quantity_vertices_dict(vertex_quantity_dict)
206
207
208        if verbose: print 'Domain: Done'
209
210
211
212
213    def set_default_order(self, n):
214        """Set default (spatial) order to either 1 or 2
215        """
216
217        msg = 'Default order must be either 1 or 2. I got %s' %n
218        assert n in [1,2], msg
219
220        self.default_order = n
221        self.order = self.default_order
222
223
224    #Public interface to Domain
225    def get_conserved_quantities(self, vol_id, vertex=None, edge=None):
226        """Get conserved quantities at volume vol_id
227
228        If vertex is specified use it as index for vertex values
229        If edge is specified use it as index for edge values
230        If neither are specified use centroid values
231        If both are specified an exeception is raised
232
233        Return value: Vector of length == number_of_conserved quantities
234
235        """
236
237        from Numeric import zeros, Float
238
239        if not (vertex is None or edge is None):
240            msg = 'Values for both vertex and edge was specified.'
241            msg += 'Only one (or none) is allowed.'
242            raise msg
243
244        q = zeros( len(self.conserved_quantities), Float)
245
246        for i, name in enumerate(self.conserved_quantities):
247            Q = self.quantities[name]
248            if vertex is not None:
249                q[i] = Q.vertex_values[vol_id, vertex]
250            elif edge is not None:
251                q[i] = Q.edge_values[vol_id, edge]
252            else:
253                q[i] = Q.centroid_values[vol_id]
254
255        return q
256
257    def set_time(self, time=0.0):
258        """Set the model time (seconds)"""
259
260        self.time = time
261
262    def set_quantity_vertices_dict(self, quantity_dict):
263        """Set values for named quantities.
264        The index is the quantity
265
266        name: Name of quantity
267        X: Compatible list, Numeric array, const or function (see below)
268
269        The values will be stored in elements following their
270        internal ordering.
271
272        """
273        for key in quantity_dict.keys():
274            self.set_quantity(key, quantity_dict[key], location='vertices')
275
276
277    def set_quantity(self, name, *args, **kwargs):
278        """Set values for named quantity
279
280
281        One keyword argument is documented here:
282        expression = None, # Arbitrary expression
283
284        expression:
285          Arbitrary expression involving quantity names
286
287        See Quantity.set_values for further documentation.
288        """
289
290        #FIXME (Ole): Allow new quantities here
291        #from quantity import Quantity, Conserved_quantity
292        #Create appropriate quantity object
293        ##if name in self.conserved_quantities:
294        ##    self.quantities[name] = Conserved_quantity(self)
295        ##else:
296        ##    self.quantities[name] = Quantity(self)
297
298
299        #Do the expression stuff
300        if kwargs.has_key('expression'):
301            expression = kwargs['expression']
302            del kwargs['expression']
303
304            Q = self.create_quantity_from_expression(expression)
305            kwargs['quantity'] = Q
306
307        #Assign values
308        self.quantities[name].set_values(*args, **kwargs)
309
310
311    def get_quantity(self, name, location='vertices', indices = None):
312        """Get quantity object.
313
314        name: Name of quantity
315
316        See methods inside the quantity object for more options
317        """
318
319        return self.quantities[name] #.get_values( location, indices = indices)
320
321
322    def get_quantity_object(self, name):
323        """Get object for named quantity
324
325        name: Name of quantity
326
327        FIXME: Obsolete
328        """
329
330        print 'get_quantity_object has been deprecated. Please use get_quantity'
331        return self.quantities[name]
332
333
334    def create_quantity_from_expression(self, expression):
335        """Create new quantity from other quantities using arbitrary expression
336
337        Combine existing quantities in domain using expression and return
338        result as a new quantity.
339
340        Note, the new quantity could e.g. be used in set_quantity
341
342        Valid expressions are limited to operators defined in class Quantity
343
344        Example:
345
346
347        """
348
349        from anuga.abstract_2d_finite_volumes.util import apply_expression_to_dictionary
350        return apply_expression_to_dictionary(expression, self.quantities)
351
352
353
354
355    def set_boundary(self, boundary_map):
356        """Associate boundary objects with tagged boundary segments.
357
358        Input boundary_map is a dictionary of boundary objects keyed
359        by symbolic tags to matched against tags in the internal dictionary
360        self.boundary.
361
362        As result one pointer to a boundary object is stored for each vertex
363        in the list self.boundary_objects.
364        More entries may point to the same boundary object
365
366        Schematically the mapping is from two dictionaries to one list
367        where the index is used as pointer to the boundary_values arrays
368        within each quantity.
369
370        self.boundary:          (vol_id, edge_id): tag
371        boundary_map (input):   tag: boundary_object
372        ----------------------------------------------
373        self.boundary_objects:  ((vol_id, edge_id), boundary_object)
374
375
376        Pre-condition:
377          self.boundary has been built.
378
379        Post-condition:
380          self.boundary_objects is built
381
382        If a tag from the domain doesn't appear in the input dictionary an
383        exception is raised.
384        However, if a tag is not used to the domain, no error is thrown.
385        FIXME: This would lead to implementation of a
386        default boundary condition
387
388        Note: If a segment is listed in the boundary dictionary and if it is
389        not None, it *will* become a boundary -
390        even if there is a neighbouring triangle.
391        This would be the case for internal boundaries
392
393        Boundary objects that are None will be skipped.
394
395        FIXME: If set_boundary is called multiple times and if Boundary
396        object is changed into None, the neighbour structure will not be
397        restored!!!
398
399
400        """
401
402        self.boundary_objects = []
403        self.boundary_map = boundary_map  #Store for use with eg. boundary_stats.
404
405        #FIXME: Try to remove the sorting and fix test_mesh.py
406        x = self.boundary.keys()
407        x.sort()
408
409        #Loop through edges that lie on the boundary and associate them with
410        #callable boundary objects depending on their tags
411        for k, (vol_id, edge_id) in enumerate(x):
412            tag = self.boundary[ (vol_id, edge_id) ]
413
414            if boundary_map.has_key(tag):
415                B = boundary_map[tag]  #Get callable boundary object
416
417                if B is not None:
418                    self.boundary_objects.append( ((vol_id, edge_id), B) )
419                    self.neighbours[vol_id, edge_id] = -len(self.boundary_objects)
420                else:
421                    pass
422                    #FIXME: Check and perhaps fix neighbour structure
423
424
425            else:
426                msg = 'ERROR (domain.py): Tag "%s" has not been ' %tag
427                msg += 'bound to a boundary object.\n'
428                msg += 'All boundary tags defined in domain must appear '
429                msg += 'in the supplied dictionary.\n'
430                msg += 'The tags are: %s' %self.get_boundary_tags()
431                raise msg
432
433
434    def set_region(self, *args, **kwargs):
435        """
436        This method is used to set quantities based on a regional tag.
437       
438        It is most often called with the following parameters;
439        (self, tag, quantity, X, location='vertices')
440        tag: the name of the regional tag used to specify the region
441        quantity: Name of quantity to change
442        X: const or function - how the quantity is changed
443        location: Where values are to be stored.
444            Permissible options are: vertices, centroid and unique vertices
445
446        A callable region class or a list of callable region classes
447        can also be passed into this function.
448        """
449        #print "*args", args
450        #print "**kwargs", kwargs
451        if len(args) == 1:
452            self._set_region(*args, **kwargs)
453        else:
454            #Assume it is arguments for the region.set_region function
455            func = region_set_region(*args, **kwargs)
456            self._set_region(func)
457           
458       
459    def _set_region(self, functions):
460        # The order of functions in the list is used.
461        if type(functions) not in [types.ListType,types.TupleType]:
462            functions = [functions]
463        for function in functions:
464            for tag in self.tagged_elements.keys():
465                function(tag, self.tagged_elements[tag], self)
466
467
468    #MISC
469    def check_integrity(self):
470        Mesh.check_integrity(self)
471
472        for quantity in self.conserved_quantities:
473            msg = 'Conserved quantities must be a subset of all quantities'
474            assert quantity in self.quantities, msg
475
476        ##assert hasattr(self, 'boundary_objects')
477
478    def write_time(self):
479        print self.timestepping_statistics()
480
481        #Old version
482        #if self.min_timestep == self.max_timestep:
483        #    print 'Time = %.4f, delta t = %.8f, steps=%d (%d)'\
484        #          %(self.time, self.min_timestep, self.number_of_steps,
485        #            self.number_of_first_order_steps)
486        #elif self.min_timestep > self.max_timestep:
487        #    print 'Time = %.4f, steps=%d (%d)'\
488        #          %(self.time, self.number_of_steps,
489        #            self.number_of_first_order_steps)
490        #else:
491        #    print 'Time = %.4f, delta t in [%.8f, %.8f], steps=%d (%d)'\
492        #          %(self.time, self.min_timestep,
493        #            self.max_timestep, self.number_of_steps,
494        #            self.number_of_first_order_steps)
495
496    def timestepping_statistics(self):
497        """Return string with time stepping statistics for printing or logging
498        """
499
500        msg = ''
501        if self.min_timestep == self.max_timestep:
502            msg += 'Time = %.4f, delta t = %.8f, steps=%d (%d)'\
503                   %(self.time, self.min_timestep, self.number_of_steps,
504                     self.number_of_first_order_steps)
505        elif self.min_timestep > self.max_timestep:
506            msg += 'Time = %.4f, steps=%d (%d)'\
507                   %(self.time, self.number_of_steps,
508                     self.number_of_first_order_steps)
509        else:
510            msg += 'Time = %.4f, delta t in [%.8f, %.8f], steps=%d (%d)'\
511                   %(self.time, self.min_timestep,
512                     self.max_timestep, self.number_of_steps,
513                     self.number_of_first_order_steps)
514
515        return msg
516
517
518    def write_boundary_statistics(self, quantities = None, tags = None):
519        print self.boundary_statistics(quantities, tags)
520
521    def boundary_statistics(self, quantities = None, tags = None):
522        """Output statistics about boundary forcing at each timestep
523
524
525        Input:
526          quantities: either None, a string or a list of strings naming the quantities to be reported
527          tags:       either None, a string or a list of strings naming the tags to be reported
528
529
530        Example output:
531        Tag 'wall':
532            stage in [2, 5.5]
533            xmomentum in []
534            ymomentum in []
535        Tag 'ocean'
536
537
538        If quantities are specified only report on those. Otherwise take all conserved quantities.
539        If tags are specified only report on those, otherwise take all tags.
540
541        """
542
543        #Input checks
544        import types, string
545
546        if quantities is None:
547            quantities = self.conserved_quantities
548        elif type(quantities) == types.StringType:
549            quantities = [quantities] #Turn it into a list
550
551        msg = 'Keyword argument quantities must be either None, '
552        msg += 'string or list. I got %s' %str(quantities)
553        assert type(quantities) == types.ListType, msg
554
555
556        if tags is None:
557            tags = self.get_boundary_tags()
558        elif type(tags) == types.StringType:
559            tags = [tags] #Turn it into a list
560
561        msg = 'Keyword argument tags must be either None, '
562        msg += 'string or list. I got %s' %str(tags)
563        assert type(tags) == types.ListType, msg
564
565        #Determine width of longest quantity name (for cosmetic purposes)
566        maxwidth = 0
567        for name in quantities:
568            w = len(name)
569            if w > maxwidth:
570                maxwidth = w
571
572        #Output stats
573        msg = 'Boundary values at time %.4f:\n' %self.time
574        for tag in tags:
575            msg += '    %s:\n' %tag
576
577            for name in quantities:
578                q = self.quantities[name]
579
580                #Find range of boundary values for tag and q
581                maxval = minval = None
582                for i, ((vol_id, edge_id), B) in\
583                        enumerate(self.boundary_objects):
584                    if self.boundary[(vol_id, edge_id)] == tag:
585                        v = q.boundary_values[i]
586                        if minval is None or v < minval: minval = v
587                        if maxval is None or v > maxval: maxval = v
588
589                if minval is None or maxval is None:
590                    msg += '        Sorry no information available about' +\
591                           ' tag %s and quantity %s\n' %(tag, name)
592                else:
593                    msg += '        %s in [%12.8f, %12.8f]\n'\
594                           %(string.ljust(name, maxwidth), minval, maxval)
595
596
597        return msg
598
599
600    def get_name(self):
601        return self.filename
602
603    def set_name(self, name):
604        self.filename = name
605
606    def get_datadir(self):
607        return self.datadir
608
609    def set_datadir(self, name):
610        self.datadir = name
611
612
613
614    #def set_defaults(self):
615    #    """Set default values for uninitialised quantities.
616    #    Should be overridden or specialised by specific modules
617    #    """#
618    #
619    #    for name in self.conserved_quantities + self.other_quantities:
620    #        self.set_quantity(name, 0.0)
621
622
623    ###########################
624    #Main components of evolve
625
626    def evolve(self,
627               yieldstep = None,
628               finaltime = None,
629               duration = None,
630               skip_initial_step = False):
631        """Evolve model through time starting from self.starttime.
632
633
634        yieldstep: Interval between yields where results are stored,
635                   statistics written and domain inspected or
636                   possibly modified. If omitted the internal predefined
637                   max timestep is used.
638                   Internally, smaller timesteps may be taken.
639
640        duration: Duration of simulation
641
642        finaltime: Time where simulation should end
643
644        If both duration and finaltime are given an exception is thrown.
645
646
647        skip_initial_step: Boolean flag that decides whether the first
648        yield step is skipped or not. This is useful for example to avoid
649        duplicate steps when multiple evolve processes are dove tailed.
650
651
652        Evolve is implemented as a generator and is to be called as such, e.g.
653
654        for t in domain.evolve(yieldstep, finaltime):
655            <Do something with domain and t>
656
657
658        All times are given in seconds
659
660        """
661
662        from anuga.config import min_timestep, max_timestep, epsilon
663
664        #FIXME: Maybe lump into a larger check prior to evolving
665        msg = 'Boundary tags must be bound to boundary objects before evolving system, '
666        msg += 'e.g. using the method set_boundary.\n'
667        msg += 'This system has the boundary tags %s '\
668               %self.get_boundary_tags()
669        assert hasattr(self, 'boundary_objects'), msg
670
671        ##self.set_defaults()
672
673        if yieldstep is None:
674            yieldstep = max_timestep
675        else:
676            yieldstep = float(yieldstep)
677
678        self.order = self.default_order
679
680
681        if finaltime is not None and duration is not None:
682            print 'F', finaltime, duration
683            msg = 'Only one of finaltime and duration may be specified'
684            raise msg
685        else:
686            if finaltime is not None:
687                self.finaltime = float(finaltime)
688            if duration is not None:
689                self.finaltime = self.starttime + float(duration)
690
691
692
693
694        self.yieldtime = 0.0 #Time between 'yields'
695
696        #Initialise interval of timestep sizes (for reporting only)
697        self.min_timestep = max_timestep
698        self.max_timestep = min_timestep
699        self.number_of_steps = 0
700        self.number_of_first_order_steps = 0
701
702        #update ghosts
703        self.update_ghosts()
704
705        #Initial update of vertex and edge values
706        self.distribute_to_vertices_and_edges()
707
708        #Initial update boundary values
709        self.update_boundary()
710
711        #Or maybe restore from latest checkpoint
712        if self.checkpoint is True:
713            self.goto_latest_checkpoint()
714
715        if skip_initial_step is False:
716            yield(self.time)  #Yield initial values
717
718        while True:
719
720            #Compute fluxes across each element edge
721            self.compute_fluxes()
722
723            #Update timestep to fit yieldstep and finaltime
724            self.update_timestep(yieldstep, finaltime)
725
726            #Update conserved quantities
727            self.update_conserved_quantities()
728
729            #update ghosts
730            self.update_ghosts()
731
732            #Update vertex and edge values
733            self.distribute_to_vertices_and_edges()
734
735            #Update boundary values
736            self.update_boundary()
737
738            #Update time
739            self.time += self.timestep
740            self.yieldtime += self.timestep
741            self.number_of_steps += 1
742            if self.order == 1:
743                self.number_of_first_order_steps += 1
744
745            #Yield results
746            if finaltime is not None and self.time >= finaltime:
747
748                if self.time > finaltime:
749                    #FIXME (Ole, 30 April 2006): Do we need this check?
750                    print 'WARNING (domain.py): time overshot finaltime. Contact Ole.Nielsen@ga.gov.au'
751                    self.time = finaltime
752
753                # Yield final time and stop
754                yield(self.time)
755                break
756
757
758            if self.yieldtime >= yieldstep:
759                # Yield (intermediate) time and allow inspection of domain
760
761                if self.checkpoint is True:
762                    self.store_checkpoint()
763                    self.delete_old_checkpoints()
764
765                #Pass control on to outer loop for more specific actions
766                yield(self.time)
767
768                # Reinitialise
769                self.yieldtime = 0.0
770                self.min_timestep = max_timestep
771                self.max_timestep = min_timestep
772                self.number_of_steps = 0
773                self.number_of_first_order_steps = 0
774
775
776    def evolve_to_end(self, finaltime = 1.0):
777        """Iterate evolve all the way to the end
778        """
779
780        for _ in self.evolve(yieldstep=None, finaltime=finaltime):
781            pass
782
783
784
785    def update_boundary(self):
786        """Go through list of boundary objects and update boundary values
787        for all conserved quantities on boundary.
788        """
789
790        #FIXME: Update only those that change (if that can be worked out)
791        #FIXME: Boundary objects should not include ghost nodes.
792        for i, ((vol_id, edge_id), B) in enumerate(self.boundary_objects):
793            if B is None:
794                print 'WARNING: Ignored boundary segment %d (None)'
795            else:
796                q = B.evaluate(vol_id, edge_id)
797
798                for j, name in enumerate(self.conserved_quantities):
799                    Q = self.quantities[name]
800                    Q.boundary_values[i] = q[j]
801
802
803    def compute_fluxes(self):
804        msg = 'Method compute_fluxes must be overridden by Domain subclass'
805        raise msg
806
807
808    def update_timestep(self, yieldstep, finaltime):
809
810        from anuga.config import min_timestep, max_timestep
811
812        # self.timestep is calculated from speed of characteristics
813        # Apply CFL condition here
814        timestep = min(self.CFL*self.timestep,max_timestep)
815
816        #Record maximal and minimal values of timestep for reporting
817        self.max_timestep = max(timestep, self.max_timestep)
818        self.min_timestep = min(timestep, self.min_timestep)
819
820        #Protect against degenerate time steps
821        if timestep < min_timestep:
822
823            #Number of consecutive small steps taken b4 taking action
824            self.smallsteps += 1
825
826            if self.smallsteps > self.max_smallsteps:
827                self.smallsteps = 0 #Reset
828
829                if self.order == 1:
830                    msg = 'WARNING: Too small timestep %.16f reached '\
831                          %timestep
832                    msg += 'even after %d steps of 1 order scheme'\
833                           %self.max_smallsteps
834                    print msg
835                    timestep = min_timestep  #Try enforcing min_step
836
837                    #raise msg
838                else:
839                    #Try to overcome situation by switching to 1 order
840                    self.order = 1
841
842        else:
843            self.smallsteps = 0
844            if self.order == 1 and self.default_order == 2:
845                self.order = 2
846
847
848        #Ensure that final time is not exceeded
849        if finaltime is not None and self.time + timestep > finaltime:
850            timestep = finaltime-self.time
851
852        #Ensure that model time is aligned with yieldsteps
853        if self.yieldtime + timestep > yieldstep:
854            timestep = yieldstep-self.yieldtime
855
856        self.timestep = timestep
857
858
859
860    def compute_forcing_terms(self):
861        """If there are any forcing functions driving the system
862        they should be defined in Domain subclass and appended to
863        the list self.forcing_terms
864        """
865
866        for f in self.forcing_terms:
867            f(self)
868
869
870
871    def update_conserved_quantities(self):
872        """Update vectors of conserved quantities using previously
873        computed fluxes specified forcing functions.
874        """
875
876        from Numeric import ones, sum, equal, Float
877
878        N = self.number_of_elements
879        d = len(self.conserved_quantities)
880
881        timestep = self.timestep
882
883        #Compute forcing terms
884        self.compute_forcing_terms()
885
886        #Update conserved_quantities
887        for name in self.conserved_quantities:
888            Q = self.quantities[name]
889            Q.update(timestep)
890
891            #Clean up
892            #Note that Q.explicit_update is reset by compute_fluxes
893
894            #MH090605 commented out the following since semi_implicit_update is now re-initialized
895            #at the end of the _update function in quantity_ext.c (This is called by the
896            #preceeding Q.update(timestep) statement above).
897            #For run_profile.py with N=128, the time of update_conserved_quantities is cut from 14.00 secs
898            #to 8.35 secs
899
900            #Q.semi_implicit_update[:] = 0.0
901
902    def update_ghosts(self):
903        pass
904
905    def distribute_to_vertices_and_edges(self):
906        """Extrapolate conserved quantities from centroid to
907        vertices and edge-midpoints for each volume
908
909        Default implementation is straight first order,
910        i.e. constant values throughout each element and
911        no reference to non-conserved quantities.
912        """
913
914        for name in self.conserved_quantities:
915            Q = self.quantities[name]
916            if self.order == 1:
917                Q.extrapolate_first_order()
918            elif self.order == 2:
919                Q.extrapolate_second_order()
920                Q.limit()
921            else:
922                raise 'Unknown order'
923            Q.interpolate_from_vertices_to_edges()
924
925
926    def centroid_norm(self, quantity, normfunc):
927        """Calculate the norm of the centroid values
928        of a specific quantity, using normfunc.
929
930        normfunc should take a list to a float.
931
932        common normfuncs are provided in the module utilities.norms
933        """
934        return normfunc(self.quantities[quantity].centroid_values)
935
936
937
938##############################################
939#Initialise module
940
941#Optimisation with psyco
942from anuga.config import use_psyco
943if use_psyco:
944    try:
945        import psyco
946    except:
947        import os
948        if os.name == 'posix' and os.uname()[4] == 'x86_64':
949            pass
950            #Psyco isn't supported on 64 bit systems, but it doesn't matter
951        else:
952            msg = 'WARNING: psyco (speedup) could not import'+\
953                  ', you may want to consider installing it'
954            print msg
955    else:
956        psyco.bind(Domain.update_boundary)
957        #psyco.bind(Domain.update_timestep)     #Not worth it
958        psyco.bind(Domain.update_conserved_quantities)
959        psyco.bind(Domain.distribute_to_vertices_and_edges)
960
961
962if __name__ == "__main__":
963    pass
Note: See TracBrowser for help on using the repository browser.