source: trunk/anuga_work/development/sudi/sw_1d/shock_detector/shm/domain.py @ 7909

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

Add 1d codes.

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