source: anuga_work/development/sudi/sw_1d/periodic_waves/johns/domain.py @ 7837

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

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

File size: 44.6 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'
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    def evolve_one_rk3_step(self, yieldstep, finaltime):
946        """
947        One 3rd order RK timestep
948        Q^(1) = 3/4 Q^n + 1/4 E(h)^2 Q^n  (at time t^n + h/2)
949        Q^{n+1} = 1/3 Q^n + 2/3 E(h) Q^(1) (at time t^{n+1})
950        """
951
952        # Save initial initial conserved quantities values
953        self.backup_conserved_quantities()           
954
955        initial_time = self.time
956       
957        #--------------------------------------
958        # First euler step
959        #--------------------------------------
960
961        # Compute fluxes across each element edge
962        self.compute_fluxes()
963
964        # Update timestep to fit yieldstep and finaltime
965        self.update_timestep(yieldstep, finaltime)
966
967        # Update conserved quantities
968        self.update_conserved_quantities()
969
970        # Update ghosts
971        self.update_ghosts()
972
973        # Update vertex and edge values
974        self.distribute_to_vertices_and_edges()
975
976        # Update boundary values
977        self.update_boundary()
978
979        # Update time
980        self.time += self.timestep
981
982        #------------------------------------
983        # Second Euler step
984        #------------------------------------
985           
986        # Compute fluxes across each element edge
987        self.compute_fluxes()
988
989        # Update conserved quantities
990        self.update_conserved_quantities()
991
992        #------------------------------------
993        #Combine steps to obtain intermediate
994        #solution at time t^n + 0.5 h
995        #------------------------------------
996
997        # Combine steps
998        self.saxpy_conserved_quantities(0.25, 0.75)
999 
1000        # Update ghosts
1001        self.update_ghosts()
1002
1003        # Update vertex and edge values
1004        self.distribute_to_vertices_and_edges()
1005
1006        # Update boundary values
1007        self.update_boundary()
1008
1009        # Set substep time
1010        self.time = initial_time + self.timestep*0.5
1011
1012        #------------------------------------
1013        # Third Euler step
1014        #------------------------------------
1015           
1016        # Compute fluxes across each element edge
1017        self.compute_fluxes()
1018
1019        # Update conserved quantities
1020        self.update_conserved_quantities()
1021
1022        #------------------------------------
1023        # Combine final and initial values
1024        # and cleanup
1025        #------------------------------------
1026       
1027        # Combine steps
1028        self.saxpy_conserved_quantities(2.0/3.0, 1.0/3.0)
1029 
1030        # Update ghosts
1031        self.update_ghosts()
1032
1033        # Update vertex and edge values
1034        self.distribute_to_vertices_and_edges()
1035
1036        # Update boundary values
1037        self.update_boundary()
1038
1039        # Set new time
1040        self.time = initial_time + self.timestep       
1041       
1042
1043    def backup_conserved_quantities(self):
1044        N = len(self) # Number_of_triangles
1045
1046        # Backup conserved_quantities centroid values
1047        for name in self.conserved_quantities:
1048            Q = self.quantities[name]
1049            Q.backup_centroid_values()       
1050
1051    def saxpy_conserved_quantities(self,a,b):
1052        N = len(self) #number_of_triangles
1053
1054        # Backup conserved_quantities centroid values
1055        for name in self.conserved_quantities:
1056            Q = self.quantities[name]
1057            Q.saxpy_centroid_values(a,b)       
1058   
1059    def solve_inhomogenous_second_order(self,yieldstep, finaltime):
1060
1061        #Update timestep to fit yieldstep and finaltime
1062        self.update_timestep(yieldstep, finaltime)
1063        #Compute forcing terms
1064        self.compute_forcing_terms()
1065        #Update conserved quantities
1066        self.update_conserved_quantities(0.5*self.timestep)
1067        #Update vertex and edge values
1068        self.distribute_to_vertices_and_edges()
1069        #Update boundary values
1070        self.update_boundary()
1071
1072    def  solve_homogenous_second_order(self,yieldstep,finaltime):
1073        """Use Shu Second order timestepping to update
1074        conserved quantities
1075
1076        q^{n+1/2} = q^{n}+0.5*dt*F^{n}
1077        q^{n+1} = q^{n}+dt*F^{n+1/2}
1078        """
1079        import copy
1080        from Numeric import zeros,Float
1081
1082        N = self.number_of_elements
1083
1084        self.compute_fluxes()
1085        #Update timestep to fit yieldstep and finaltime
1086        self.update_timestep(yieldstep, finaltime)
1087        #Compute forcing terms
1088        #NOT NEEDED FOR 2ND ORDER STRANG SPLIITING
1089        #ADDING THIS WILL NEED TO REMOVE ZEROING IN COMPUTE_FORCING
1090        #self.compute_forcing_terms()
1091
1092        QC = zeros((N,len(self.conserved_quantities)),Float)
1093        QF = zeros((N,len(self.conserved_quantities)),Float)
1094
1095        i = 0
1096        for name in self.conserved_quantities:
1097            Q = self.quantities[name]
1098            #Store the centroid values at time t^n
1099            QC[:,i] = copy.copy(Q.centroid_values)
1100            QF[:,i] = copy.copy(Q.explicit_update)
1101            #Update conserved quantities
1102            Q.update(self.timestep)
1103            i+=1
1104           
1105        #Update vertex and edge values
1106        self.distribute_to_vertices_and_edges()
1107        #Update boundary values
1108        self.update_boundary()
1109       
1110        self.compute_fluxes()
1111        self.update_timestep(yieldstep, finaltime)
1112        #Compute forcing terms
1113        #NOT NEEDED FOR 2ND ORDER STRANG SPLIITING
1114        #self.compute_forcing_terms()
1115
1116        i = 0
1117        for name in self.conserved_quantities:
1118            Q = self.quantities[name]
1119            Q.centroid_values = QC[:,i]
1120            Q.explicit_update = 0.5*(Q.explicit_update+QF[:,i]) 
1121            #Update conserved quantities
1122            Q.update(self.timestep)
1123            i+=1
1124       
1125        #Update vertex and edge values
1126        self.distribute_to_vertices_and_edges()
1127        #Update boundary values
1128        self.update_boundary()
1129
1130    def  solve_homogenous_second_order_harten(self,yieldstep,finaltime):
1131        """Use Harten Second order timestepping to update
1132        conserved quantities
1133
1134        q^{n+1/2} = q^{n}+0.5*dt*F^{n}
1135        q^{n+1} = q^{n}+dt*F^{n+1/2}
1136        """
1137        import copy
1138        from Numeric import zeros,Float
1139
1140        N = self.number_of_elements
1141
1142        self.compute_fluxes()
1143        #Update timestep to fit yieldstep and finaltime
1144        self.update_timestep(yieldstep, finaltime)
1145        #Compute forcing terms
1146        #NOT NEEDED FOR 2ND ORDER STRANG SPLIITING
1147        #ADDING THIS WILL NEED TO REMOVE ZEROING IN COMPUTE_FORCING
1148        #self.compute_forcing_terms()
1149
1150        QC = zeros((N,len(self.conserved_quantities)),Float)
1151
1152        i = 0
1153        for name in self.conserved_quantities:
1154            Q = self.quantities[name]
1155            #Store the centroid values at time t^n
1156            QC[:,i] = copy.copy(Q.centroid_values)
1157            #Update conserved quantities
1158            Q.update(0.5*self.timestep)
1159            i+=1
1160           
1161        #Update vertex and edge values
1162        self.distribute_to_vertices_and_edges()
1163        #Update boundary values
1164        self.update_boundary()
1165       
1166        self.compute_fluxes()
1167        self.update_timestep(yieldstep, finaltime)
1168        #Compute forcing terms
1169        #NOT NEEDED FOR 2ND ORDER STRANG SPLIITING
1170        #self.compute_forcing_terms()
1171
1172        i = 0
1173        for name in self.conserved_quantities:
1174            Q = self.quantities[name]
1175            Q.centroid_values = QC[:,i] 
1176            #Update conserved quantities
1177            Q.update(self.timestep)
1178            i+=1
1179       
1180        #Update vertex and edge values
1181        self.distribute_to_vertices_and_edges()
1182        #Update boundary values
1183        self.update_boundary()
1184       
1185    def distribute_to_vertices_and_edges(self):
1186        """Extrapolate conserved quantities from centroid to
1187        vertices and edge-midpoints for each volume
1188
1189        Default implementation is straight first order,
1190        i.e. constant values throughout each element and
1191        no reference to non-conserved quantities.
1192        """
1193
1194        for name in self.conserved_quantities:
1195            Q = self.quantities[name]
1196            if self.order == 1:
1197                Q.extrapolate_first_order()
1198            elif self.order == 2:
1199                Q.extrapolate_second_order()
1200                #Q.limit()
1201            else:
1202                raise 'Unknown order'
1203
1204
1205    def update_ghosts(self):
1206        #pass
1207        from parameter import prescribe_at_O_JOHNS, h_0, L
1208        from config import g
1209        from math import sqrt
1210        def f_JOHNS(t):
1211            timing = t*sqrt(g*h_0)/L
1212            w, u = prescribe_at_O_JOHNS(timing)
1213            wO = w*h_0
1214            uO = u*sqrt(g*h_0)
1215            zO = -h_0
1216            hO = wO - zO
1217            pO = uO * hO
1218            #[    'stage', 'xmomentum', 'elevation', 'height', 'velocity']
1219            return [wO,  pO,   zO,      hO,      uO]
1220        tempor = f_JOHNS(self.time)
1221       
1222        self.quantities['stage'].centroid_values[0]     = tempor[0]
1223        self.quantities['xmomentum'].centroid_values[0] = tempor[1]
1224        self.quantities['elevation'].centroid_values[0] = tempor[2]
1225        self.quantities['height'].centroid_values[0]    = tempor[3]
1226        self.quantities['velocity'].centroid_values[0]  = tempor[4]
1227       
1228       
1229    def update_boundary(self):
1230        """Go through list of boundary objects and update boundary values
1231        for all conserved quantities on boundary.
1232        """
1233       
1234        for i, ((vol_id, vertex_id), B) in enumerate(self.boundary_objects):
1235            q = B.evaluate(vol_id, vertex_id)
1236            for j, name in enumerate(self.evolved_quantities):
1237                #print 'name %s j = %f \n'%(name,j)
1238                Q = self.quantities[name]
1239                Q.boundary_values[i] = q[j]
1240               
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,1] = VV[0,0] = BV[0]
1248            VV[n-1,1] = BV[1]
1249       
1250        """
1251        #Below is just for checking/testing
1252        self.quantities['stage'].vertex_values[0,1]     = self.quantities['stage'].vertex_values[0,0]
1253        self.quantities['xmomentum'].vertex_values[0,1] = self.quantities['xmomentum'].vertex_values[0,0]
1254        self.quantities['elevation'].vertex_values[0,1] = self.quantities['elevation'].vertex_values[0,0]
1255        self.quantities['height'].vertex_values[0,1]    = self.quantities['height'].vertex_values[0,0]
1256        self.quantities['velocity'].vertex_values[0,1]  = self.quantities['velocity'].vertex_values[0,0]
1257        """         
1258   
1259    def update_timestep(self, yieldstep, finaltime):
1260
1261        from config import min_timestep, max_timestep
1262
1263        # self.timestep is calculated from speed of characteristics
1264        # Apply CFL condition here
1265        timestep = min(self.CFL*self.flux_timestep, max_timestep)
1266
1267        #Record maximal and minimal values of timestep for reporting
1268        self.max_timestep = max(timestep, self.max_timestep)
1269        self.min_timestep = min(timestep, self.min_timestep)
1270
1271        #Protect against degenerate time steps
1272        if timestep < min_timestep:
1273
1274            #Number of consecutive small steps taken b4 taking action
1275            self.smallsteps += 1
1276
1277            if self.smallsteps > self.max_smallsteps:
1278                self.smallsteps = 0 #Reset
1279
1280                if self.order == 1:
1281                    msg = 'WARNING: Too small timestep %.16f reached '\
1282                          %timestep
1283                    msg += 'even after %d steps of 1 order scheme'\
1284                           %self.max_smallsteps
1285                    print msg
1286                    timestep = min_timestep  #Try enforcing min_step
1287
1288                    #raise msg
1289                else:
1290                    #Try to overcome situation by switching to 1 order
1291                    print "changing Order 1"
1292                    self.order = 1
1293
1294        else:
1295            self.smallsteps = 0
1296            if self.order == 1 and self.default_order == 2:
1297                self.order = 2
1298
1299
1300        #Ensure that final time is not exceeded
1301        if finaltime is not None and self.time + timestep > finaltime:
1302            timestep = finaltime-self.time
1303
1304        #Ensure that model time is aligned with yieldsteps
1305        if self.yieldtime + timestep > yieldstep:
1306            timestep = yieldstep-self.yieldtime
1307
1308        self.timestep = timestep
1309
1310    def update_extrema(self):
1311        pass
1312
1313    def compute_forcing_terms(self):
1314        """If there are any forcing functions driving the system
1315        they should be defined in Domain subclass and appended to
1316        the list self.forcing_terms
1317        """
1318        #Clears explicit_update needed for second order method
1319        if self.time_order == 2:
1320            for name in self.conserved_quantities:
1321                Q = self.quantities[name]
1322                Q.explicit_update[:] = 0.0
1323
1324        for f in self.forcing_terms:
1325            f(self)
1326           
1327
1328    def update_derived_quantites(self):
1329        pass
1330   
1331    #def update_conserved_quantities(self):
1332    def update_conserved_quantities(self):
1333        """Update vectors of conserved quantities using previously
1334        computed fluxes specified forcing functions.
1335        """
1336
1337        from Numeric import ones, sum, equal, Float
1338
1339        N = self.number_of_elements
1340        d = len(self.conserved_quantities)
1341
1342        timestep = self.timestep
1343
1344        #Compute forcing terms
1345        self.compute_forcing_terms()
1346
1347        #Update conserved_quantities
1348        for name in self.conserved_quantities:
1349            Q = self.quantities[name]
1350            Q.update(timestep)
1351
1352
1353
1354if __name__ == "__main__":
1355
1356    points1 = [0.0, 1.0, 2.0, 3.0]
1357    D1 = Domain(points1)
1358
1359    print D1.get_coordinate(0)
1360    print D1.get_coordinate(0,1)
1361    print 'Number of Elements = ',D1.number_of_elements
1362
1363    try:
1364        print D1.get_coordinate(3)
1365    except:
1366        pass
1367    else:
1368        msg =  'Should have raised an out of bounds exception'
1369        raise msg
Note: See TracBrowser for help on using the repository browser.