source: anuga_work/development/anuga_1d/quantity.py @ 7793

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

Setting up test function for Quantity

File size: 35.5 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
18
19
20class Quantity:
21
22   
23    def __init__(self, domain, vertex_values=None):
24        #Initialise Quantity using optional vertex values.
25       
26        from domain import Domain
27        from Numeric import array, zeros, Float
28
29        msg = 'First argument in Quantity.__init__ '
30        msg += 'must be of class Domain (or a subclass thereof)'
31        assert isinstance(domain, Domain), msg
32
33        if vertex_values is None:
34            N = domain.number_of_elements
35            self.vertex_values = zeros((N, 2), Float)
36        else:
37            self.vertex_values = array(vertex_values, Float)
38
39            N, V = self.vertex_values.shape
40            assert V == 2,\
41                   'Two vertex values per element must be specified'
42
43
44            msg = 'Number of vertex values (%d) must be consistent with'\
45                  %N
46            msg += 'number of elements in specified domain (%d).'\
47                   %domain.number_of_elements
48
49            assert N == domain.number_of_elements, msg
50
51        self.domain = domain
52
53        #Allocate space for other quantities
54        self.centroid_values = zeros(N, Float)
55        self.centroid_backup_values = zeros(N, Float)
56        #self.edge_values = zeros((N, 2), Float)
57        #edge values are values of the ends of each interval
58     
59        #Intialise centroid values
60        self.interpolate()
61
62
63        from Numeric import zeros, Float
64       
65        #Allocate space for boundary values
66        #L = len(domain.boundary)
67        self.boundary_values = zeros(2, 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 = zeros(N, Float )
74        self.semi_implicit_update = zeros(N, Float )
75
76        self.gradients = zeros(N, Float)
77        self.qmax = zeros(self.centroid_values.shape, Float)
78        self.qmin = zeros(self.centroid_values.shape, 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        import Numeric
218        x = Numeric.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    def set_values(self, X, location='vertices'):
274        """Set values for quantity
275
276        X: Compatible list, Numeric array (see below), constant or function
277        location: Where values are to be stored.
278                  Permissible options are: vertices, centroid
279                  Default is "vertices"
280
281        In case of location == 'centroid' the dimension values must
282        be a list of a Numerical array of length N, N being the number
283        of elements in the mesh. Otherwise it must be of dimension Nx2
284
285        The values will be stored in elements following their
286        internal ordering.
287
288        If values are described a function, it will be evaluated at specified points
289
290        If selected location is vertices, values for centroid and edges
291        will be assigned interpolated values.
292        In any other case, only values for the specified locations
293        will be assigned and the others will be left undefined.
294        """
295
296        if location not in ['vertices', 'centroids']:
297            msg = 'Invalid location: %s, (possible choices vertices, centroids)' %location
298            raise msg
299
300        if X is None:
301            msg = 'Given values are None'
302            raise msg
303
304        import types
305
306        if callable(X):
307            #Use function specific method
308            self.set_function_values(X, location)
309         
310        elif type(X) in [types.FloatType, types.IntType, types.LongType]:
311            if location == 'centroids':
312                self.centroid_values[:] = X
313            else:
314                self.vertex_values[:] = X
315
316        else:
317            #Use array specific method
318            self.set_array_values(X, location)
319
320        if location == 'vertices':
321            #Intialise centroid and edge values
322            self.interpolate()
323
324
325
326
327
328    def set_function_values(self, f, location='vertices'):
329        """Set values for quantity using specified function
330
331        f: x -> z Function where x and z are arrays
332        location: Where values are to be stored.
333                  Permissible options are: vertices, centroid
334                  Default is "vertices"
335        """
336       
337        if location == 'centroids':
338         
339            P = self.domain.centroids
340            self.set_values(f(P), location)
341        else:
342            #Vertices
343           
344            P = self.domain.get_vertices()
345           
346            for i in range(2):
347               
348                self.vertex_values[:,i] = f(P[:,i])
349               
350    def set_array_values(self, values, location='vertices'):
351        """Set values for quantity
352
353        values: Numeric array
354        location: Where values are to be stored.
355                  Permissible options are: vertices, centroid, edges
356                  Default is "vertices"
357
358        In case of location == 'centroid' the dimension values must
359        be a list of a Numerical array of length N, N being the number
360        of elements in the mesh. Otherwise it must be of dimension Nx2
361
362        The values will be stored in elements following their
363        internal ordering.
364
365        If selected location is vertices, values for centroid
366        will be assigned interpolated values.
367        In any other case, only values for the specified locations
368        will be assigned and the others will be left undefined.
369        """
370
371        from Numeric import array, Float
372
373        values = array(values).astype(Float)
374
375        N = self.centroid_values.shape[0]
376
377        msg = 'Number of values must match number of elements'
378        assert values.shape[0] == N, msg
379
380        if location == 'centroids':
381            assert len(values.shape) == 1, 'Values array must be 1d'
382            self.centroid_values = values
383        #elif location == 'edges':
384        #    assert len(values.shape) == 2, 'Values array must be 2d'
385        #    msg = 'Array must be N x 2'
386        #    self.edge_values = values
387        else:
388            assert len(values.shape) == 2, 'Values array must be 2d'
389            msg = 'Array must be N x 2'
390            assert values.shape[1] == 2, msg
391
392            self.vertex_values[:] = values
393
394
395    def get_values(self, location='vertices', indices = None):
396        """get values for quantity
397
398        return X, Compatible list, Numeric array (see below)
399        location: Where values are to be stored.
400                  Permissible options are: vertices, centroid
401                  and unique vertices. Default is 'vertices'
402
403        In case of location == 'centroids' the dimension values must
404        be a list of a Numerical array of length N, N being the number
405        of elements. Otherwise it must be of dimension Nx3
406
407        The returned values with be a list the length of indices
408        (N if indices = None).  Each value will be a list of the three
409        vertex values for this quantity.
410
411        Indices is the set of element ids that the operation applies to.
412
413        """
414        from Numeric import take
415
416        if location not in ['vertices', 'centroids', 'unique vertices']:
417            msg = 'Invalid location: %s' %location
418            raise msg
419
420        import types, Numeric
421        assert type(indices) in [types.ListType, types.NoneType,
422                                 Numeric.ArrayType],\
423                                 'Indices must be a list or None'
424
425        if location == 'centroids':
426            if (indices ==  None):
427                indices = range(len(self))
428            return take(self.centroid_values,indices)
429        elif location == 'unique vertices':
430            if (indices ==  None):
431                indices=range(self.domain.coordinates.shape[0])
432            vert_values = []
433            #Go through list of unique vertices
434            for unique_vert_id in indices:
435                cells = self.domain.vertexlist[unique_vert_id]
436
437                #In case there are unused points
438                if cells is None:
439                    msg = 'Unique vertex not associated with cells'
440                    raise msg
441
442                # Go through all cells, vertex pairs
443                # Average the values
444                sum = 0
445                for cell_id, vertex_id in cells:
446                    sum += self.vertex_values[cell_id, vertex_id]
447                vert_values.append(sum/len(cells))
448            return Numeric.array(vert_values)
449        else:
450            if (indices ==  None):
451                indices = range(len(self))
452            return take(self.vertex_values,indices)
453
454
455    def get_vertex_values(self,
456                          x=True,
457                          smooth = None,
458                          precision = None,
459                          reduction = None):
460        """Return vertex values like an OBJ format
461
462        The vertex values are returned as one sequence in the 1D float array A.
463        If requested the coordinates will be returned in 1D arrays X.
464
465        The connectivity is represented as an integer array, V, of dimension
466        M x 2, where M is the number of volumes. Each row has two indices
467        into the X, A arrays defining the element.
468
469        if smooth is True, vertex values corresponding to one common
470        coordinate set will be smoothed according to the given
471        reduction operator. In this case vertex coordinates will be
472        de-duplicated.
473
474        If no smoothings is required, vertex coordinates and values will
475        be aggregated as a concatenation of values at
476        vertices 0, vertices 1
477
478
479        Calling convention
480        if x is True:
481           X,A,V = get_vertex_values
482        else:
483           A,V = get_vertex_values
484
485        """
486
487        from Numeric import concatenate, zeros, Float, Int, array, reshape
488
489
490        if smooth is None:
491            smooth = self.domain.smooth
492
493        if precision is None:
494            precision = Float
495
496        if reduction is None:
497            reduction = self.domain.reduction
498
499        #Create connectivity
500
501        if smooth == True:
502
503            V = self.domain.get_vertices()
504            N = len(self.domain.vertexlist)
505            #N = len(self.domain.vertices)
506            A = zeros(N, precision)
507
508            #Smoothing loop
509            for k in range(N):
510                L = self.domain.vertexlist[k]
511                #L = self.domain.vertices[k]
512
513                #Go through all triangle, vertex pairs
514                #contributing to vertex k and register vertex value
515
516                if L is None: continue #In case there are unused points
517
518                contributions = []
519                for volume_id, vertex_id in L:
520                    v = self.vertex_values[volume_id, vertex_id]
521                    contributions.append(v)
522
523                A[k] = reduction(contributions)
524
525            if x is True:
526                 #X = self.domain.coordinates[:,0].astype(precision)
527                 X = self.domain.coordinates[:].astype(precision)
528                 #Y = self.domain.coordinates[:,1].astype(precision)
529
530                 #return X, Y, A, V
531                 return X, A, V
532           
533            #else:
534            return A, V
535        else:
536            #Don't smooth
537            #obj machinery moved to general_mesh
538
539            # Create a V like [[0 1 2], [3 4 5]....[3*m-2 3*m-1 3*m]]
540            # These vert_id's will relate to the verts created below
541            #m = len(self.domain)  #Number of volumes
542            #M = 3*m        #Total number of unique vertices
543            #V = reshape(array(range(M)).astype(Int), (m,3))
544
545            #V = self.domain.get_triangles(obj=True)
546            V = self.domain.get_vertices
547            #FIXME use get_vertices, when ready
548
549            A = self.vertex_values.flat
550
551            #Do vertex coordinates
552            if x is True:
553                X = self.domain.get_vertex_coordinates()
554
555                #X = C[:,0:6:2].copy()
556                #Y = C[:,1:6:2].copy()
557
558                return X.flat, A, V
559            else:
560                return A, V
561
562    def get_integral(self):
563        """Compute the integral of quantity across entire domain
564        """
565        integral = 0
566        for k in range(self.domain.number_of_elements):
567            area = self.domain.areas[k]
568            qc = self.centroid_values[k]
569            integral += qc*area
570
571        return integral
572
573
574    def update(self, timestep):
575        """Update centroid values based on values stored in
576        explicit_update and semi_implicit_update as well as given timestep
577        """
578
579        from Numeric import sum, equal, ones, Float
580
581        N = self.centroid_values.shape[0]
582
583        #Explicit updates
584        self.centroid_values += timestep*self.explicit_update
585       
586        #Semi implicit updates
587        denominator = ones(N, Float)-timestep*self.semi_implicit_update
588
589        if sum(equal(denominator, 0.0)) > 0.0:
590            msg = 'Zero division in semi implicit update. Call Stephen :-)'
591            raise msg
592        else:
593            #Update conserved_quantities from semi implicit updates
594            self.centroid_values /= denominator
595
596
597    def compute_gradients(self):
598        """Compute gradients of piecewise linear function defined by centroids of
599        neighbouring volumes.
600        """
601
602
603        from Numeric import array, zeros, Float
604
605        N = self.centroid_values.shape[0]
606
607
608        G = self.gradients
609        Q = self.centroid_values
610        X = self.domain.centroids
611
612        for k in range(N):
613
614            # first and last elements have boundaries
615
616            if k == 0:
617
618                #Get data
619                k0 = k
620                k1 = k+1
621                k2 = k+2
622
623                q0 = Q[k0]
624                q1 = Q[k1]
625                q2 = Q[k2]
626
627                x0 = X[k0] #V0 centroid
628                x1 = X[k1] #V1 centroid
629                x2 = X[k2]
630
631                #Gradient
632                #G[k] = (q1 - q0)/(x1 - x0)
633               
634                G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0)
635                G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1)
636
637            elif k == N-1:
638
639                #Get data
640                k0 = k
641                k1 = k-1
642                k2 = k-2
643
644                q0 = Q[k0]
645                q1 = Q[k1]
646                q2 = Q[k2]
647
648                x0 = X[k0] #V0 centroid
649                x1 = X[k1] #V1 centroid
650                x2 = X[k2]
651
652                #Gradient
653                #G[k] = (q1 - q0)/(x1 - x0)
654               
655                G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0)
656                G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1)
657
658##                q0 = Q[k0]
659##                q1 = Q[k1]
660##
661##                x0 = X[k0] #V0 centroid
662##                x1 = X[k1] #V1 centroid
663##
664##                #Gradient
665##                G[k] = (q1 - q0)/(x1 - x0)
666
667            else:
668                #Interior Volume (2 neighbours)
669
670                #Get data
671                k0 = k-1
672                k2 = k+1
673
674                q0 = Q[k0]
675                q1 = Q[k]
676                q2 = Q[k2]
677
678                x0 = X[k0] #V0 centroid
679                x1 = X[k]  #V1 centroid (Self)
680                x2 = X[k2] #V2 centroid
681
682                #Gradient
683                #G[k] = (q2-q0)/(x2-x0)
684                G[k] = ((q0-q1)/(x0-x1)*(x2-x1) - (q2-q1)/(x2-x1)*(x0-x1))/(x2-x0)
685
686
687    def compute_minmod_gradients(self):
688        """Compute gradients of piecewise linear function defined by centroids of
689        neighbouring volumes.
690        """
691
692        #print 'compute_minmod_gradients'
693       
694        from Numeric import array, zeros, Float,sign
695       
696        def xmin(a,b):
697            return 0.5*(sign(a)+sign(b))*min(abs(a),abs(b))
698
699        def xmic(t,a,b):
700            return xmin(t*xmin(a,b), 0.50*(a+b) )
701
702
703
704        N = self.centroid_values.shape[0]
705
706
707        G = self.gradients
708        Q = self.centroid_values
709        X = self.domain.centroids
710
711        for k in range(N):
712
713            # first and last elements have boundaries
714
715            if k == 0:
716
717                #Get data
718                k0 = k
719                k1 = k+1
720                k2 = k+2
721
722                q0 = Q[k0]
723                q1 = Q[k1]
724                q2 = Q[k2]
725
726                x0 = X[k0] #V0 centroid
727                x1 = X[k1] #V1 centroid
728                x2 = X[k2]
729
730                #Gradient
731                #G[k] = (q1 - q0)/(x1 - x0)
732               
733                G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0)
734                G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1)
735
736            elif k == N-1:
737
738                #Get data
739                k0 = k
740                k1 = k-1
741                k2 = k-2
742
743                q0 = Q[k0]
744                q1 = Q[k1]
745                q2 = Q[k2]
746
747                x0 = X[k0] #V0 centroid
748                x1 = X[k1] #V1 centroid
749                x2 = X[k2]
750
751                #Gradient
752                #G[k] = (q1 - q0)/(x1 - x0)
753               
754                G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0)
755                G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1)
756
757##                #Get data
758##                k0 = k
759##                k1 = k-1
760##
761##                q0 = Q[k0]
762##                q1 = Q[k1]
763##
764##                x0 = X[k0] #V0 centroid
765##                x1 = X[k1] #V1 centroid
766##
767##                #Gradient
768##                G[k] = (q1 - q0)/(x1 - x0)
769
770            elif (self.domain.wet_nodes[k,0] == 2) & (self.domain.wet_nodes[k,1] == 2):
771                G[k] = 0.0
772
773            else:
774                #Interior Volume (2 neighbours)
775
776                #Get data
777                k0 = k-1
778                k2 = k+1
779
780                q0 = Q[k0]
781                q1 = Q[k]
782                q2 = Q[k2]
783
784                x0 = X[k0] #V0 centroid
785                x1 = X[k]  #V1 centroid (Self)
786                x2 = X[k2] #V2 centroid
787
788                # assuming uniform grid
789                d1 = (q1 - q0)/(x1-x0)
790                d2 = (q2 - q1)/(x2-x1)
791
792                #Gradient
793                #G[k] = (d1+d2)*0.5
794                #G[k] = (d1*(x2-x1) - d2*(x0-x1))/(x2-x0)               
795                G[k] = xmic( self.domain.beta, d1, d2 )
796       
797
798    def extrapolate_first_order(self):
799        """Extrapolate conserved quantities from centroid to
800        vertices for each volume using
801        first order scheme.
802        """
803
804        qc = self.centroid_values
805        qv = self.vertex_values
806
807        for i in range(2):
808            qv[:,i] = qc
809
810
811    def extrapolate_second_order(self):
812        """Extrapolate conserved quantities from centroid to
813        vertices for each volume using
814        second order scheme.
815        """
816        if self.domain.limiter == "pyvolution":
817            #Z = self.gradients
818            #print "gradients 1",Z
819            self.compute_gradients()
820            #print "gradients 2",Z
821
822            #Z = self.gradients
823            #print "gradients 1",Z
824            #self.compute_minmod_gradients()
825            #print "gradients 2", Z
826
827            G = self.gradients
828            V = self.domain.vertices
829            qc = self.centroid_values
830            qv = self.vertex_values       
831
832            #Check each triangle
833            for k in range(self.domain.number_of_elements):
834                #Centroid coordinates
835                x = self.domain.centroids[k]
836
837                #vertex coordinates
838                x0, x1 = V[k,:]
839
840                #Extrapolate
841                qv[k,0] = qc[k] + G[k]*(x0-x)
842                qv[k,1] = qc[k] + G[k]*(x1-x)
843            self.limit_pyvolution()
844        elif self.domain.limiter == "minmod_steve":
845            self.limit_minmod()
846        else:
847            self.limit_range()
848       
849       
850
851    def limit_minmod(self):
852        #Z = self.gradients
853        #print "gradients 1",Z
854        self.compute_minmod_gradients()
855        #print "gradients 2", Z
856
857        G = self.gradients
858        V = self.domain.vertices
859        qc = self.centroid_values
860        qv = self.vertex_values       
861       
862        #Check each triangle
863        for k in range(self.domain.number_of_elements):
864            #Centroid coordinates
865            x = self.domain.centroids[k]
866
867            #vertex coordinates
868            x0, x1 = V[k,:]
869
870            #Extrapolate
871            qv[k,0] = qc[k] + G[k]*(x0-x)
872            qv[k,1] = qc[k] + G[k]*(x1-x)
873
874 
875    def limit_pyvolution(self):
876        """
877        Limit slopes for each volume to eliminate artificial variance
878        introduced by e.g. second order extrapolator
879
880        This is an unsophisticated limiter as it does not take into
881        account dependencies among quantities.
882
883        precondition:
884        vertex values are estimated from gradient
885        postcondition:
886        vertex values are updated
887        """
888        from Numeric import zeros, Float
889
890        N = self.domain.number_of_elements
891        beta = self.domain.beta
892        #beta = 0.8
893
894        qc = self.centroid_values
895        qv = self.vertex_values
896
897        #Find min and max of this and neighbour's centroid values
898        qmax = self.qmax
899        qmin = self.qmin
900
901        for k in range(N):
902            qmax[k] = qmin[k] = qc[k]
903            for i in range(2):
904                n = self.domain.neighbours[k,i]
905                if n >= 0:
906                    qn = qc[n] #Neighbour's centroid value
907
908                    qmin[k] = min(qmin[k], qn)
909                    qmax[k] = max(qmax[k], qn)
910
911
912        #Diffences between centroids and maxima/minima
913        dqmax = qmax - qc
914        dqmin = qmin - qc
915
916        #Deltas between vertex and centroid values
917        dq = zeros(qv.shape, Float)
918        for i in range(2):
919            dq[:,i] = qv[:,i] - qc
920
921        #Phi limiter
922        for k in range(N):
923
924            #Find the gradient limiter (phi) across vertices
925            phi = 1.0
926            for i in range(2):
927                r = 1.0
928                if (dq[k,i] > 0): r = dqmax[k]/dq[k,i]
929                if (dq[k,i] < 0): r = dqmin[k]/dq[k,i]
930
931                phi = min( min(r*beta, 1), phi )
932
933            #Then update using phi limiter
934            for i in range(2):
935                qv[k,i] = qc[k] + phi*dq[k,i]
936
937    def limit_range(self):
938        import sys
939        from Numeric import zeros, Float
940        from util import minmod, minmod_kurganov, maxmod, vanleer, vanalbada
941        limiter = self.domain.limiter
942        #print limiter
943       
944        #print 'limit_range'
945        N = self.domain.number_of_elements
946        qc = self.centroid_values
947        qv = self.vertex_values
948        C = self.domain.centroids
949        X = self.domain.vertices
950        beta_p = zeros(N,Float)
951        beta_m = zeros(N,Float)
952        beta_x = zeros(N,Float)
953       
954        for k in range(N):
955       
956            n0 = self.domain.neighbours[k,0]
957            n1 = self.domain.neighbours[k,1]
958           
959            if ( n0 >= 0) & (n1 >= 0):
960                #SLOPE DERIVATIVE LIMIT
961                beta_p[k] = (qc[k]-qc[k-1])/(C[k]-C[k-1])
962                beta_m[k] = (qc[k+1]-qc[k])/(C[k+1]-C[k])
963                beta_x[k] = (qc[k+1]-qc[k-1])/(C[k+1]-C[k-1])
964               
965        dq = zeros(qv.shape, Float)
966        for i in range(2):
967            dq[:,i] =self.domain.vertices[:,i]-self.domain.centroids
968           
969        #Phi limiter
970        for k in range(N):
971            n0 = self.domain.neighbours[k,0]
972            n1 = self.domain.neighbours[k,1]
973            if n0 < 0:
974                phi = (qc[k+1] - qc[k])/(C[k+1] - C[k])
975            elif n1 < 0:
976                phi = (qc[k] - qc[k-1])/(C[k] - C[k-1])
977            #elif (self.domain.wet_nodes[k,0] == 2) & (self.domain.wet_nodes[k,1] == 2):
978            #    phi = 0.0
979            else:
980                if limiter == "minmod":
981                    phi = minmod(beta_p[k],beta_m[k])
982
983                elif limiter == "minmod_kurganov":#Change this
984                    # Also known as monotonized central difference limiter
985                    # if theta = 2.0
986                    theta = 2.0 
987                    phi = minmod_kurganov(theta*beta_p[k],theta*beta_m[k],beta_x[k])
988                elif limiter == "superbee":
989                    slope1 = minmod(beta_m[k],2.0*beta_p[k])
990                    slope2 = minmod(2.0*beta_m[k],beta_p[k])
991                    phi = maxmod(slope1,slope2)
992
993                elif limiter == "vanleer":
994                    phi = vanleer(beta_p[k],beta_m[k])
995
996                elif limiter == "vanalbada":
997                    phi = vanalbada(beta_m[k],beta_p[k])
998           
999            for i in range(2):
1000                qv[k,i] = qc[k] + phi*dq[k,i]
1001
1002    def limit_steve_slope(self):   
1003
1004        import sys
1005        from Numeric import zeros, Float
1006        from util import minmod, minmod_kurganov, maxmod, vanleer
1007
1008        N = self.domain.number_of_elements
1009        limiter = self.domain.limiter
1010        limiter_type = self.domain.limiter_type
1011           
1012        qc = self.centroid_values
1013        qv = self.vertex_values
1014
1015        #Find min and max of this and neighbour's centroid values
1016        beta_p = zeros(N,Float)
1017        beta_m = zeros(N,Float)
1018        beta_x = zeros(N,Float)
1019        C = self.domain.centroids
1020        X = self.domain.vertices
1021
1022        for k in range(N):
1023       
1024            n0 = self.domain.neighbours[k,0]
1025            n1 = self.domain.neighbours[k,1]
1026           
1027            if (n0 >= 0) & (n1 >= 0):
1028                # Check denominator not zero
1029                if (qc[k+1]-qc[k]) == 0.0:
1030                    beta_p[k] = float(sys.maxint)
1031                    beta_m[k] = float(sys.maxint)
1032                else:
1033                    #STEVE LIMIT
1034                    beta_p[k] = (qc[k]-qc[k-1])/(qc[k+1]-qc[k])
1035                    beta_m[k] = (qc[k+2]-qc[k+1])/(qc[k+1]-qc[k])
1036
1037        #Deltas between vertex and centroid values
1038        dq = zeros(qv.shape, Float)
1039        for i in range(2):
1040            dq[:,i] =self.domain.vertices[:,i]-self.domain.centroids
1041           
1042        #Phi limiter
1043        for k in range(N):
1044               
1045            phi = 0.0
1046            if limiter == "flux_minmod":
1047                #FLUX MINMOD
1048                phi = minmod_kurganov(1.0,beta_m[k],beta_p[k])
1049            elif limiter == "flux_superbee":
1050                #FLUX SUPERBEE
1051                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
1052            elif limiter == "flux_muscl":
1053                #FLUX MUSCL
1054                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])))
1055            elif limiter == "flux_vanleer":
1056                #FLUX VAN LEER
1057                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
1058               
1059                #Then update using phi limiter
1060                n = self.domain.neighbours[k,1]
1061                if n>=0:
1062                    #qv[k,0] = qc[k] - 0.5*phi*(qc[k+1]-qc[k])
1063                    #qv[k,1] = qc[k] + 0.5*phi*(qc[k+1]-qc[k])
1064                    qv[k,0] = qc[k] + 0.5*phi*(qv[k,0]-qc[k])
1065                    qv[k,1] = qc[k] + 0.5*phi*(qv[k,1]-qc[k])
1066                else:
1067                    qv[k,i] = qc[k]
1068
1069    def backup_centroid_values(self):
1070        # Call correct module function
1071        # (either from this module or C-extension)
1072        #backup_centroid_values(self)
1073
1074        self.centroid_backup_values[:] = (self.centroid_values).astype('f')
1075
1076    def saxpy_centroid_values(self,a,b):
1077        # Call correct module function
1078        # (either from this module or C-extension)
1079        self.centroid_values[:] = (a*self.centroid_values + b*self.centroid_backup_values).astype('f')
1080       
1081
1082
1083def newLinePlot(title='Simple Plot'):
1084    import pylab as g
1085    g.ion()
1086    g.hold(False)
1087    g.title(title)
1088    g.xlabel('x')
1089    g.ylabel('y')
1090   
1091
1092def linePlot(x,y):
1093    import pylab as g
1094    g.plot(x.flat,y.flat)
1095
1096
1097def closePlots():
1098    import pylab as g
1099    g.close('all')
1100   
1101if __name__ == "__main__":
1102    #from domain import Domain
1103    from shallow_water_domain import Domain     
1104    from Numeric import arange
1105   
1106    points1 = [0.0, 1.0, 2.0, 3.0]
1107    vertex_values = [[1.0,2.0],[4.0,5.0],[-1.0,2.0]]
1108
1109    D1 = Domain(points1)
1110
1111    Q1 = Quantity(D1, vertex_values)
1112
1113    print Q1.vertex_values
1114    print Q1.centroid_values
1115
1116    new_vertex_values = [[2.0,1.0],[3.0,4.0],[-2.0,4.0]]
1117
1118    Q1.set_values(new_vertex_values)
1119
1120    print Q1.vertex_values
1121    print Q1.centroid_values
1122
1123    new_centroid_values = [20,30,40]
1124    Q1.set_values(new_centroid_values,'centroids')
1125
1126    print Q1.vertex_values
1127    print Q1.centroid_values
1128
1129    class FunClass:
1130        def __init__(self,value):
1131            self.value = value
1132
1133        def __call__(self,x):
1134            return self.value*(x**2)
1135
1136
1137    fun = FunClass(1.0)
1138    Q1.set_values(fun,'vertices')
1139
1140    print Q1.vertex_values
1141    print Q1.centroid_values
1142
1143    Xc = Q1.domain.vertices
1144    Qc = Q1.vertex_values
1145    print Xc
1146    print Qc
1147
1148    Qc[1,0] = 3
1149
1150    Q1.extrapolate_second_order()
1151    #Q1.limit_minmod()
1152
1153    newLinePlot('plots')
1154    linePlot(Xc,Qc)
1155    raw_input('press return')
1156
1157    points2 = arange(10)
1158    D2 = Domain(points2)
1159
1160    Q2 = Quantity(D2)
1161    Q2.set_values(fun,'vertices')
1162    Xc = Q2.domain.vertices
1163    Qc = Q2.vertex_values
1164    linePlot(Xc,Qc)
1165    raw_input('press return')
1166
1167   
1168    Q2.extrapolate_second_order()
1169    #Q2.limit_minmod()
1170    Xc = Q2.domain.vertices
1171    Qc = Q2.vertex_values
1172    print Q2.centroid_values
1173    print Qc
1174    linePlot(Xc,Qc)
1175    raw_input('press return')
1176
1177
1178    for i in range(10):
1179        import pylab as g
1180        g.hold(True)
1181        fun = FunClass(i/10.0)
1182        Q2.set_values(fun,'centroids')
1183        Q2.extrapolate_second_order()
1184        #Q2.limit_minmod()
1185        Qc = Q2.vertex_values
1186        linePlot(Xc,Qc)
1187        raw_input('press return')
1188
1189    raw_input('press return to quit')
1190closePlots()
Note: See TracBrowser for help on using the repository browser.