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

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

Moving 2010 project

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