source: inundation/pyvolution/domain.py @ 2319

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

Comments

File size: 25.1 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
369        Input:
370          quantities: either None, a string or a list of strings naming the quantities to be reported
371          tags:       either None, a string or a list of strings naming the tags to be reported         
372
373
374        Example output:
375        Tag 'wall':
376            stage in [2, 5.5]
377            xmomentum in []
378            ymomentum in []
379        Tag 'ocean'
380
381
382        If quantities are specified only report on those. Otherwise take all conserved quantities.
383        If tags are specified only report on those, otherwise take all tags.
384
385        """
386
387        #Input checks
388        import types, string
389
390        if quantities is None:
391            quantities = self.conserved_quantities
392        elif type(quantities) == types.StringType:
393            quantities = [quantities] #Turn it into a list
394
395        msg = 'Keyword argument quantities must be either None, '
396        msg += 'string or list. I got %s' %str(quantities)
397        assert type(quantities) == types.ListType, msg
398
399
400        if tags is None:
401            tags = self.get_boundary_tags()
402        elif type(tags) == types.StringType:
403            tags = [tags] #Turn it into a list
404
405        msg = 'Keyword argument tags must be either None, '
406        msg += 'string or list. I got %s' %str(tags)
407        assert type(tags) == types.ListType, msg
408
409        #Determine width of longest quantity name (for cosmetic purposes)
410        maxwidth = 0
411        for name in quantities:
412            w = len(name)
413            if w > maxwidth:
414                maxwidth = w
415
416        #Output stats
417        msg = 'Boundary values at time %.4f:\n' %self.time
418        for tag in tags:
419            msg += '    %s:\n' %tag
420
421            for name in quantities:
422                q = self.quantities[name]
423
424                #Find range of boundary values for tag and q
425                maxval = minval = None
426                for i, ((vol_id, edge_id), B) in\
427                        enumerate(self.boundary_objects):
428                    if self.boundary[(vol_id, edge_id)] == tag:
429                        v = q.boundary_values[i]
430                        if minval is None or v < minval: minval = v
431                        if maxval is None or v > maxval: maxval = v
432
433                if minval is None or maxval is None:
434                    msg += '        Sorry no information available about' +\
435                           ' tag %s and quantity %s\n' %(tag, name)
436                else:
437                    msg += '        %s in [%12.8f, %12.8f]\n'\
438                           %(string.ljust(name, maxwidth), minval, maxval)
439
440
441        return msg
442
443
444    def get_name(self):
445        return self.filename
446
447    def set_name(self, name):
448        self.filename = name
449
450    def get_datadir(self):
451        return self.datadir
452
453    def set_datadir(self, name):
454        self.datadir = name
455
456
457
458    #def set_defaults(self):
459    #    """Set default values for uninitialised quantities.
460    #    Should be overridden or specialised by specific modules
461    #    """#
462    #
463    #    for name in self.conserved_quantities + self.other_quantities:
464    #        self.set_quantity(name, 0.0)
465
466
467    ###########################
468    #Main components of evolve
469
470    def evolve(self, yieldstep = None, finaltime = None,
471               skip_initial_step = False):
472        """Evolve model from time=0.0 to finaltime yielding results
473        every yieldstep.
474
475        Internally, smaller timesteps may be taken.
476
477        Evolve is implemented as a generator and is to be called as such, e.g.
478
479        for t in domain.evolve(timestep, yieldstep, finaltime):
480            <Do something with domain and t>
481
482        """
483
484        from config import min_timestep, max_timestep, epsilon
485
486        #FIXME: Maybe lump into a larger check prior to evolving
487        msg = 'Boundary tags must be bound to boundary objects before evolving system, '
488        msg += 'e.g. using the method set_boundary.\n'
489        msg += 'This system has the boundary tags %s ' %self.get_boundary_tags()
490        assert hasattr(self, 'boundary_objects'), msg
491
492        ##self.set_defaults()
493
494        if yieldstep is None:
495            yieldstep = max_timestep
496        else:
497            yieldstep = float(yieldstep)
498
499        self.order = self.default_order
500
501
502        self.yieldtime = 0.0 #Time between 'yields'
503
504        #Initialise interval of timestep sizes (for reporting only)
505        self.min_timestep = max_timestep
506        self.max_timestep = min_timestep
507        self.finaltime = finaltime
508        self.number_of_steps = 0
509        self.number_of_first_order_steps = 0
510
511        #update ghosts
512        self.update_ghosts()
513
514        #Initial update of vertex and edge values
515        self.distribute_to_vertices_and_edges()
516
517        #Initial update boundary values
518        self.update_boundary()
519
520        #Or maybe restore from latest checkpoint
521        if self.checkpoint is True:
522            self.goto_latest_checkpoint()
523
524        if skip_initial_step is False:
525            yield(self.time)  #Yield initial values
526
527        while True:
528
529            #Compute fluxes across each element edge
530            self.compute_fluxes()
531
532            #Update timestep to fit yieldstep and finaltime
533            self.update_timestep(yieldstep, finaltime)
534
535            #Update conserved quantities
536            self.update_conserved_quantities()
537
538            #update ghosts
539            self.update_ghosts()
540
541            #Update vertex and edge values
542            self.distribute_to_vertices_and_edges()
543
544            #Update boundary values
545            self.update_boundary()
546
547            #Update time
548            self.time += self.timestep
549            self.yieldtime += self.timestep
550            self.number_of_steps += 1
551            if self.order == 1:
552                self.number_of_first_order_steps += 1
553
554            #Yield results
555            if finaltime is not None and abs(self.time - finaltime) < epsilon:
556
557                #FIXME: There is a rare situation where the
558                #final time step is stored twice. Can we make a test?
559
560                # Yield final time and stop
561                yield(self.time)
562                break
563
564
565            if abs(self.yieldtime - yieldstep) < epsilon:
566                # Yield (intermediate) time and allow inspection of domain
567
568                if self.checkpoint is True:
569                    self.store_checkpoint()
570                    self.delete_old_checkpoints()
571
572                #Pass control on to outer loop for more specific actions
573                yield(self.time)
574
575                # Reinitialise
576                self.yieldtime = 0.0
577                self.min_timestep = max_timestep
578                self.max_timestep = min_timestep
579                self.number_of_steps = 0
580                self.number_of_first_order_steps = 0
581
582
583    def evolve_to_end(self, finaltime = 1.0):
584        """Iterate evolve all the way to the end
585        """
586
587        for _ in self.evolve(yieldstep=None, finaltime=finaltime):
588            pass
589
590
591
592    def update_boundary(self):
593        """Go through list of boundary objects and update boundary values
594        for all conserved quantities on boundary.
595        """
596
597        #FIXME: Update only those that change (if that can be worked out)
598        #FIXME: Boundary objects should not include ghost nodes.
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
732
733
734##############################################
735#Initialise module
736
737#Optimisation with psyco
738from config import use_psyco
739if use_psyco:
740    try:
741        import psyco
742    except:
743        import os
744        if os.name == 'posix' and os.uname()[4] == 'x86_64':
745            pass
746            #Psyco isn't supported on 64 bit systems, but it doesn't matter
747        else:
748            msg = 'WARNING: psyco (speedup) could not import'+\
749                  ', you may want to consider installing it'
750            print msg
751    else:
752        psyco.bind(Domain.update_boundary)
753        #psyco.bind(Domain.update_timestep)     #Not worth it
754        psyco.bind(Domain.update_conserved_quantities)
755        psyco.bind(Domain.distribute_to_vertices_and_edges)
756
757
758if __name__ == "__main__":
759    pass
Note: See TracBrowser for help on using the repository browser.