source: inundation/pyvolution/domain.py @ 3354

Last change on this file since 3354 was 3085, checked in by ole, 19 years ago

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