source: inundation-numpy-branch/pyvolution/domain.py @ 2608

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

Played with custom exceptions for ANUGA

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