source: trunk/anuga_work/development/2010-projects/anuga_1d/base/quantity.py @ 7884

Last change on this file since 7884 was 7884, checked in by steve, 12 years ago

Moving 2010 project

File size: 35.8 KB
Line 
1"""
2Class Quantity - Implements values at each 1d element
3
4To create:
5
6   Quantity(domain, vertex_values)
7
8   domain: Associated domain structure. Required.
9
10   vertex_values: N x 2 array of values at each vertex for each element.
11                  Default None
12
13   If vertex_values are None Create array of zeros compatible with domain.
14   Otherwise check that it is compatible with dimenions of domain.
15   Otherwise raise an exception
16"""
17
18import numpy
19
20
21
22class Quantity:
23
24   
25    def __init__(self, domain, vertex_values=None):
26        #Initialise Quantity using optional vertex values.
27       
28        from anuga_1d.base.generic_domain import Generic_domain
29
30        msg = 'First argument in Quantity.__init__ '
31        msg += 'must be of class Generic_domain (or a subclass thereof)'
32        assert isinstance(domain, Generic_domain), msg
33
34        if vertex_values is None:
35            N = domain.number_of_elements
36            self.vertex_values = numpy.zeros((N, 2), numpy.float)
37        else:
38            self.vertex_values = numpy.array(vertex_values, numpy.float)
39
40            N, V = self.vertex_values.shape
41            assert V == 2,\
42                   'Two vertex values per element must be specified'
43
44
45            msg = 'Number of vertex values (%d) must be consistent with'\
46                  %N
47            msg += 'number of elements in specified domain (%d).'\
48                   %domain.number_of_elements
49
50            assert N == domain.number_of_elements, msg
51
52        self.domain = domain
53
54        #Allocate space for other quantities
55        self.centroid_values = numpy.zeros(N, numpy.float)
56        self.centroid_backup_values = numpy.zeros(N, numpy.float)
57        #self.edge_values = numpy.zeros((N, 2), numpy.float)
58        #edge values are values of the ends of each interval
59     
60        #Intialise centroid values
61        self.interpolate()
62
63
64       
65        #Allocate space for boundary values
66        #L = len(domain.boundary)
67        self.boundary_values = numpy.zeros(2, numpy.float) #assumes no parrellism
68
69        #Allocate space for updates of conserved quantities by
70        #flux calculations and forcing functions
71       
72        self.N = domain.number_of_elements
73        N = self.N
74        self.explicit_update = numpy.zeros(N, numpy.float )
75        self.semi_implicit_update = numpy.zeros(N, numpy.float )
76
77        self.gradients = numpy.zeros(N, numpy.float)
78        self.qmax = numpy.zeros(self.centroid_values.shape, numpy.float)
79        self.qmin = numpy.zeros(self.centroid_values.shape, numpy.float)
80
81        #These are taken from domain but can be set for each quantity
82        # if so desired
83        self.beta = 0.0
84        self.limiter = 'vanleer'
85
86
87        self.beta_p = numpy.zeros(N,numpy.float)
88        self.beta_m = numpy.zeros(N,numpy.float)
89        self.beta_x = numpy.zeros(N,numpy.float)
90
91
92        self.dx  = numpy.zeros((N,2), numpy.float)
93        self.phi = numpy.zeros(N, numpy.float)
94
95
96    def __len__(self):
97        """
98        Returns number of intervals.
99        """
100        return self.centroid_values.shape[0]
101
102    def __neg__(self):
103        """Negate all values in this quantity giving meaning to the
104        expression -Q where Q is an instance of class Quantity
105        """
106
107        Q = Quantity(self.domain)
108        Q.set_values_from_numeric(-self.vertex_values)
109        return Q
110
111
112
113    def __add__(self, other):
114        """Add to self anything that could populate a quantity
115
116        E.g other can be a constant, an array, a function, another quantity
117        (except for a filename or points, attributes (for now))
118        - see set_values_from_numeric for details
119        """
120
121        Q = Quantity(self.domain)
122        Q.set_values_from_numeric(other)
123
124        result = Quantity(self.domain)
125        result.set_values_from_numeric(self.vertex_values + Q.vertex_values)
126        return result
127
128    def __radd__(self, other):
129        """Handle cases like 7+Q, where Q is an instance of class Quantity
130        """
131
132        return self + other
133
134    def __sub__(self, other):
135        return self + -other            # Invoke self.__neg__()
136
137    def __mul__(self, other):
138        """Multiply self with anything that could populate a quantity
139
140        E.g other can be a constant, an array, a function, another quantity
141        (except for a filename or points, attributes (for now))
142        - see set_values_from_numeric for details
143        """
144
145        if isinstance(other, Quantity):
146            Q = other
147        else:
148            Q = Quantity(self.domain)
149            Q.set_values_from_numeric(other)
150
151        result = Quantity(self.domain)
152
153        # The product of vertex_values, edge_values and centroid_values
154        # are calculated and assigned directly without using
155        # set_values_from_numeric (which calls interpolate). Otherwise
156        # centroid values wouldn't be products from q1 and q2
157        result.vertex_values = self.vertex_values * Q.vertex_values
158        result.centroid_values = self.centroid_values * Q.centroid_values
159
160        return result
161
162    def __rmul__(self, other):
163        """Handle cases like 3*Q, where Q is an instance of class Quantity
164        """
165
166        return self * other
167
168    def __div__(self, other):
169        """Divide self with anything that could populate a quantity
170
171        E.g other can be a constant, an array, a function, another quantity
172        (except for a filename or points, attributes (for now))
173        - see set_values_from_numeric for details
174
175        Zero division is dealt with by adding an epsilon to the divisore
176        FIXME (Ole): Replace this with native INF once we migrate to NumPy
177        """
178
179        if isinstance(other, Quantity):
180            Q = other
181        else:
182            Q = Quantity(self.domain)
183            Q.set_values_from_numeric(other)
184
185        result = Quantity(self.domain)
186
187        # The quotient of vertex_values, edge_values and centroid_values
188        # are calculated and assigned directly without using
189        # set_values_from_numeric (which calls interpolate). Otherwise
190        # centroid values wouldn't be quotient of q1 and q2
191        result.vertex_values = self.vertex_values/(Q.vertex_values + epsilon)
192        result.centroid_values = self.centroid_values/(Q.centroid_values + epsilon)
193
194        return result
195
196    def __rdiv__(self, other):
197        """Handle cases like 3/Q, where Q is an instance of class Quantity
198        """
199
200        return self / other
201
202    def __pow__(self, other):
203        """Raise quantity to (numerical) power
204
205        As with __mul__ vertex values are processed entry by entry
206        while centroid and edge values are re-interpolated.
207
208        Example using __pow__:
209          Q = (Q1**2 + Q2**2)**0.5
210        """
211
212        if isinstance(other, Quantity):
213            Q = other
214        else:
215            Q = Quantity(self.domain)
216            Q.set_values_from_numeric(other)
217
218        result = Quantity(self.domain)
219
220        # The power of vertex_values, edge_values and centroid_values
221        # are calculated and assigned directly without using
222        # set_values_from_numeric (which calls interpolate). Otherwise
223        # centroid values wouldn't be correct
224        result.vertex_values = self.vertex_values ** other
225        result.centroid_values = self.centroid_values ** other
226
227        return result
228
229    def set_values_from_numeric(self, numeric):
230       
231        x = numpy.array([1.0,2.0])
232        y = [1.0,2.0]
233
234        if type(numeric) == type(y):
235            self.set_values_from_array(numeric)
236        elif type(numeric) == type(x):
237            self.set_values_from_array(numeric)
238        elif callable(numeric):
239            self.set_values_from_function(numeric)
240        elif isinstance(numeric, Quantity):
241            self.set_values_from_quantity(numeric)
242        else:   # see if it's coercible to a float (float, int or long, etc)
243            try:
244                numeric = float(numeric)
245            except ValueError:
246                msg = ("Illegal type for variable 'numeric': %s" % type(numeric))
247                raise Exception(msg)
248            self.set_values_from_constant(numeric)
249
250    def set_values_from_constant(self,numeric):
251
252        self.vertex_values[:,:]   = numeric
253        self.centroid_values[:,] = numeric
254
255
256    def set_values_from_array(self,numeric):
257
258        self.vertex_values[:,:]   = numeric
259        self.interpolate()
260
261
262    def set_values_from_quantity(self,quantity):
263
264        self.vertex_values[:,:]   = quantity.vertex_values
265        self.centroid_values[:,] = quantity.centroid_values
266
267    def set_values_from_function(self,function):
268
269        self.vertex_values[:,:]   = map(function, self.domain.vertices)
270        self.centroid_values[:,] = map(function, self.domain.centroids)
271
272
273    def interpolate(self):
274        """
275        Compute interpolated values at centroid
276        Pre-condition: vertex_values have been set
277        """
278
279        N = self.vertex_values.shape[0]
280        for i in range(N):
281            v0 = self.vertex_values[i, 0]
282            v1 = self.vertex_values[i, 1]
283
284            self.centroid_values[i] = (v0 + v1)/2.0
285
286
287
288
289    def set_values(self, X, location='vertices'):
290        """Set values for quantity
291
292        X: Compatible list, Numeric array (see below), constant or function
293        location: Where values are to be stored.
294                  Permissible options are: vertices, centroid
295                  Default is "vertices"
296
297        In case of location == 'centroid' the dimension values must
298        be a list of a Numerical array of length N, N being the number
299        of elements in the mesh. Otherwise it must be of dimension Nx2
300
301        The values will be stored in elements following their
302        internal ordering.
303
304        If values are described a function, it will be evaluated at specified points
305
306        If selected location is vertices, values for centroid and edges
307        will be assigned interpolated values.
308        In any other case, only values for the specified locations
309        will be assigned and the others will be left undefined.
310        """
311
312        if location not in ['vertices', 'centroids', 'unique_vertices']:
313            msg = 'Invalid location: %s, (possible choices vertices, centroids, unique_vertices)' %location
314            raise Exception(msg)
315
316        if X is None:
317            msg = 'Given values are None'
318            raise Exception(msg)
319
320        import types
321
322        if callable(X):
323            #Use function specific method
324            self.set_function_values(X, location)
325         
326        elif type(X) in [types.FloatType, types.IntType, types.LongType]:
327           
328            self.centroid_values[:,] = float(X)
329            self.vertex_values[:,:] = float(X)
330
331        elif isinstance(X, Quantity):
332            self.set_array_values(X.vertex_values, location)
333
334        else:
335            #Use array specific method
336            self.set_array_values(X, location)
337
338        if location == 'vertices' or location == 'unique_vertices':
339            #Intialise centroid
340            self.interpolate()
341
342        if location == 'centroid':
343            self.extrapolate_first_order()
344
345
346
347
348
349    def set_function_values(self, f, location='vertices'):
350        """Set values for quantity using specified function
351
352        f: x -> z Function where x and z are arrays
353        location: Where values are to be stored.
354                  Permissible options are: vertices, centroid
355                  Default is "vertices"
356        """
357       
358        if location == 'centroids':
359         
360            P = self.domain.centroids
361            self.set_values(f(P), location)
362        else:
363            #Vertices
364           
365            P = self.domain.get_vertices()
366           
367            for i in range(2):
368               
369                self.vertex_values[:,i] = f(P[:,i])
370               
371    def set_array_values(self, values, location='vertices'):
372        """Set values for quantity
373
374        values: Numeric array
375        location: Where values are to be stored.
376                  Permissible options are: vertices, centroid, unique_vertices
377                  Default is "vertices"
378
379        In case of location == 'centroid' the dimension values must
380        be a list of a Numerical array of length N, N being the number
381        of elements in the mesh. If location == 'unique_vertices' the
382        dimension values must be a list or a Numeric array of length N+1.
383        Otherwise it must be of dimension Nx2
384
385        The values will be stored in elements following their
386        internal ordering.
387
388        If selected location is vertices, values for centroid
389        will be assigned interpolated values.
390        In any other case, only values for the specified locations
391        will be assigned and the others will be left undefined.
392        """
393
394
395        values = numpy.array(values).astype(numpy.float)
396
397        N = self.centroid_values.shape[0]
398
399
400        if location == 'centroids':
401            msg = 'Number of values must match number of elements'
402            assert values.shape[0] == N, msg
403            assert len(values.shape) == 1, 'Values array must be 1d'
404
405            for i in range(N):
406                self.centroid_values[i] = values[i]
407
408            self.vertex_values[:,0] = values
409            self.vertex_values[:,1] = values
410 
411        if location == 'vertices':
412            msg = 'Number of values must match number of elements'
413            assert values.shape[0] == N, msg
414            assert len(values.shape) == 2, 'Values array must be 2d'
415            msg = 'Array must be N x 2'
416            assert values.shape[1] == 2, msg
417
418            self.vertex_values[:,:] = values[:,:]
419
420        if location == 'unique_vertices':
421            msg = 'Number of values must match number of elements +1'
422            assert values.shape[0] == N+1, msg
423            assert len(values.shape) == 1, 'Values array must be 1d'
424
425            self.vertex_values[:,0] = values[0:-1]
426            self.vertex_values[:,1] = values[1:N+1] 
427
428
429           
430
431
432    def get_values(self, location='vertices', indices = None):
433        """get values for quantity
434
435        return X, Compatible list, Numeric array (see below)
436        location: Where values are to be stored.
437                  Permissible options are: vertices, centroid
438                  and unique vertices. Default is 'vertices'
439
440        In case of location == 'centroids' the dimension values must
441        be a list of a Numerical array of length N, N being the number
442        of elements. Otherwise it must be of dimension Nx3
443
444        The returned values with be a list the length of indices
445        (N if indices = None).  Each value will be a list of the three
446        vertex values for this quantity.
447
448        Indices is the set of element ids that the operation applies to.
449
450        """
451
452
453        if location not in ['vertices', 'centroids', 'unique vertices']:
454            msg = 'Invalid location: %s' %location
455            raise Exception(msg)
456
457        import types, numpy
458        assert type(indices) in [types.ListType, types.NoneType,
459                                 numpy.ArrayType],\
460                                 'Indices must be a list or None'
461
462        if location == 'centroids':
463            if (indices ==  None):
464                indices = range(len(self))
465            return take(self.centroid_values,indices)
466        elif location == 'unique vertices':
467            if (indices ==  None):
468                indices=range(self.domain.coordinates.shape[0])
469            vert_values = []
470            #Go through list of unique vertices
471            for unique_vert_id in indices:
472                cells = self.domain.vertexlist[unique_vert_id]
473
474                #In case there are unused points
475                if cells is None:
476                    msg = 'Unique vertex not associated with cells'
477                    raise Exception(msg)
478
479                # Go through all cells, vertex pairs
480                # Average the values
481                sum = 0
482                for cell_id, vertex_id in cells:
483                    sum += self.vertex_values[cell_id, vertex_id]
484                vert_values.append(sum/len(cells))
485            return numpy.array(vert_values)
486        else:
487            if (indices ==  None):
488                indices = range(len(self))
489            return take(self.vertex_values,indices)
490
491
492    def get_vertex_values(self,
493                          x=True,
494                          smooth = None,
495                          precision = None,
496                          reduction = None):
497        """Return vertex values like an OBJ format
498
499        The vertex values are returned as one sequence in the 1D float array A.
500        If requested the coordinates will be returned in 1D arrays X.
501
502        The connectivity is represented as an integer array, V, of dimension
503        M x 2, where M is the number of volumes. Each row has two indices
504        into the X, A arrays defining the element.
505
506        if smooth is True, vertex values corresponding to one common
507        coordinate set will be smoothed according to the given
508        reduction operator. In this case vertex coordinates will be
509        de-duplicated.
510
511        If no smoothings is required, vertex coordinates and values will
512        be aggregated as a concatenation of values at
513        vertices 0, vertices 1
514
515
516        Calling convention
517        if x is True:
518           X,A,V = get_vertex_values
519        else:
520           A,V = get_vertex_values
521
522        """
523
524
525
526
527        if smooth is None:
528            smooth = self.domain.smooth
529
530        if precision is None:
531            precision = numpy.float
532
533        if reduction is None:
534            reduction = self.domain.reduction
535
536        #Create connectivity
537
538        if smooth == True:
539
540            V = self.domain.get_vertices()
541            N = len(self.domain.vertexlist)
542            #N = len(self.domain.vertices)
543            A = numpy.zeros(N, precision)
544
545            #Smoothing loop
546            for k in range(N):
547                L = self.domain.vertexlist[k]
548                #L = self.domain.vertices[k]
549
550                #Go through all triangle, vertex pairs
551                #contributing to vertex k and register vertex value
552
553                if L is None: continue #In case there are unused points
554
555                contributions = []
556                for volume_id, vertex_id in L:
557                    v = self.vertex_values[volume_id, vertex_id]
558                    contributions.append(v)
559
560                A[k] = reduction(contributions)
561
562            if x is True:
563                 #X = self.domain.coordinates[:,0].astype(precision)
564                 X = self.domain.coordinates[:].astype(precision)
565                 #Y = self.domain.coordinates[:,1].astype(precision)
566
567                 #return X, Y, A, V
568                 return X, A, V
569           
570            #else:
571            return A, V
572        else:
573            #Don't smooth
574            #obj machinery moved to general_mesh
575
576            # Create a V like [[0 1 2], [3 4 5]....[3*m-2 3*m-1 3*m]]
577            # These vert_id's will relate to the verts created below
578            #m = len(self.domain)  #Number of volumes
579            #M = 3*m        #Total number of unique vertices
580            #V = reshape(array(range(M)).astype(Int), (m,3))
581
582            #V = self.domain.get_triangles(obj=True)
583            V = self.domain.get_vertices
584            #FIXME use get_vertices, when ready
585
586            A = self.vertex_values.flat
587
588            #Do vertex coordinates
589            if x is True:
590                X = self.domain.get_vertex_coordinates()
591
592                #X = C[:,0:6:2].copy()
593                #Y = C[:,1:6:2].copy()
594
595                return X.flat, A, V
596            else:
597                return A, V
598
599    def get_integral(self):
600        """Compute the integral of quantity across entire domain
601        """
602        integral = 0
603        for k in range(self.domain.number_of_elements):
604            area = self.domain.areas[k]
605            qc = self.centroid_values[k]
606            integral += qc*area
607
608        return integral
609
610    def get_beta(self,beta):
611        """Get limiting parameter
612        """
613       
614        return self.beta
615
616    def set_beta(self,beta):
617        """Set limiting parameter
618        """
619
620        #Probably should test that it is not too large
621        self.beta = beta
622
623
624    def get_limiter(self):
625        return self.limiter
626
627    def set_limiter(self,limiter):
628
629        possible_limiters = \
630        ['pyvolution', 'minmod_steve', 'minmod', 'minmod_kurganov', 'superbee', 'vanleer', 'vanalbada']
631
632        if limiter in possible_limiters:
633            self.limiter = limiter
634            return
635
636        msg = '%s is an incorrect limiter type.\n'% limiter
637        msg += 'Possible types are: '+ ", ".join(["%s" % el for el in possible_limiters])
638        raise Exception, msg
639
640    def update(self, timestep):
641        """Update centroid values based on values stored in
642        explicit_update and semi_implicit_update as well as given timestep
643        """
644
645
646        #Explicit updates
647        self.centroid_values += timestep*self.explicit_update
648       
649        #Semi implicit updates
650        denominator = 1.0-timestep*self.semi_implicit_update
651
652#        if sum(numpy.equal(denominator, 0.0)) > 0.0:
653#            msg = 'Zero division in semi implicit update. Call Stephen :-)'
654#            raise Exception(msg)
655#        else:
656#            #Update conserved_quantities from semi implicit updates
657#            self.centroid_values /= denominator
658#           
659
660        #Update conserved_quantities from semi implicit updates
661        self.centroid_values /= denominator
662
663
664    def compute_gradients(self):
665        """Compute gradients of piecewise linear function defined by centroids of
666        neighbouring volumes.
667        """
668
669
670        N = self.centroid_values.shape[0]
671
672
673        G = self.gradients
674        Q = self.centroid_values
675        X = self.domain.centroids
676
677        # first element
678
679        k = 0
680
681        #Get data
682        k0 = k
683        k1 = k+1
684        k2 = k+2
685
686        q0 = Q[k0]
687        q1 = Q[k1]
688        q2 = Q[k2]
689
690        x0 = X[k0] #V0 centroid
691        x1 = X[k1] #V1 centroid
692        x2 = X[k2]
693
694        #Gradient
695        #G[k] = (q1 - q0)/(x1 - x0)
696
697        G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0)
698        G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1)
699
700        #last element
701        k = N-1
702
703
704        k0 = k
705        k1 = k-1
706        k2 = k-2
707
708        q0 = Q[k0]
709        q1 = Q[k1]
710        q2 = Q[k2]
711
712        x0 = X[k0] #V0 centroid
713        x1 = X[k1] #V1 centroid
714        x2 = X[k2]
715
716        #Gradient
717        #G[k] = (q1 - q0)/(x1 - x0)
718
719        G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0)
720        G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1)
721
722
723        #Interior Volume (2 neighbours)
724
725
726        q0 = Q[0:-2]
727        q1 = Q[1:-1]
728        q2 = Q[2:]
729
730        x0 = X[0:-2] #V0 centroid
731        x1 = X[1:-1]  #V1 centroid (Self)
732        x2 = X[2:] #V2 centroid
733
734        #Gradient
735        #G[k] = (q2-q0)/(x2-x0)
736        G[1:-1] = ((q0-q1)/(x0-x1)*(x2-x1) - (q2-q1)/(x2-x1)*(x0-x1))/(x2-x0)
737
738
739
740    def compute_minmod_gradients(self):
741        """Compute gradients of piecewise linear function defined by centroids of
742        neighbouring volumes.
743        """
744
745        #print 'compute_minmod_gradients'
746        from numpy import sign
747
748       
749        def xmin(a,b):
750            from numpy import sign, minimum
751            return 0.5*(sign(a)+sign(b))*minimum(abs(a),abs(b))
752
753        def xmic(t,a,b):
754            return xmin(t*xmin(a,b), 0.50*(a+b) )
755
756
757
758        N = self.centroid_values.shape[0]
759
760
761        G = self.gradients
762        Q = self.centroid_values
763        X = self.domain.centroids
764
765        #-----------------
766        #first element
767        #-----------------
768        k = 0
769
770        k0 = k
771        k1 = k+1
772        k2 = k+2
773
774        q0 = Q[k0]
775        q1 = Q[k1]
776        q2 = Q[k2]
777
778        x0 = X[k0] #V0 centroid
779        x1 = X[k1] #V1 centroid
780        x2 = X[k2]
781
782        #Gradient
783        #G[k] = (q1 - q0)/(x1 - x0)
784
785        G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0)
786        G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1)
787
788        #-------------------
789        # Last element
790        #-------------------
791        k = N-1
792
793        k0 = k
794        k1 = k-1
795        k2 = k-2
796
797        q0 = Q[k0]
798        q1 = Q[k1]
799        q2 = Q[k2]
800
801        x0 = X[k0] #V0 centroid
802        x1 = X[k1] #V1 centroid
803        x2 = X[k2]
804
805        #Gradient
806        #G[k] = (q1 - q0)/(x1 - x0)
807
808        G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0)
809        G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1)
810
811
812
813        #------------------------------
814        #Interior Volume (2 neighbours)
815        #------------------------------
816
817        q0 = Q[0:-2]
818        q1 = Q[1:-1]
819        q2 = Q[2:]
820
821        x0 = X[0:-2] #V0 centroid
822        x1 = X[1:-1]  #V1 centroid (Self)
823        x2 = X[2:] #V2 centroid
824
825        # assuming uniform grid
826        d1 = (q1 - q0)/(x1-x0)
827        d2 = (q2 - q1)/(x2-x1)
828
829        #Gradient
830        G[1:-1] = xmic( self.beta, d1, d2 )     
831
832       
833
834    def extrapolate_first_order(self):
835        """Extrapolate conserved quantities from centroid to
836        vertices for each volume using
837        first order scheme.
838        """
839
840        qc = self.centroid_values
841        qv = self.vertex_values
842
843        for i in range(2):
844            qv[:,i] = qc
845
846
847    def extrapolate_second_order(self):
848        """Extrapolate conserved quantities from centroid to
849        vertices for each volume using
850        second order scheme.
851        """
852
853        if self.limiter == "pyvolution":
854            self.limit_pyvolution()
855
856        elif self.limiter == "minmod_steve":
857            self.limit_minmod()
858
859        else:
860            self.limit_range()
861       
862
863
864
865
866    def find_qmax_qmin(self):
867        """ Find min and max of this and neighbour's centroid values"""
868
869        from numpy import maximum, minimum
870
871        qmax = self.qmax
872        qmin = self.qmin
873
874        qc = self.centroid_values
875       
876        qmax[:] = qc
877        qmin[:] = qc
878       
879        # check left
880        qmax[1:] = maximum(qmax[1:], qc[0:-1])
881        qmin[1:] = minimum(qmin[1:], qc[0:-1])
882       
883        # check right
884        qmax[0:-1] = maximum(qmax[0:-1], qc[1:])
885        qmin[0:-1] = minimum(qmin[0:-1], qc[1:])
886
887
888
889#        for k in range(N):
890#            qmax[k] = qmin[k] = qc[k]
891#            for i in range(2):
892#                n = self.domain.neighbours[k,i]
893#                if n >= 0:
894#                    qn = qc[n] #Neighbour's centroid value
895#
896#                    qmin[k] = min(qmin[k], qn)
897#                    qmax[k] = max(qmax[k], qn)
898
899
900
901    def limit_minmod(self):
902        #Z = self.gradients
903        #print "gradients 1",Z
904        self.compute_minmod_gradients()
905        #print "gradients 2", Z
906
907        G = self.gradients
908        V = self.domain.vertices
909        qc = self.centroid_values
910        qv = self.vertex_values       
911
912        x = self.domain.centroids
913
914        x0 = V[:,0]
915        x1 = V[:,1]
916
917        #Extrapolate
918        qv[:,0] = qc + G*(x0-x)
919        qv[:,1] = qc + G*(x1-x)
920
921#        #Check each triangle
922#        for k in range(self.domain.number_of_elements):
923#            #Centroid coordinates
924#            x = self.domain.centroids[k]
925#
926#            #vertex coordinates
927#            x0, x1 = V[k,:]
928#
929#            #Extrapolate
930#            qv[k,0] = qc[k] + G[k]*(x0-x)
931#            qv[k,1] = qc[k] + G[k]*(x1-x)
932
933 
934    def limit_pyvolution(self):
935        """
936        Limit slopes for each volume to eliminate artificial variance
937        introduced by e.g. second order extrapolator
938
939        This is an unsophisticated limiter as it does not take into
940        account dependencies among quantities.
941
942        precondition:
943        vertex values are estimated from gradient
944        postcondition:
945        vertex values are updated
946        """
947
948
949        N = self.domain.number_of_elements
950        beta = self.domain.beta
951        #beta = 0.8
952
953        self.compute_gradients()
954
955
956        G = self.gradients
957        V = self.domain.vertices
958        C = self.domain.centroids
959        qc = self.centroid_values
960        qv = self.vertex_values
961
962        V0 = V[:,0]
963        V1 = V[:,1]
964
965        # Extrapolate to vertices
966        qv[:,0] = qc + G*(V0-C)
967        qv[:,1] = qc + G*(V1-C)
968
969
970        # Find max and min values
971        self.find_qmax_qmin()
972
973        qmax = self.qmax
974        qmin = self.qmin
975
976        #Diffences between centroids and maxima/minima
977        dqmax = qmax - qc
978        dqmin = qmin - qc
979
980        #Deltas between vertex and centroid values
981        dq = numpy.zeros(qv.shape, numpy.float)
982
983        dq[:,0] = qv[:,0] - qc
984        dq[:,1] = qv[:,1] - qc
985
986        phi = numpy.ones(qc.shape, numpy.float)
987
988        r0 = numpy.where(dq[:,0]>0.0,dqmax/dq[:,0],1.0)
989        r0 = numpy.where(dq[:,0]<0.0,dqmin/dq[:,0],r0)
990
991        r1 = numpy.where(dq[:,1]>0.0,dqmax/dq[:,1],1.0)
992        r1 = numpy.where(dq[:,1]<0.0,dqmin/dq[:,1],r1)
993
994        phi = numpy.min(r0*beta,numpy.min(r1*beta,1.0))
995
996        qv[:,0] = qc + phi*dq[:,0]
997        qv[:,1] = qc + phi*dq[:,1]
998
999#        #Phi limiter
1000#        for k in range(N):
1001#
1002#            #Find the gradient limiter (phi) across vertices
1003#            phi = 1.0
1004#            for i in range(2):
1005#                r = 1.0
1006#                if (dq[k,i] > 0): r = dqmax[k]/dq[k,i]
1007#                if (dq[k,i] < 0): r = dqmin[k]/dq[k,i]
1008#
1009#                phi = min( min(r*beta, 1), phi )
1010#
1011#            #Then update using phi limiter
1012#            for i in range(2):
1013#                qv[k,i] = qc[k] + phi*dq[k,i]
1014
1015    def limit_range(self):
1016        import sys
1017
1018        from limiters_python import minmod, minmod_kurganov, minmod_kurganov_old, maxmod, vanleer, vanalbada
1019
1020        limiter = self.get_limiter()
1021        #print limiter
1022       
1023        #print 'limit_range'
1024        N = self.N
1025        qc = self.centroid_values
1026        qv = self.vertex_values
1027        xc = self.domain.centroids
1028        x0 = self.domain.vertices[:,0]
1029        x1 = self.domain.vertices[:,1]
1030
1031        beta_p = self.beta_p
1032        beta_m = self.beta_m
1033        beta_x = self.beta_x
1034        phi = self.phi
1035        dx  = self.dx
1036
1037
1038        beta_p[1:]  = (qc[1:]-qc[:-1])/(xc[1:]-xc[:-1])
1039        beta_m[:-1] = beta_p[1:]
1040        beta_x[1:-1] = (qc[2:]-qc[:-2])/(xc[2:]-xc[:-2])
1041
1042        dx[:,0] = x0 - xc
1043        dx[:,1] = x1 - xc
1044
1045        phi[0] = (qc[1] - qc[0])/(xc[1] - xc[0])
1046        phi[N-1] = (qc[N-1] - qc[N-2])/(xc[N-1] - xc[N-2])
1047
1048
1049        if limiter == "minmod":
1050            phi[1:-1] = minmod(beta_p[1:-1],beta_m[1:-1])
1051
1052        elif limiter == "vanleer":
1053            phi[1:-1] = vanleer(beta_p[1:-1],beta_m[1:-1])
1054           
1055        elif limiter == "vanalbada":
1056            phi[1:-1] = vanalbada(beta_p[1:-1],beta_m[1:-1])
1057
1058        elif limiter == "minmod_kurganov":
1059            theta = self.beta
1060            phi[1:-1] = minmod_kurganov(theta*beta_p[1:-1],theta*beta_m[1:-1], beta_x[1:-1])
1061
1062        elif limiter == "superbee":
1063            slope1 = minmod(beta_m[1:-1],2.0*beta_p[1:-1])
1064            slope2 = minmod(2.0*beta_m[1:-1],beta_p[1:-1])
1065            phi[1:-1] = maxmod(slope1,slope2)
1066
1067        else:
1068            msg = 'Unknown limiter'
1069            raise Exception, msg
1070
1071
1072       
1073        qv[:,0] = qc + phi*dx[:,0]
1074        qv[:,1] = qc + phi*dx[:,1]
1075
1076
1077
1078
1079    def limit_steve_slope(self):   
1080
1081        import sys
1082
1083        from util import minmod, minmod_kurganov, maxmod, vanleer
1084
1085        N = self.domain.number_of_elements
1086        limiter = self.domain.limiter
1087        limiter_type = self.domain.limiter_type
1088           
1089        qc = self.centroid_values
1090        qv = self.vertex_values
1091
1092        #Find min and max of this and neighbour's centroid values
1093        beta_p = numpy.zeros(N,numpy.float)
1094        beta_m = numpy.zeros(N,numpy.float)
1095        beta_x = numpy.zeros(N,numpy.float)
1096        C = self.domain.centroids
1097        X = self.domain.vertices
1098
1099        for k in range(N):
1100       
1101            n0 = self.domain.neighbours[k,0]
1102            n1 = self.domain.neighbours[k,1]
1103           
1104            if (n0 >= 0) & (n1 >= 0):
1105                # Check denominator not zero
1106                if (qc[k+1]-qc[k]) == 0.0:
1107                    beta_p[k] = float(sys.maxint)
1108                    beta_m[k] = float(sys.maxint)
1109                else:
1110                    #STEVE LIMIT
1111                    beta_p[k] = (qc[k]-qc[k-1])/(qc[k+1]-qc[k])
1112                    beta_m[k] = (qc[k+2]-qc[k+1])/(qc[k+1]-qc[k])
1113
1114        #Deltas between vertex and centroid values
1115        dq = numpy.zeros(qv.shape, numpy.float)
1116        for i in range(2):
1117            dq[:,i] =self.domain.vertices[:,i]-self.domain.centroids
1118           
1119        #Phi limiter
1120        for k in range(N):
1121               
1122            phi = 0.0
1123            if limiter == "flux_minmod":
1124                #FLUX MINMOD
1125                phi = minmod_kurganov(1.0,beta_m[k],beta_p[k])
1126            elif limiter == "flux_superbee":
1127                #FLUX SUPERBEE
1128                phi = max(0.0,min(1.0,2.0*beta_m[k]),min(2.0,beta_m[k]))+max(0.0,min(1.0,2.0*beta_p[k]),min(2.0,beta_p[k]))-1.0
1129            elif limiter == "flux_muscl":
1130                #FLUX MUSCL
1131                phi = max(0.0,min(2.0,2.0*beta_m[k],2.0*beta_p[k],0.5*(beta_m[k]+beta_p[k])))
1132            elif limiter == "flux_vanleer":
1133                #FLUX VAN LEER
1134                phi = (beta_m[k]+abs(beta_m[k]))/(1.0+abs(beta_m[k]))+(beta_p[k]+abs(beta_p[k]))/(1.0+abs(beta_p[k]))-1.0
1135               
1136                #Then update using phi limiter
1137                n = self.domain.neighbours[k,1]
1138                if n>=0:
1139                    #qv[k,0] = qc[k] - 0.5*phi*(qc[k+1]-qc[k])
1140                    #qv[k,1] = qc[k] + 0.5*phi*(qc[k+1]-qc[k])
1141                    qv[k,0] = qc[k] + 0.5*phi*(qv[k,0]-qc[k])
1142                    qv[k,1] = qc[k] + 0.5*phi*(qv[k,1]-qc[k])
1143                else:
1144                    qv[k,i] = qc[k]
1145
1146    def backup_centroid_values(self):
1147        # Call correct module function
1148        # (either from this module or C-extension)
1149        #backup_centroid_values(self)
1150
1151        self.centroid_backup_values[:,] = (self.centroid_values).astype('f')
1152
1153    def saxpy_centroid_values(self,a,b):
1154        # Call correct module function
1155        # (either from this module or C-extension)
1156        self.centroid_values[:,] = (a*self.centroid_values + b*self.centroid_backup_values).astype('f')
1157       
1158
1159
1160
1161   
1162if __name__ == "__main__":
1163    #from domain import Domain
1164    from generic_domain import Generic_domain as Domain
1165
1166
1167    def newLinePlot(title='Simple Plot'):
1168        import pylab as g
1169        g.ion()
1170        g.hold(False)
1171        g.title(title)
1172        g.xlabel('x')
1173        g.ylabel('y')
1174
1175
1176    def linePlot(x,y):
1177        import pylab as g
1178        g.plot(x.flat,y.flat)
1179
1180
1181    def closePlots():
1182        import pylab as g
1183        g.close('all')
1184
1185       
1186   
1187    points1 = [0.0, 1.0, 2.0, 3.0]
1188    vertex_values = [[1.0,2.0],[4.0,5.0],[-1.0,2.0]]
1189
1190    D1 = Domain(points1)
1191
1192    Q1 = Quantity(D1, vertex_values)
1193
1194    print Q1.vertex_values
1195    print Q1.centroid_values
1196
1197    new_vertex_values = [[2.0,1.0],[3.0,4.0],[-2.0,4.0]]
1198
1199    Q1.set_values(new_vertex_values)
1200
1201    print Q1.vertex_values
1202    print Q1.centroid_values
1203
1204    new_centroid_values = [20,30,40]
1205    Q1.set_values(new_centroid_values,'centroids')
1206
1207    print Q1.vertex_values
1208    print Q1.centroid_values
1209
1210    class FunClass:
1211        def __init__(self,value):
1212            self.value = value
1213
1214        def __call__(self,x):
1215            return self.value*(x**2)
1216
1217
1218    fun = FunClass(1.0)
1219    Q1.set_values(fun,'vertices')
1220
1221    print Q1.vertex_values
1222    print Q1.centroid_values
1223
1224    Xc = Q1.domain.vertices
1225    Qc = Q1.vertex_values
1226    print Xc
1227    print Qc
1228
1229    Qc[1,0] = 3
1230
1231    Q1.extrapolate_second_order()
1232    #Q1.limit_minmod()
1233
1234    newLinePlot('plots')
1235    linePlot(Xc,Qc)
1236    raw_input('press return')
1237
1238    points2 = numpy.arange(10)
1239    D2 = Domain(points2)
1240
1241    Q2 = Quantity(D2)
1242    Q2.set_values(fun,'vertices')
1243    Xc = Q2.domain.vertices
1244    Qc = Q2.vertex_values
1245    linePlot(Xc,Qc)
1246    raw_input('press return')
1247
1248   
1249    Q2.extrapolate_second_order()
1250    #Q2.limit_minmod()
1251    Xc = Q2.domain.vertices
1252    Qc = Q2.vertex_values
1253    print Q2.centroid_values
1254    print Qc
1255    linePlot(Xc,Qc)
1256    raw_input('press return')
1257
1258
1259    for i in range(10):
1260        import pylab as g
1261        g.hold(True)
1262        fun = FunClass(i/10.0)
1263        Q2.set_values(fun,'centroids')
1264        Q2.extrapolate_second_order()
1265        #Q2.limit_minmod()
1266        Qc = Q2.vertex_values
1267        linePlot(Xc,Qc)
1268        raw_input('press return')
1269
1270    raw_input('press return to quit')
1271
1272    closePlots()
Note: See TracBrowser for help on using the repository browser.