source: inundation/pyvolution/domain.py @ 1740

Last change on this file since 1740 was 1697, checked in by steve, 19 years ago

Found problem with File_Boundary as used in validation test LWRU2. Have setup new BC Transmissive_Momentum_Set _Stage

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