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

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

Played with parallel example

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