source: development/pyvolution-1d/domain.py @ 3322

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