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

Last change on this file since 7837 was 7837, checked in by mungkasi, 12 years ago

Again, adding some codes for 1d problems on debris avalanche and periodic waves.

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