# source:anuga_work/development/shallow_water_1d/domain_t2.py@4960

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