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

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