source: anuga_work/development/2010-projects/anuga_1d/base/generic_domain.py @ 7860

Last change on this file since 7860 was 7860, checked in by steve, 12 years ago

Continuing to numpy the for loops

File size: 45.9 KB
Line 
1"""Class Domain - 1D 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
10import numpy
11from generic_boundary_conditions import *
12
13
14class Generic_domain:
15
16    def __init__(self,
17                 coordinates,
18                 boundary = None,
19                 conserved_quantities = None,
20                 evolved_quantities = None,
21                 other_quantities = None,
22                 tagged_elements = None):
23        """
24        Build 1D elements from x coordinates
25        """
26
27        from config import timestepping_method
28        from config import CFL
29       
30        #Store Points
31        self.coordinates = numpy.array(coordinates, numpy.float)
32
33
34        #Register number of Elements
35        self.number_of_elements = N = len(self.coordinates)-1
36
37        self.beta = 1.0
38        self.set_limiter("minmod_kurganov")
39        self.set_CFL(CFL)
40        self.set_timestepping_method(timestepping_method)
41
42
43        # work array for special vertices (like wet or dry cells
44        self.special_vertices = numpy.zeros((N,2), numpy.int) # should this be here
45
46        #Allocate space for neighbour and boundary structures
47        self.neighbours = numpy.zeros((N, 2), numpy.int)
48        #self.neighbour_edges = numpy.zeros((N, 2), numpy.int)
49        self.neighbour_vertices   = numpy.zeros((N, 2), numpy.int)
50        self.number_of_boundaries = numpy.zeros(N, numpy.int)
51        self.surrogate_neighbours = numpy.zeros((N, 2), numpy.int)
52       
53        #Allocate space for geometric quantities
54        self.vertices  = numpy.zeros((N, 2), numpy.float)
55        self.centroids = numpy.zeros(N, numpy.float)
56        self.areas     = numpy.zeros(N, numpy.float)
57
58        self.max_speed_array = numpy.zeros(N, numpy.float)
59     
60
61        self.normals = numpy.zeros((N, 2), numpy.float)
62
63
64        xl = self.coordinates[0:-1]
65        xr = self.coordinates[1: ]
66
67        self.vertices[:,0] = xl
68        self.vertices[:,1] = xr
69
70        self.centroids[:] = (xl+xr)/2.0
71
72        msg = 'Coordinates should be ordered, smallest to largest'
73        assert (xr>xl).all(), msg
74
75        #The normal vectors
76        # - point outward from each edge
77        # - are orthogonal to the edge
78        # - have unit length
79        # - Are enumerated by left vertex then right vertex normals
80
81
82        self.normals[:,0] = -1.0
83        self.normals[:,1] =  1.0
84
85        self.areas[:] = (xr-xl)
86
87        #Initialise Neighbours (-1 means that it is a boundary neighbour)
88
89        self.neighbours[:, :] = -1
90
91        #Initialise edge ids of neighbours
92        self.neighbour_vertices[:, :] = -1
93
94
95#        for i in range(N):
96#            xl = self.coordinates[i]
97#            xr = self.coordinates[i+1]
98#            self.vertices[i,0] = xl
99#            self.vertices[i,1] = xr
100#
101#            centroid = (xl+xr)/2.0
102#            self.centroids[i] = centroid
103#
104#            msg = 'Coordinates should be ordered, smallest to largest'
105#            assert xr>xl, msg
106#
107#            #The normal vectors
108#            # - point outward from each edge
109#            # - are orthogonal to the edge
110#            # - have unit length
111#            # - Are enumerated by left vertex then right vertex normals
112#
113#            nl = -1.0
114#            nr =  1.0
115#            self.normals[i,:] = [nl, nr]
116#
117#            self.areas[i] = (xr-xl)
118#
119## #         print 'N', N
120## #         print 'Centroid', self.centroids
121## #         print 'Areas', self.areas
122## #         print 'Vertex_Coordinates', self.vertices
123#
124#            #Initialise Neighbours (-1 means that it is a boundary neighbour)
125#            self.neighbours[i, :] = [-1, -1]
126#            #Initialise edge ids of neighbours
127#            #Initialise vertex ids of neighbours
128#            #In case of boundaries this slot is not used
129#            #self.neighbour_edges[i, :] = [-1, -1]
130#            self.neighbour_vertices[i, :] = [-1, -1]
131
132        self.build_vertexlist()
133
134        #Build neighbour structure
135        self.build_neighbour_structure()
136
137        #Build surrogate neighbour structure
138        self.build_surrogate_neighbour_structure()         
139
140        #Build boundary dictionary mapping (id, vertex) to symbolic tags
141        self.build_boundary_dictionary(boundary)
142
143        #Build tagged element  dictionary mapping (tag) to array of elements
144        self.build_tagged_elements_dictionary(tagged_elements)
145       
146        from quantity import Quantity
147        #from quantity_domain import Quantity, Conserved_quantity
148       
149        #List of quantity names entering
150        #the conservation equations
151        #(Must be a subset of quantities)
152        if conserved_quantities is None:
153            self.conserved_quantities = []
154        else:
155            self.conserved_quantities = conserved_quantities
156
157        if evolved_quantities is None:
158            self.evolved_quantities = self.conserved_quantities
159        else:
160            self.evolved_quantities = evolved_quantities
161
162        if other_quantities is None:
163            self.other_quantities = []
164        else:
165            self.other_quantities = other_quantities
166
167
168        #Build dictionary of Quantity instances keyed by quantity names
169        self.quantities = {}
170
171        #print self.conserved_quantities
172        #print self.evolved_quantities
173       
174
175        #FIXME: remove later - maybe OK, though....
176        for name in self.evolved_quantities:
177            self.quantities[name] = Quantity(self)
178        for name in self.other_quantities:
179            self.quantities[name] = Quantity(self)
180
181        #Create an empty list for explicit forcing terms
182        self.forcing_terms = []
183
184        #Defaults
185        from config import max_smallsteps, beta_w, beta_h, epsilon, CFL
186        self.beta_w = beta_w
187        self.beta_h = beta_h
188        self.epsilon = epsilon
189
190        #FIXME: Maybe have separate orders for h-limiter and w-limiter?
191        #Or maybe get rid of order altogether and use beta_w and beta_h
192        self.default_order = 1
193        self.order = self.default_order
194
195        self.default_time_order = 1
196        self.time_order = self.default_time_order
197       
198        self.smallsteps = 0
199        self.max_smallsteps = max_smallsteps
200        self.number_of_steps = 0
201        self.number_of_first_order_steps = 0
202
203        #Model time
204        self.time = 0.0
205        self.finaltime = None
206        self.min_timestep = self.max_timestep = 0.0
207        self.starttime = 0 #Physical starttime if any (0 is 1 Jan 1970 00:00:00)
208        #Checkpointing and storage
209        from config import default_datadir
210        self.set_datadir(default_datadir)
211        self.filename = 'domain'
212        self.checkpoint = False
213
214    def __len__(self):
215        return self.number_of_elements
216
217    def build_vertexlist(self):
218        """Build vertexlist index by vertex ids and for each entry (point id)
219        build a list of (triangles, vertex_id) pairs that use the point
220        as vertex.
221
222        Preconditions:
223          self.coordinates and self.triangles are defined
224
225        Postcondition:
226          self.vertexlist is built
227        """
228
229        vertexlist = [None]*len(self.coordinates)
230        for i in range(self.number_of_elements):
231
232            #a = self.triangles[i, 0]
233            #b = self.triangles[i, 1]
234            #c = self.triangles[i, 2]
235            a = i
236            b = i + 1
237
238            #Register the vertices v as lists of
239            #(triangle_id, vertex_id) tuples associated with them
240            #This is used for smoothing
241            #for vertex_id, v in enumerate([a,b,c]):
242            for vertex_id, v in enumerate([a,b]):
243                if vertexlist[v] is None:
244                    vertexlist[v] = []
245
246                vertexlist[v].append( (i, vertex_id) )
247
248        self.vertexlist = vertexlist
249
250       
251    def build_neighbour_structure(self):
252        """Update all registered triangles to point to their neighbours.
253
254        Also, keep a tally of the number of boundaries for each triangle
255
256        Postconditions:
257          neighbours and neighbour_edges is populated
258          neighbours and neighbour_vertices is populated
259          number_of_boundaries integer array is defined.
260        """
261
262        #Step 1:
263        #Build dictionary mapping from segments (2-tuple of points)
264        #to left hand side edge (facing neighbouring triangle)
265
266        N = self.number_of_elements
267        neighbourdict = {}
268        #l_edge = 0
269        #r_edge = 1
270        l_vertex = 0
271        r_vertex = 1
272        for i in range(N):
273
274            #Register all segments as keys mapping to current triangle
275            #and segment id
276            #a = self.triangles[i, 0]
277            #b = self.triangles[i, 1]
278            #c = self.triangles[i, 2]
279            a = self.vertices[i,0]
280            b = self.vertices[i,1]
281           
282            """
283            if neighbourdict.has_key((a,b)):
284                    msg = "Edge 2 of triangle %d is duplicating edge %d of triangle %d.\n" %(i,neighbourdict[a,b][1],neighbourdict[a,b][0])
285                    raise msg
286            if neighbourdict.has_key((b,c)):
287                    msg = "Edge 0 of triangle %d is duplicating edge %d of triangle %d.\n" %(i,neighbourdict[b,c][1],neighbourdict[b,c][0])
288                    raise msg
289            if neighbourdict.has_key((c,a)):
290                    msg = "Edge 1 of triangle %d is duplicating edge %d of triangle %d.\n" %(i,neighbourdict[c,a][1],neighbourdict[c,a][0])
291                    raise msg
292                    """
293            #neighbourdict[a,b] = (i, 2) #(id, edge)
294            #neighbourdict[b,c] = (i, 0) #(id, edge)
295            #neighbourdict[c,a] = (i, 1) #(id, edge)
296            #neighbourdict[a,b] = (i, 1) #(id, edge)
297            #neighbourdict[b,a] = (i, 0) #(id, edge)
298            #neighbourdict[a,l_edge] = (i, 0) #(id, edge)
299            #neighbourdict[b,r_edge] = (i, 1) #(id, edge)
300            neighbourdict[a,l_vertex] = (i, 0) #(id, vertex)
301            neighbourdict[b,r_vertex] = (i, 1) #(id, vertex)
302
303
304        #Step 2:
305        #Go through triangles again, but this time
306        #reverse direction of segments and lookup neighbours.
307        for i in range(N):
308            #a = self.triangles[i, 0]
309            #b = self.triangles[i, 1]
310            #c = self.triangles[i, 2]
311           
312            a = self.vertices[i,0]
313            b = self.vertices[i,1]
314           
315            #self.number_of_boundaries[i] = 3
316            self.number_of_boundaries[i] = 2
317            #if neighbourdict.has_key((b,l_edge)):
318            if neighbourdict.has_key((b,l_vertex)):
319                #self.neighbours[i, 1] = neighbourdict[b,l_edge][0]
320                #self.neighbour_edges[i, 1] = neighbourdict[b,l_edge][1]
321                self.neighbours[i, 1] = neighbourdict[b,l_vertex][0]
322                self.neighbour_vertices[i, 1] = neighbourdict[b,l_vertex][1]
323                self.number_of_boundaries[i] -= 1
324
325            #if neighbourdict.has_key((a,r_edge)):
326            if neighbourdict.has_key((a,r_vertex)):
327                #self.neighbours[i, 0] = neighbourdict[a,r_edge][0]
328                #self.neighbour_edges[i, 0] = neighbourdict[a,r_edge][1]
329                self.neighbours[i, 0] = neighbourdict[a,r_vertex][0]
330                self.neighbour_vertices[i, 0] = neighbourdict[a,r_vertex][1]
331                self.number_of_boundaries[i] -= 1
332               
333            #if neighbourdict.has_key((b,a)):
334            #    self.neighbours[i, 1] = neighbourdict[b,a][0]
335            #    self.neighbour_edges[i, 1] = neighbourdict[b,a][1]
336            #    self.number_of_boundaries[i] -= 1
337
338            #if neighbourdict.has_key((c,b)):
339            #    self.neighbours[i, 0] = neighbourdict[c,b][0]
340            #    self.neighbour_edges[i, 0] = neighbourdict[c,b][1]
341            #    self.number_of_boundaries[i] -= 1
342
343            #if neighbourdict.has_key((a,b)):
344            #    self.neighbours[i, 0] = neighbourdict[a,b][0]
345            #    self.neighbour_edges[i, 0] = neighbourdict[a,b][1]
346            #    self.number_of_boundaries[i] -= 1
347
348    def build_surrogate_neighbour_structure(self):
349        """Build structure where each triangle edge points to its neighbours
350        if they exist. Otherwise point to the triangle itself.
351
352        The surrogate neighbour structure is useful for computing gradients
353        based on centroid values of neighbours.
354
355        Precondition: Neighbour structure is defined
356        Postcondition:
357          Surrogate neighbour structure is defined:
358          surrogate_neighbours: i0, i1, i2 where all i_k >= 0 point to
359          triangles.
360
361        """
362
363        N = self.number_of_elements
364        for i in range(N):
365            #Find all neighbouring volumes that are not boundaries
366            #for k in range(3):
367            for k in range(2):   
368                if self.neighbours[i, k] < 0:
369                    self.surrogate_neighbours[i, k] = i #Point this triangle
370                else:
371                    self.surrogate_neighbours[i, k] = self.neighbours[i, k]
372
373    def build_boundary_dictionary(self, boundary = None):
374        """Build or check the dictionary of boundary tags.
375         self.boundary is a dictionary of tags,
376         keyed by volume id and edge:
377         { (id, edge): tag, ... }
378
379         Postconditions:
380            self.boundary is defined.
381        """
382
383        from config import default_boundary_tag
384
385        if boundary is None:
386            boundary = {}
387            for vol_id in range(self.number_of_elements):
388                #for edge_id in range(0, 3):
389                #for edge_id in range(0, 2):
390                for vertex_id in range(0, 2):   
391                    #if self.neighbours[vol_id, edge_id] < 0:
392                    if self.neighbours[vol_id, vertex_id] < 0:   
393                        #boundary[(vol_id, edge_id)] = default_boundary_tag
394                        boundary[(vol_id, vertex_id)] = default_boundary_tag
395        else:
396            #Check that all keys in given boundary exist
397            #for vol_id, edge_id in boundary.keys():
398            for vol_id, vertex_id in boundary.keys():   
399                #msg = 'Segment (%d, %d) does not exist' %(vol_id, edge_id)
400                msg = 'Segment (%d, %d) does not exist' %(vol_id, vertex_id)
401                a, b = self.neighbours.shape
402                #assert vol_id < a and edge_id < b, msg
403                assert vol_id < a and vertex_id < b, msg
404
405                #FIXME: This assert violates internal boundaries (delete it)
406                #msg = 'Segment (%d, %d) is not a boundary' %(vol_id, edge_id)
407                #assert self.neighbours[vol_id, edge_id] < 0, msg
408
409            #Check that all boundary segments are assigned a tag
410            for vol_id in range(self.number_of_elements):
411                #for edge_id in range(0, 3):
412                #for edge_id in range(0, 2):
413                for vertex_id in range(0, 2):   
414                    #if self.neighbours[vol_id, edge_id] < 0:
415                    if self.neighbours[vol_id, vertex_id] < 0:   
416                        #if not boundary.has_key( (vol_id, edge_id) ):
417                        if not boundary.has_key( (vol_id, vertex_id) ):   
418                            msg = 'WARNING: Given boundary does not contain '
419                            #msg += 'tags for edge (%d, %d). '\
420                            #       %(vol_id, edge_id)
421                            msg += 'tags for vertex (%d, %d). '\
422                                   %(vol_id, vertex_id)
423                            msg += 'Assigning default tag (%s).'\
424                                   %default_boundary_tag
425
426                            #FIXME: Print only as per verbosity
427                            #print msg
428
429                            #FIXME: Make this situation an error in the future
430                            #and make another function which will
431                            #enable default boundary-tags where
432                            #tags a not specified
433                            #boundary[ (vol_id, edge_id) ] =\
434                            boundary[ (vol_id, vertex_id) ] =\
435                                      default_boundary_tag
436
437
438
439        self.boundary = boundary
440
441    def build_tagged_elements_dictionary(self, tagged_elements = None):
442        """Build the dictionary of element tags.
443         self.tagged_elements is a dictionary of element arrays,
444         keyed by tag:
445         { (tag): [e1, e2, e3..] }
446
447         Postconditions:
448            self.element_tag is defined
449        """
450
451        if tagged_elements is None:
452            tagged_elements = {}
453        else:
454            #Check that all keys in given boundary exist
455            for tag in tagged_elements.keys():
456                tagged_elements[tag] = array(tagged_elements[tag]).astype(numpy.int)
457
458                msg = 'Not all elements exist. '
459                assert max(tagged_elements[tag]) < self.number_of_elements, msg
460        #print "tagged_elements", tagged_elements
461        self.tagged_elements = tagged_elements
462
463
464    def set_quantities_to_be_stored(self, q):
465        """Specify which quantities will be stored in the sww file.
466
467        q must be either:
468          - the name of a quantity
469          - a list of quantity names
470          - None
471
472        In the two first cases, the named quantities will be stored at each
473        yieldstep
474        (This is in addition to the quantities elevation and friction) 
475        If q is None, storage will be switched off altogether.
476        """
477
478
479        if q is None:
480            self.quantities_to_be_stored = []           
481            self.store = False
482            return
483
484        if isinstance(q, basestring):
485            q = [q] # Turn argument into a list
486
487        #Check correcness   
488        for quantity_name in q:
489            msg = 'Quantity %s is not a valid conserved quantity' %quantity_name
490            assert quantity_name in self.conserved_quantities, msg
491       
492        self.quantities_to_be_stored = q
493       
494
495
496
497
498    def get_boundary_tags(self):
499        """Return list of available boundary tags
500        """
501
502        tags = {}
503        for v in self.boundary.values():
504            tags[v] = 1
505
506        return tags.keys()
507       
508    def get_vertex_coordinates(self, obj = False):
509        """Return all vertex coordinates.
510        Return all vertex coordinates for all triangles as an Nx6 array
511        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)
512
513        if obj is True, the x/y pairs are returned in a 3*N x 2 array.
514        FIXME, we might make that the default.
515        FIXME Maybe use keyword: continuous = False for this condition?
516
517       
518        """
519
520        if obj is True:
521       
522            #V = self.vertex_coordinates
523            V = self.vertices
524            #return concatenate( (V[:,0:2], V[:,2:4], V[:,4:6]), axis=0)
525
526            N = V.shape[0]
527            #return reshape(V, (3*N, 2))
528            return numpy.reshape(V, (N, 2))
529        else:   
530            #return self.vertex_coordinates
531            return self.vertices
532   
533    def get_conserved_quantities(self, vol_id, vertex=None):#, edge=None):
534        """Get conserved quantities at volume vol_id
535
536        If vertex is specified use it as index for vertex values
537        If edge is specified use it as index for edge values
538        If neither are specified use centroid values
539        If both are specified an exeception is raised
540
541        Return value: Vector of length == number_of_conserved quantities
542
543        """
544
545        #if not (vertex is None):# or edge is None):
546        #    msg = 'Values for both vertex and edge was specified.'
547        #    msg += 'Only one (or none) is allowed.'
548       #     raise msg
549
550        q = numpy.zeros( len(self.conserved_quantities), numpy.float)
551
552        for i, name in enumerate(self.conserved_quantities):
553            Q = self.quantities[name]
554            if vertex is not None:
555                q[i] = Q.vertex_values[vol_id, vertex]
556            #elif edge is not None:
557            #    q[i] = Q.edge_values[vol_id, edge]
558            else:
559                q[i] = Q.centroid_values[vol_id]
560
561        return q
562
563
564    def get_evolved_quantities(self, vol_id, vertex=None):#, edge=None):
565        """Get evolved quantities at volume vol_id
566
567        If vertex is specified use it as index for vertex values
568        If edge is specified use it as index for edge values
569        If neither are specified use centroid values
570        If both are specified an exeception is raised
571
572        Return value: Vector of length == number_of_evolved quantities
573
574        """
575
576
577
578        #if not (vertex is None):# or edge is None):
579        #    msg = 'Values for both vertex and edge was specified.'
580        #    msg += 'Only one (or none) is allowed.'
581       #     raise msg
582
583        q = numpy.zeros( len(self.evolved_quantities), numpy.float)
584
585        for i, name in enumerate(self.evolved_quantities):
586            Q = self.quantities[name]
587            if vertex is not None:
588                q[i] = Q.vertex_values[vol_id, vertex]
589            #elif edge is not None:
590            #    q[i] = Q.edge_values[vol_id, edge]
591            else:
592                q[i] = Q.centroid_values[vol_id]
593
594        return q
595       
596
597    def get_centroids(self):
598        """Return all coordinates of centroids
599        Return x coordinate of centroid for each element as a N array
600        """
601
602        return self.centroids
603
604    def get_vertices(self):
605        """Return all coordinates of centroids
606        Return x coordinate of centroid for each element as a N array
607        """
608
609        return self.vertices
610
611    def get_coordinate(self, elem_id, vertex=None):
612        """Return coordinate of centroid,
613        or left or right vertex.
614        Left vertex (vertex=0). Right vertex (vertex=1)
615        """
616
617        if vertex is None:
618            return self.centroids[elem_id]
619        else:
620            return self.vertices[elem_id,vertex]
621
622    def get_area(self, elem_id=None):
623        """Return area of element id
624        """
625
626        if elem_id is None:
627            return sum(self.areas)
628        else:
629            return self.areas[elem_id]
630
631
632    def get_quantity(self, name, location='vertices', indices = None):
633        """Get values for named quantity
634
635        name: Name of quantity
636
637        In case of location == 'centroids' the dimension values must
638        be a list of a numpyal array of length N, N being the number
639        of elements. Otherwise it must be of dimension Nx3.
640
641        Indices is the set of element ids that the operation applies to.
642
643        The values will be stored in elements following their
644        internal ordering.
645        """
646
647        return self.quantities[name].get_values( location, indices = indices)
648
649    def get_centroid_coordinates(self):
650        """Return all centroid coordinates.
651        Return all centroid coordinates for all triangles as an Nx2 array
652        (ordered as x0, y0 for each triangle)
653        """
654        return self.centroids
655
656
657    def get_timestepping_method(self):
658        return self.timestepping_method
659
660    def set_timestepping_method(self,timestepping_method):
661       
662        if timestepping_method in ['euler', 'rk2', 'rk3']:
663            self.timestepping_method = timestepping_method
664            return
665
666        if timestepping_method in [1, 2, 3]:
667            self.timestepping_method = ['euler', 'rk2', 'rk3'][timestepping_method-1]
668            return
669
670        msg = '%s is an incorrect timestepping type'% timestepping_method
671        raise Exception, msg
672
673
674    def set_quantity(self, name, *args, **kwargs):
675        """Set values for named quantity
676
677
678        One keyword argument is documented here:
679        expression = None, # Arbitrary expression
680
681        expression:
682          Arbitrary expression involving quantity names
683
684        See Quantity.set_values for further documentation.
685        """
686
687        #FIXME (Ole): Allow new quantities here
688        #from quantity import Quantity, Conserved_quantity
689        #Create appropriate quantity object
690        # #if name in self.conserved_quantities:
691        # #    self.quantities[name] = Conserved_quantity(self)
692        # #else:
693        # #    self.quantities[name] = Quantity(self)
694
695
696        #Do the expression stuff
697        if kwargs.has_key('expression'):
698            expression = kwargs['expression']
699            del kwargs['expression']
700
701            Q = self.create_quantity_from_expression(expression)
702            kwargs['quantity'] = Q
703       
704        #Assign values
705        self.quantities[name].set_values(*args, **kwargs)
706
707    def set_boundary(self, boundary_map):
708        """Associate boundary objects with tagged boundary segments.
709
710        Input boundary_map is a dictionary of boundary objects keyed
711        by symbolic tags to matched against tags in the internal dictionary
712        self.boundary.
713
714        As result one pointer to a boundary object is stored for each vertex
715        in the list self.boundary_objects.
716        More entries may point to the same boundary object
717
718        Schematically the mapping is from two dictionaries to one list
719        where the index is used as pointer to the boundary_values arrays
720        within each quantity.
721
722        self.boundary:          (vol_id, edge_id): tag
723        boundary_map (input):   tag: boundary_object
724        ----------------------------------------------
725        self.boundary_objects:  ((vol_id, edge_id), boundary_object)
726
727
728        Pre-condition:
729          self.boundary has been built.
730
731        Post-condition:
732          self.boundary_objects is built
733
734        If a tag from the domain doesn't appear in the input dictionary an
735        exception is raised.
736        However, if a tag is not used to the domain, no error is thrown.
737        FIXME: This would lead to implementation of a
738        default boundary condition
739
740        Note: If a segment is listed in the boundary dictionary and if it is
741        not None, it *will* become a boundary -
742        even if there is a neighbouring triangle.
743        This would be the case for internal boundaries
744
745        Boundary objects that are None will be skipped.
746
747        FIXME: If set_boundary is called multiple times and if Boundary
748        object is changed into None, the neighbour structure will not be
749        restored!!!
750        """
751
752        self.boundary_objects = []
753        self.boundary_map = boundary_map  #Store for use with eg. boundary_stats.
754
755        #FIXME: Try to remove the sorting and fix test_mesh.py
756        x = self.boundary.keys()
757        x.sort()
758
759        #Loop through edges that lie on the boundary and associate them with
760        #callable boundary objects depending on their tags
761        #for k, (vol_id, edge_id) in enumerate(x):
762        for k, (vol_id, vertex_id) in enumerate(x):
763            #tag = self.boundary[ (vol_id, edge_id) ]
764            tag = self.boundary[ (vol_id, vertex_id) ]
765
766            if boundary_map.has_key(tag):
767                B = boundary_map[tag]  #Get callable boundary object
768
769                if B is not None:
770                    #self.boundary_objects.append( ((vol_id, edge_id), B) )
771                    #self.neighbours[vol_id, edge_id] = -len(self.boundary_objects)
772                    self.boundary_objects.append( ((vol_id, vertex_id), B) )
773                    self.neighbours[vol_id, vertex_id] = -len(self.boundary_objects)
774                else:
775                    pass
776                    #FIXME: Check and perhaps fix neighbour structure
777
778            else:
779                msg = 'ERROR (domain.py): Tag "%s" has not been ' %tag
780                msg += 'bound to a boundary object.\n'
781                msg += 'All boundary tags defined in domain must appear '
782                msg += 'in the supplied dictionary.\n'
783                msg += 'The tags are: %s' %self.get_boundary_tags()
784                raise Exception, msg
785
786
787
788    def check_integrity(self):
789        #Mesh.check_integrity(self)
790
791        #print self.quantities
792        #print self.conserved_quantities
793       
794        for quantity in self.conserved_quantities:
795            msg = 'Conserved quantities must be a subset of all quantities'
796            assert quantity in self.quantities, msg
797
798        for quantity in self.evolved_quantities:
799            msg = 'Evolved quantities must be a subset of all quantities'
800            assert quantity in self.quantities, msg           
801
802        # #assert hasattr(self, 'boundary_objects')
803
804    def write_time(self):
805        print self.timestepping_statistics()
806
807    def timestepping_statistics(self):
808        """Return string with time stepping statistics for printing or logging
809        """
810
811        msg = ''
812        if self.min_timestep == self.max_timestep:
813            msg += 'Time = %.4f, delta t = %.8f, steps=%d (%d)'\
814                   %(self.time, self.min_timestep, self.number_of_steps,
815                     self.number_of_first_order_steps)
816        elif self.min_timestep > self.max_timestep:
817            msg += 'Time = %.4f, steps=%d (%d)'\
818                   %(self.time, self.number_of_steps,
819                     self.number_of_first_order_steps)
820        else:
821            msg += 'Time = %.4f, delta t in [%.8f, %.8f], steps=%d (%d)'\
822                   %(self.time, self.min_timestep,
823                     self.max_timestep, self.number_of_steps,
824                     self.number_of_first_order_steps)
825
826        return msg
827
828    def get_name(self):
829        return self.filename
830
831    def set_name(self, name):
832        self.filename = name
833
834    def get_datadir(self):
835        return self.datadir
836
837    def set_datadir(self, name):
838        self.datadir = name
839
840    def set_CFL(self, cfl):
841        if cfl > 1.0:
842            print 'WARNING: Setting CFL condition to %f which is greater than 1' % cfl
843        self.CFL = cfl
844
845    def get_CFL(self):
846        return self.CFL
847   
848    def set_filename(self, name):
849        self.filename = name
850
851    def get_filename(self):
852        return self.filename
853
854
855    def get_spatial_order(self):
856        return self.order
857
858    def set_spatial_order(self, order):
859
860        possible_orders = [1,2]
861
862        if order in possible_orders:
863            self.order = order
864            return
865
866        msg = '%s is an incorrect spatial order.\n'% limiter
867        msg += 'Possible orders are: '+ ", ".join(["%s" % el for el in possible_orders])
868        raise Exception, msg
869   
870
871    def get_limiter(self):
872        return self.limiter
873
874    def set_limiter(self,limiter):
875
876        possible_limiters = \
877        ['pyvolution', 'minmod_steve', 'minmod', 'minmod_kurganov', 'superbee', 'vanleer', 'vanalbada']
878
879        if limiter in possible_limiters:
880            self.limiter = limiter
881            return
882
883        msg = '%s is an incorrect limiter type.\n'% limiter
884        msg += 'Possible types are: '+ ", ".join(["%s" % el for el in possible_limiters])
885        raise Exception, msg
886
887       
888    #--------------------------
889    # Main components of evolve
890    #--------------------------   
891
892    def evolve(self, yieldstep = None,
893               finaltime = None,
894               duration = None,
895               skip_initial_step = False):
896        """Evolve model through time starting from self.starttime.
897
898
899        yieldstep: numpy.interval between yields where results are stored,
900                   statistics written and domain inspected or
901                   possibly modified. If omitted the internal predefined
902                   max timestep is used.
903                   numpy.internally, smaller timesteps may be taken.
904
905        duration: Duration of simulation
906
907        finaltime: Time where simulation should end. This is currently
908        relative time.  So it's the same as duration.
909
910        If both duration and finaltime are given an exception is thrown.
911
912
913        skip_initial_step: Boolean flag that decides whether the first
914        yield step is skipped or not. This is useful for example to avoid
915        duplicate steps when multiple evolve processes are dove tailed.
916
917
918        Evolve is implemented as a generator and is to be called as such, e.g.
919
920        for t in domain.evolve(yieldstep, finaltime):
921            <Do something with domain and t>
922
923
924        All times are given in seconds
925
926        """
927
928        from config import min_timestep, max_timestep, epsilon
929
930        # FIXME: Maybe lump into a larger check prior to evolving
931        msg = 'Boundary tags must be bound to boundary objects before '
932        msg += 'evolving system, '
933        msg += 'e.g. using the method set_boundary.\n'
934        msg += 'This system has the boundary tags %s '\
935               %self.get_boundary_tags()
936        assert hasattr(self, 'boundary_objects'), msg
937
938
939        if yieldstep is None:
940            yieldstep = max_timestep
941        else:
942            yieldstep = numpy.float(yieldstep)
943
944        self._order_ = self.default_order
945
946
947        if finaltime is not None and duration is not None:
948            # print 'F', finaltime, duration
949            msg = 'Only one of finaltime and duration may be specified'
950            raise msg
951        else:
952            if finaltime is not None:
953                self.finaltime = numpy.float(finaltime)
954            if duration is not None:
955                self.finaltime = self.starttime + numpy.float(duration)
956
957
958
959        N = len(self) # Number of triangles
960        self.yieldtime = 0.0 # Track time between 'yields'
961
962        # Initialise interval of timestep sizes (for reporting only)
963        self.min_timestep = max_timestep
964        self.max_timestep = min_timestep
965        self.number_of_steps = 0
966        self.number_of_first_order_steps = 0
967
968
969        # Update ghosts
970        self.update_ghosts()
971
972        # Initial update of vertex and edge values
973        self.distribute_to_vertices_and_edges()
974
975        # Update extrema if necessary (for reporting)
976        self.update_extrema()
977       
978        # Initial update boundary values
979        self.update_boundary()
980
981        # Or maybe restore from latest checkpoint
982        if self.checkpoint is True:
983            self.goto_latest_checkpoint()
984
985        if skip_initial_step is False:
986            yield(self.time)  # Yield initial values
987
988        while True:
989
990            # Evolve One Step, using appropriate timestepping method
991            if self.get_timestepping_method() == 'euler':
992                self.evolve_one_euler_step(yieldstep,finaltime)
993               
994            elif self.get_timestepping_method() == 'rk2':
995                self.evolve_one_rk2_step(yieldstep,finaltime)
996
997            elif self.get_timestepping_method() == 'rk3':
998                self.evolve_one_rk3_step(yieldstep,finaltime)               
999
1000           
1001            # Update extrema if necessary (for reporting)
1002            self.update_extrema()           
1003
1004
1005           
1006            self.yieldtime += self.timestep
1007            self.number_of_steps += 1
1008            if self._order_ == 1:
1009                self.number_of_first_order_steps += 1
1010
1011
1012            # Yield results
1013            if finaltime is not None and self.time >= finaltime-epsilon:
1014
1015                if self.time > finaltime:
1016                    # FIXME (Ole, 30 April 2006): Do we need this check?
1017                    # Probably not (Ole, 18 September 2008). Now changed to
1018                    # Exception
1019                    msg = 'WARNING (domain.py): time overshot finaltime. '
1020                    msg += 'Contact Ole.Nielsen@ga.gov.au'
1021                    raise Exception, msg
1022                   
1023
1024                # Yield final time and stop
1025                self.time = finaltime
1026                yield(self.time)
1027                break
1028
1029            if self.yieldtime >= yieldstep:
1030                # Yield (intermediate) time and allow inspection of domain
1031
1032                if self.checkpoint is True:
1033                    self.store_checkpoint()
1034                    self.delete_old_checkpoints()
1035
1036                # Pass control on to outer loop for more specific actions
1037
1038                yield(self.time)
1039
1040                # Reinitialise
1041                self.yieldtime = 0.0
1042                self.min_timestep = max_timestep
1043                self.max_timestep = min_timestep
1044                self.number_of_steps = 0
1045                self.number_of_first_order_steps = 0
1046                #self.max_speed_array = 0.0
1047
1048
1049    def evolve_one_euler_step(self, yieldstep, finaltime):
1050        """
1051        One Euler Time Step
1052        Q^{n+1} = E(h) Q^n
1053        """
1054
1055
1056        # Compute fluxes across each element edge
1057        self.compute_fluxes()
1058
1059        # Update timestep to fit yieldstep and finaltime
1060        self.update_timestep(yieldstep, finaltime)
1061
1062        # Update conserved quantities
1063        self.update_conserved_quantities()
1064
1065        # Update ghosts
1066        self.update_ghosts()
1067
1068        # Update vertex and edge values
1069        self.distribute_to_vertices_and_edges()
1070
1071        # Update boundary values
1072        self.update_boundary()
1073
1074        # Update time
1075        self.time += self.timestep
1076
1077       
1078
1079
1080    def evolve_one_rk2_step(self, yieldstep, finaltime):
1081        """
1082        One 2nd order RK timestep
1083        Q^{n+1} = 0.5 Q^n + 0.5 E(h)^2 Q^n
1084        """
1085
1086        # Save initial conserved quantities values
1087        self.backup_conserved_quantities()           
1088
1089        #--------------------------------------
1090        # First euler step
1091        #--------------------------------------
1092
1093        # Compute fluxes across each element edge
1094        self.compute_fluxes()
1095
1096        # Update timestep to fit yieldstep and finaltime
1097        self.update_timestep(yieldstep, finaltime)
1098
1099        # Update conserved quantities
1100        self.update_conserved_quantities()
1101
1102        # Update ghosts
1103        self.update_ghosts()
1104
1105        # Update vertex and edge values
1106        self.distribute_to_vertices_and_edges()
1107
1108        # Update boundary values
1109        self.update_boundary()
1110
1111        # Update time
1112        self.time += self.timestep
1113
1114        #------------------------------------
1115        # Second Euler step
1116        #------------------------------------
1117           
1118        # Compute fluxes across each element edge
1119        self.compute_fluxes()
1120
1121        # Update conserved quantities
1122        self.update_conserved_quantities()
1123
1124        #------------------------------------
1125        # Combine initial and final values
1126        # of conserved quantities and cleanup
1127        #------------------------------------
1128       
1129        # Combine steps
1130        self.saxpy_conserved_quantities(0.5, 0.5)
1131
1132        #-----------------------------------
1133        # clean up vertex values
1134        #-----------------------------------
1135 
1136        # Update ghosts
1137        self.update_ghosts()
1138
1139        # Update vertex and edge values
1140        self.distribute_to_vertices_and_edges()
1141
1142        # Update boundary values
1143        self.update_boundary()
1144
1145
1146
1147    def evolve_one_rk3_step(self, yieldstep, finaltime):
1148        """
1149        One 3rd order RK timestep
1150        Q^(1) = 3/4 Q^n + 1/4 E(h)^2 Q^n  (at time t^n + h/2)
1151        Q^{n+1} = 1/3 Q^n + 2/3 E(h) Q^(1) (at time t^{n+1})
1152        """
1153
1154        # Save initial initial conserved quantities values
1155        self.backup_conserved_quantities()           
1156
1157        initial_time = self.time
1158       
1159        #--------------------------------------
1160        # First euler step
1161        #--------------------------------------
1162
1163        # Compute fluxes across each element edge
1164        self.compute_fluxes()
1165
1166        # Update timestep to fit yieldstep and finaltime
1167        self.update_timestep(yieldstep, finaltime)
1168
1169        # Update conserved quantities
1170        self.update_conserved_quantities()
1171
1172        # Update ghosts
1173        self.update_ghosts()
1174
1175        # Update vertex and edge values
1176        self.distribute_to_vertices_and_edges()
1177
1178        # Update boundary values
1179        self.update_boundary()
1180
1181        # Update time
1182        self.time += self.timestep
1183
1184        #------------------------------------
1185        # Second Euler step
1186        #------------------------------------
1187           
1188        # Compute fluxes across each element edge
1189        self.compute_fluxes()
1190
1191        # Update conserved quantities
1192        self.update_conserved_quantities()
1193
1194        #------------------------------------
1195        #Combine steps to obtain intermediate
1196        #solution at time t^n + 0.5 h
1197        #------------------------------------
1198
1199        # Combine steps
1200        self.saxpy_conserved_quantities(0.25, 0.75)
1201 
1202        # Update ghosts
1203        self.update_ghosts()
1204
1205        # Update vertex and edge values
1206        self.distribute_to_vertices_and_edges()
1207
1208        # Update boundary values
1209        self.update_boundary()
1210
1211        # Set substep time
1212        self.time = initial_time + self.timestep*0.5
1213
1214        #------------------------------------
1215        # Third Euler step
1216        #------------------------------------
1217           
1218        # Compute fluxes across each element edge
1219        self.compute_fluxes()
1220
1221        # Update conserved quantities
1222        self.update_conserved_quantities()
1223
1224        #------------------------------------
1225        # Combine final and initial values
1226        # and cleanup
1227        #------------------------------------
1228       
1229        # Combine steps
1230        self.saxpy_conserved_quantities(2.0/3.0, 1.0/3.0)
1231 
1232        # Update ghosts
1233        self.update_ghosts()
1234
1235        # Update vertex and edge values
1236        self.distribute_to_vertices_and_edges()
1237
1238        # Update boundary values
1239        self.update_boundary()
1240
1241        # Set new time
1242        self.time = initial_time + self.timestep       
1243       
1244
1245    def backup_conserved_quantities(self):
1246        N = len(self) # Number_of_triangles
1247
1248        # Backup conserved_quantities centroid values
1249        for name in self.conserved_quantities:
1250            Q = self.quantities[name]
1251            Q.backup_centroid_values()       
1252
1253    def saxpy_conserved_quantities(self,a,b):
1254        N = len(self) #number_of_triangles
1255
1256        # Backup conserved_quantities centroid values
1257        for name in self.conserved_quantities:
1258            Q = self.quantities[name]
1259            Q.saxpy_centroid_values(a,b)       
1260
1261
1262
1263       
1264    def distribute_to_vertices_and_edges(self):
1265        """Extrapolate conserved quantities from centroid to
1266        vertices and edge-midpoints for each volume
1267
1268        Default implementation is straight first order,
1269        i.e. constant values throughout each element and
1270        no reference to non-conserved quantities.
1271        """
1272
1273        for name in self.conserved_quantities:
1274            Q = self.quantities[name]
1275            if self.order == 1:
1276                Q.extrapolate_first_order()
1277            elif self.order == 2:
1278                Q.extrapolate_second_order()
1279                #Q.limit()
1280            else:
1281                raise 'Unknown order'
1282            #Q.interpolate_from_vertices_to_edges()
1283
1284
1285    def update_ghosts(self):
1286        pass
1287   
1288    def update_boundary(self):
1289        """Go through list of boundary objects and update boundary values
1290        for all conserved quantities on boundary.
1291        """
1292
1293        #FIXME: Update only those that change (if that can be worked out)
1294        #FIXME: Boundary objects should not include ghost nodes.
1295        #for i, ((vol_id, edge_id), B) in enumerate(self.boundary_objects):
1296        #    q = B.evaluate(vol_id, edge_id)
1297        for i, ((vol_id, vertex_id), B) in enumerate(self.boundary_objects):
1298            q = B.evaluate(vol_id, vertex_id)
1299            #print 'q ',q
1300            for j, name in enumerate(self.evolved_quantities):
1301                #print 'name %s j = %f \n'%(name,j)
1302                Q = self.quantities[name]
1303
1304                Q.boundary_values[i] = q[j]
1305                #print 'Q=',Q
1306   
1307    def update_timestep(self, yieldstep, finaltime):
1308
1309        from config import min_timestep, max_timestep
1310
1311        # self.timestep is calculated from speed of characteristics
1312        # Apply CFL condition here
1313        timestep = min(self.CFL*self.flux_timestep, max_timestep)
1314
1315        #Record maximal and minimal values of timestep for reporting
1316        self.max_timestep = max(timestep, self.max_timestep)
1317        self.min_timestep = min(timestep, self.min_timestep)
1318
1319        #Protect against degenerate time steps
1320        if timestep < min_timestep:
1321
1322            #Number of consecutive small steps taken b4 taking action
1323            self.smallsteps += 1
1324
1325            if self.smallsteps > self.max_smallsteps:
1326                self.smallsteps = 0 #Reset
1327
1328                if self.order == 1:
1329                    msg = 'WARNING: Too small timestep %.16f reached '\
1330                          %timestep
1331                    msg += 'even after %d steps of 1 order scheme'\
1332                           %self.max_smallsteps
1333                    print msg
1334                    timestep = min_timestep  #Try enforcing min_step
1335
1336                    #raise msg
1337                else:
1338                    #Try to overcome situation by switching to 1 order
1339                    print "changing Order 1"
1340                    self.order = 1
1341
1342        else:
1343            self.smallsteps = 0
1344            if self.order == 1 and self.default_order == 2:
1345                self.order = 2
1346
1347
1348        #Ensure that final time is not exceeded
1349        if finaltime is not None and self.time + timestep > finaltime:
1350            timestep = finaltime-self.time
1351
1352        #Ensure that model time is aligned with yieldsteps
1353        if self.yieldtime + timestep > yieldstep:
1354            timestep = yieldstep-self.yieldtime
1355
1356        self.timestep = timestep
1357
1358    def update_extrema(self):
1359        pass
1360
1361    def compute_forcing_terms(self):
1362        """If there are any forcing functions driving the system
1363        they should be defined in Domain subclass and appended to
1364        the list self.forcing_terms
1365        """
1366        #Clears explicit_update needed for second order method
1367        if self.time_order == 2:
1368            for name in self.conserved_quantities:
1369                Q = self.quantities[name]
1370                Q.explicit_update[:] = 0.0
1371
1372        for f in self.forcing_terms:
1373            f(self)
1374           
1375
1376    def update_derived_quantites(self):
1377        pass
1378   
1379    #def update_conserved_quantities(self):
1380    def update_conserved_quantities(self):
1381        """Update vectors of conserved quantities using previously
1382        computed fluxes specified forcing functions.
1383        """
1384
1385
1386        timestep = self.timestep
1387
1388        #Compute forcing terms
1389        self.compute_forcing_terms()
1390
1391        #Update conserved_quantities
1392        for name in self.conserved_quantities:
1393            Q = self.quantities[name]
1394            Q.update(timestep)
1395
1396
1397
1398if __name__ == "__main__":
1399
1400    points1 = [0.0, 1.0, 2.0, 3.0]
1401    D1 = Domain(points1)
1402
1403    print D1.get_coordinate(0)
1404    print D1.get_coordinate(0,1)
1405    print 'Number of Elements = ',D1.number_of_elements
1406
1407    try:
1408        print D1.get_coordinate(3)
1409    except:
1410        pass
1411    else:
1412        msg =  'Should have raised an out of bounds exception'
1413        raise msg
1414
1415    #points2 = [0.0, 1.0, 2.0, 3.0, 2.5]
1416    #D2 = Domain(points2)
Note: See TracBrowser for help on using the repository browser.