# source:anuga_work/development/shallow_water_1d_old/domain_order2.py@5533

Last change on this file since 5533 was 4959, checked in by steve, 16 years ago

Copied John Jakeman's and Sudi Mungkasi's 1d shallow water code

File size: 36.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 *
10from coordinate_transforms.geo_reference import Geo_reference
11
12class Domain:
13
14    def __init__(self, coordinates, boundary = None,
15                 conserved_quantities = None, other_quantities = None,
16                 tagged_elements = None, geo_reference = None):
17        """
18        Build 1D elements from x coordinates
19        """
20
21        from Numeric import array, zeros, Float, Int
22
23        #Store Points
24        self.coordinates = array(coordinates)
25
26        if geo_reference is None:
27            self.geo_reference = 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        if self.time_order == 2:
763            self.compute_timestep()
764
765        while True:
766            if self.time_order == 1:
767                #Compute fluxes across each element edge
768                self.compute_fluxes()
769                #Update timestep to fit yieldstep and finaltime
770                self.update_timestep(yieldstep, finaltime)
771                #Compute forcing terms
772                self.compute_forcing_terms()
773                #Update conserved quantities
774                self.update_conserved_quantities(self.timestep)
775                #update ghosts
776                #self.update_ghosts()
777                #Update vertex and edge values
778                self.distribute_to_vertices_and_edges()
779                #Update boundary values
780                self.update_boundary()
781
782            elif self.time_order == 2:
783                if self.flux_first:
784                    #Solve homogeneous operator for half a timestep
785                    self.compute_fluxes()
786                    self.update_timestep(yieldstep, finaltime)
787                    self.update_conserved_quantities(0.5*self.timestep)
788                    self.distribute_to_vertices_and_edges()
789                    self.update_boundary()
790
791                    #Solve inhomogeneous operator for full timestep
792                    self.compute_forcing_terms()
793                    self.update_conserved_quantities(self.timestep)
794                    self.distribute_to_vertices_and_edges()
795                    self.update_boundary()
796
797                    #Solve homogeneous operator for half a timestep
798                    self.compute_fluxes()
799                    #self.update_timestep(yieldstep, finaltime)
800                    self.update_conserved_quantities(0.5*self.timestep)
801                    self.distribute_to_vertices_and_edges()
802                    self.update_boundary()
803
804                elif not self.flux_first:
805                    #Solve inhomogeneous operator for half a timestep
806                    #self.compute_timestep()
807                    self.update_timestep(yieldstep, finaltime) #note adding this line improves performance
808                    #TEMP FIX: Need to develop better way to calculate timestep
809                    #This uses previous timestep as new timestep
810                    self.compute_forcing_terms()
811                    self.update_conserved_quantities(0.5*self.timestep)
812                    self.distribute_to_vertices_and_edges()
813                    self.update_boundary()
814
815                    #Solve homogeneous operator for full timestep
816                    #Compute fluxes across each element edge
817                    self.compute_fluxes()
818                    #Update timestep to fit yieldstep and finaltime
819                    self.update_timestep(yieldstep, finaltime)
820                    #Update conserved quantities
821                    self.update_conserved_quantities(self.timestep)
822                    #Update vertex and edge values
823                    self.distribute_to_vertices_and_edges()
824                    #Update boundary values
825                    self.update_boundary()
826
827
828                    ##Solve inhomogeneous operator for half a timestep
829                    self.compute_forcing_terms()
830                    self.update_conserved_quantities(0.5*self.timestep)
831                    self.distribute_to_vertices_and_edges()
832                    self.update_boundary()
833
834
835            #print 'timestep', self.timestep
836
837            #Update time
838            self.time += self.timestep
839            self.yieldtime += self.timestep
840            self.number_of_steps += 1
841            if self.order == 1:
842                self.number_of_first_order_steps += 1
843
844            #Yield results
845            if finaltime is not None and abs(self.time - finaltime) < epsilon:
846
847                #FIXME: There is a rare situation where the
848                #final time step is stored twice. Can we make a test?
849
850                # Yield final time and stop
851                yield(self.time)
852                break
853
854
855            if abs(self.yieldtime - yieldstep) < epsilon:
856                # Yield (intermediate) time and allow inspection of domain
857
858                if self.checkpoint is True:
859                    self.store_checkpoint()
860                    self.delete_old_checkpoints()
861
862                #Pass control on to outer loop for more specific actions
863                yield(self.time)
864
865                # Reinitialise
866                self.yieldtime = 0.0
867                self.min_timestep = max_timestep
868                self.max_timestep = min_timestep
869                self.number_of_steps = 0
870                self.number_of_first_order_steps = 0
871
872
873    def distribute_to_vertices_and_edges(self):
874        """Extrapolate conserved quantities from centroid to
875        vertices and edge-midpoints for each volume
876
877        Default implementation is straight first order,
878        i.e. constant values throughout each element and
879        no reference to non-conserved quantities.
880        """
881
882        for name in self.conserved_quantities:
883            Q = self.quantities[name]
884            if self.order == 1:
885                Q.extrapolate_first_order()
886            elif self.order == 2:
887                Q.extrapolate_second_order()
888                #Q.limit()
889            else:
890                raise 'Unknown order'
891            #Q.interpolate_from_vertices_to_edges()
892
893
894    def update_boundary(self):
895        """Go through list of boundary objects and update boundary values
896        for all conserved quantities on boundary.
897        """
898
899        #FIXME: Update only those that change (if that can be worked out)
900        #FIXME: Boundary objects should not include ghost nodes.
901        #for i, ((vol_id, edge_id), B) in enumerate(self.boundary_objects):
902        #    q = B.evaluate(vol_id, edge_id)
903        for i, ((vol_id, vertex_id), B) in enumerate(self.boundary_objects):
904            q = B.evaluate(vol_id, vertex_id)
905
906            for j, name in enumerate(self.conserved_quantities):
907                Q = self.quantities[name]
908                Q.boundary_values[i] = q[j]
909
910    def update_timestep(self, yieldstep, finaltime):
911
912        from config import min_timestep
913
914        # self.timestep is calculated from speed of characteristics
915        # Apply CFL condition here
916        timestep = self.CFL*self.timestep
917
918        #Record maximal and minimal values of timestep for reporting
919        self.max_timestep = max(timestep, self.max_timestep)
920        self.min_timestep = min(timestep, self.min_timestep)
921
922        #Protect against degenerate time steps
923        if timestep < min_timestep:
924
925            #Number of consecutive small steps taken b4 taking action
926            self.smallsteps += 1
927
928            if self.smallsteps > self.max_smallsteps:
929                self.smallsteps = 0 #Reset
930
931                if self.order == 1:
932                    msg = 'WARNING: Too small timestep %.16f reached '\
933                          %timestep
934                    msg += 'even after %d steps of 1 order scheme'\
935                           %self.max_smallsteps
936                    print msg
937                    timestep = min_timestep  #Try enforcing min_step
938
939                    #raise msg
940                else:
941                    #Try to overcome situation by switching to 1 order
942                    print "changing Order 1"
943                    self.order = 1
944
945        else:
946            self.smallsteps = 0
947            if self.order == 1 and self.default_order == 2:
948                self.order = 2
949
950
951        #Ensure that final time is not exceeded
952        if finaltime is not None and self.time + timestep > finaltime:
953            timestep = finaltime-self.time
954
955        #Ensure that model time is aligned with yieldsteps
956        if self.yieldtime + timestep > yieldstep:
957            timestep = yieldstep-self.yieldtime
958
959        self.timestep = timestep
960
961
962    def compute_forcing_terms(self):
963        """If there are any forcing functions driving the system
964        they should be defined in Domain subclass and appended to
965        the list self.forcing_terms
966        """
967        #Clears explicit_update needed for second order method
968        if self.time_order == 2:
969            for name in self.conserved_quantities:
970                Q = self.quantities[name]
971                Q.explicit_update[:] = 0.0
972
973        for f in self.forcing_terms:
974            f(self)
975
976
977    def update_derived_quantites(self):
978        pass
979
980    #def update_conserved_quantities(self):
981    def update_conserved_quantities(self,timestep):
982        """Update vectors of conserved quantities using previously
983        computed fluxes specified forcing functions.
984        """
985
986        from Numeric import ones, sum, equal, Float
987
988        N = self.number_of_elements
989        d = len(self.conserved_quantities)
990
991        #timestep = self.timestep
992
993        #Compute forcing terms
994        #self.compute_forcing_terms()
995
996        #Update conserved_quantities
997        for name in self.conserved_quantities:
998            Q = self.quantities[name]
999            Q.update(timestep)
1000
1001            #Clean up
1002            #Note that Q.explicit_update is reset by compute_fluxes
1003
1004            #MH090605 commented out the following since semi_implicit_update is now re-initialized
1005            #at the end of the _update function in quantity_ext.c (This is called by the
1006            #preceeding Q.update(timestep) statement above).
1007            #For run_profile.py with N=128, the time of update_conserved_quantities is cut from 14.00 secs
1008            #to 8.35 secs
1009
1010            #Q.semi_implicit_update[:] = 0.0
1011
1012if __name__ == "__main__":
1013
1014    points1 = [0.0, 1.0, 2.0, 3.0]
1015    D1 = Domain(points1)
1016
1017    print D1.get_coordinate(0)
1018    print D1.get_coordinate(0,1)
1019    print 'Number of Elements = ',D1.number_of_elements
1020
1021    try:
1022        print D1.get_coordinate(3)
1023    except:
1024        pass
1025    else:
1026        msg =  'Should have raised an out of bounds exception'
1027        raise msg
1028
1029    #points2 = [0.0, 1.0, 2.0, 3.0, 2.5]
1030    #D2 = Domain(points2)
1031
Note: See TracBrowser for help on using the repository browser.