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

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

Adding in a few new files

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