source: inundation/ga/storm_surge/pyvolution/domain.py @ 1454

Last change on this file since 1454 was 1402, checked in by steve, 20 years ago

Periodic boundary working via ghosts

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