source: inundation/pyvolution/domain.py @ 2793

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

Renaming and comments

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