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

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

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

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