source: inundation/pyvolution/domain.py @ 2851

Last change on this file since 2851 was 2851, checked in by ole, 18 years ago

Another stab at removing potential for duplicate timesteps

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