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

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

Almost got the unit test working for 1d quantity
File size: 36.3 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
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 msg
302
303        if X is None:
304            msg = 'Given values are None'
305            raise 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        from Numeric import array, Float
382
383        values = array(values).astype(Float)
384
385        N = self.centroid_values.shape[0]
386
387
388        if location == 'centroids':
389            msg = 'Number of values must match number of elements'
390            assert values.shape[0] == N, msg
391            assert len(values.shape) == 1, 'Values array must be 1d'
392
393            for i in range(N):
394                self.centroid_values[i] = values[i]
395
396            self.vertex_values[:,0] = values
397            self.vertex_values[:,1] = values
398 
399        if location == 'vertices':
400            msg = 'Number of values must match number of elements'
401            assert values.shape[0] == N, msg
402            assert len(values.shape) == 2, 'Values array must be 2d'
403            msg = 'Array must be N x 2'
404            assert values.shape[1] == 2, msg
405
406            self.vertex_values[:,:] = values[:,:]
407
408        if location == 'unique_vertices':
409            msg = 'Number of values must match number of elements +1'
410            assert values.shape[0] == N+1, msg
411            assert len(values.shape) == 1, 'Values array must be 1d'
412
413            self.vertex_values[:,0] = values[0:-1]
414            self.vertex_values[:,1] = values[1:N+1] 
415
416
417           
418
419
420    def get_values(self, location='vertices', indices = None):
421        """get values for quantity
422
423        return X, Compatible list, Numeric array (see below)
424        location: Where values are to be stored.
425                  Permissible options are: vertices, centroid
426                  and unique vertices. Default is 'vertices'
427
428        In case of location == 'centroids' the dimension values must
429        be a list of a Numerical array of length N, N being the number
430        of elements. Otherwise it must be of dimension Nx3
431
432        The returned values with be a list the length of indices
433        (N if indices = None).  Each value will be a list of the three
434        vertex values for this quantity.
435
436        Indices is the set of element ids that the operation applies to.
437
438        """
439        from Numeric import take
440
441        if location not in ['vertices', 'centroids', 'unique vertices']:
442            msg = 'Invalid location: %s' %location
443            raise msg
444
445        import types, Numeric
446        assert type(indices) in [types.ListType, types.NoneType,
447                                 Numeric.ArrayType],\
448                                 'Indices must be a list or None'
449
450        if location == 'centroids':
451            if (indices ==  None):
452                indices = range(len(self))
453            return take(self.centroid_values,indices)
454        elif location == 'unique vertices':
455            if (indices ==  None):
456                indices=range(self.domain.coordinates.shape[0])
457            vert_values = []
458            #Go through list of unique vertices
459            for unique_vert_id in indices:
460                cells = self.domain.vertexlist[unique_vert_id]
461
462                #In case there are unused points
463                if cells is None:
464                    msg = 'Unique vertex not associated with cells'
465                    raise msg
466
467                # Go through all cells, vertex pairs
468                # Average the values
469                sum = 0
470                for cell_id, vertex_id in cells:
471                    sum += self.vertex_values[cell_id, vertex_id]
472                vert_values.append(sum/len(cells))
473            return Numeric.array(vert_values)
474        else:
475            if (indices ==  None):
476                indices = range(len(self))
477            return take(self.vertex_values,indices)
478
479
480    def get_vertex_values(self,
481                          x=True,
482                          smooth = None,
483                          precision = None,
484                          reduction = None):
485        """Return vertex values like an OBJ format
486
487        The vertex values are returned as one sequence in the 1D float array A.
488        If requested the coordinates will be returned in 1D arrays X.
489
490        The connectivity is represented as an integer array, V, of dimension
491        M x 2, where M is the number of volumes. Each row has two indices
492        into the X, A arrays defining the element.
493
494        if smooth is True, vertex values corresponding to one common
495        coordinate set will be smoothed according to the given
496        reduction operator. In this case vertex coordinates will be
497        de-duplicated.
498
499        If no smoothings is required, vertex coordinates and values will
500        be aggregated as a concatenation of values at
501        vertices 0, vertices 1
502
503
504        Calling convention
505        if x is True:
506           X,A,V = get_vertex_values
507        else:
508           A,V = get_vertex_values
509
510        """
511
512        from Numeric import concatenate, zeros, Float, Int, array, reshape
513
514
515        if smooth is None:
516            smooth = self.domain.smooth
517
518        if precision is None:
519            precision = Float
520
521        if reduction is None:
522            reduction = self.domain.reduction
523
524        #Create connectivity
525
526        if smooth == True:
527
528            V = self.domain.get_vertices()
529            N = len(self.domain.vertexlist)
530            #N = len(self.domain.vertices)
531            A = zeros(N, precision)
532
533            #Smoothing loop
534            for k in range(N):
535                L = self.domain.vertexlist[k]
536                #L = self.domain.vertices[k]
537
538                #Go through all triangle, vertex pairs
539                #contributing to vertex k and register vertex value
540
541                if L is None: continue #In case there are unused points
542
543                contributions = []
544                for volume_id, vertex_id in L:
545                    v = self.vertex_values[volume_id, vertex_id]
546                    contributions.append(v)
547
548                A[k] = reduction(contributions)
549
550            if x is True:
551                 #X = self.domain.coordinates[:,0].astype(precision)
552                 X = self.domain.coordinates[:].astype(precision)
553                 #Y = self.domain.coordinates[:,1].astype(precision)
554
555                 #return X, Y, A, V
556                 return X, A, V
557           
558            #else:
559            return A, V
560        else:
561            #Don't smooth
562            #obj machinery moved to general_mesh
563
564            # Create a V like [[0 1 2], [3 4 5]....[3*m-2 3*m-1 3*m]]
565            # These vert_id's will relate to the verts created below
566            #m = len(self.domain)  #Number of volumes
567            #M = 3*m        #Total number of unique vertices
568            #V = reshape(array(range(M)).astype(Int), (m,3))
569
570            #V = self.domain.get_triangles(obj=True)
571            V = self.domain.get_vertices
572            #FIXME use get_vertices, when ready
573
574            A = self.vertex_values.flat
575
576            #Do vertex coordinates
577            if x is True:
578                X = self.domain.get_vertex_coordinates()
579
580                #X = C[:,0:6:2].copy()
581                #Y = C[:,1:6:2].copy()
582
583                return X.flat, A, V
584            else:
585                return A, V
586
587    def get_integral(self):
588        """Compute the integral of quantity across entire domain
589        """
590        integral = 0
591        for k in range(self.domain.number_of_elements):
592            area = self.domain.areas[k]
593            qc = self.centroid_values[k]
594            integral += qc*area
595
596        return integral
597
598
599    def update(self, timestep):
600        """Update centroid values based on values stored in
601        explicit_update and semi_implicit_update as well as given timestep
602        """
603
604        from Numeric import sum, equal, ones, Float
605
606        N = self.centroid_values.shape[0]
607
608        #Explicit updates
609        self.centroid_values += timestep*self.explicit_update
610       
611        #Semi implicit updates
612        denominator = ones(N, Float)-timestep*self.semi_implicit_update
613
614        if sum(equal(denominator, 0.0)) > 0.0:
615            msg = 'Zero division in semi implicit update. Call Stephen :-)'
616            raise msg
617        else:
618            #Update conserved_quantities from semi implicit updates
619            self.centroid_values /= denominator
620
621
622    def compute_gradients(self):
623        """Compute gradients of piecewise linear function defined by centroids of
624        neighbouring volumes.
625        """
626
627
628        from Numeric import array, zeros, Float
629
630        N = self.centroid_values.shape[0]
631
632
633        G = self.gradients
634        Q = self.centroid_values
635        X = self.domain.centroids
636
637        for k in range(N):
638
639            # first and last elements have boundaries
640
641            if k == 0:
642
643                #Get data
644                k0 = k
645                k1 = k+1
646                k2 = k+2
647
648                q0 = Q[k0]
649                q1 = Q[k1]
650                q2 = Q[k2]
651
652                x0 = X[k0] #V0 centroid
653                x1 = X[k1] #V1 centroid
654                x2 = X[k2]
655
656                #Gradient
657                #G[k] = (q1 - q0)/(x1 - x0)
658               
659                G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0)
660                G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1)
661
662            elif k == N-1:
663
664                #Get data
665                k0 = k
666                k1 = k-1
667                k2 = k-2
668
669                q0 = Q[k0]
670                q1 = Q[k1]
671                q2 = Q[k2]
672
673                x0 = X[k0] #V0 centroid
674                x1 = X[k1] #V1 centroid
675                x2 = X[k2]
676
677                #Gradient
678                #G[k] = (q1 - q0)/(x1 - x0)
679               
680                G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0)
681                G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1)
682
683##                q0 = Q[k0]
684##                q1 = Q[k1]
685##
686##                x0 = X[k0] #V0 centroid
687##                x1 = X[k1] #V1 centroid
688##
689##                #Gradient
690##                G[k] = (q1 - q0)/(x1 - x0)
691
692            else:
693                #Interior Volume (2 neighbours)
694
695                #Get data
696                k0 = k-1
697                k2 = k+1
698
699                q0 = Q[k0]
700                q1 = Q[k]
701                q2 = Q[k2]
702
703                x0 = X[k0] #V0 centroid
704                x1 = X[k]  #V1 centroid (Self)
705                x2 = X[k2] #V2 centroid
706
707                #Gradient
708                #G[k] = (q2-q0)/(x2-x0)
709                G[k] = ((q0-q1)/(x0-x1)*(x2-x1) - (q2-q1)/(x2-x1)*(x0-x1))/(x2-x0)
710
711
712    def compute_minmod_gradients(self):
713        """Compute gradients of piecewise linear function defined by centroids of
714        neighbouring volumes.
715        """
716
717        #print 'compute_minmod_gradients'
718       
719        from Numeric import array, zeros, Float,sign
720       
721        def xmin(a,b):
722            return 0.5*(sign(a)+sign(b))*min(abs(a),abs(b))
723
724        def xmic(t,a,b):
725            return xmin(t*xmin(a,b), 0.50*(a+b) )
726
727
728
729        N = self.centroid_values.shape[0]
730
731
732        G = self.gradients
733        Q = self.centroid_values
734        X = self.domain.centroids
735
736        for k in range(N):
737
738            # first and last elements have boundaries
739
740            if k == 0:
741
742                #Get data
743                k0 = k
744                k1 = k+1
745                k2 = k+2
746
747                q0 = Q[k0]
748                q1 = Q[k1]
749                q2 = Q[k2]
750
751                x0 = X[k0] #V0 centroid
752                x1 = X[k1] #V1 centroid
753                x2 = X[k2]
754
755                #Gradient
756                #G[k] = (q1 - q0)/(x1 - x0)
757               
758                G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0)
759                G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1)
760
761            elif k == N-1:
762
763                #Get data
764                k0 = k
765                k1 = k-1
766                k2 = k-2
767
768                q0 = Q[k0]
769                q1 = Q[k1]
770                q2 = Q[k2]
771
772                x0 = X[k0] #V0 centroid
773                x1 = X[k1] #V1 centroid
774                x2 = X[k2]
775
776                #Gradient
777                #G[k] = (q1 - q0)/(x1 - x0)
778               
779                G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0)
780                G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1)
781
782##                #Get data
783##                k0 = k
784##                k1 = k-1
785##
786##                q0 = Q[k0]
787##                q1 = Q[k1]
788##
789##                x0 = X[k0] #V0 centroid
790##                x1 = X[k1] #V1 centroid
791##
792##                #Gradient
793##                G[k] = (q1 - q0)/(x1 - x0)
794
795            elif (self.domain.wet_nodes[k,0] == 2) & (self.domain.wet_nodes[k,1] == 2):
796                G[k] = 0.0
797
798            else:
799                #Interior Volume (2 neighbours)
800
801                #Get data
802                k0 = k-1
803                k2 = k+1
804
805                q0 = Q[k0]
806                q1 = Q[k]
807                q2 = Q[k2]
808
809                x0 = X[k0] #V0 centroid
810                x1 = X[k]  #V1 centroid (Self)
811                x2 = X[k2] #V2 centroid
812
813                # assuming uniform grid
814                d1 = (q1 - q0)/(x1-x0)
815                d2 = (q2 - q1)/(x2-x1)
816
817                #Gradient
818                #G[k] = (d1+d2)*0.5
819                #G[k] = (d1*(x2-x1) - d2*(x0-x1))/(x2-x0)               
820                G[k] = xmic( self.domain.beta, d1, d2 )
821       
822
823    def extrapolate_first_order(self):
824        """Extrapolate conserved quantities from centroid to
825        vertices for each volume using
826        first order scheme.
827        """
828
829        qc = self.centroid_values
830        qv = self.vertex_values
831
832        for i in range(2):
833            qv[:,i] = qc
834
835
836    def extrapolate_second_order(self):
837        """Extrapolate conserved quantities from centroid to
838        vertices for each volume using
839        second order scheme.
840        """
841        if self.domain.limiter == "pyvolution":
842            #Z = self.gradients
843            #print "gradients 1",Z
844            self.compute_gradients()
845            #print "gradients 2",Z
846
847            #Z = self.gradients
848            #print "gradients 1",Z
849            #self.compute_minmod_gradients()
850            #print "gradients 2", Z
851
852            G = self.gradients
853            V = self.domain.vertices
854            qc = self.centroid_values
855            qv = self.vertex_values       
856
857            #Check each triangle
858            for k in range(self.domain.number_of_elements):
859                #Centroid coordinates
860                x = self.domain.centroids[k]
861
862                #vertex coordinates
863                x0, x1 = V[k,:]
864
865                #Extrapolate
866                qv[k,0] = qc[k] + G[k]*(x0-x)
867                qv[k,1] = qc[k] + G[k]*(x1-x)
868            self.limit_pyvolution()
869        elif self.domain.limiter == "minmod_steve":
870            self.limit_minmod()
871        else:
872            self.limit_range()
873       
874       
875
876    def limit_minmod(self):
877        #Z = self.gradients
878        #print "gradients 1",Z
879        self.compute_minmod_gradients()
880        #print "gradients 2", Z
881
882        G = self.gradients
883        V = self.domain.vertices
884        qc = self.centroid_values
885        qv = self.vertex_values       
886       
887        #Check each triangle
888        for k in range(self.domain.number_of_elements):
889            #Centroid coordinates
890            x = self.domain.centroids[k]
891
892            #vertex coordinates
893            x0, x1 = V[k,:]
894
895            #Extrapolate
896            qv[k,0] = qc[k] + G[k]*(x0-x)
897            qv[k,1] = qc[k] + G[k]*(x1-x)
898
899 
900    def limit_pyvolution(self):
901        """
902        Limit slopes for each volume to eliminate artificial variance
903        introduced by e.g. second order extrapolator
904
905        This is an unsophisticated limiter as it does not take into
906        account dependencies among quantities.
907
908        precondition:
909        vertex values are estimated from gradient
910        postcondition:
911        vertex values are updated
912        """
913        from Numeric import zeros, Float
914
915        N = self.domain.number_of_elements
916        beta = self.domain.beta
917        #beta = 0.8
918
919        qc = self.centroid_values
920        qv = self.vertex_values
921
922        #Find min and max of this and neighbour's centroid values
923        qmax = self.qmax
924        qmin = self.qmin
925
926        for k in range(N):
927            qmax[k] = qmin[k] = qc[k]
928            for i in range(2):
929                n = self.domain.neighbours[k,i]
930                if n >= 0:
931                    qn = qc[n] #Neighbour's centroid value
932
933                    qmin[k] = min(qmin[k], qn)
934                    qmax[k] = max(qmax[k], qn)
935
936
937        #Diffences between centroids and maxima/minima
938        dqmax = qmax - qc
939        dqmin = qmin - qc
940
941        #Deltas between vertex and centroid values
942        dq = zeros(qv.shape, Float)
943        for i in range(2):
944            dq[:,i] = qv[:,i] - qc
945
946        #Phi limiter
947        for k in range(N):
948
949            #Find the gradient limiter (phi) across vertices
950            phi = 1.0
951            for i in range(2):
952                r = 1.0
953                if (dq[k,i] > 0): r = dqmax[k]/dq[k,i]
954                if (dq[k,i] < 0): r = dqmin[k]/dq[k,i]
955
956                phi = min( min(r*beta, 1), phi )
957
958            #Then update using phi limiter
959            for i in range(2):
960                qv[k,i] = qc[k] + phi*dq[k,i]
961
962    def limit_range(self):
963        import sys
964        from Numeric import zeros, Float
965        from util import minmod, minmod_kurganov, maxmod, vanleer, vanalbada
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 = zeros(N,Float)
976        beta_m = zeros(N,Float)
977        beta_x = zeros(N,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 = zeros(qv.shape, 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        from Numeric import zeros, Float
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 = zeros(N,Float)
1042        beta_m = zeros(N,Float)
1043        beta_x = zeros(N,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 = zeros(qv.shape, 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 shallow_water_domain import Domain     
1129    from Numeric import arange
1130   
1131    points1 = [0.0, 1.0, 2.0, 3.0]
1132    vertex_values = [[1.0,2.0],[4.0,5.0],[-1.0,2.0]]
1133
1134    D1 = Domain(points1)
1135
1136    Q1 = Quantity(D1, vertex_values)
1137
1138    print Q1.vertex_values
1139    print Q1.centroid_values
1140
1141    new_vertex_values = [[2.0,1.0],[3.0,4.0],[-2.0,4.0]]
1142
1143    Q1.set_values(new_vertex_values)
1144
1145    print Q1.vertex_values
1146    print Q1.centroid_values
1147
1148    new_centroid_values = [20,30,40]
1149    Q1.set_values(new_centroid_values,'centroids')
1150
1151    print Q1.vertex_values
1152    print Q1.centroid_values
1153
1154    class FunClass:
1155        def __init__(self,value):
1156            self.value = value
1157
1158        def __call__(self,x):
1159            return self.value*(x**2)
1160
1161
1162    fun = FunClass(1.0)
1163    Q1.set_values(fun,'vertices')
1164
1165    print Q1.vertex_values
1166    print Q1.centroid_values
1167
1168    Xc = Q1.domain.vertices
1169    Qc = Q1.vertex_values
1170    print Xc
1171    print Qc
1172
1173    Qc[1,0] = 3
1174
1175    Q1.extrapolate_second_order()
1176    #Q1.limit_minmod()
1177
1178    newLinePlot('plots')
1179    linePlot(Xc,Qc)
1180    raw_input('press return')
1181
1182    points2 = arange(10)
1183    D2 = Domain(points2)
1184
1185    Q2 = Quantity(D2)
1186    Q2.set_values(fun,'vertices')
1187    Xc = Q2.domain.vertices
1188    Qc = Q2.vertex_values
1189    linePlot(Xc,Qc)
1190    raw_input('press return')
1191
1192   
1193    Q2.extrapolate_second_order()
1194    #Q2.limit_minmod()
1195    Xc = Q2.domain.vertices
1196    Qc = Q2.vertex_values
1197    print Q2.centroid_values
1198    print Qc
1199    linePlot(Xc,Qc)
1200    raw_input('press return')
1201
1202
1203    for i in range(10):
1204        import pylab as g
1205        g.hold(True)
1206        fun = FunClass(i/10.0)
1207        Q2.set_values(fun,'centroids')
1208        Q2.extrapolate_second_order()
1209        #Q2.limit_minmod()
1210        Qc = Q2.vertex_values
1211        linePlot(Xc,Qc)
1212        raw_input('press return')
1213
1214    raw_input('press return to quit')
1215closePlots()
Note: See TracBrowser for help on using the repository browser.