source: trunk/anuga_work/development/sudi/sw_1d/avalanche/B_momentum/domain.py @ 7915

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