source: inundation/pyvolution/domain.py @ 1916

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

Implemented operator overloading for (pow) in class Quantity and tested.
Wrote create_quantity_from_expression and test.

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