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

Last change on this file since 7855 was 7855, checked in by steve, 14 years ago

Changing for loop to numpy.where

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