source: trunk/anuga_work/development/sudi/sw_1d/periodic_waves/cg/domain.py @ 7922

Last change on this file since 7922 was 7922, checked in by mungkasi, 14 years ago

Modifying codes so that arrays are compatible with numpy instead of Numeric.

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