source: anuga_work/development/anuga_1d/domain.py @ 5742

Last change on this file since 5742 was 5742, checked in by steve, 16 years ago

Found error in get_python_double

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