source: development/pyvolution-1d/quantity.py @ 3425

Last change on this file since 3425 was 3425, checked in by steve, 18 years ago
File size: 24.1 KB
Line 
1"""Class Quantity - Implements values at each 1d element
2
3To create:
4
5   Quantity(domain, vertex_values)
6
7   domain: Associated domain structure. Required.
8
9   vertex_values: N x 2 array of values at each vertex for each element.
10                  Default None
11
12   If vertex_values are None Create array of zeros compatible with domain.
13   Otherwise check that it is compatible with dimenions of domain.
14   Otherwise raise an exception
15"""
16
17
18class Quantity:
19
20    def __init__(self, domain, vertex_values=None):
21
22        from domain import Domain
23        from Numeric import array, zeros, Float
24
25        msg = 'First argument in Quantity.__init__ '
26        msg += 'must be of class Domain (or a subclass thereof)'
27        assert isinstance(domain, Domain), msg
28
29        if vertex_values is None:
30            N = domain.number_of_elements
31            self.vertex_values = zeros((N, 2), Float)
32        else:
33            self.vertex_values = array(vertex_values, Float)
34
35            N, V = self.vertex_values.shape
36            assert V == 2,\
37                   'Two vertex values per element must be specified'
38
39
40            msg = 'Number of vertex values (%d) must be consistent with'\
41                  %N
42            msg += 'number of elements in specified domain (%d).'\
43                   %domain.number_of_elements
44
45            assert N == domain.number_of_elements, msg
46
47        self.domain = domain
48
49        #Allocate space for other quantities
50        self.centroid_values = zeros(N, Float)
51        #self.edge_values = zeros((N, 2), Float)
52        #edge values are values of the ends of each interval
53       
54        #does oe dimension need edge values???
55
56        #Intialise centroid values
57        self.interpolate()
58
59    #Methods for operator overloading
60    def __len__(self):
61        return self.centroid_values.shape[0]
62   
63    def interpolate(self):
64        """Compute interpolated values at centroid
65        Pre-condition: vertex_values have been set
66        """
67
68        N = self.vertex_values.shape[0]
69        for i in range(N):
70            v0 = self.vertex_values[i, 0]
71            v1 = self.vertex_values[i, 1]
72
73            self.centroid_values[i] = (v0 + v1)/2
74
75    def set_values(self, X, location='vertices'):
76        """Set values for quantity
77
78        X: Compatible list, Numeric array (see below), constant or function
79        location: Where values are to be stored.
80                  Permissible options are: vertices, centroid
81                  Default is "vertices"
82
83        In case of location == 'centroid' the dimension values must
84        be a list of a Numerical array of length N, N being the number
85        of elements in the mesh. Otherwise it must be of dimension Nx3
86
87        The values will be stored in elements following their
88        internal ordering.
89
90        If values are described a function, it will be evaluated at specified points
91
92        If selected location is vertices, values for centroid and edges
93        will be assigned interpolated values.
94        In any other case, only values for the specified locations
95        will be assigned and the others will be left undefined.
96        """
97
98        if location not in ['vertices', 'centroids']:#, 'edges']:
99            msg = 'Invalid location: %s, (possible choices vertices, centroids)' %location
100            raise msg
101
102        if X is None:
103            msg = 'Given values are None'
104            raise msg
105
106        import types
107
108        if callable(X):
109            #Use function specific method
110            self.set_function_values(X, location)
111        elif type(X) in [types.FloatType, types.IntType, types.LongType]:
112            if location == 'centroids':
113                self.centroid_values[:] = X
114            else:
115                self.vertex_values[:] = X
116
117        else:
118            #Use array specific method
119            self.set_array_values(X, location)
120
121        if location == 'vertices':
122            #Intialise centroid and edge values
123            self.interpolate()
124
125
126
127
128    def set_function_values(self, f, location='vertices'):
129        """Set values for quantity using specified function
130
131        f: x -> z Function where x and z are arrays
132        location: Where values are to be stored.
133                  Permissible options are: vertices, centroid
134                  Default is "vertices"
135        """
136       
137        if location == 'centroids':
138           
139            P = self.domain.centroids
140            self.set_values(f(P), location)
141        else:
142            #Vertices
143            P = self.domain.get_vertices()
144           
145            for i in range(2):
146                self.vertex_values[:,i] = f(P[:,i])       
147
148    def set_array_values(self, values, location='vertices'):
149        """Set values for quantity
150
151        values: Numeric array
152        location: Where values are to be stored.
153                  Permissible options are: vertices, centroid, edges
154                  Default is "vertices"
155
156        In case of location == 'centroid' the dimension values must
157        be a list of a Numerical array of length N, N being the number
158        of elements in the mesh. Otherwise it must be of dimension Nx2
159
160        The values will be stored in elements following their
161        internal ordering.
162
163        If selected location is vertices, values for centroid
164        will be assigned interpolated values.
165        In any other case, only values for the specified locations
166        will be assigned and the others will be left undefined.
167        """
168
169        from Numeric import array, Float
170
171        values = array(values).astype(Float)
172
173        N = self.centroid_values.shape[0]
174
175        msg = 'Number of values must match number of elements'
176        assert values.shape[0] == N, msg
177
178        if location == 'centroids':
179            assert len(values.shape) == 1, 'Values array must be 1d'
180            self.centroid_values = values
181        #elif location == 'edges':
182        #    assert len(values.shape) == 2, 'Values array must be 2d'
183        #    msg = 'Array must be N x 2'
184        #    self.edge_values = values
185        else:
186            assert len(values.shape) == 2, 'Values array must be 2d'
187            msg = 'Array must be N x 2'
188            assert values.shape[1] == 2, msg
189
190            self.vertex_values = values
191
192
193    def get_values(self, location='vertices', indices = None):
194        """get values for quantity
195
196        return X, Compatible list, Numeric array (see below)
197        location: Where values are to be stored.
198                  Permissible options are: vertices, edges, centroid
199                  and unique vertices. Default is 'vertices'
200
201        In case of location == 'centroids' the dimension values must
202        be a list of a Numerical array of length N, N being the number
203        of elements. Otherwise it must be of dimension Nx3
204
205        The returned values with be a list the length of indices
206        (N if indices = None).  Each value will be a list of the three
207        vertex values for this quantity.
208
209        Indices is the set of element ids that the operation applies to.
210
211        """
212        from Numeric import take
213
214        #if location not in ['vertices', 'centroids', 'edges', 'unique vertices']:
215        if location not in ['vertices', 'centroids', 'unique vertices']:
216            msg = 'Invalid location: %s' %location
217            raise msg
218
219        import types, Numeric
220        assert type(indices) in [types.ListType, types.NoneType,
221                                 Numeric.ArrayType],\
222                                 'Indices must be a list or None'
223
224        if location == 'centroids':
225            if (indices ==  None):
226                indices = range(len(self))
227            return take(self.centroid_values,indices)
228        #elif location == 'edges':
229        #    if (indices ==  None):
230        #        indices = range(len(self))
231        #    return take(self.edge_values,indices)
232        elif location == 'unique vertices':
233            if (indices ==  None):
234                indices=range(self.domain.coordinates.shape[0])
235            vert_values = []
236            #Go through list of unique vertices
237            for unique_vert_id in indices:
238                triangles = self.domain.vertexlist[unique_vert_id]
239
240                #In case there are unused points
241                if triangles is None:
242                    msg = 'Unique vertex not associated with triangles'
243                    raise msg
244
245                # Go through all triangle, vertex pairs
246                # Average the values
247                sum = 0
248                for triangle_id, vertex_id in triangles:
249                    sum += self.vertex_values[triangle_id, vertex_id]
250                vert_values.append(sum/len(triangles))
251            return Numeric.array(vert_values)
252        else:
253            if (indices ==  None):
254                indices = range(len(self))
255            return take(self.vertex_values,indices)
256
257    #Method for outputting model results
258    #FIXME: Split up into geometric and numeric stuff.
259    #FIXME: Geometric (X,Y,V) should live in mesh.py
260    #FIXME: STill remember to move XY to mesh
261    def get_vertex_values(self,
262                          #xy=True,
263                          x=True,
264                          smooth = None,
265                          precision = None,
266                          reduction = None):
267        """Return vertex values like an OBJ format
268
269        The vertex values are returned as one sequence in the 1D float array A.
270        If requested the coordinates will be returned in 1D arrays X and Y.
271
272        The connectivity is represented as an integer array, V, of dimension
273        M x 3, where M is the number of volumes. Each row has three indices
274        into the X, Y, A arrays defining the triangle.
275
276        if smooth is True, vertex values corresponding to one common
277        coordinate set will be smoothed according to the given
278        reduction operator. In this case vertex coordinates will be
279        de-duplicated.
280
281        If no smoothings is required, vertex coordinates and values will
282        be aggregated as a concatenation of values at
283        vertices 0, vertices 1 and vertices 2
284
285
286        Calling convention
287        if xy is True:
288           X,Y,A,V = get_vertex_values
289        else:
290           A,V = get_vertex_values
291
292        """
293
294        from Numeric import concatenate, zeros, Float, Int, array, reshape
295
296
297        if smooth is None:
298            smooth = self.domain.smooth
299
300        if precision is None:
301            precision = Float
302
303        if reduction is None:
304            reduction = self.domain.reduction
305
306        #Create connectivity
307
308        if smooth == True:
309
310            V = self.domain.get_vertices()
311            N = len(self.domain.vertexlist)
312            #N = len(self.domain.vertices)
313            A = zeros(N, precision)
314
315            #Smoothing loop
316            for k in range(N):
317                L = self.domain.vertexlist[k]
318                #L = self.domain.vertices[k]
319
320                #Go through all triangle, vertex pairs
321                #contributing to vertex k and register vertex value
322
323                if L is None: continue #In case there are unused points
324
325                contributions = []
326                for volume_id, vertex_id in L:
327                    v = self.vertex_values[volume_id, vertex_id]
328                    contributions.append(v)
329
330                A[k] = reduction(contributions)
331
332
333            #if xy is True:
334            if x is True:
335                 #X = self.domain.coordinates[:,0].astype(precision)
336                 X = self.domain.coordinates[:].astype(precision)
337                 #Y = self.domain.coordinates[:,1].astype(precision)
338
339                 #return X, Y, A, V
340                 return X, A, V
341           
342            #else:
343            return A, V
344        else:
345            #Don't smooth
346            #obj machinery moved to general_mesh
347
348            # Create a V like [[0 1 2], [3 4 5]....[3*m-2 3*m-1 3*m]]
349            # These vert_id's will relate to the verts created below
350            #m = len(self.domain)  #Number of volumes
351            #M = 3*m        #Total number of unique vertices
352            #V = reshape(array(range(M)).astype(Int), (m,3))
353
354            #V = self.domain.get_triangles(obj=True)
355            V = self.domain.get_vertices
356            #FIXME use get_vertices, when ready
357
358            A = self.vertex_values.flat
359
360            #Do vertex coordinates
361            #if xy is True:
362            if x is True:
363                C = self.domain.get_vertex_coordinates()
364
365                X = C[:,0:6:2].copy()
366                Y = C[:,1:6:2].copy()
367
368                return X.flat, Y.flat, A, V
369            else:
370                return A, V
371
372
373class Conserved_quantity(Quantity):
374    """Class conserved quantity adds to Quantity:
375
376    storage and method for updating, and
377    methods for extrapolation from centropid to vertices inluding
378    gradients and limiters
379    """
380
381    def __init__(self, domain, vertex_values=None):
382        Quantity.__init__(self, domain, vertex_values)
383
384        from Numeric import zeros, Float
385       
386        #Allocate space for boundary values
387        #L = len(domain.boundary)
388        self.boundary_values = zeros(2, Float) #assumes no parrellism
389
390        #Allocate space for updates of conserved quantities by
391        #flux calculations and forcing functions
392       
393        N = domain.number_of_elements
394        self.explicit_update = zeros(N, Float )
395        self.semi_implicit_update = zeros(N, Float )
396
397        self.gradients = zeros(N, Float)
398        self.qmax = zeros(self.centroid_values.shape, Float)
399        self.qmin = zeros(self.centroid_values.shape, Float)
400
401        self.beta = domain.beta
402
403
404    def update(self, timestep):
405        """Update centroid values based on values stored in
406        explicit_update and semi_implicit_update as well as given timestep
407        """
408
409        from Numeric import sum, equal, ones, Float
410
411        N = self.centroid_values.shape[0]
412
413        #Explicit updates
414        self.centroid_values += timestep*self.explicit_update
415
416##         #Semi implicit updates
417##         denominator = ones(N, Float)-timestep*self.semi_implicit_update
418
419##         if sum(equal(denominator, 0.0)) > 0.0:
420##             msg = 'Zero division in semi implicit update. Call Stephen :-)'
421##             raise msg
422##         else:
423##             #Update conserved_quantities from semi implicit updates
424##             self.centroid_values /= denominator
425
426
427    def compute_gradients(self):
428        """Compute gradients of piecewise linear function defined by centroids of
429        neighbouring volumes.
430        """
431
432
433        from Numeric import array, zeros, Float
434
435        N = self.centroid_values.shape[0]
436
437
438        G = self.gradients
439        Q = self.centroid_values
440        X = self.domain.centroids
441
442        for k in range(N):
443
444            # first and last elements have boundaries
445
446            if k == 0:
447
448                #Get data
449                k0 = k
450                k1 = k+1
451
452                q0 = Q[k0]
453                q1 = Q[k1]
454
455                x0 = X[k0] #V0 centroid
456                x1 = X[k1] #V1 centroid
457
458                #Gradient
459                G[k] = (q1 - q0)/(x1 - x0)
460
461            elif k == N-1:
462
463                #Get data
464                k0 = k
465                k1 = k-1
466
467                q0 = Q[k0]
468                q1 = Q[k1]
469
470                x0 = X[k0] #V0 centroid
471                x1 = X[k1] #V1 centroid
472
473                #Gradient
474                G[k] = (q1 - q0)/(x1 - x0)
475
476            else:
477                #Interior Volume (2 neighbours)
478
479                #Get data
480                k0 = k-1
481                k2 = k+1
482
483                q0 = Q[k0]
484                q1 = Q[k]
485                q2 = Q[k2]
486
487                x0 = X[k0] #V0 centroid
488                x1 = X[k]  #V1 centroid (Self)
489                x2 = X[k2] #V2 centroid
490
491                #Gradient
492                #G[k] = (q2-q0)/(x2-x0)
493                G[k] = ((q0-q1)/(x0-x1)*(x2-x1) - (q2-q1)/(x2-x1)*(x0-x1))/(x2-x0)
494
495
496    def compute_minmod_gradients(self):
497        """Compute gradients of piecewise linear function defined by centroids of
498        neighbouring volumes.
499        """
500
501        from Numeric import array, zeros, Float,sign
502       
503        def xmin(a,b):
504            return 0.5*(sign(a)+sign(b))*min(abs(a),abs(b))
505
506        def xmic(t,a,b):
507            return xmin(t*xmin(a,b), 0.50*(a+b) )
508
509
510
511        N = self.centroid_values.shape[0]
512
513
514        G = self.gradients
515        Q = self.centroid_values
516        X = self.domain.centroids
517
518        for k in range(N):
519
520            # first and last elements have boundaries
521
522            if k == 0:
523
524                #Get data
525                k0 = k
526                k1 = k+1
527
528                q0 = Q[k0]
529                q1 = Q[k1]
530
531                x0 = X[k0] #V0 centroid
532                x1 = X[k1] #V1 centroid
533
534                #Gradient
535                G[k] = (q1 - q0)/(x1 - x0)
536
537            elif k == N-1:
538
539                #Get data
540                k0 = k
541                k1 = k-1
542
543                q0 = Q[k0]
544                q1 = Q[k1]
545
546                x0 = X[k0] #V0 centroid
547                x1 = X[k1] #V1 centroid
548
549                #Gradient
550                G[k] = (q1 - q0)/(x1 - x0)
551
552            else:
553                #Interior Volume (2 neighbours)
554
555                #Get data
556                k0 = k-1
557                k2 = k+1
558
559                q0 = Q[k0]
560                q1 = Q[k]
561                q2 = Q[k2]
562
563                x0 = X[k0] #V0 centroid
564                x1 = X[k]  #V1 centroid (Self)
565                x2 = X[k2] #V2 centroid
566
567                # assuming uniform grid
568                d1 = (q1 - q0)/(x1-x0)
569                d2 = (q2 - q1)/(x2-x1)
570
571                #Gradient
572                #G[k] = (d1+d2)*0.5
573                #G[k] = (d1*(x2-x1) - d2*(x0-x1))/(x2-x0)               
574                G[k] = xmic( self.domain.beta, d1, d2 )
575       
576
577    def extrapolate_first_order(self):
578        """Extrapolate conserved quantities from centroid to
579        vertices for each volume using
580        first order scheme.
581        """
582
583        qc = self.centroid_values
584        qv = self.vertex_values
585
586        for i in range(2):
587            qv[:,i] = qc
588
589
590    def extrapolate_second_order(self):
591        """Extrapolate conserved quantities from centroid to
592        vertices for each volume using
593        second order scheme.
594        """
595
596        Z = self.gradients
597        #print "gradients 1",Z
598       
599        self.compute_minmod_gradients()
600        #print "gradients 2", Z
601
602        G = self.gradients
603        V = self.domain.vertices
604        qc = self.centroid_values
605        qv = self.vertex_values       
606       
607        #Check each triangle
608        for k in range(self.domain.number_of_elements):
609            #Centroid coordinates
610            x = self.domain.centroids[k]
611
612            #vertex coordinates
613            x0, x1 = V[k,:]
614
615            #Extrapolate
616            qv[k,0] = qc[k] + G[k]*(x0-x)
617            qv[k,1] = qc[k] + G[k]*(x1-x)
618
619    """    def limit(self):
620        Limit slopes for each volume to eliminate artificial variance
621        introduced by e.g. second order extrapolator
622
623        This is an unsophisticated limiter as it does not take into
624        account dependencies among quantities.
625
626        precondition:
627        vertex values are estimated from gradient
628        postcondition:
629        vertex values are updated
630       
631        print "in Q.limit"
632        from Numeric import zeros, Float
633
634        N = self.domain.number_of_elements
635        beta = self.beta
636
637        qc = self.centroid_values
638        qv = self.vertex_values
639
640        #Find min and max of this and neighbour's centroid values
641        qmax = self.qmax
642        qmin = self.qmin
643
644        for k in range(N):
645            qmax[k] = qmin[k] = qc[k]
646
647            for i in [-1,1]:
648                n = k+i
649                if (n >= 0) & (n <= N-1):
650                    qn = qc[n] #Neighbour's centroid value
651
652                    qmin[k] = min(qmin[k], qn)
653                    qmax[k] = max(qmax[k], qn)
654
655        #Phi limiter
656        for k in range(N):
657
658            #Diffences between centroids and maxima/minima
659            dqmax = qmax[k] - qc[k]
660            dqmin = qmin[k] - qc[k]
661
662            #Deltas between vertex and centroid values
663            dq = [0.0, 0.0]
664            for i in range(2):
665                dq[i] = qv[k,i] - qc[k]
666
667            #Find the gradient limiter (phi) across vertices
668            phi = 1.0
669            for i in range(2):
670                r = 1.0
671                if (dq[i] > 0): r = dqmax/dq[i]
672                if (dq[i] < 0): r = dqmin/dq[i]
673
674                phi = min( min(r*beta, 1), phi )
675
676            #Then update using phi limiter
677            for i in range(2):
678                qv[k,i] = qc[k] + phi*dq[i]
679    """
680
681    def limit(self):
682        """Limit slopes for each volume to eliminate artificial variance
683        introduced by e.g. second order extrapolator
684
685        This is an unsophisticated limiter as it does not take into
686        account dependencies among quantities.
687
688        precondition:
689        vertex values are estimated from gradient
690        postcondition:
691        vertex values are updated
692        """
693
694        from Numeric import zeros, Float
695
696        N = self.domain.number_of_elements
697        #beta = self.beta
698        beta = 0.8
699
700        qc = self.centroid_values
701        qv = self.vertex_values
702
703        #Find min and max of this and neighbour's centroid values
704        qmax = self.qmax
705        qmin = self.qmin
706
707        for k in range(N):
708            qmax[k] = qmin[k] = qc[k]
709            for i in range(2):
710                n = self.domain.neighbours[k,i]
711                if n >= 0:
712                    qn = qc[n] #Neighbour's centroid value
713
714                    qmin[k] = min(qmin[k], qn)
715                    qmax[k] = max(qmax[k], qn)
716
717
718        #Diffences between centroids and maxima/minima
719        dqmax = qmax - qc
720        dqmin = qmin - qc
721
722        #Deltas between vertex and centroid values
723        dq = zeros(qv.shape, Float)
724        for i in range(2):
725            dq[:,i] = qv[:,i] - qc
726
727        #Phi limiter
728        for k in range(N):
729
730            #Find the gradient limiter (phi) across vertices
731            phi = 1.0
732            for i in range(2):
733                r = 1.0
734                if (dq[k,i] > 0): r = dqmax[k]/dq[k,i]
735                if (dq[k,i] < 0): r = dqmin[k]/dq[k,i]
736
737                phi = min( min(r*beta, 1), phi )
738
739            #Then update using phi limiter
740            for i in range(2):
741                qv[k,i] = qc[k] + phi*dq[k,i]
742
743       
744
745def newLinePlot(title='Simple Plot'):
746    import Gnuplot
747    g = Gnuplot.Gnuplot()
748    g.title(title)
749    g('set data style linespoints') 
750    g.xlabel('x')
751    g.ylabel('y')
752    return g
753
754def linePlot(g,x,y):
755    import Gnuplot
756    g.plot(Gnuplot.PlotItems.Data(x.flat,y.flat))
757
758
759
760
761
762if __name__ == "__main__":
763    from domain import Domain
764    from Numeric import arange
765   
766    points1 = [0.0, 1.0, 2.0, 3.0]
767    vertex_values = [[1.0,2.0],[4.0,5.0],[-1.0,2.0]]
768
769    D1 = Domain(points1)
770
771    Q1 = Conserved_quantity(D1, vertex_values)
772
773    print Q1.vertex_values
774    print Q1.centroid_values
775
776    new_vertex_values = [[2.0,1.0],[3.0,4.0],[-2.0,4.0]]
777
778    Q1.set_values(new_vertex_values)
779
780    print Q1.vertex_values
781    print Q1.centroid_values
782
783    new_centroid_values = [20,30,40]
784    Q1.set_values(new_centroid_values,'centroids')
785
786    print Q1.vertex_values
787    print Q1.centroid_values
788
789    class FunClass:
790        def __init__(self,value):
791            self.value = value
792
793        def __call__(self,x):
794            return self.value*(x**2)
795
796
797    fun = FunClass(1.0)
798    Q1.set_values(fun,'vertices')
799
800    print Q1.vertex_values
801    print Q1.centroid_values
802
803    Xc = Q1.domain.vertices
804    Qc = Q1.vertex_values
805    print Xc
806    print Qc
807
808    Qc[1,0] = 3
809
810    Q1.beta = 1.0
811    Q1.extrapolate_second_order()
812    Q1.limit()
813
814    g1 = newLinePlot('plot 1')
815    linePlot(g1,Xc,Qc)
816
817    points2 = arange(10)
818    D2 = Domain(points2)
819
820    Q2 = Conserved_quantity(D2)
821    Q2.set_values(fun,'vertices')
822    Xc = Q2.domain.vertices
823    Qc = Q2.vertex_values
824
825    g2 = newLinePlot('plot 2')
826    linePlot(g2,Xc,Qc)
827
828
829   
830    Q2.extrapolate_second_order()
831    Q2.limit()
832    Xc = Q2.domain.vertices
833    Qc = Q2.vertex_values
834
835    print Q2.centroid_values
836    print Qc
837    raw_input('press_return')
838   
839    g3 = newLinePlot('plot 3')
840    linePlot(g3,Xc,Qc)
841    raw_input('press return')
842
843
844    for i in range(10):
845        fun = FunClass(i/10.0)
846        Q2.set_values(fun,'vertices')
847        Qc = Q2.vertex_values
848        linePlot(g3,Xc,Qc)
849
850    raw_input('press return')
851
Note: See TracBrowser for help on using the repository browser.