source: anuga_work/development/pipeflow/domain.py @ 7777

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

Copying 1d code over to pipeflow directory

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 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):
590        """Return area of element id
591        """
592
593        return self.areas[elem_id]
594
595    def get_quantity(self, name, location='vertices', indices = None):
596        """Get values for named quantity
597
598        name: Name of quantity
599
600        In case of location == 'centroids' the dimension values must
601        be a list of a numpyal array of length N, N being the number
602        of elements. Otherwise it must be of dimension Nx3.
603
604        Indices is the set of element ids that the operation applies to.
605
606        The values will be stored in elements following their
607        internal ordering.
608        """
609
610        return self.quantities[name].get_values( location, indices = indices)
611
612    def get_centroid_coordinates(self):
613        """Return all centroid coordinates.
614        Return all centroid coordinates for all triangles as an Nx2 array
615        (ordered as x0, y0 for each triangle)
616        """
617        return self.centroids
618
619
620    def get_timestepping_method(self):
621        return self.timestepping_method
622
623    def set_timestepping_method(self,timestepping_method):
624       
625        if timestepping_method in ['euler', 'rk2', 'rk3']:
626            self.timestepping_method = timestepping_method
627            return
628
629        msg = '%s is an incorrect timestepping type'% timestepping_method
630        raise Exception, msg
631
632
633    def set_quantity(self, name, *args, **kwargs):
634        """Set values for named quantity
635
636
637        One keyword argument is documented here:
638        expression = None, # Arbitrary expression
639
640        expression:
641          Arbitrary expression involving quantity names
642
643        See Quantity.set_values for further documentation.
644        """
645
646        #FIXME (Ole): Allow new quantities here
647        #from quantity import Quantity, Conserved_quantity
648        #Create appropriate quantity object
649        # #if name in self.conserved_quantities:
650        # #    self.quantities[name] = Conserved_quantity(self)
651        # #else:
652        # #    self.quantities[name] = Quantity(self)
653
654
655        #Do the expression stuff
656        if kwargs.has_key('expression'):
657            expression = kwargs['expression']
658            del kwargs['expression']
659
660            Q = self.create_quantity_from_expression(expression)
661            kwargs['quantity'] = Q
662       
663        #Assign values
664        self.quantities[name].set_values(*args, **kwargs)
665
666    def set_boundary(self, boundary_map):
667        """Associate boundary objects with tagged boundary segments.
668
669        Input boundary_map is a dictionary of boundary objects keyed
670        by symbolic tags to matched against tags in the internal dictionary
671        self.boundary.
672
673        As result one pointer to a boundary object is stored for each vertex
674        in the list self.boundary_objects.
675        More entries may point to the same boundary object
676
677        Schematically the mapping is from two dictionaries to one list
678        where the index is used as pointer to the boundary_values arrays
679        within each quantity.
680
681        self.boundary:          (vol_id, edge_id): tag
682        boundary_map (input):   tag: boundary_object
683        ----------------------------------------------
684        self.boundary_objects:  ((vol_id, edge_id), boundary_object)
685
686
687        Pre-condition:
688          self.boundary has been built.
689
690        Post-condition:
691          self.boundary_objects is built
692
693        If a tag from the domain doesn't appear in the input dictionary an
694        exception is raised.
695        However, if a tag is not used to the domain, no error is thrown.
696        FIXME: This would lead to implementation of a
697        default boundary condition
698
699        Note: If a segment is listed in the boundary dictionary and if it is
700        not None, it *will* become a boundary -
701        even if there is a neighbouring triangle.
702        This would be the case for internal boundaries
703
704        Boundary objects that are None will be skipped.
705
706        FIXME: If set_boundary is called multiple times and if Boundary
707        object is changed into None, the neighbour structure will not be
708        restored!!!
709        """
710
711        self.boundary_objects = []
712        self.boundary_map = boundary_map  #Store for use with eg. boundary_stats.
713
714        #FIXME: Try to remove the sorting and fix test_mesh.py
715        x = self.boundary.keys()
716        x.sort()
717
718        #Loop through edges that lie on the boundary and associate them with
719        #callable boundary objects depending on their tags
720        #for k, (vol_id, edge_id) in enumerate(x):
721        for k, (vol_id, vertex_id) in enumerate(x):
722            #tag = self.boundary[ (vol_id, edge_id) ]
723            tag = self.boundary[ (vol_id, vertex_id) ]
724
725            if boundary_map.has_key(tag):
726                B = boundary_map[tag]  #Get callable boundary object
727
728                if B is not None:
729                    #self.boundary_objects.append( ((vol_id, edge_id), B) )
730                    #self.neighbours[vol_id, edge_id] = -len(self.boundary_objects)
731                    self.boundary_objects.append( ((vol_id, vertex_id), B) )
732                    self.neighbours[vol_id, vertex_id] = -len(self.boundary_objects)
733                else:
734                    pass
735                    #FIXME: Check and perhaps fix neighbour structure
736
737            else:
738                msg = 'ERROR (domain.py): Tag "%s" has not been ' %tag
739                msg += 'bound to a boundary object.\n'
740                msg += 'All boundary tags defined in domain must appear '
741                msg += 'in the supplied dictionary.\n'
742                msg += 'The tags are: %s' %self.get_boundary_tags()
743                raise msg
744
745
746
747    def check_integrity(self):
748        #Mesh.check_integrity(self)
749
750        #print self.quantities
751        #print self.conserved_quantities
752       
753        for quantity in self.conserved_quantities:
754            msg = 'Conserved quantities must be a subset of all quantities'
755            assert quantity in self.quantities, msg
756
757        for quantity in self.evolved_quantities:
758            msg = 'Evolved quantities must be a subset of all quantities'
759            assert quantity in self.quantities, msg           
760
761        # #assert hasattr(self, 'boundary_objects')
762
763    def write_time(self):
764        print self.timestepping_statistics()
765
766    def timestepping_statistics(self):
767        """Return string with time stepping statistics for printing or logging
768        """
769
770        msg = ''
771        if self.min_timestep == self.max_timestep:
772            msg += 'Time = %.4f, delta t = %.8f, steps=%d (%d)'\
773                   %(self.time, self.min_timestep, self.number_of_steps,
774                     self.number_of_first_order_steps)
775        elif self.min_timestep > self.max_timestep:
776            msg += 'Time = %.4f, steps=%d (%d)'\
777                   %(self.time, self.number_of_steps,
778                     self.number_of_first_order_steps)
779        else:
780            msg += 'Time = %.4f, delta t in [%.8f, %.8f], steps=%d (%d)'\
781                   %(self.time, self.min_timestep,
782                     self.max_timestep, self.number_of_steps,
783                     self.number_of_first_order_steps)
784
785        return msg
786
787    def get_name(self):
788        return self.filename
789
790    def set_name(self, name):
791        self.filename = name
792
793    def get_datadir(self):
794        return self.datadir
795
796    def set_datadir(self, name):
797        self.datadir = name
798
799    def set_CFL(self, cfl):
800        if cfl > 1.0:
801            print 'WARNING: Setting CFL condition to %f which is greater than 1' % cfl
802        self.CFL = cfl
803
804    def get_CFL(self):
805        return self.CFL
806   
807    def set_filename(self, name):
808        self.filename = name
809
810    def get_filename(self):
811        return self.filename
812
813    def get_limiter(self):
814        return self.limiter
815
816    def set_limiter(self,limiter):
817
818        possible_limiters = \
819        ['pyvolution', 'minmod_steve', 'minmod', 'minmod_kurganov', 'superbee', 'vanleer', 'vanalbada']
820
821        if limiter in possible_limiters:
822            self.limiter = limiter
823            return
824
825        msg = '%s is an incorrect limiter type.\n'% limiter
826        msg += 'Possible types are: '+ ", ".join(["%s" % el for el in possible_limiters])
827        raise Exception, msg
828
829       
830    #--------------------------
831    # Main components of evolve
832    #--------------------------   
833
834    def evolve(self, yieldstep = None,
835               finaltime = None,
836               duration = None,
837               skip_initial_step = False):
838        """Evolve model through time starting from self.starttime.
839
840
841        yieldstep: numpy.interval between yields where results are stored,
842                   statistics written and domain inspected or
843                   possibly modified. If omitted the internal predefined
844                   max timestep is used.
845                   numpy.internally, smaller timesteps may be taken.
846
847        duration: Duration of simulation
848
849        finaltime: Time where simulation should end. This is currently
850        relative time.  So it's the same as duration.
851
852        If both duration and finaltime are given an exception is thrown.
853
854
855        skip_initial_step: Boolean flag that decides whether the first
856        yield step is skipped or not. This is useful for example to avoid
857        duplicate steps when multiple evolve processes are dove tailed.
858
859
860        Evolve is implemented as a generator and is to be called as such, e.g.
861
862        for t in domain.evolve(yieldstep, finaltime):
863            <Do something with domain and t>
864
865
866        All times are given in seconds
867
868        """
869
870        from config import min_timestep, max_timestep, epsilon
871
872        # FIXME: Maybe lump into a larger check prior to evolving
873        msg = 'Boundary tags must be bound to boundary objects before '
874        msg += 'evolving system, '
875        msg += 'e.g. using the method set_boundary.\n'
876        msg += 'This system has the boundary tags %s '\
877               %self.get_boundary_tags()
878        assert hasattr(self, 'boundary_objects'), msg
879
880
881        if yieldstep is None:
882            yieldstep = max_timestep
883        else:
884            yieldstep = numpy.float(yieldstep)
885
886        self._order_ = self.default_order
887
888
889        if finaltime is not None and duration is not None:
890            # print 'F', finaltime, duration
891            msg = 'Only one of finaltime and duration may be specified'
892            raise msg
893        else:
894            if finaltime is not None:
895                self.finaltime = numpy.float(finaltime)
896            if duration is not None:
897                self.finaltime = self.starttime + numpy.float(duration)
898
899
900
901        N = len(self) # Number of triangles
902        self.yieldtime = 0.0 # Track time between 'yields'
903
904        # Initialise interval of timestep sizes (for reporting only)
905        self.min_timestep = max_timestep
906        self.max_timestep = min_timestep
907        self.number_of_steps = 0
908        self.number_of_first_order_steps = 0
909
910
911        # Update ghosts
912        self.update_ghosts()
913
914        # Initial update of vertex and edge values
915        self.distribute_to_vertices_and_edges()
916
917        # Update extrema if necessary (for reporting)
918        self.update_extrema()
919       
920        # Initial update boundary values
921        self.update_boundary()
922
923        # Or maybe restore from latest checkpoint
924        if self.checkpoint is True:
925            self.goto_latest_checkpoint()
926
927        if skip_initial_step is False:
928            yield(self.time)  # Yield initial values
929
930        while True:
931
932            # Evolve One Step, using appropriate timestepping method
933            if self.get_timestepping_method() == 'euler':
934                self.evolve_one_euler_step(yieldstep,finaltime)
935               
936            elif self.get_timestepping_method() == 'rk2':
937                self.evolve_one_rk2_step(yieldstep,finaltime)
938
939            elif self.get_timestepping_method() == 'rk3':
940                self.evolve_one_rk3_step(yieldstep,finaltime)               
941
942           
943            # Update extrema if necessary (for reporting)
944            self.update_extrema()           
945
946
947           
948            self.yieldtime += self.timestep
949            self.number_of_steps += 1
950            if self._order_ == 1:
951                self.number_of_first_order_steps += 1
952
953
954            # Yield results
955            if finaltime is not None and self.time >= finaltime-epsilon:
956
957                if self.time > finaltime:
958                    # FIXME (Ole, 30 April 2006): Do we need this check?
959                    # Probably not (Ole, 18 September 2008). Now changed to
960                    # Exception
961                    msg = 'WARNING (domain.py): time overshot finaltime. '
962                    msg += 'Contact Ole.Nielsen@ga.gov.au'
963                    raise Exception, msg
964                   
965
966                # Yield final time and stop
967                self.time = finaltime
968                yield(self.time)
969                break
970
971            if self.yieldtime >= yieldstep:
972                # Yield (intermediate) time and allow inspection of domain
973
974                if self.checkpoint is True:
975                    self.store_checkpoint()
976                    self.delete_old_checkpoints()
977
978                # Pass control on to outer loop for more specific actions
979
980                yield(self.time)
981
982                # Reinitialise
983                self.yieldtime = 0.0
984                self.min_timestep = max_timestep
985                self.max_timestep = min_timestep
986                self.number_of_steps = 0
987                self.number_of_first_order_steps = 0
988                #self.max_speed_array = 0.0
989
990
991    def evolve_one_euler_step(self, yieldstep, finaltime):
992        """
993        One Euler Time Step
994        Q^{n+1} = E(h) Q^n
995        """
996
997
998        # Compute fluxes across each element edge
999        self.compute_fluxes()
1000
1001        # Update timestep to fit yieldstep and finaltime
1002        self.update_timestep(yieldstep, finaltime)
1003
1004        # Update conserved quantities
1005        self.update_conserved_quantities()
1006
1007        # Update ghosts
1008        self.update_ghosts()
1009
1010        # Update vertex and edge values
1011        self.distribute_to_vertices_and_edges()
1012
1013        # Update boundary values
1014        self.update_boundary()
1015
1016        # Update time
1017        self.time += self.timestep
1018
1019       
1020
1021
1022    def evolve_one_rk2_step(self, yieldstep, finaltime):
1023        """
1024        One 2nd order RK timestep
1025        Q^{n+1} = 0.5 Q^n + 0.5 E(h)^2 Q^n
1026        """
1027
1028        # Save initial conserved quantities values
1029        self.backup_conserved_quantities()           
1030
1031        #--------------------------------------
1032        # First euler step
1033        #--------------------------------------
1034
1035        # Compute fluxes across each element edge
1036        self.compute_fluxes()
1037
1038        # Update timestep to fit yieldstep and finaltime
1039        self.update_timestep(yieldstep, finaltime)
1040
1041        # Update conserved quantities
1042        self.update_conserved_quantities()
1043
1044        # Update ghosts
1045        self.update_ghosts()
1046
1047        # Update vertex and edge values
1048        self.distribute_to_vertices_and_edges()
1049
1050        # Update boundary values
1051        self.update_boundary()
1052
1053        # Update time
1054        self.time += self.timestep
1055
1056        #------------------------------------
1057        # Second Euler step
1058        #------------------------------------
1059           
1060        # Compute fluxes across each element edge
1061        self.compute_fluxes()
1062
1063        # Update conserved quantities
1064        self.update_conserved_quantities()
1065
1066        #------------------------------------
1067        # Combine initial and final values
1068        # of conserved quantities and cleanup
1069        #------------------------------------
1070       
1071        # Combine steps
1072        self.saxpy_conserved_quantities(0.5, 0.5)
1073
1074        #-----------------------------------
1075        # clean up vertex values
1076        #-----------------------------------
1077 
1078        # Update ghosts
1079        self.update_ghosts()
1080
1081        # Update vertex and edge values
1082        self.distribute_to_vertices_and_edges()
1083
1084        # Update boundary values
1085        self.update_boundary()
1086
1087
1088
1089    def evolve_one_rk3_step(self, yieldstep, finaltime):
1090        """
1091        One 3rd order RK timestep
1092        Q^(1) = 3/4 Q^n + 1/4 E(h)^2 Q^n  (at time t^n + h/2)
1093        Q^{n+1} = 1/3 Q^n + 2/3 E(h) Q^(1) (at time t^{n+1})
1094        """
1095
1096        # Save initial initial conserved quantities values
1097        self.backup_conserved_quantities()           
1098
1099        initial_time = self.time
1100       
1101        #--------------------------------------
1102        # First euler step
1103        #--------------------------------------
1104
1105        # Compute fluxes across each element edge
1106        self.compute_fluxes()
1107
1108        # Update timestep to fit yieldstep and finaltime
1109        self.update_timestep(yieldstep, finaltime)
1110
1111        # Update conserved quantities
1112        self.update_conserved_quantities()
1113
1114        # Update ghosts
1115        self.update_ghosts()
1116
1117        # Update vertex and edge values
1118        self.distribute_to_vertices_and_edges()
1119
1120        # Update boundary values
1121        self.update_boundary()
1122
1123        # Update time
1124        self.time += self.timestep
1125
1126        #------------------------------------
1127        # Second Euler step
1128        #------------------------------------
1129           
1130        # Compute fluxes across each element edge
1131        self.compute_fluxes()
1132
1133        # Update conserved quantities
1134        self.update_conserved_quantities()
1135
1136        #------------------------------------
1137        #Combine steps to obtain intermediate
1138        #solution at time t^n + 0.5 h
1139        #------------------------------------
1140
1141        # Combine steps
1142        self.saxpy_conserved_quantities(0.25, 0.75)
1143 
1144        # Update ghosts
1145        self.update_ghosts()
1146
1147        # Update vertex and edge values
1148        self.distribute_to_vertices_and_edges()
1149
1150        # Update boundary values
1151        self.update_boundary()
1152
1153        # Set substep time
1154        self.time = initial_time + self.timestep*0.5
1155
1156        #------------------------------------
1157        # Third Euler step
1158        #------------------------------------
1159           
1160        # Compute fluxes across each element edge
1161        self.compute_fluxes()
1162
1163        # Update conserved quantities
1164        self.update_conserved_quantities()
1165
1166        #------------------------------------
1167        # Combine final and initial values
1168        # and cleanup
1169        #------------------------------------
1170       
1171        # Combine steps
1172        self.saxpy_conserved_quantities(2.0/3.0, 1.0/3.0)
1173 
1174        # Update ghosts
1175        self.update_ghosts()
1176
1177        # Update vertex and edge values
1178        self.distribute_to_vertices_and_edges()
1179
1180        # Update boundary values
1181        self.update_boundary()
1182
1183        # Set new time
1184        self.time = initial_time + self.timestep       
1185       
1186
1187    def backup_conserved_quantities(self):
1188        N = len(self) # Number_of_triangles
1189
1190        # Backup conserved_quantities centroid values
1191        for name in self.conserved_quantities:
1192            Q = self.quantities[name]
1193            Q.backup_centroid_values()       
1194
1195    def saxpy_conserved_quantities(self,a,b):
1196        N = len(self) #number_of_triangles
1197
1198        # Backup conserved_quantities centroid values
1199        for name in self.conserved_quantities:
1200            Q = self.quantities[name]
1201            Q.saxpy_centroid_values(a,b)       
1202
1203
1204
1205       
1206    def distribute_to_vertices_and_edges(self):
1207        """Extrapolate conserved quantities from centroid to
1208        vertices and edge-midpoints for each volume
1209
1210        Default implementation is straight first order,
1211        i.e. constant values throughout each element and
1212        no reference to non-conserved quantities.
1213        """
1214
1215        for name in self.conserved_quantities:
1216            Q = self.quantities[name]
1217            if self.order == 1:
1218                Q.extrapolate_first_order()
1219            elif self.order == 2:
1220                Q.extrapolate_second_order()
1221                #Q.limit()
1222            else:
1223                raise 'Unknown order'
1224            #Q.interpolate_from_vertices_to_edges()
1225
1226
1227    def update_ghosts(self):
1228        pass
1229   
1230    def update_boundary(self):
1231        """Go through list of boundary objects and update boundary values
1232        for all conserved quantities on boundary.
1233        """
1234
1235        #FIXME: Update only those that change (if that can be worked out)
1236        #FIXME: Boundary objects should not include ghost nodes.
1237        #for i, ((vol_id, edge_id), B) in enumerate(self.boundary_objects):
1238        #    q = B.evaluate(vol_id, edge_id)
1239        for i, ((vol_id, vertex_id), B) in enumerate(self.boundary_objects):
1240            q = B.evaluate(vol_id, vertex_id)
1241            #print 'q ',q
1242            for j, name in enumerate(self.evolved_quantities):
1243                #print 'name %s j = %f \n'%(name,j)
1244                Q = self.quantities[name]
1245                Q.boundary_values[i] = q[j]
1246                #print 'Q=',Q
1247   
1248    def update_timestep(self, yieldstep, finaltime):
1249
1250        from config import min_timestep, max_timestep
1251
1252        # self.timestep is calculated from speed of characteristics
1253        # Apply CFL condition here
1254        timestep = min(self.CFL*self.flux_timestep, max_timestep)
1255
1256        #Record maximal and minimal values of timestep for reporting
1257        self.max_timestep = max(timestep, self.max_timestep)
1258        self.min_timestep = min(timestep, self.min_timestep)
1259
1260        #Protect against degenerate time steps
1261        if timestep < min_timestep:
1262
1263            #Number of consecutive small steps taken b4 taking action
1264            self.smallsteps += 1
1265
1266            if self.smallsteps > self.max_smallsteps:
1267                self.smallsteps = 0 #Reset
1268
1269                if self.order == 1:
1270                    msg = 'WARNING: Too small timestep %.16f reached '\
1271                          %timestep
1272                    msg += 'even after %d steps of 1 order scheme'\
1273                           %self.max_smallsteps
1274                    print msg
1275                    timestep = min_timestep  #Try enforcing min_step
1276
1277                    #raise msg
1278                else:
1279                    #Try to overcome situation by switching to 1 order
1280                    print "changing Order 1"
1281                    self.order = 1
1282
1283        else:
1284            self.smallsteps = 0
1285            if self.order == 1 and self.default_order == 2:
1286                self.order = 2
1287
1288
1289        #Ensure that final time is not exceeded
1290        if finaltime is not None and self.time + timestep > finaltime:
1291            timestep = finaltime-self.time
1292
1293        #Ensure that model time is aligned with yieldsteps
1294        if self.yieldtime + timestep > yieldstep:
1295            timestep = yieldstep-self.yieldtime
1296
1297        self.timestep = timestep
1298
1299    def update_extrema(self):
1300        pass
1301
1302    def compute_forcing_terms(self):
1303        """If there are any forcing functions driving the system
1304        they should be defined in Domain subclass and appended to
1305        the list self.forcing_terms
1306        """
1307        #Clears explicit_update needed for second order method
1308        if self.time_order == 2:
1309            for name in self.conserved_quantities:
1310                Q = self.quantities[name]
1311                Q.explicit_update[:] = 0.0
1312
1313        for f in self.forcing_terms:
1314            f(self)
1315           
1316
1317    def update_derived_quantites(self):
1318        pass
1319   
1320    #def update_conserved_quantities(self):
1321    def update_conserved_quantities(self):
1322        """Update vectors of conserved quantities using previously
1323        computed fluxes specified forcing functions.
1324        """
1325
1326
1327        timestep = self.timestep
1328
1329        #Compute forcing terms
1330        self.compute_forcing_terms()
1331
1332        #Update conserved_quantities
1333        for name in self.conserved_quantities:
1334            Q = self.quantities[name]
1335            Q.update(timestep)
1336
1337
1338
1339if __name__ == "__main__":
1340
1341    points1 = [0.0, 1.0, 2.0, 3.0]
1342    D1 = Domain(points1)
1343
1344    print D1.get_coordinate(0)
1345    print D1.get_coordinate(0,1)
1346    print 'Number of Elements = ',D1.number_of_elements
1347
1348    try:
1349        print D1.get_coordinate(3)
1350    except:
1351        pass
1352    else:
1353        msg =  'Should have raised an out of bounds exception'
1354        raise msg
1355
1356    #points2 = [0.0, 1.0, 2.0, 3.0, 2.5]
1357    #D2 = Domain(points2)
Note: See TracBrowser for help on using the repository browser.