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

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

Updating domain etc to take conserved, evolved and other quantitiies

File size: 52.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"""
9
10from generic_boundary_conditions import *
11
12
13class Domain:
14
15    def __init__(self,
16                 coordinates,
17                 boundary = None,
18                 conserved_quantities = None,
19                 evolved_quantities = None,
20                 other_quantities = None,
21                 tagged_elements = None):
22        """
23        Build 1D elements from x coordinates
24        """
25
26        from Numeric import array, zeros, Float, Int
27
28        from config import timestepping_method
29        from config import CFL
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
39        self.set_limiter("minmod_kurganov")
40        self.set_CFL(CFL)
41        self.set_timestepping_method(timestepping_method)
42
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
57        self.max_speed_array = zeros(N, Float)
58     
59
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           
86# #         print 'N', N
87# #         print 'Centroid', self.centroids
88# #         print 'Areas', self.areas
89# #         print 'Vertex_Coordinates', self.vertices
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
125        if evolved_quantities is None:
126            self.evolved_quantities = self.conserved_quantities
127        else:
128            self.evolved_quantities = evolved_quantities
129
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
139        #print self.conserved_quantities
140        #print self.evolved_quantities
141       
142
143        #FIXME: remove later - maybe OK, though....
144        for name in self.evolved_quantities:
145            self.quantities[name] = Quantity(self)
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
178        self.set_datadir(default_datadir)
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
433    def get_boundary_tags(self):
434        """Return list of available boundary tags
435        """
436
437        tags = {}
438        for v in self.boundary.values():
439            tags[v] = 1
440
441        return tags.keys()
442       
443    def get_vertex_coordinates(self, obj = False):
444        """Return all vertex coordinates.
445        Return all vertex coordinates for all triangles as an Nx6 array
446        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)
447
448        if obj is True, the x/y pairs are returned in a 3*N x 2 array.
449        FIXME, we might make that the default.
450        FIXME Maybe use keyword: continuous = False for this condition?
451
452       
453        """
454
455        if obj is True:
456            from Numeric import concatenate, reshape
457            #V = self.vertex_coordinates
458            V = self.vertices
459            #return concatenate( (V[:,0:2], V[:,2:4], V[:,4:6]), axis=0)
460
461            N = V.shape[0]
462            #return reshape(V, (3*N, 2))
463            return reshape(V, (N, 2))
464        else:   
465            #return self.vertex_coordinates
466            return self.vertices
467   
468    def get_conserved_quantities(self, vol_id, vertex=None):#, edge=None):
469        """Get conserved quantities at volume vol_id
470
471        If vertex is specified use it as index for vertex values
472        If edge is specified use it as index for edge values
473        If neither are specified use centroid values
474        If both are specified an exeception is raised
475
476        Return value: Vector of length == number_of_conserved quantities
477
478        """
479
480        from Numeric import zeros, Float
481
482        #if not (vertex is None):# or edge is None):
483        #    msg = 'Values for both vertex and edge was specified.'
484        #    msg += 'Only one (or none) is allowed.'
485       #     raise msg
486
487        q = zeros( len(self.conserved_quantities), Float)
488
489        for i, name in enumerate(self.conserved_quantities):
490            Q = self.quantities[name]
491            if vertex is not None:
492                q[i] = Q.vertex_values[vol_id, vertex]
493            #elif edge is not None:
494            #    q[i] = Q.edge_values[vol_id, edge]
495            else:
496                q[i] = Q.centroid_values[vol_id]
497
498        return q
499
500
501    def get_evolved_quantities(self, vol_id, vertex=None):#, edge=None):
502        """Get evolved quantities at volume vol_id
503
504        If vertex is specified use it as index for vertex values
505        If edge is specified use it as index for edge values
506        If neither are specified use centroid values
507        If both are specified an exeception is raised
508
509        Return value: Vector of length == number_of_evolved quantities
510
511        """
512
513        from Numeric import zeros, Float
514
515        #if not (vertex is None):# or edge is None):
516        #    msg = 'Values for both vertex and edge was specified.'
517        #    msg += 'Only one (or none) is allowed.'
518       #     raise msg
519
520        q = zeros( len(self.evolved_quantities), Float)
521
522        for i, name in enumerate(self.evolved_quantities):
523            Q = self.quantities[name]
524            if vertex is not None:
525                q[i] = Q.vertex_values[vol_id, vertex]
526            #elif edge is not None:
527            #    q[i] = Q.edge_values[vol_id, edge]
528            else:
529                q[i] = Q.centroid_values[vol_id]
530
531        return q
532       
533
534    def get_centroids(self):
535        """Return all coordinates of centroids
536        Return x coordinate of centroid for each element as a N array
537        """
538
539        return self.centroids
540
541    def get_vertices(self):
542        """Return all coordinates of centroids
543        Return x coordinate of centroid for each element as a N array
544        """
545
546        return self.vertices
547
548    def get_coordinate(self, elem_id, vertex=None):
549        """Return coordinate of centroid,
550        or left or right vertex.
551        Left vertex (vertex=0). Right vertex (vertex=1)
552        """
553
554        if vertex is None:
555            return self.centroids[elem_id]
556        else:
557            return self.vertices[elem_id,vertex]
558
559    def get_area(self, elem_id):
560        """Return area of element id
561        """
562
563        return self.areas[elem_id]
564
565    def get_quantity(self, name, location='vertices', indices = None):
566        """Get values for named quantity
567
568        name: Name of quantity
569
570        In case of location == 'centroids' the dimension values must
571        be a list of a Numerical array of length N, N being the number
572        of elements. Otherwise it must be of dimension Nx3.
573
574        Indices is the set of element ids that the operation applies to.
575
576        The values will be stored in elements following their
577        internal ordering.
578        """
579
580        return self.quantities[name].get_values( location, indices = indices)
581
582    def get_centroid_coordinates(self):
583        """Return all centroid coordinates.
584        Return all centroid coordinates for all triangles as an Nx2 array
585        (ordered as x0, y0 for each triangle)
586        """
587        return self.centroids
588
589
590    def get_timestepping_method(self):
591        return self.timestepping_method
592
593    def set_timestepping_method(self,timestepping_method):
594       
595        if timestepping_method in ['euler', 'rk2', 'rk3']:
596            self.timestepping_method = timestepping_method
597            return
598
599        msg = '%s is an incorrect timestepping type'% timestepping_method
600        raise Exception, msg
601
602
603    def set_quantity(self, name, *args, **kwargs):
604        """Set values for named quantity
605
606
607        One keyword argument is documented here:
608        expression = None, # Arbitrary expression
609
610        expression:
611          Arbitrary expression involving quantity names
612
613        See Quantity.set_values for further documentation.
614        """
615
616        #FIXME (Ole): Allow new quantities here
617        #from quantity import Quantity, Conserved_quantity
618        #Create appropriate quantity object
619        # #if name in self.conserved_quantities:
620        # #    self.quantities[name] = Conserved_quantity(self)
621        # #else:
622        # #    self.quantities[name] = Quantity(self)
623
624
625        #Do the expression stuff
626        if kwargs.has_key('expression'):
627            expression = kwargs['expression']
628            del kwargs['expression']
629
630            Q = self.create_quantity_from_expression(expression)
631            kwargs['quantity'] = Q
632       
633        #Assign values
634        self.quantities[name].set_values(*args, **kwargs)
635
636    def set_boundary(self, boundary_map):
637        """Associate boundary objects with tagged boundary segments.
638
639        Input boundary_map is a dictionary of boundary objects keyed
640        by symbolic tags to matched against tags in the internal dictionary
641        self.boundary.
642
643        As result one pointer to a boundary object is stored for each vertex
644        in the list self.boundary_objects.
645        More entries may point to the same boundary object
646
647        Schematically the mapping is from two dictionaries to one list
648        where the index is used as pointer to the boundary_values arrays
649        within each quantity.
650
651        self.boundary:          (vol_id, edge_id): tag
652        boundary_map (input):   tag: boundary_object
653        ----------------------------------------------
654        self.boundary_objects:  ((vol_id, edge_id), boundary_object)
655
656
657        Pre-condition:
658          self.boundary has been built.
659
660        Post-condition:
661          self.boundary_objects is built
662
663        If a tag from the domain doesn't appear in the input dictionary an
664        exception is raised.
665        However, if a tag is not used to the domain, no error is thrown.
666        FIXME: This would lead to implementation of a
667        default boundary condition
668
669        Note: If a segment is listed in the boundary dictionary and if it is
670        not None, it *will* become a boundary -
671        even if there is a neighbouring triangle.
672        This would be the case for internal boundaries
673
674        Boundary objects that are None will be skipped.
675
676        FIXME: If set_boundary is called multiple times and if Boundary
677        object is changed into None, the neighbour structure will not be
678        restored!!!
679        """
680
681        self.boundary_objects = []
682        self.boundary_map = boundary_map  #Store for use with eg. boundary_stats.
683
684        #FIXME: Try to remove the sorting and fix test_mesh.py
685        x = self.boundary.keys()
686        x.sort()
687
688        #Loop through edges that lie on the boundary and associate them with
689        #callable boundary objects depending on their tags
690        #for k, (vol_id, edge_id) in enumerate(x):
691        for k, (vol_id, vertex_id) in enumerate(x):
692            #tag = self.boundary[ (vol_id, edge_id) ]
693            tag = self.boundary[ (vol_id, vertex_id) ]
694
695            if boundary_map.has_key(tag):
696                B = boundary_map[tag]  #Get callable boundary object
697
698                if B is not None:
699                    #self.boundary_objects.append( ((vol_id, edge_id), B) )
700                    #self.neighbours[vol_id, edge_id] = -len(self.boundary_objects)
701                    self.boundary_objects.append( ((vol_id, vertex_id), B) )
702                    self.neighbours[vol_id, vertex_id] = -len(self.boundary_objects)
703                else:
704                    pass
705                    #FIXME: Check and perhaps fix neighbour structure
706
707            else:
708                msg = 'ERROR (domain.py): Tag "%s" has not been ' %tag
709                msg += 'bound to a boundary object.\n'
710                msg += 'All boundary tags defined in domain must appear '
711                msg += 'in the supplied dictionary.\n'
712                msg += 'The tags are: %s' %self.get_boundary_tags()
713                raise msg
714
715
716
717    def check_integrity(self):
718        #Mesh.check_integrity(self)
719
720        #print self.quantities
721        #print self.conserved_quantities
722       
723        for quantity in self.conserved_quantities:
724            msg = 'Conserved quantities must be a subset of all quantities'
725            assert quantity in self.quantities, msg
726
727        for quantity in self.evolved_quantities:
728            msg = 'Evolved quantities must be a subset of all quantities'
729            assert quantity in self.quantities, msg           
730
731        # #assert hasattr(self, 'boundary_objects')
732
733    def write_time(self):
734        print self.timestepping_statistics()
735
736    def timestepping_statistics(self):
737        """Return string with time stepping statistics for printing or logging
738        """
739
740        msg = ''
741        if self.min_timestep == self.max_timestep:
742            msg += 'Time = %.4f, delta t = %.8f, steps=%d (%d)'\
743                   %(self.time, self.min_timestep, self.number_of_steps,
744                     self.number_of_first_order_steps)
745        elif self.min_timestep > self.max_timestep:
746            msg += 'Time = %.4f, steps=%d (%d)'\
747                   %(self.time, self.number_of_steps,
748                     self.number_of_first_order_steps)
749        else:
750            msg += 'Time = %.4f, delta t in [%.8f, %.8f], steps=%d (%d)'\
751                   %(self.time, self.min_timestep,
752                     self.max_timestep, self.number_of_steps,
753                     self.number_of_first_order_steps)
754
755        return msg
756
757    def get_name(self):
758        return self.filename
759
760    def set_name(self, name):
761        self.filename = name
762
763    def get_datadir(self):
764        return self.datadir
765
766    def set_datadir(self, name):
767        self.datadir = name
768
769    def set_CFL(self, cfl):
770        if cfl > 1.0:
771            print 'WARNING: Setting CFL condition to %f which is greater than 1' % cfl
772        self.CFL = cfl
773
774    def get_CFL(self):
775        return self.CFL
776   
777    def set_filename(self, name):
778        self.filename = name
779
780    def get_filename(self):
781        return self.filename
782
783    def get_limiter(self):
784        return self.limiter
785
786    def set_limiter(self,limiter):
787
788        possible_limiters = \
789        ['pyvolution', 'minmod_steve', 'minmod', 'minmod_kurganov', 'superbee', 'vanleer', 'vanalbada']
790
791        if limiter in possible_limiters:
792            self.limiter = limiter
793            return
794
795        msg = '%s is an incorrect limiter type.\n'% limiter
796        msg += 'Possible types are: '+ ", ".join(["%s" % el for el in possible_limiters])
797        raise Exception, msg
798
799       
800    #--------------------------
801    # Main components of evolve
802    #--------------------------   
803
804    def evolve(self, yieldstep = None,
805               finaltime = None,
806               duration = None,
807               skip_initial_step = False):
808        """Evolve model through time starting from self.starttime.
809
810
811        yieldstep: Interval between yields where results are stored,
812                   statistics written and domain inspected or
813                   possibly modified. If omitted the internal predefined
814                   max timestep is used.
815                   Internally, smaller timesteps may be taken.
816
817        duration: Duration of simulation
818
819        finaltime: Time where simulation should end. This is currently
820        relative time.  So it's the same as duration.
821
822        If both duration and finaltime are given an exception is thrown.
823
824
825        skip_initial_step: Boolean flag that decides whether the first
826        yield step is skipped or not. This is useful for example to avoid
827        duplicate steps when multiple evolve processes are dove tailed.
828
829
830        Evolve is implemented as a generator and is to be called as such, e.g.
831
832        for t in domain.evolve(yieldstep, finaltime):
833            <Do something with domain and t>
834
835
836        All times are given in seconds
837
838        """
839
840        from config import min_timestep, max_timestep, epsilon
841
842        # FIXME: Maybe lump into a larger check prior to evolving
843        msg = 'Boundary tags must be bound to boundary objects before '
844        msg += 'evolving system, '
845        msg += 'e.g. using the method set_boundary.\n'
846        msg += 'This system has the boundary tags %s '\
847               %self.get_boundary_tags()
848        assert hasattr(self, 'boundary_objects'), msg
849
850
851        if yieldstep is None:
852            yieldstep = max_timestep
853        else:
854            yieldstep = float(yieldstep)
855
856        self._order_ = self.default_order
857
858
859        if finaltime is not None and duration is not None:
860            # print 'F', finaltime, duration
861            msg = 'Only one of finaltime and duration may be specified'
862            raise msg
863        else:
864            if finaltime is not None:
865                self.finaltime = float(finaltime)
866            if duration is not None:
867                self.finaltime = self.starttime + float(duration)
868
869
870
871        N = len(self) # Number of triangles
872        self.yieldtime = 0.0 # Track time between 'yields'
873
874        # Initialise interval of timestep sizes (for reporting only)
875        self.min_timestep = max_timestep
876        self.max_timestep = min_timestep
877        self.number_of_steps = 0
878        self.number_of_first_order_steps = 0
879
880
881        # Update ghosts
882        self.update_ghosts()
883
884        # Initial update of vertex and edge values
885        self.distribute_to_vertices_and_edges()
886
887        # Update extrema if necessary (for reporting)
888        self.update_extrema()
889       
890        # Initial update boundary values
891        self.update_boundary()
892
893        # Or maybe restore from latest checkpoint
894        if self.checkpoint is True:
895            self.goto_latest_checkpoint()
896
897        if skip_initial_step is False:
898            yield(self.time)  # Yield initial values
899
900        while True:
901
902            # Evolve One Step, using appropriate timestepping method
903            if self.get_timestepping_method() == 'euler':
904                self.evolve_one_euler_step(yieldstep,finaltime)
905               
906            elif self.get_timestepping_method() == 'rk2':
907                self.evolve_one_rk2_step(yieldstep,finaltime)
908
909            elif self.get_timestepping_method() == 'rk3':
910                self.evolve_one_rk3_step(yieldstep,finaltime)               
911           
912            # Update extrema if necessary (for reporting)
913            self.update_extrema()           
914
915
916            self.yieldtime += self.timestep
917            self.number_of_steps += 1
918            if self._order_ == 1:
919                self.number_of_first_order_steps += 1
920
921            # Yield results
922            if finaltime is not None and self.time >= finaltime-epsilon:
923
924                if self.time > finaltime:
925                    # FIXME (Ole, 30 April 2006): Do we need this check?
926                    # Probably not (Ole, 18 September 2008). Now changed to
927                    # Exception
928                    msg = 'WARNING (domain.py): time overshot finaltime. '
929                    msg += 'Contact Ole.Nielsen@ga.gov.au'
930                    raise Exception, msg
931                   
932
933                # Yield final time and stop
934                self.time = finaltime
935                yield(self.time)
936                break
937
938
939            if self.yieldtime >= yieldstep:
940                # Yield (intermediate) time and allow inspection of domain
941
942                if self.checkpoint is True:
943                    self.store_checkpoint()
944                    self.delete_old_checkpoints()
945
946                # Pass control on to outer loop for more specific actions
947                yield(self.time)
948
949                # Reinitialise
950                self.yieldtime = 0.0
951                self.min_timestep = max_timestep
952                self.max_timestep = min_timestep
953                self.number_of_steps = 0
954                self.number_of_first_order_steps = 0
955                self.max_speed_array = 0.0
956
957
958    def evolve_one_euler_step(self, yieldstep, finaltime):
959        """
960        One Euler Time Step
961        Q^{n+1} = E(h) Q^n
962        """
963
964        # Compute fluxes across each element edge
965        self.compute_fluxes()
966
967        # Update timestep to fit yieldstep and finaltime
968        self.update_timestep(yieldstep, finaltime)
969
970        # Update conserved quantities
971        self.update_conserved_quantities()
972
973        # Update ghosts
974        self.update_ghosts()
975
976        # Update vertex and edge values
977        self.distribute_to_vertices_and_edges()
978
979        # Update boundary values
980        self.update_boundary()
981
982        # Update time
983        self.time += self.timestep
984
985       
986
987
988    def evolve_one_rk2_step(self, yieldstep, finaltime):
989        """
990        One 2nd order RK timestep
991        Q^{n+1} = 0.5 Q^n + 0.5 E(h)^2 Q^n
992        """
993
994        # Save initial conserved quantities values
995        self.backup_conserved_quantities()           
996
997        #--------------------------------------
998        # First euler step
999        #--------------------------------------
1000
1001        # Compute fluxes across each element edge
1002        self.compute_fluxes()
1003
1004        # Update timestep to fit yieldstep and finaltime
1005        self.update_timestep(yieldstep, finaltime)
1006
1007        # Update conserved quantities
1008        self.update_conserved_quantities()
1009
1010        # Update ghosts
1011        self.update_ghosts()
1012
1013        # Update vertex and edge values
1014        self.distribute_to_vertices_and_edges()
1015
1016        # Update boundary values
1017        self.update_boundary()
1018
1019        # Update time
1020        self.time += self.timestep
1021
1022        #------------------------------------
1023        # Second Euler step
1024        #------------------------------------
1025           
1026        # Compute fluxes across each element edge
1027        self.compute_fluxes()
1028
1029        # Update conserved quantities
1030        self.update_conserved_quantities()
1031
1032        #------------------------------------
1033        # Combine initial and final values
1034        # of conserved quantities and cleanup
1035        #------------------------------------
1036       
1037        # Combine steps
1038        self.saxpy_conserved_quantities(0.5, 0.5)
1039
1040        #-----------------------------------
1041        # clean up vertex values
1042        #-----------------------------------
1043 
1044        # Update ghosts
1045        self.update_ghosts()
1046
1047        # Update vertex and edge values
1048        self.distribute_to_vertices_and_edges()
1049
1050        # Update boundary values
1051        self.update_boundary()
1052
1053
1054
1055    def evolve_one_rk3_step(self, yieldstep, finaltime):
1056        """
1057        One 3rd order RK timestep
1058        Q^(1) = 3/4 Q^n + 1/4 E(h)^2 Q^n  (at time t^n + h/2)
1059        Q^{n+1} = 1/3 Q^n + 2/3 E(h) Q^(1) (at time t^{n+1})
1060        """
1061
1062        # Save initial initial conserved quantities values
1063        self.backup_conserved_quantities()           
1064
1065        initial_time = self.time
1066       
1067        #--------------------------------------
1068        # First euler step
1069        #--------------------------------------
1070
1071        # Compute fluxes across each element edge
1072        self.compute_fluxes()
1073
1074        # Update timestep to fit yieldstep and finaltime
1075        self.update_timestep(yieldstep, finaltime)
1076
1077        # Update conserved quantities
1078        self.update_conserved_quantities()
1079
1080        # Update ghosts
1081        self.update_ghosts()
1082
1083        # Update vertex and edge values
1084        self.distribute_to_vertices_and_edges()
1085
1086        # Update boundary values
1087        self.update_boundary()
1088
1089        # Update time
1090        self.time += self.timestep
1091
1092        #------------------------------------
1093        # Second Euler step
1094        #------------------------------------
1095           
1096        # Compute fluxes across each element edge
1097        self.compute_fluxes()
1098
1099        # Update conserved quantities
1100        self.update_conserved_quantities()
1101
1102        #------------------------------------
1103        #Combine steps to obtain intermediate
1104        #solution at time t^n + 0.5 h
1105        #------------------------------------
1106
1107        # Combine steps
1108        self.saxpy_conserved_quantities(0.25, 0.75)
1109 
1110        # Update ghosts
1111        self.update_ghosts()
1112
1113        # Update vertex and edge values
1114        self.distribute_to_vertices_and_edges()
1115
1116        # Update boundary values
1117        self.update_boundary()
1118
1119        # Set substep time
1120        self.time = initial_time + self.timestep*0.5
1121
1122        #------------------------------------
1123        # Third Euler step
1124        #------------------------------------
1125           
1126        # Compute fluxes across each element edge
1127        self.compute_fluxes()
1128
1129        # Update conserved quantities
1130        self.update_conserved_quantities()
1131
1132        #------------------------------------
1133        # Combine final and initial values
1134        # and cleanup
1135        #------------------------------------
1136       
1137        # Combine steps
1138        self.saxpy_conserved_quantities(2.0/3.0, 1.0/3.0)
1139 
1140        # Update ghosts
1141        self.update_ghosts()
1142
1143        # Update vertex and edge values
1144        self.distribute_to_vertices_and_edges()
1145
1146        # Update boundary values
1147        self.update_boundary()
1148
1149        # Set new time
1150        self.time = initial_time + self.timestep       
1151       
1152
1153    def backup_conserved_quantities(self):
1154        N = len(self) # Number_of_triangles
1155
1156        # Backup conserved_quantities centroid values
1157        for name in self.conserved_quantities:
1158            Q = self.quantities[name]
1159            Q.backup_centroid_values()       
1160
1161    def saxpy_conserved_quantities(self,a,b):
1162        N = len(self) #number_of_triangles
1163
1164        # Backup conserved_quantities centroid values
1165        for name in self.conserved_quantities:
1166            Q = self.quantities[name]
1167            Q.saxpy_centroid_values(a,b)       
1168
1169
1170    #==============================
1171    # John Jakeman's old evolve code
1172    #=============================
1173
1174    def evolve_john(self, yieldstep = None, finaltime = None,
1175               skip_initial_step = False):
1176        """Evolve model from time=0.0 to finaltime yielding results
1177        every yieldstep.
1178
1179        Internally, smaller timesteps may be taken.
1180
1181        Evolve is implemented as a generator and is to be called as such, e.g.
1182
1183        for t in domain.evolve(timestep, yieldstep, finaltime):
1184            <Do something with domain and t>
1185
1186        """
1187
1188        from config import min_timestep, max_timestep, epsilon
1189
1190        #FIXME: Maybe lump into a larger check prior to evolving
1191        msg = 'Boundary tags must be bound to boundary objects before evolving system, '
1192        msg += 'e.g. using the method set_boundary.\n'
1193        msg += 'This system has the boundary tags %s ' %self.get_boundary_tags()
1194        assert hasattr(self, 'boundary_objects'), msg
1195
1196        # #self.set_defaults()
1197
1198        if yieldstep is None:
1199            yieldstep = max_timestep
1200        else:
1201            yieldstep = float(yieldstep)
1202
1203        self.order = self.default_order
1204       
1205        self.time_order = self.default_time_order
1206
1207        self.yieldtime = 0.0 #Time between 'yields'
1208
1209        #Initialise interval of timestep sizes (for reporting only)
1210        # SEEMS WIERD
1211        self.min_timestep = max_timestep
1212        self.max_timestep = min_timestep
1213        self.finaltime = finaltime
1214        self.number_of_steps = 0
1215        self.number_of_first_order_steps = 0
1216       
1217        #update ghosts
1218        self.update_ghosts()
1219   
1220        #Initial update of vertex and edge values
1221        self.distribute_to_vertices_and_edges()
1222
1223        #Initial update boundary values
1224        self.update_boundary()
1225
1226        #Or maybe restore from latest checkpoint
1227        if self.checkpoint is True:
1228            self.goto_latest_checkpoint()
1229
1230        if skip_initial_step is False:
1231            yield(self.time)  #Yield initial values   
1232
1233        while True:
1234            if self.time_order == 1:
1235                #Compute fluxes across each element edge
1236                self.compute_fluxes()
1237                #Update timestep to fit yieldstep and finaltime
1238                self.update_timestep(yieldstep, finaltime)
1239                #Compute forcing terms
1240                self.compute_forcing_terms()
1241                #Update conserved quantities
1242                self.update_conserved_quantities(self.timestep)
1243                #update ghosts
1244                #self.update_ghosts()
1245                #Update vertex and edge values
1246                self.distribute_to_vertices_and_edges()               
1247                #Update boundary values
1248                self.update_boundary()
1249
1250            elif self.time_order == 2:
1251               
1252                self.compute_timestep() #self.compute_fluxes() !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1253
1254                #Solve inhomogeneous operator for half a timestep
1255                self.solve_inhomogenous_second_order(yieldstep, finaltime)
1256               
1257                #Solve homogeneous operator for full timestep using
1258                #Harten second order timestepping
1259                self.solve_homogenous_second_order(yieldstep,finaltime)
1260
1261                #Solve inhomogeneous operator for half a timestep
1262                self.solve_inhomogenous_second_order(yieldstep, finaltime)
1263               
1264            #Update time
1265            self.time += self.timestep
1266            self.yieldtime += self.timestep
1267            self.number_of_steps += 1
1268            if self.order == 1:
1269                self.number_of_first_order_steps += 1
1270
1271            #Yield results
1272            if finaltime is not None and abs(self.time - finaltime) < epsilon:
1273
1274                #FIXME: There is a rare situation where the
1275                #final time step is stored twice. Can we make a test?
1276
1277                # Yield final time and stop
1278                yield(self.time)
1279                break
1280
1281
1282            if abs(self.yieldtime - yieldstep) < epsilon:
1283                # Yield (intermediate) time and allow inspection of domain
1284
1285                if self.checkpoint is True:
1286                    self.store_checkpoint()
1287                    self.delete_old_checkpoints()
1288
1289                #Pass control on to outer loop for more specific actions
1290                yield(self.time)
1291
1292                # Reinitialise
1293                self.yieldtime = 0.0
1294                self.min_timestep = max_timestep
1295                self.max_timestep = min_timestep
1296                self.number_of_steps = 0
1297                self.number_of_first_order_steps = 0
1298
1299    def solve_inhomogenous_second_order(self,yieldstep, finaltime):
1300
1301        #Update timestep to fit yieldstep and finaltime
1302        self.update_timestep(yieldstep, finaltime)
1303        #Compute forcing terms
1304        self.compute_forcing_terms()
1305        #Update conserved quantities
1306        self.update_conserved_quantities(0.5*self.timestep)
1307        #Update vertex and edge values
1308        self.distribute_to_vertices_and_edges()
1309        #Update boundary values
1310        self.update_boundary()
1311
1312    def  solve_homogenous_second_order(self,yieldstep,finaltime):
1313        """Use Shu Second order timestepping to update
1314        conserved quantities
1315
1316        q^{n+1/2} = q^{n}+0.5*dt*F^{n}
1317        q^{n+1} = q^{n}+dt*F^{n+1/2}
1318        """
1319        import copy
1320        from Numeric import zeros,Float
1321
1322        N = self.number_of_elements
1323
1324        self.compute_fluxes()
1325        #Update timestep to fit yieldstep and finaltime
1326        self.update_timestep(yieldstep, finaltime)
1327        #Compute forcing terms
1328        #NOT NEEDED FOR 2ND ORDER STRANG SPLIITING
1329        #ADDING THIS WILL NEED TO REMOVE ZEROING IN COMPUTE_FORCING
1330        #self.compute_forcing_terms()
1331
1332        QC = zeros((N,len(self.conserved_quantities)),Float)
1333        QF = zeros((N,len(self.conserved_quantities)),Float)
1334
1335        i = 0
1336        for name in self.conserved_quantities:
1337            Q = self.quantities[name]
1338            #Store the centroid values at time t^n
1339            QC[:,i] = copy.copy(Q.centroid_values)
1340            QF[:,i] = copy.copy(Q.explicit_update)
1341            #Update conserved quantities
1342            Q.update(self.timestep)
1343            i+=1
1344           
1345        #Update vertex and edge values
1346        self.distribute_to_vertices_and_edges()
1347        #Update boundary values
1348        self.update_boundary()
1349       
1350        self.compute_fluxes()
1351        self.update_timestep(yieldstep, finaltime)
1352        #Compute forcing terms
1353        #NOT NEEDED FOR 2ND ORDER STRANG SPLIITING
1354        #self.compute_forcing_terms()
1355
1356        i = 0
1357        for name in self.conserved_quantities:
1358            Q = self.quantities[name]
1359            Q.centroid_values = QC[:,i]
1360            Q.explicit_update = 0.5*(Q.explicit_update+QF[:,i]) 
1361            #Update conserved quantities
1362            Q.update(self.timestep)
1363            i+=1
1364       
1365        #Update vertex and edge values
1366        self.distribute_to_vertices_and_edges()
1367        #Update boundary values
1368        self.update_boundary()
1369
1370    def  solve_homogenous_second_order_harten(self,yieldstep,finaltime):
1371        """Use Harten Second order timestepping to update
1372        conserved quantities
1373
1374        q^{n+1/2} = q^{n}+0.5*dt*F^{n}
1375        q^{n+1} = q^{n}+dt*F^{n+1/2}
1376        """
1377        import copy
1378        from Numeric import zeros,Float
1379
1380        N = self.number_of_elements
1381
1382        self.compute_fluxes()
1383        #Update timestep to fit yieldstep and finaltime
1384        self.update_timestep(yieldstep, finaltime)
1385        #Compute forcing terms
1386        #NOT NEEDED FOR 2ND ORDER STRANG SPLIITING
1387        #ADDING THIS WILL NEED TO REMOVE ZEROING IN COMPUTE_FORCING
1388        #self.compute_forcing_terms()
1389
1390        QC = zeros((N,len(self.conserved_quantities)),Float)
1391
1392        i = 0
1393        for name in self.conserved_quantities:
1394            Q = self.quantities[name]
1395            #Store the centroid values at time t^n
1396            QC[:,i] = copy.copy(Q.centroid_values)
1397            #Update conserved quantities
1398            Q.update(0.5*self.timestep)
1399            i+=1
1400           
1401        #Update vertex and edge values
1402        self.distribute_to_vertices_and_edges()
1403        #Update boundary values
1404        self.update_boundary()
1405       
1406        self.compute_fluxes()
1407        self.update_timestep(yieldstep, finaltime)
1408        #Compute forcing terms
1409        #NOT NEEDED FOR 2ND ORDER STRANG SPLIITING
1410        #self.compute_forcing_terms()
1411
1412        i = 0
1413        for name in self.conserved_quantities:
1414            Q = self.quantities[name]
1415            Q.centroid_values = QC[:,i] 
1416            #Update conserved quantities
1417            Q.update(self.timestep)
1418            i+=1
1419       
1420        #Update vertex and edge values
1421        self.distribute_to_vertices_and_edges()
1422        #Update boundary values
1423        self.update_boundary()
1424       
1425    def distribute_to_vertices_and_edges(self):
1426        """Extrapolate conserved quantities from centroid to
1427        vertices and edge-midpoints for each volume
1428
1429        Default implementation is straight first order,
1430        i.e. constant values throughout each element and
1431        no reference to non-conserved quantities.
1432        """
1433
1434        for name in self.conserved_quantities:
1435            Q = self.quantities[name]
1436            if self.order == 1:
1437                Q.extrapolate_first_order()
1438            elif self.order == 2:
1439                Q.extrapolate_second_order()
1440                #Q.limit()
1441            else:
1442                raise 'Unknown order'
1443            #Q.interpolate_from_vertices_to_edges()
1444
1445
1446    def update_ghosts(self):
1447        pass
1448   
1449    def update_boundary(self):
1450        """Go through list of boundary objects and update boundary values
1451        for all conserved quantities on boundary.
1452        """
1453
1454        #FIXME: Update only those that change (if that can be worked out)
1455        #FIXME: Boundary objects should not include ghost nodes.
1456        #for i, ((vol_id, edge_id), B) in enumerate(self.boundary_objects):
1457        #    q = B.evaluate(vol_id, edge_id)
1458        for i, ((vol_id, vertex_id), B) in enumerate(self.boundary_objects):
1459            q = B.evaluate(vol_id, vertex_id)
1460            #print 'q ',q
1461            for j, name in enumerate(self.evolved_quantities):
1462                #print 'name %s j = %f \n'%(name,j)
1463                Q = self.quantities[name]
1464                Q.boundary_values[i] = q[j]
1465   
1466    def update_timestep(self, yieldstep, finaltime):
1467
1468        from config import min_timestep, max_timestep
1469
1470        # self.timestep is calculated from speed of characteristics
1471        # Apply CFL condition here
1472        timestep = min(self.CFL*self.flux_timestep, max_timestep)
1473
1474        #Record maximal and minimal values of timestep for reporting
1475        self.max_timestep = max(timestep, self.max_timestep)
1476        self.min_timestep = min(timestep, self.min_timestep)
1477
1478        #Protect against degenerate time steps
1479        if timestep < min_timestep:
1480
1481            #Number of consecutive small steps taken b4 taking action
1482            self.smallsteps += 1
1483
1484            if self.smallsteps > self.max_smallsteps:
1485                self.smallsteps = 0 #Reset
1486
1487                if self.order == 1:
1488                    msg = 'WARNING: Too small timestep %.16f reached '\
1489                          %timestep
1490                    msg += 'even after %d steps of 1 order scheme'\
1491                           %self.max_smallsteps
1492                    print msg
1493                    timestep = min_timestep  #Try enforcing min_step
1494
1495                    #raise msg
1496                else:
1497                    #Try to overcome situation by switching to 1 order
1498                    print "changing Order 1"
1499                    self.order = 1
1500
1501        else:
1502            self.smallsteps = 0
1503            if self.order == 1 and self.default_order == 2:
1504                self.order = 2
1505
1506
1507        #Ensure that final time is not exceeded
1508        if finaltime is not None and self.time + timestep > finaltime:
1509            timestep = finaltime-self.time
1510
1511        #Ensure that model time is aligned with yieldsteps
1512        if self.yieldtime + timestep > yieldstep:
1513            timestep = yieldstep-self.yieldtime
1514
1515        self.timestep = timestep
1516
1517    def update_extrema(self):
1518        pass
1519
1520    def compute_forcing_terms(self):
1521        """If there are any forcing functions driving the system
1522        they should be defined in Domain subclass and appended to
1523        the list self.forcing_terms
1524        """
1525        #Clears explicit_update needed for second order method
1526        if self.time_order == 2:
1527            for name in self.conserved_quantities:
1528                Q = self.quantities[name]
1529                Q.explicit_update[:] = 0.0
1530
1531        for f in self.forcing_terms:
1532            f(self)
1533           
1534
1535    def update_derived_quantites(self):
1536        pass
1537   
1538    #def update_conserved_quantities(self):
1539    def update_conserved_quantities(self):
1540        """Update vectors of conserved quantities using previously
1541        computed fluxes specified forcing functions.
1542        """
1543
1544        from Numeric import ones, sum, equal, Float
1545
1546        N = self.number_of_elements
1547        d = len(self.conserved_quantities)
1548
1549        timestep = self.timestep
1550
1551        #Compute forcing terms
1552        self.compute_forcing_terms()
1553
1554        #Update conserved_quantities
1555        for name in self.conserved_quantities:
1556            Q = self.quantities[name]
1557            Q.update(timestep)
1558
1559
1560
1561if __name__ == "__main__":
1562
1563    points1 = [0.0, 1.0, 2.0, 3.0]
1564    D1 = Domain(points1)
1565
1566    print D1.get_coordinate(0)
1567    print D1.get_coordinate(0,1)
1568    print 'Number of Elements = ',D1.number_of_elements
1569
1570    try:
1571        print D1.get_coordinate(3)
1572    except:
1573        pass
1574    else:
1575        msg =  'Should have raised an out of bounds exception'
1576        raise msg
1577
1578    #points2 = [0.0, 1.0, 2.0, 3.0, 2.5]
1579    #D2 = Domain(points2)
Note: See TracBrowser for help on using the repository browser.