source: anuga_work/development/pipeflow/generic_domain.py @ 7823

Last change on this file since 7823 was 7823, checked in by steve, 14 years ago

Adding old files as well as workinng files from anuga_1d

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