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

Last change on this file since 7793 was 7793, checked in by steve, 14 years ago

Setting up test function for Quantity

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