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

Last change on this file since 6454 was 6453, checked in by steve, 15 years ago

Added Padarn's (2008/2009 summer student) files

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