source: inundation/pyvolution/domain.py @ 2267

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