source: inundation/pyvolution/domain.py @ 3447

Last change on this file since 3447 was 3447, checked in by duncan, 18 years ago

Making set_region easier to use.

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