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

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

Working version of comp_flux_ext

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