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

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

Changed name of generic folder

File size: 36.1 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
76        self.gradients = numpy.zeros(N, numpy.float)
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
100    def __add__(self, other):
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
115    def __radd__(self, other):
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
610        #Semi implicit updates
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
632        G = self.gradients
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
731        G = self.gradients
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        if self.domain.limiter == "pyvolution":
841            #Z = self.gradients
842            #print "gradients 1",Z
844            #print "gradients 2",Z
845
846            #Z = self.gradients
847            #print "gradients 1",Z
849            #print "gradients 2", Z
850
851            G = self.gradients
852            V = self.domain.vertices
853            qc = self.centroid_values
854            qv = self.vertex_values
855
856            #Check each triangle
857            for k in range(self.domain.number_of_elements):
858                #Centroid coordinates
859                x = self.domain.centroids[k]
860
861                #vertex coordinates
862                x0, x1 = V[k,:]
863
864                #Extrapolate
865                qv[k,0] = qc[k] + G[k]*(x0-x)
866                qv[k,1] = qc[k] + G[k]*(x1-x)
867            self.limit_pyvolution()
868
869        elif self.domain.limiter == "minmod_steve":
870            self.limit_minmod()
871        else:
872            self.limit_range()
873
874
875
876    def limit_minmod(self):
877        #Z = self.gradients
878        #print "gradients 1",Z
880        #print "gradients 2", Z
881
882        G = self.gradients
883        V = self.domain.vertices
884        qc = self.centroid_values
885        qv = self.vertex_values
886
887        #Check each triangle
888        for k in range(self.domain.number_of_elements):
889            #Centroid coordinates
890            x = self.domain.centroids[k]
891
892            #vertex coordinates
893            x0, x1 = V[k,:]
894
895            #Extrapolate
896            qv[k,0] = qc[k] + G[k]*(x0-x)
897            qv[k,1] = qc[k] + G[k]*(x1-x)
898
899
900    def limit_pyvolution(self):
901        """
902        Limit slopes for each volume to eliminate artificial variance
903        introduced by e.g. second order extrapolator
904
905        This is an unsophisticated limiter as it does not take into
906        account dependencies among quantities.
907
908        precondition:
909        vertex values are estimated from gradient
910        postcondition:
911        vertex values are updated
912        """
913
914
915        N = self.domain.number_of_elements
916        beta = self.domain.beta
917        #beta = 0.8
918
919        qc = self.centroid_values
920        qv = self.vertex_values
921
922        #Find min and max of this and neighbour's centroid values
923        qmax = self.qmax
924        qmin = self.qmin
925
926        for k in range(N):
927            qmax[k] = qmin[k] = qc[k]
928            for i in range(2):
929                n = self.domain.neighbours[k,i]
930                if n >= 0:
931                    qn = qc[n] #Neighbour's centroid value
932
933                    qmin[k] = min(qmin[k], qn)
934                    qmax[k] = max(qmax[k], qn)
935
936
937        #Diffences between centroids and maxima/minima
938        dqmax = qmax - qc
939        dqmin = qmin - qc
940
941        #Deltas between vertex and centroid values
942        dq = numpy.zeros(qv.shape, numpy.float)
943        for i in range(2):
944            dq[:,i] = qv[:,i] - qc
945
946        #Phi limiter
947        for k in range(N):
948
949            #Find the gradient limiter (phi) across vertices
950            phi = 1.0
951            for i in range(2):
952                r = 1.0
953                if (dq[k,i] > 0): r = dqmax[k]/dq[k,i]
954                if (dq[k,i] < 0): r = dqmin[k]/dq[k,i]
955
956                phi = min( min(r*beta, 1), phi )
957
958            #Then update using phi limiter
959            for i in range(2):
960                qv[k,i] = qc[k] + phi*dq[k,i]
961
962    def limit_range(self):
963        import sys
964
965        from util import minmod, minmod_kurganov, maxmod, vanleer, vanalbada
966
967        limiter = self.domain.limiter
968        #print limiter
969
970        #print 'limit_range'
971        N = self.domain.number_of_elements
972        qc = self.centroid_values
973        qv = self.vertex_values
974        C = self.domain.centroids
975        X = self.domain.vertices
976        beta_p = numpy.zeros(N,numpy.float)
977        beta_m = numpy.zeros(N,numpy.float)
978        beta_x = numpy.zeros(N,numpy.float)
979
980        for k in range(N):
981
982            n0 = self.domain.neighbours[k,0]
983            n1 = self.domain.neighbours[k,1]
984
985            if ( n0 >= 0) & (n1 >= 0):
986                #SLOPE DERIVATIVE LIMIT
987                beta_p[k] = (qc[k]-qc[k-1])/(C[k]-C[k-1])
988                beta_m[k] = (qc[k+1]-qc[k])/(C[k+1]-C[k])
989                beta_x[k] = (qc[k+1]-qc[k-1])/(C[k+1]-C[k-1])
990
991        dq = numpy.zeros(qv.shape, numpy.float)
992        for i in range(2):
993            dq[:,i] =self.domain.vertices[:,i]-self.domain.centroids
994
995        #Phi limiter
996        for k in range(N):
997            n0 = self.domain.neighbours[k,0]
998            n1 = self.domain.neighbours[k,1]
999            if n0 < 0:
1000                phi = (qc[k+1] - qc[k])/(C[k+1] - C[k])
1001            elif n1 < 0:
1002                phi = (qc[k] - qc[k-1])/(C[k] - C[k-1])
1003            #elif (self.domain.wet_nodes[k,0] == 2) & (self.domain.wet_nodes[k,1] == 2):
1004            #    phi = 0.0
1005            else:
1006                if limiter == "minmod":
1007                    phi = minmod(beta_p[k],beta_m[k])
1008
1009                elif limiter == "minmod_kurganov":#Change this
1010                    # Also known as monotonized central difference limiter
1011                    # if theta = 2.0
1012                    theta = 2.0
1013                    phi = minmod_kurganov(theta*beta_p[k],theta*beta_m[k],beta_x[k])
1014                elif limiter == "superbee":
1015                    slope1 = minmod(beta_m[k],2.0*beta_p[k])
1016                    slope2 = minmod(2.0*beta_m[k],beta_p[k])
1017                    phi = maxmod(slope1,slope2)
1018
1019                elif limiter == "vanleer":
1020                    phi = vanleer(beta_p[k],beta_m[k])
1021
1022                elif limiter == "vanalbada":
1023                    phi = vanalbada(beta_m[k],beta_p[k])
1024
1025            for i in range(2):
1026                qv[k,i] = qc[k] + phi*dq[k,i]
1027
1028    def limit_steve_slope(self):
1029
1030        import sys
1031
1032        from util import minmod, minmod_kurganov, maxmod, vanleer
1033
1034        N = self.domain.number_of_elements
1035        limiter = self.domain.limiter
1036        limiter_type = self.domain.limiter_type
1037
1038        qc = self.centroid_values
1039        qv = self.vertex_values
1040
1041        #Find min and max of this and neighbour's centroid values
1042        beta_p = numpy.zeros(N,numpy.float)
1043        beta_m = numpy.zeros(N,numpy.float)
1044        beta_x = numpy.zeros(N,numpy.float)
1045        C = self.domain.centroids
1046        X = self.domain.vertices
1047
1048        for k in range(N):
1049
1050            n0 = self.domain.neighbours[k,0]
1051            n1 = self.domain.neighbours[k,1]
1052
1053            if (n0 >= 0) & (n1 >= 0):
1054                # Check denominator not zero
1055                if (qc[k+1]-qc[k]) == 0.0:
1056                    beta_p[k] = float(sys.maxint)
1057                    beta_m[k] = float(sys.maxint)
1058                else:
1059                    #STEVE LIMIT
1060                    beta_p[k] = (qc[k]-qc[k-1])/(qc[k+1]-qc[k])
1061                    beta_m[k] = (qc[k+2]-qc[k+1])/(qc[k+1]-qc[k])
1062
1063        #Deltas between vertex and centroid values
1064        dq = numpy.zeros(qv.shape, numpy.float)
1065        for i in range(2):
1066            dq[:,i] =self.domain.vertices[:,i]-self.domain.centroids
1067
1068        #Phi limiter
1069        for k in range(N):
1070
1071            phi = 0.0
1072            if limiter == "flux_minmod":
1073                #FLUX MINMOD
1074                phi = minmod_kurganov(1.0,beta_m[k],beta_p[k])
1075            elif limiter == "flux_superbee":
1076                #FLUX SUPERBEE
1077                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
1078            elif limiter == "flux_muscl":
1079                #FLUX MUSCL
1080                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])))
1081            elif limiter == "flux_vanleer":
1082                #FLUX VAN LEER
1083                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
1084
1085                #Then update using phi limiter
1086                n = self.domain.neighbours[k,1]
1087                if n>=0:
1088                    #qv[k,0] = qc[k] - 0.5*phi*(qc[k+1]-qc[k])
1089                    #qv[k,1] = qc[k] + 0.5*phi*(qc[k+1]-qc[k])
1090                    qv[k,0] = qc[k] + 0.5*phi*(qv[k,0]-qc[k])
1091                    qv[k,1] = qc[k] + 0.5*phi*(qv[k,1]-qc[k])
1092                else:
1093                    qv[k,i] = qc[k]
1094
1095    def backup_centroid_values(self):
1096        # Call correct module function
1097        # (either from this module or C-extension)
1098        #backup_centroid_values(self)
1099
1100        self.centroid_backup_values[:,] = (self.centroid_values).astype('f')
1101
1102    def saxpy_centroid_values(self,a,b):
1103        # Call correct module function
1104        # (either from this module or C-extension)
1105        self.centroid_values[:,] = (a*self.centroid_values + b*self.centroid_backup_values).astype('f')
1106
1107
1108
1109def newLinePlot(title='Simple Plot'):
1110    import pylab as g
1111    g.ion()
1112    g.hold(False)
1113    g.title(title)
1114    g.xlabel('x')
1115    g.ylabel('y')
1116
1117
1118def linePlot(x,y):
1119    import pylab as g
1120    g.plot(x.flat,y.flat)
1121
1122
1123def closePlots():
1124    import pylab as g
1125    g.close('all')
1126
1127if __name__ == "__main__":
1128    #from domain import Domain
1129    from generic_domain import Generic_domain as Domain
1130
1131    points1 = [0.0, 1.0, 2.0, 3.0]
1132    vertex_values = [[1.0,2.0],[4.0,5.0],[-1.0,2.0]]
1133
1134    D1 = Domain(points1)
1135
1136    Q1 = Quantity(D1, vertex_values)
1137
1138    print Q1.vertex_values
1139    print Q1.centroid_values
1140
1141    new_vertex_values = [[2.0,1.0],[3.0,4.0],[-2.0,4.0]]
1142
1143    Q1.set_values(new_vertex_values)
1144
1145    print Q1.vertex_values
1146    print Q1.centroid_values
1147
1148    new_centroid_values = [20,30,40]
1149    Q1.set_values(new_centroid_values,'centroids')
1150
1151    print Q1.vertex_values
1152    print Q1.centroid_values
1153
1154    class FunClass:
1155        def __init__(self,value):
1156            self.value = value
1157
1158        def __call__(self,x):
1159            return self.value*(x**2)
1160
1161
1162    fun = FunClass(1.0)
1163    Q1.set_values(fun,'vertices')
1164
1165    print Q1.vertex_values
1166    print Q1.centroid_values
1167
1168    Xc = Q1.domain.vertices
1169    Qc = Q1.vertex_values
1170    print Xc
1171    print Qc
1172
1173    Qc[1,0] = 3
1174
1175    Q1.extrapolate_second_order()
1176    #Q1.limit_minmod()
1177
1178    newLinePlot('plots')
1179    linePlot(Xc,Qc)
1180    raw_input('press return')
1181
1182    points2 = numpy.arange(10)
1183    D2 = Domain(points2)
1184
1185    Q2 = Quantity(D2)
1186    Q2.set_values(fun,'vertices')
1187    Xc = Q2.domain.vertices
1188    Qc = Q2.vertex_values
1189    linePlot(Xc,Qc)
1190    raw_input('press return')
1191
1192
1193    Q2.extrapolate_second_order()
1194    #Q2.limit_minmod()
1195    Xc = Q2.domain.vertices
1196    Qc = Q2.vertex_values
1197    print Q2.centroid_values
1198    print Qc
1199    linePlot(Xc,Qc)
1200    raw_input('press return')
1201
1202
1203    for i in range(10):
1204        import pylab as g
1205        g.hold(True)
1206        fun = FunClass(i/10.0)
1207        Q2.set_values(fun,'centroids')
1208        Q2.extrapolate_second_order()
1209        #Q2.limit_minmod()
1210        Qc = Q2.vertex_values
1211        linePlot(Xc,Qc)
1212        raw_input('press return')
1213