source: anuga_work/development/flow_1d/generic_1d/quantity.py @ 7830

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

Moving channel code to numpy

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