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

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

Refactored references to domain.filename away.
Use

domain.set_name()
domain.get_name()

instead.

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