# source:anuga_work/development/2010-projects/anuga_1d/base/quantity.py@7852

Last change on this file since 7852 was 7852, checked in by steve, 13 years ago

Moving calculation of limiters to numpy calculations

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