source: trunk/anuga_work/development/shallow_water_1d_old/domain_adjust.py @ 7884

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

Copied John Jakeman's and Sudi Mungkasi's 1d shallow water code
File size: 39.4 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 *
10from 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, geo_reference = None):
17        """
18        Build 1D elements from x coordinates
19        """
20
21        from Numeric import array, zeros, Float, Int
22       
23        #Store Points
24        self.coordinates = array(coordinates)
25
26        if geo_reference is None:
27            self.geo_reference = Geo_reference() #Use defaults
28        else:
29            self.geo_reference = geo_reference
30
31        #Register number of Elements
32        self.number_of_elements = N = len(self.coordinates)-1
33
34        self.beta = 1.0
35        self.limiter = "minmod_kurganov"
36        self.split = False
37        self.wet_nodes = zeros((N,2), Int) # should this be here
38
39        #Allocate space for neighbour and boundary structures
40        self.neighbours = zeros((N, 2), Int)
41        #self.neighbour_edges = zeros((N, 2), Int)
42        self.neighbour_vertices = zeros((N, 2), Int)
43        self.number_of_boundaries = zeros(N, Int)
44        self.surrogate_neighbours = zeros((N, 2), Int)
45       
46        #Allocate space for geometric quantities
47        self.vertices  = zeros((N, 2), Float)
48        self.centroids = zeros(N, Float)
49        self.areas     = zeros(N, Float)
50
51        self.normals = zeros((N, 2), Float)
52       
53        for i in range(N):
54            xl = self.coordinates[i]
55            xr = self.coordinates[i+1]
56            self.vertices[i,0] = xl
57            self.vertices[i,1] = xr
58
59            centroid = (xl+xr)/2
60            self.centroids[i] = centroid
61
62            msg = 'Coordinates should be ordered, smallest to largest'
63            assert xr>xl, msg
64           
65            #The normal vectors
66            # - point outward from each edge
67            # - are orthogonal to the edge
68            # - have unit length
69            # - Are enumerated by left vertex then right vertex normals
70
71            nl = -1.0
72            nr =  1.0
73            self.normals[i,:] = [nl, nr]
74
75            self.areas[i] = (xr-xl)
76           
77##         print 'N', N
78##         print 'Centroid', self.centroids
79##         print 'Areas', self.areas
80##         print 'Vertex_Coordinates', self.vertices
81           
82            #Initialise Neighbours (-1 means that it is a boundary neighbour)
83            self.neighbours[i, :] = [-1, -1]
84            #Initialise edge ids of neighbours
85            #Initialise vertex ids of neighbours
86            #In case of boundaries this slot is not used
87            #self.neighbour_edges[i, :] = [-1, -1]
88            self.neighbour_vertices[i, :] = [-1, -1]
89
90        self.build_vertexlist()
91
92        #Build neighbour structure
93        self.build_neighbour_structure()
94
95        #Build surrogate neighbour structure
96        self.build_surrogate_neighbour_structure()         
97
98        #Build boundary dictionary mapping (id, edge) to symbolic tags
99        #Build boundary dictionary mapping (id, vertex) to symbolic tags
100        self.build_boundary_dictionary(boundary)
101
102        #Build tagged element  dictionary mapping (tag) to array of elements
103        self.build_tagged_elements_dictionary(tagged_elements)
104       
105        from quantity_adjust import Quantity, Conserved_quantity
106       
107        #List of quantity names entering
108        #the conservation equations
109        #(Must be a subset of quantities)
110        if conserved_quantities is None:
111            self.conserved_quantities = []
112        else:
113            self.conserved_quantities = conserved_quantities
114
115        if other_quantities is None:
116            self.other_quantities = []
117        else:
118            self.other_quantities = other_quantities
119
120
121        #Build dictionary of Quantity instances keyed by quantity names
122        self.quantities = {}
123
124        #FIXME: remove later - maybe OK, though....
125        for name in self.conserved_quantities:
126            self.quantities[name] = Conserved_quantity(self)
127        for name in self.other_quantities:
128            self.quantities[name] = Quantity(self)
129
130        #Create an empty list for explicit forcing terms
131        self.forcing_terms = []
132
133        #Defaults
134        from config import max_smallsteps, beta_w, beta_h, epsilon, CFL
135        self.beta_w = beta_w
136        self.beta_h = beta_h
137        self.epsilon = epsilon
138
139        #FIXME: Maybe have separate orders for h-limiter and w-limiter?
140        #Or maybe get rid of order altogether and use beta_w and beta_h
141        self.default_order = 1
142        self.order = self.default_order
143
144        self.default_time_order = 1
145        self.time_order = self.default_time_order
146       
147        self.smallsteps = 0
148        self.max_smallsteps = max_smallsteps
149        self.number_of_steps = 0
150        self.number_of_first_order_steps = 0
151        self.CFL = CFL
152
153        #Model time
154        self.time = 0.0
155        self.finaltime = None
156        self.min_timestep = self.max_timestep = 0.0
157        self.starttime = 0 #Physical starttime if any (0 is 1 Jan 1970 00:00:00)
158        #Checkpointing and storage
159        from config import default_datadir
160        self.datadir = default_datadir
161        self.filename = 'domain'
162        self.checkpoint = False
163
164    def __len__(self):
165        return self.number_of_elements
166
167    def build_vertexlist(self):
168        """Build vertexlist index by vertex ids and for each entry (point id)
169        build a list of (triangles, vertex_id) pairs that use the point
170        as vertex.
171
172        Preconditions:
173          self.coordinates and self.triangles are defined
174
175        Postcondition:
176          self.vertexlist is built
177        """
178        from Numeric import array
179
180        vertexlist = [None]*len(self.coordinates)
181        for i in range(self.number_of_elements):
182
183            #a = self.triangles[i, 0]
184            #b = self.triangles[i, 1]
185            #c = self.triangles[i, 2]
186            a = i
187            b = i + 1
188
189            #Register the vertices v as lists of
190            #(triangle_id, vertex_id) tuples associated with them
191            #This is used for smoothing
192            #for vertex_id, v in enumerate([a,b,c]):
193            for vertex_id, v in enumerate([a,b]):
194                if vertexlist[v] is None:
195                    vertexlist[v] = []
196
197                vertexlist[v].append( (i, vertex_id) )
198
199        self.vertexlist = vertexlist
200
201       
202    def build_neighbour_structure(self):
203        """Update all registered triangles to point to their neighbours.
204
205        Also, keep a tally of the number of boundaries for each triangle
206
207        Postconditions:
208          neighbours and neighbour_edges is populated
209          neighbours and neighbour_vertices is populated
210          number_of_boundaries integer array is defined.
211        """
212
213        #Step 1:
214        #Build dictionary mapping from segments (2-tuple of points)
215        #to left hand side edge (facing neighbouring triangle)
216
217        N = self.number_of_elements
218        neighbourdict = {}
219        #l_edge = 0
220        #r_edge = 1
221        l_vertex = 0
222        r_vertex = 1
223        for i in range(N):
224
225            #Register all segments as keys mapping to current triangle
226            #and segment id
227            #a = self.triangles[i, 0]
228            #b = self.triangles[i, 1]
229            #c = self.triangles[i, 2]
230            a = self.vertices[i,0]
231            b = self.vertices[i,1]
232           
233            """
234            if neighbourdict.has_key((a,b)):
235                    msg = "Edge 2 of triangle %d is duplicating edge %d of triangle %d.\n" %(i,neighbourdict[a,b][1],neighbourdict[a,b][0])
236                    raise msg
237            if neighbourdict.has_key((b,c)):
238                    msg = "Edge 0 of triangle %d is duplicating edge %d of triangle %d.\n" %(i,neighbourdict[b,c][1],neighbourdict[b,c][0])
239                    raise msg
240            if neighbourdict.has_key((c,a)):
241                    msg = "Edge 1 of triangle %d is duplicating edge %d of triangle %d.\n" %(i,neighbourdict[c,a][1],neighbourdict[c,a][0])
242                    raise msg
243                    """
244            #neighbourdict[a,b] = (i, 2) #(id, edge)
245            #neighbourdict[b,c] = (i, 0) #(id, edge)
246            #neighbourdict[c,a] = (i, 1) #(id, edge)
247            #neighbourdict[a,b] = (i, 1) #(id, edge)
248            #neighbourdict[b,a] = (i, 0) #(id, edge)
249            #neighbourdict[a,l_edge] = (i, 0) #(id, edge)
250            #neighbourdict[b,r_edge] = (i, 1) #(id, edge)
251            neighbourdict[a,l_vertex] = (i, 0) #(id, vertex)
252            neighbourdict[b,r_vertex] = (i, 1) #(id, vertex)
253
254
255        #Step 2:
256        #Go through triangles again, but this time
257        #reverse direction of segments and lookup neighbours.
258        for i in range(N):
259            #a = self.triangles[i, 0]
260            #b = self.triangles[i, 1]
261            #c = self.triangles[i, 2]
262           
263            a = self.vertices[i,0]
264            b = self.vertices[i,1]
265           
266            #self.number_of_boundaries[i] = 3
267            self.number_of_boundaries[i] = 2
268            #if neighbourdict.has_key((b,l_edge)):
269            if neighbourdict.has_key((b,l_vertex)):
270                #self.neighbours[i, 1] = neighbourdict[b,l_edge][0]
271                #self.neighbour_edges[i, 1] = neighbourdict[b,l_edge][1]
272                self.neighbours[i, 1] = neighbourdict[b,l_vertex][0]
273                self.neighbour_vertices[i, 1] = neighbourdict[b,l_vertex][1]
274                self.number_of_boundaries[i] -= 1
275
276            #if neighbourdict.has_key((a,r_edge)):
277            if neighbourdict.has_key((a,r_vertex)):
278                #self.neighbours[i, 0] = neighbourdict[a,r_edge][0]
279                #self.neighbour_edges[i, 0] = neighbourdict[a,r_edge][1]
280                self.neighbours[i, 0] = neighbourdict[a,r_vertex][0]
281                self.neighbour_vertices[i, 0] = neighbourdict[a,r_vertex][1]
282                self.number_of_boundaries[i] -= 1
283               
284            #if neighbourdict.has_key((b,a)):
285            #    self.neighbours[i, 1] = neighbourdict[b,a][0]
286            #    self.neighbour_edges[i, 1] = neighbourdict[b,a][1]
287            #    self.number_of_boundaries[i] -= 1
288
289            #if neighbourdict.has_key((c,b)):
290            #    self.neighbours[i, 0] = neighbourdict[c,b][0]
291            #    self.neighbour_edges[i, 0] = neighbourdict[c,b][1]
292            #    self.number_of_boundaries[i] -= 1
293
294            #if neighbourdict.has_key((a,b)):
295            #    self.neighbours[i, 0] = neighbourdict[a,b][0]
296            #    self.neighbour_edges[i, 0] = neighbourdict[a,b][1]
297            #    self.number_of_boundaries[i] -= 1
298
299    def build_surrogate_neighbour_structure(self):
300        """Build structure where each triangle edge points to its neighbours
301        if they exist. Otherwise point to the triangle itself.
302
303        The surrogate neighbour structure is useful for computing gradients
304        based on centroid values of neighbours.
305
306        Precondition: Neighbour structure is defined
307        Postcondition:
308          Surrogate neighbour structure is defined:
309          surrogate_neighbours: i0, i1, i2 where all i_k >= 0 point to
310          triangles.
311
312        """
313
314        N = self.number_of_elements
315        for i in range(N):
316            #Find all neighbouring volumes that are not boundaries
317            #for k in range(3):
318            for k in range(2):   
319                if self.neighbours[i, k] < 0:
320                    self.surrogate_neighbours[i, k] = i #Point this triangle
321                else:
322                    self.surrogate_neighbours[i, k] = self.neighbours[i, k]
323
324    def build_boundary_dictionary(self, boundary = None):
325        """Build or check the dictionary of boundary tags.
326         self.boundary is a dictionary of tags,
327         keyed by volume id and edge:
328         { (id, edge): tag, ... }
329
330         Postconditions:
331            self.boundary is defined.
332        """
333
334        from config import default_boundary_tag
335
336        if boundary is None:
337            boundary = {}
338            for vol_id in range(self.number_of_elements):
339                #for edge_id in range(0, 3):
340                #for edge_id in range(0, 2):
341                for vertex_id in range(0, 2):   
342                    #if self.neighbours[vol_id, edge_id] < 0:
343                    if self.neighbours[vol_id, vertex_id] < 0:   
344                        #boundary[(vol_id, edge_id)] = default_boundary_tag
345                        boundary[(vol_id, vertex_id)] = default_boundary_tag
346        else:
347            #Check that all keys in given boundary exist
348            #for vol_id, edge_id in boundary.keys():
349            for vol_id, vertex_id in boundary.keys():   
350                #msg = 'Segment (%d, %d) does not exist' %(vol_id, edge_id)
351                msg = 'Segment (%d, %d) does not exist' %(vol_id, vertex_id)
352                a, b = self.neighbours.shape
353                #assert vol_id < a and edge_id < b, msg
354                assert vol_id < a and vertex_id < b, msg
355
356                #FIXME: This assert violates internal boundaries (delete it)
357                #msg = 'Segment (%d, %d) is not a boundary' %(vol_id, edge_id)
358                #assert self.neighbours[vol_id, edge_id] < 0, msg
359
360            #Check that all boundary segments are assigned a tag
361            for vol_id in range(self.number_of_elements):
362                #for edge_id in range(0, 3):
363                #for edge_id in range(0, 2):
364                for vertex_id in range(0, 2):   
365                    #if self.neighbours[vol_id, edge_id] < 0:
366                    if self.neighbours[vol_id, vertex_id] < 0:   
367                        #if not boundary.has_key( (vol_id, edge_id) ):
368                        if not boundary.has_key( (vol_id, vertex_id) ):   
369                            msg = 'WARNING: Given boundary does not contain '
370                            #msg += 'tags for edge (%d, %d). '\
371                            #       %(vol_id, edge_id)
372                            msg += 'tags for vertex (%d, %d). '\
373                                   %(vol_id, vertex_id)
374                            msg += 'Assigning default tag (%s).'\
375                                   %default_boundary_tag
376
377                            #FIXME: Print only as per verbosity
378                            #print msg
379
380                            #FIXME: Make this situation an error in the future
381                            #and make another function which will
382                            #enable default boundary-tags where
383                            #tags a not specified
384                            #boundary[ (vol_id, edge_id) ] =\
385                            boundary[ (vol_id, vertex_id) ] =\
386                                      default_boundary_tag
387
388
389
390        self.boundary = boundary
391
392    def build_tagged_elements_dictionary(self, tagged_elements = None):
393        """Build the dictionary of element tags.
394         self.tagged_elements is a dictionary of element arrays,
395         keyed by tag:
396         { (tag): [e1, e2, e3..] }
397
398         Postconditions:
399            self.element_tag is defined
400        """
401        from Numeric import array, Int
402
403        if tagged_elements is None:
404            tagged_elements = {}
405        else:
406            #Check that all keys in given boundary exist
407            for tag in tagged_elements.keys():
408                tagged_elements[tag] = array(tagged_elements[tag]).astype(Int)
409
410                msg = 'Not all elements exist. '
411                assert max(tagged_elements[tag]) < self.number_of_elements, msg
412        #print "tagged_elements", tagged_elements
413        self.tagged_elements = tagged_elements
414
415    def get_boundary_tags(self):
416        """Return list of available boundary tags
417        """
418
419        tags = {}
420        for v in self.boundary.values():
421            tags[v] = 1
422
423        return tags.keys()
424       
425    def get_vertex_coordinates(self, obj = False):
426        """Return all vertex coordinates.
427        Return all vertex coordinates for all triangles as an Nx6 array
428        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)
429
430        if obj is True, the x/y pairs are returned in a 3*N x 2 array.
431        FIXME, we might make that the default.
432        FIXME Maybe use keyword: continuous = False for this condition?
433
434       
435        """
436
437        if obj is True:
438            from Numeric import concatenate, reshape
439            #V = self.vertex_coordinates
440            V = self.vertices
441            #return concatenate( (V[:,0:2], V[:,2:4], V[:,4:6]), axis=0)
442
443            N = V.shape[0]
444            #return reshape(V, (3*N, 2))
445            return reshape(V, (N, 2))
446        else:   
447            #return self.vertex_coordinates
448            return self.vertices
449   
450    def get_conserved_quantities(self, vol_id, vertex=None):#, edge=None):
451        """Get conserved quantities at volume vol_id
452
453        If vertex is specified use it as index for vertex values
454        If edge is specified use it as index for edge values
455        If neither are specified use centroid values
456        If both are specified an exeception is raised
457
458        Return value: Vector of length == number_of_conserved quantities
459
460        """
461
462        from Numeric import zeros, Float
463
464        #if not (vertex is None):# or edge is None):
465        #    msg = 'Values for both vertex and edge was specified.'
466        #    msg += 'Only one (or none) is allowed.'
467       #     raise msg
468
469        q = zeros( len(self.conserved_quantities), Float)
470
471        for i, name in enumerate(self.conserved_quantities):
472            Q = self.quantities[name]
473            if vertex is not None:
474                q[i] = Q.vertex_values[vol_id, vertex]
475            #elif edge is not None:
476            #    q[i] = Q.edge_values[vol_id, edge]
477            else:
478                q[i] = Q.centroid_values[vol_id]
479
480        return q
481       
482
483    def get_centroids(self):
484        """Return all coordinates of centroids
485        Return x coordinate of centroid for each element as a N array
486        """
487
488        return self.centroids
489
490    def get_vertices(self):
491        """Return all coordinates of centroids
492        Return x coordinate of centroid for each element as a N array
493        """
494
495        return self.vertices
496
497    def get_coordinate(self, elem_id, vertex=None):
498        """Return coordinate of centroid,
499        or left or right vertex.
500        Left vertex (vertex=0). Right vertex (vertex=1)
501        """
502
503        if vertex is None:
504            return self.centroids[elem_id]
505        else:
506            return self.vertices[elem_id,vertex]
507
508    def get_area(self, elem_id):
509        """Return area of element id
510        """
511
512        return self.areas[elem_id]
513
514    def get_quantity(self, name, location='vertices', indices = None):
515        """Get values for named quantity
516
517        name: Name of quantity
518
519        In case of location == 'centroids' the dimension values must
520        be a list of a Numerical array of length N, N being the number
521        of elements. Otherwise it must be of dimension Nx3.
522
523        Indices is the set of element ids that the operation applies to.
524
525        The values will be stored in elements following their
526        internal ordering.
527        """
528
529        return self.quantities[name].get_values( location, indices = indices)
530
531    def get_centroid_coordinates(self):
532        """Return all centroid coordinates.
533        Return all centroid coordinates for all triangles as an Nx2 array
534        (ordered as x0, y0 for each triangle)
535        """
536        return self.centroids
537
538    def set_quantity(self, name, *args, **kwargs):
539        """Set values for named quantity
540
541
542        One keyword argument is documented here:
543        expression = None, # Arbitrary expression
544
545        expression:
546          Arbitrary expression involving quantity names
547
548        See Quantity.set_values for further documentation.
549        """
550
551        #FIXME (Ole): Allow new quantities here
552        #from quantity import Quantity, Conserved_quantity
553        #Create appropriate quantity object
554        ##if name in self.conserved_quantities:
555        ##    self.quantities[name] = Conserved_quantity(self)
556        ##else:
557        ##    self.quantities[name] = Quantity(self)
558
559
560        #Do the expression stuff
561        if kwargs.has_key('expression'):
562            expression = kwargs['expression']
563            del kwargs['expression']
564
565            Q = self.create_quantity_from_expression(expression)
566            kwargs['quantity'] = Q
567       
568        #Assign values
569        self.quantities[name].set_values(*args, **kwargs)
570
571    def set_boundary(self, boundary_map):
572        """Associate boundary objects with tagged boundary segments.
573
574        Input boundary_map is a dictionary of boundary objects keyed
575        by symbolic tags to matched against tags in the internal dictionary
576        self.boundary.
577
578        As result one pointer to a boundary object is stored for each vertex
579        in the list self.boundary_objects.
580        More entries may point to the same boundary object
581
582        Schematically the mapping is from two dictionaries to one list
583        where the index is used as pointer to the boundary_values arrays
584        within each quantity.
585
586        self.boundary:          (vol_id, edge_id): tag
587        boundary_map (input):   tag: boundary_object
588        ----------------------------------------------
589        self.boundary_objects:  ((vol_id, edge_id), boundary_object)
590
591
592        Pre-condition:
593          self.boundary has been built.
594
595        Post-condition:
596          self.boundary_objects is built
597
598        If a tag from the domain doesn't appear in the input dictionary an
599        exception is raised.
600        However, if a tag is not used to the domain, no error is thrown.
601        FIXME: This would lead to implementation of a
602        default boundary condition
603
604        Note: If a segment is listed in the boundary dictionary and if it is
605        not None, it *will* become a boundary -
606        even if there is a neighbouring triangle.
607        This would be the case for internal boundaries
608
609        Boundary objects that are None will be skipped.
610
611        FIXME: If set_boundary is called multiple times and if Boundary
612        object is changed into None, the neighbour structure will not be
613        restored!!!
614        """
615
616        self.boundary_objects = []
617       
618
619
620
621       
622        self.boundary_map = boundary_map  #Store for use with eg. boundary_stats.
623
624        #FIXME: Try to remove the sorting and fix test_mesh.py
625        x = self.boundary.keys()
626        x.sort()
627
628        #Loop through edges that lie on the boundary and associate them with
629        #callable boundary objects depending on their tags
630        #for k, (vol_id, edge_id) in enumerate(x):
631        for k, (vol_id, vertex_id) in enumerate(x):
632            #tag = self.boundary[ (vol_id, edge_id) ]
633            tag = self.boundary[ (vol_id, vertex_id) ]
634
635            if boundary_map.has_key(tag):
636                B = boundary_map[tag]  #Get callable boundary object
637
638                if B is not None:
639                    #self.boundary_objects.append( ((vol_id, edge_id), B) )
640                    #self.neighbours[vol_id, edge_id] = -len(self.boundary_objects)
641                    self.boundary_objects.append( ((vol_id, vertex_id), B) )
642                    self.neighbours[vol_id, vertex_id] = -len(self.boundary_objects)
643                else:
644                    pass
645                    #FIXME: Check and perhaps fix neighbour structure
646
647            else:
648                msg = 'ERROR (domain.py): Tag "%s" has not been ' %tag
649                msg += 'bound to a boundary object.\n'
650                msg += 'All boundary tags defined in domain must appear '
651                msg += 'in the supplied dictionary.\n'
652                msg += 'The tags are: %s' %self.get_boundary_tags()
653                raise msg
654
655
656
657    def check_integrity(self):
658        #Mesh.check_integrity(self)
659       
660        for quantity in self.conserved_quantities:
661            msg = 'Conserved quantities must be a subset of all quantities'
662            assert quantity in self.quantities, msg
663
664        ##assert hasattr(self, 'boundary_objects')
665
666    def write_time(self):
667        print self.timestepping_statistics()
668
669    def timestepping_statistics(self):
670        """Return string with time stepping statistics for printing or logging
671        """
672
673        msg = ''
674        if self.min_timestep == self.max_timestep:
675            msg += 'Time = %.4f, delta t = %.8f, steps=%d (%d)'\
676                   %(self.time, self.min_timestep, self.number_of_steps,
677                     self.number_of_first_order_steps)
678        elif self.min_timestep > self.max_timestep:
679            msg += 'Time = %.4f, steps=%d (%d)'\
680                   %(self.time, self.number_of_steps,
681                     self.number_of_first_order_steps)
682        else:
683            msg += 'Time = %.4f, delta t in [%.8f, %.8f], steps=%d (%d)'\
684                   %(self.time, self.min_timestep,
685                     self.max_timestep, self.number_of_steps,
686                     self.number_of_first_order_steps)
687
688        return msg
689
690    def get_name(self):
691        return self.filename
692
693    def set_name(self, name):
694        self.filename = name
695
696    def get_datadir(self):
697        return self.datadir
698
699    def set_datadir(self, name):
700        self.datadir = name
701
702    #Main components of evolve
703    def evolve(self, yieldstep = None, finaltime = None,
704               skip_initial_step = False):
705        """Evolve model from time=0.0 to finaltime yielding results
706        every yieldstep.
707
708        Internally, smaller timesteps may be taken.
709
710        Evolve is implemented as a generator and is to be called as such, e.g.
711
712        for t in domain.evolve(timestep, yieldstep, finaltime):
713            <Do something with domain and t>
714
715        """
716
717        from config import min_timestep, max_timestep, epsilon
718
719        #FIXME: Maybe lump into a larger check prior to evolving
720        msg = 'Boundary tags must be bound to boundary objects before evolving system, '
721        msg += 'e.g. using the method set_boundary.\n'
722        msg += 'This system has the boundary tags %s ' %self.get_boundary_tags()
723        assert hasattr(self, 'boundary_objects'), msg
724
725        ##self.set_defaults()
726
727        if yieldstep is None:
728            yieldstep = max_timestep
729        else:
730            yieldstep = float(yieldstep)
731
732        self.order = self.default_order
733        self.time_order = self.default_time_order
734
735        self.yieldtime = 0.0 #Time between 'yields'
736
737        #Initialise interval of timestep sizes (for reporting only)
738        # SEEMS WIERD
739        self.min_timestep = max_timestep
740        self.max_timestep = min_timestep
741        self.finaltime = finaltime
742        self.number_of_steps = 0
743        self.number_of_first_order_steps = 0
744       
745        #update ghosts
746        #self.update_ghosts()
747   
748        #Initial update of vertex and edge values
749        self.distribute_to_vertices_and_edges()
750
751        self.distribute_stage_to_height()
752
753        #Initial update boundary values
754        self.update_boundary()
755
756        #Or maybe restore from latest checkpoint
757        if self.checkpoint is True:
758            self.goto_latest_checkpoint()
759
760        if skip_initial_step is False:
761            yield(self.time)  #Yield initial values   
762
763        while True:
764            if self.time_order == 1:
765                #Compute fluxes across each element edge
766                self.compute_fluxes()
767                #Update timestep to fit yieldstep and finaltime
768                self.update_timestep(yieldstep, finaltime)
769                #Compute forcing terms
770                self.compute_forcing_terms()
771                #Update conserved quantities
772                self.update_conserved_quantities(self.timestep)
773
774                self.adjust_partially_submerged_cells()
775                #update ghosts
776                #self.update_ghosts()
777                #Update vertex and edge values
778                self.distribute_to_vertices_and_edges()
779                self.distribute_stage_to_height()
780                #Update boundary values
781                self.update_boundary()
782
783            elif self.time_order == 2:
784               
785                self.compute_timestep()
786
787                #Solve inhomogeneous operator for half a timestep
788                self.solve_inhomogenous_second_order(yieldstep, finaltime)
789               
790                #Solve homogeneous operator for full timestep using
791                #Harten second order timestepping
792                self.solve_homogenous_second_order(yieldstep,finaltime)
793
794                #Solve inhomogeneous operator for half a timestep
795                self.solve_inhomogenous_second_order(yieldstep, finaltime)
796               
797            #Update time
798            self.time += self.timestep
799            self.yieldtime += self.timestep
800            self.number_of_steps += 1
801            if self.order == 1:
802                self.number_of_first_order_steps += 1
803
804            #Yield results
805            if finaltime is not None and abs(self.time - finaltime) < epsilon:
806
807                #FIXME: There is a rare situation where the
808                #final time step is stored twice. Can we make a test?
809
810                # Yield final time and stop
811                yield(self.time)
812                break
813
814
815            if abs(self.yieldtime - yieldstep) < epsilon:
816                # Yield (intermediate) time and allow inspection of domain
817
818                if self.checkpoint is True:
819                    self.store_checkpoint()
820                    self.delete_old_checkpoints()
821
822                #Pass control on to outer loop for more specific actions
823                yield(self.time)
824
825                # Reinitialise
826                self.yieldtime = 0.0
827                self.min_timestep = max_timestep
828                self.max_timestep = min_timestep
829                self.number_of_steps = 0
830                self.number_of_first_order_steps = 0
831
832    def solve_inhomogenous_second_order(self,yieldstep, finaltime):
833
834        #Update timestep to fit yieldstep and finaltime
835        self.update_timestep(yieldstep, finaltime)
836        #Compute forcing terms
837        self.compute_forcing_terms()
838        #Update conserved quantities
839        self.update_conserved_quantities(0.5*self.timestep)
840        self.adjust_partially_submerged_cells()
841        #Update vertex and edge values
842        self.distribute_to_vertices_and_edges()
843        self.distribute_stage_to_height()
844        #Update boundary values
845        self.update_boundary()
846
847    def  solve_homogenous_second_order(self,yieldstep,finaltime):
848        """Use Shu Second order timestepping to update
849        conserved quantities
850
851        q^{n+1/2} = q^{n}+0.5*dt*F^{n}
852        q^{n+1} = q^{n}+dt*F^{n+1/2}
853        """
854        import copy
855        from Numeric import zeros,Float
856
857        N = self.number_of_elements
858
859        self.compute_fluxes()
860        #Update timestep to fit yieldstep and finaltime
861        self.update_timestep(yieldstep, finaltime)
862        #Compute forcing terms
863        #NOT NEEDED FOR 2ND ORDER STRANG SPLIITING
864        #ADDING THIS WILL NEED TO REMOVE ZEROING IN COMPUTE_FORCING
865        #self.compute_forcing_terms()
866
867        QC = zeros((N,len(self.conserved_quantities)),Float)
868        QF = zeros((N,len(self.conserved_quantities)),Float)
869
870        i = 0
871        for name in self.conserved_quantities:
872            Q = self.quantities[name]
873            #Store the centroid values at time t^n
874            QC[:,i] = copy.copy(Q.centroid_values)
875            QF[:,i] = copy.copy(Q.explicit_update)
876            #Update conserved quantities
877            Q.update(self.timestep)
878            i+=1
879
880        self.adjust_partially_submerged_cells()
881        #Update vertex and edge values
882        self.distribute_to_vertices_and_edges()
883        self.distribute_stage_to_height()
884        #Update boundary values
885        self.update_boundary()
886       
887        self.compute_fluxes()
888        self.update_timestep(yieldstep, finaltime)
889        #Compute forcing terms
890        #NOT NEEDED FOR 2ND ORDER STRANG SPLIITING
891        #self.compute_forcing_terms()
892
893        i = 0
894        for name in self.conserved_quantities:
895            Q = self.quantities[name]
896            Q.centroid_values = QC[:,i]
897            Q.explicit_update = 0.5*(Q.explicit_update+QF[:,i]) 
898            #Update conserved quantities
899            Q.update(self.timestep)
900            i+=1
901
902        self.adjust_partially_submerged_cells()
903       
904        #Update vertex and edge values
905        self.distribute_to_vertices_and_edges()
906        self.distribute_stage_to_height()
907        #Update boundary values
908        self.update_boundary()
909
910    def  solve_homogenous_second_order_harten(self,yieldstep,finaltime):
911        """Use Harten Second order timestepping to update
912        conserved quantities
913
914        q^{n+1/2} = q^{n}+0.5*dt*F^{n}
915        q^{n+1} = q^{n}+dt*F^{n+1/2}
916        """
917        import copy
918        from Numeric import zeros,Float
919
920        N = self.number_of_elements
921
922        self.compute_fluxes()
923        #Update timestep to fit yieldstep and finaltime
924        self.update_timestep(yieldstep, finaltime)
925        #Compute forcing terms
926        #NOT NEEDED FOR 2ND ORDER STRANG SPLIITING
927        #ADDING THIS WILL NEED TO REMOVE ZEROING IN COMPUTE_FORCING
928        #self.compute_forcing_terms()
929
930        QC = zeros((N,len(self.conserved_quantities)),Float)
931
932        i = 0
933        for name in self.conserved_quantities:
934            Q = self.quantities[name]
935            #Store the centroid values at time t^n
936            QC[:,i] = copy.copy(Q.centroid_values)
937            #Update conserved quantities
938            Q.update(0.5*self.timestep)
939            i+=1
940           
941        #Update vertex and edge values
942        self.distribute_to_vertices_and_edges()
943        #Update boundary values
944        self.update_boundary()
945       
946        self.compute_fluxes()
947        self.update_timestep(yieldstep, finaltime)
948        #Compute forcing terms
949        #NOT NEEDED FOR 2ND ORDER STRANG SPLIITING
950        #self.compute_forcing_terms()
951
952        i = 0
953        for name in self.conserved_quantities:
954            Q = self.quantities[name]
955            Q.centroid_values = QC[:,i] 
956            #Update conserved quantities
957            Q.update(self.timestep)
958            i+=1
959       
960        #Update vertex and edge values
961        self.distribute_to_vertices_and_edges()
962        #Update boundary values
963        self.update_boundary()
964
965
966    def distribute_to_vertices_and_edges(self):
967        """Extrapolate conserved quantities from centroid to
968        vertices and edge-midpoints for each volume
969
970        Default implementation is straight first order,
971        i.e. constant values throughout each element and
972        no reference to non-conserved quantities.
973        """
974
975        for name in self.conserved_quantities:
976            Q = self.quantities[name]
977            if self.order == 1:
978                Q.extrapolate_first_order()
979            elif self.order == 2:
980                Q.extrapolate_second_order()
981                Q.limit()
982            else:
983                raise 'Unknown order'
984            #Q.interpolate_from_vertices_to_edges()
985
986
987    def update_boundary(self):
988        """Go through list of boundary objects and update boundary values
989        for all conserved quantities on boundary.
990        """
991
992        #FIXME: Update only those that change (if that can be worked out)
993        #FIXME: Boundary objects should not include ghost nodes.
994        #for i, ((vol_id, edge_id), B) in enumerate(self.boundary_objects):
995        #    q = B.evaluate(vol_id, edge_id)
996        for i, ((vol_id, vertex_id), B) in enumerate(self.boundary_objects):
997            q = B.evaluate(vol_id, vertex_id)
998
999            for j, name in enumerate(self.conserved_quantities):
1000                Q = self.quantities[name]
1001                Q.boundary_values[i] = q[j]
1002
1003    def update_timestep(self, yieldstep, finaltime):
1004
1005        from config import min_timestep
1006
1007        # self.timestep is calculated from speed of characteristics
1008        # Apply CFL condition here
1009        timestep = self.CFL*self.timestep
1010
1011        #Record maximal and minimal values of timestep for reporting
1012        self.max_timestep = max(timestep, self.max_timestep)
1013        self.min_timestep = min(timestep, self.min_timestep)
1014
1015        #Protect against degenerate time steps
1016        if timestep < min_timestep:
1017
1018            #Number of consecutive small steps taken b4 taking action
1019            self.smallsteps += 1
1020
1021            if self.smallsteps > self.max_smallsteps:
1022                self.smallsteps = 0 #Reset
1023
1024                if self.order == 1:
1025                    msg = 'WARNING: Too small timestep %.16f reached '\
1026                          %timestep
1027                    msg += 'even after %d steps of 1 order scheme'\
1028                           %self.max_smallsteps
1029                    print msg
1030                    timestep = min_timestep  #Try enforcing min_step
1031
1032                    #raise msg
1033                else:
1034                    #Try to overcome situation by switching to 1 order
1035                    print "changing Order 1"
1036                    self.order = 1
1037
1038        else:
1039            self.smallsteps = 0
1040            if self.order == 1 and self.default_order == 2:
1041                self.order = 2
1042
1043
1044        #Ensure that final time is not exceeded
1045        if finaltime is not None and self.time + timestep > finaltime:
1046            timestep = finaltime-self.time
1047
1048        #Ensure that model time is aligned with yieldsteps
1049        if self.yieldtime + timestep > yieldstep:
1050            timestep = yieldstep-self.yieldtime
1051
1052        self.timestep = timestep
1053
1054
1055    def compute_forcing_terms(self):
1056        """If there are any forcing functions driving the system
1057        they should be defined in Domain subclass and appended to
1058        the list self.forcing_terms
1059        """
1060        #print "in compute_forcing_terms"
1061       
1062        #Clears explicit_update needed for second order method
1063        if self.time_order == 2:
1064            for name in self.conserved_quantities:
1065                Q = self.quantities[name]
1066                Q.explicit_update[:] = 0.0
1067
1068        for f in self.forcing_terms:
1069            f(self)
1070           
1071
1072    def update_derived_quantites(self):
1073        pass
1074   
1075    #def update_conserved_quantities(self):
1076    def update_conserved_quantities(self,timestep):
1077        """Update vectors of conserved quantities using previously
1078        computed fluxes specified forcing functions.
1079        """
1080
1081        from Numeric import ones, sum, equal, Float
1082
1083        N = self.number_of_elements
1084        d = len(self.conserved_quantities)
1085
1086        #timestep = self.timestep
1087
1088        #Compute forcing terms
1089        #self.compute_forcing_terms()
1090
1091        #print "in update conserved quantities"
1092
1093        #Update conserved_quantities
1094        for name in self.conserved_quantities:
1095            Q = self.quantities[name]
1096            #print "integral", Q.get_integral()
1097            Q.update(timestep)
1098            #print "integral", Q.get_integral()
1099
1100            #Clean up
1101            #Note that Q.explicit_update is reset by compute_fluxes
1102
1103            #MH090605 commented out the following since semi_implicit_update is now re-initialized
1104            #at the end of the _update function in quantity_ext.c (This is called by the
1105            #preceeding Q.update(timestep) statement above).
1106            #For run_profile.py with N=128, the time of update_conserved_quantities is cut from 14.00 secs
1107            #to 8.35 secs
1108
1109            #Q.semi_implicit_update[:] = 0.0
1110
1111if __name__ == "__main__":
1112
1113    points1 = [0.0, 1.0, 2.0, 3.0]
1114    D1 = Domain(points1)
1115
1116    print D1.get_coordinate(0)
1117    print D1.get_coordinate(0,1)
1118    print 'Number of Elements = ',D1.number_of_elements
1119
1120    try:
1121        print D1.get_coordinate(3)
1122    except:
1123        pass
1124    else:
1125        msg =  'Should have raised an out of bounds exception'
1126        raise msg
1127
1128    #points2 = [0.0, 1.0, 2.0, 3.0, 2.5]
1129    #D2 = Domain(points2)
1130
Note: See TracBrowser for help on using the repository browser.