source: trunk/anuga_work/development/sudi/sw_1d/periodic_waves/johns/quantity.py @ 7922

Last change on this file since 7922 was 7922, checked in by mungkasi, 14 years ago

Modifying codes so that arrays are compatible with numpy instead of Numeric.

File size: 30.8 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
18
19class Quantity:
20
21   
22    def __init__(self, domain, vertex_values=None):
23        #Initialise Quantity using optional vertex values.
24
25        from domain import Domain
26        from Numeric import Float
27        from numpy import array, zeros
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 Float
64        from numpy import zeros
65       
66        #Allocate space for boundary values
67        #L = len(domain.boundary)
68        self.boundary_values = zeros(2, Float) #assumes no parrellism
69
70        #Allocate space for updates of conserved quantities by
71        #flux calculations and forcing functions
72       
73        N = domain.number_of_elements
74        self.explicit_update = zeros(N, Float )
75        self.semi_implicit_update = zeros(N, Float )
76
77        self.gradients = zeros(N, Float)
78        self.qmax = zeros(self.centroid_values.shape, Float)
79        self.qmin = zeros(self.centroid_values.shape, Float)
80
81        self.beta = domain.beta       
82
83
84    def __len__(self):
85        """
86        Returns number of intervals.
87        """
88        return self.centroid_values.shape[0]
89   
90    def interpolate(self):
91        """
92        Compute interpolated values at centroid
93        Pre-condition: vertex_values have been set
94        """
95
96        N = self.vertex_values.shape[0]
97        for i in range(N):
98            v0 = self.vertex_values[i, 0]
99            v1 = self.vertex_values[i, 1]
100
101            self.centroid_values[i] = (v0 + v1)/2.0
102
103    def set_values(self, X, location='vertices'):
104        """Set values for quantity
105
106        X: Compatible list, Numeric array (see below), constant or function
107        location: Where values are to be stored.
108                  Permissible options are: vertices, centroid
109                  Default is "vertices"
110
111        In case of location == 'centroid' the dimension values must
112        be a list of a Numerical array of length N, N being the number
113        of elements in the mesh. Otherwise it must be of dimension Nx2
114
115        The values will be stored in elements following their
116        internal ordering.
117
118        If values are described a function, it will be evaluated at specified points
119
120        If selected location is vertices, values for centroid and edges
121        will be assigned interpolated values.
122        In any other case, only values for the specified locations
123        will be assigned and the others will be left undefined.
124        """
125
126        if location not in ['vertices', 'centroids']:
127            msg = 'Invalid location: %s, (possible choices vertices, centroids)' %location
128            raise msg
129
130        if X is None:
131            msg = 'Given values are None'
132            raise msg
133
134        import types
135
136        if callable(X):
137            #Use function specific method
138            self.set_function_values(X, location)
139        elif type(X) in [types.FloatType, types.IntType, types.LongType]:
140            if location == 'centroids':
141                self.centroid_values[:] = X
142            else:
143                self.vertex_values[:] = X
144
145        else:
146            #Use array specific method
147            self.set_array_values(X, location)
148
149        if location == 'vertices':
150            #Intialise centroid and edge values
151            self.interpolate()
152
153
154
155
156
157    def set_function_values(self, f, location='vertices'):
158        """Set values for quantity using specified function
159
160        f: x -> z Function where x and z are arrays
161        location: Where values are to be stored.
162                  Permissible options are: vertices, centroid
163                  Default is "vertices"
164        """
165       
166        if location == 'centroids':
167           
168            P = self.domain.centroids
169            self.set_values(f(P), location)
170        else:
171            #Vertices
172            P = self.domain.get_vertices()
173           
174            for i in range(2):
175                self.vertex_values[:,i] = f(P[:,i])       
176
177    def set_array_values(self, values, location='vertices'):
178        """Set values for quantity
179
180        values: Numeric array
181        location: Where values are to be stored.
182                  Permissible options are: vertices, centroid, edges
183                  Default is "vertices"
184
185        In case of location == 'centroid' the dimension values must
186        be a list of a Numerical array of length N, N being the number
187        of elements in the mesh. Otherwise it must be of dimension Nx2
188
189        The values will be stored in elements following their
190        internal ordering.
191
192        If selected location is vertices, values for centroid
193        will be assigned interpolated values.
194        In any other case, only values for the specified locations
195        will be assigned and the others will be left undefined.
196        """
197
198        from Numeric import Float
199        from numpy import array
200
201        values = array(values).astype(Float)
202
203        N = self.centroid_values.shape[0]
204
205        msg = 'Number of values must match number of elements'
206        assert values.shape[0] == N, msg
207
208        if location == 'centroids':
209            assert len(values.shape) == 1, 'Values array must be 1d'
210            self.centroid_values = values
211        #elif location == 'edges':
212        #    assert len(values.shape) == 2, 'Values array must be 2d'
213        #    msg = 'Array must be N x 2'
214        #    self.edge_values = values
215        else:
216            assert len(values.shape) == 2, 'Values array must be 2d'
217            msg = 'Array must be N x 2'
218            assert values.shape[1] == 2, msg
219
220            self.vertex_values = values
221
222
223    def get_values(self, location='vertices', indices = None):
224        """get values for quantity
225
226        return X, Compatible list, Numeric array (see below)
227        location: Where values are to be stored.
228                  Permissible options are: vertices, centroid
229                  and unique vertices. Default is 'vertices'
230
231        In case of location == 'centroids' the dimension values must
232        be a list of a Numerical array of length N, N being the number
233        of elements. Otherwise it must be of dimension Nx3
234
235        The returned values with be a list the length of indices
236        (N if indices = None).  Each value will be a list of the three
237        vertex values for this quantity.
238
239        Indices is the set of element ids that the operation applies to.
240
241        """
242        from numpy import take
243
244        if location not in ['vertices', 'centroids', 'unique vertices']:
245            msg = 'Invalid location: %s' %location
246            raise msg
247
248        import types, Numeric
249        assert type(indices) in [types.ListType, types.NoneType,
250                                 Numeric.ArrayType],\
251                                 'Indices must be a list or None'
252
253        if location == 'centroids':
254            if (indices ==  None):
255                indices = range(len(self))
256            return take(self.centroid_values,indices)
257        elif location == 'unique vertices':
258            if (indices ==  None):
259                indices=range(self.domain.coordinates.shape[0])
260            vert_values = []
261            #Go through list of unique vertices
262            for unique_vert_id in indices:
263                cells = self.domain.vertexlist[unique_vert_id]
264
265                #In case there are unused points
266                if cells is None:
267                    msg = 'Unique vertex not associated with cells'
268                    raise msg
269
270                # Go through all cells, vertex pairs
271                # Average the values
272                sum = 0
273                for cell_id, vertex_id in cells:
274                    sum += self.vertex_values[cell_id, vertex_id]
275                vert_values.append(sum/len(cells))
276            return Numeric.array(vert_values)
277        else:
278            if (indices ==  None):
279                indices = range(len(self))
280            return take(self.vertex_values,indices)
281
282
283    def get_vertex_values(self,
284                          x=True,
285                          smooth = None,
286                          precision = None,
287                          reduction = None):
288        """Return vertex values like an OBJ format
289
290        The vertex values are returned as one sequence in the 1D float array A.
291        If requested the coordinates will be returned in 1D arrays X.
292
293        The connectivity is represented as an integer array, V, of dimension
294        M x 2, where M is the number of volumes. Each row has two indices
295        into the X, A arrays defining the element.
296
297        if smooth is True, vertex values corresponding to one common
298        coordinate set will be smoothed according to the given
299        reduction operator. In this case vertex coordinates will be
300        de-duplicated.
301
302        If no smoothings is required, vertex coordinates and values will
303        be aggregated as a concatenation of values at
304        vertices 0, vertices 1
305
306
307        Calling convention
308        if x is True:
309           X,A,V = get_vertex_values
310        else:
311           A,V = get_vertex_values
312
313        """
314
315        from Numeric import Float, Int
316        from numpy import concatenate, zeros, array, reshape       
317
318
319        if smooth is None:
320            smooth = self.domain.smooth
321
322        if precision is None:
323            precision = Float
324
325        if reduction is None:
326            reduction = self.domain.reduction
327
328        #Create connectivity
329
330        if smooth == True:
331
332            V = self.domain.get_vertices()
333            N = len(self.domain.vertexlist)
334            #N = len(self.domain.vertices)
335            A = zeros(N, precision)
336
337            #Smoothing loop
338            for k in range(N):
339                L = self.domain.vertexlist[k]
340                #L = self.domain.vertices[k]
341
342                #Go through all triangle, vertex pairs
343                #contributing to vertex k and register vertex value
344
345                if L is None: continue #In case there are unused points
346
347                contributions = []
348                for volume_id, vertex_id in L:
349                    v = self.vertex_values[volume_id, vertex_id]
350                    contributions.append(v)
351
352                A[k] = reduction(contributions)
353
354            if x is True:
355                 #X = self.domain.coordinates[:,0].astype(precision)
356                 X = self.domain.coordinates[:].astype(precision)
357                 #Y = self.domain.coordinates[:,1].astype(precision)
358
359                 #return X, Y, A, V
360                 return X, A, V
361           
362            #else:
363            return A, V
364        else:
365            #Don't smooth
366            #obj machinery moved to general_mesh
367
368            # Create a V like [[0 1 2], [3 4 5]....[3*m-2 3*m-1 3*m]]
369            # These vert_id's will relate to the verts created below
370            #m = len(self.domain)  #Number of volumes
371            #M = 3*m        #Total number of unique vertices
372            #V = reshape(array(range(M)).astype(Int), (m,3))
373
374            #V = self.domain.get_triangles(obj=True)
375            V = self.domain.get_vertices
376            #FIXME use get_vertices, when ready
377
378            A = self.vertex_values.flat
379
380            #Do vertex coordinates
381            if x is True:
382                X = self.domain.get_vertex_coordinates()
383
384                #X = C[:,0:6:2].copy()
385                #Y = C[:,1:6:2].copy()
386
387                return X.flat, A, V
388            else:
389                return A, V
390
391    def get_integral(self):
392        """Compute the integral of quantity across entire domain
393        """
394        integral = 0
395        for k in range(self.domain.number_of_elements):
396            area = self.domain.areas[k]
397            qc = self.centroid_values[k]
398            integral += qc*area
399
400        return integral
401
402
403    def update(self, timestep):
404        """Update centroid values based on values stored in
405        explicit_update and semi_implicit_update as well as given timestep
406        """
407
408        from Numeric import Float
409        from numpy import sum, equal, ones
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        #print 'compute_gradient'
433
434        from Numeric import Float
435        from numpy import array, zeros
436
437        N = self.centroid_values.shape[0]
438
439
440        G = self.gradients
441        Q = self.centroid_values
442        X = self.domain.centroids
443
444        for k in range(N):
445
446            # first and last elements have boundaries
447
448            if k == 0:
449
450                #Get data
451                k0 = k
452                k1 = k+1
453                k2 = k+2
454
455                q0 = Q[k0]
456                q1 = Q[k1]
457                q2 = Q[k2]
458
459                x0 = X[k0] #V0 centroid
460                x1 = X[k1] #V1 centroid
461                x2 = X[k2]
462
463                #Gradient
464                #G[k] = (q1 - q0)/(x1 - x0)
465               
466                G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0)
467                G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1)
468
469            elif k == N-1:
470
471                #Get data
472                k0 = k
473                k1 = k-1
474                k2 = k-2
475
476                q0 = Q[k0]
477                q1 = Q[k1]
478                q2 = Q[k2]
479
480                x0 = X[k0] #V0 centroid
481                x1 = X[k1] #V1 centroid
482                x2 = X[k2]
483
484                #Gradient
485                #G[k] = (q1 - q0)/(x1 - x0)
486               
487                G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0)
488                G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1)
489
490##                q0 = Q[k0]
491##                q1 = Q[k1]
492##
493##                x0 = X[k0] #V0 centroid
494##                x1 = X[k1] #V1 centroid
495##
496##                #Gradient
497##                G[k] = (q1 - q0)/(x1 - x0)
498
499            else:
500                #Interior Volume (2 neighbours)
501
502                #Get data
503                k0 = k-1
504                k2 = k+1
505
506                q0 = Q[k0]
507                q1 = Q[k]
508                q2 = Q[k2]
509
510                x0 = X[k0] #V0 centroid
511                x1 = X[k]  #V1 centroid (Self)
512                x2 = X[k2] #V2 centroid
513
514                #Gradient
515                #G[k] = (q2-q0)/(x2-x0)
516                G[k] = ((q0-q1)/(x0-x1)*(x2-x1) - (q2-q1)/(x2-x1)*(x0-x1))/(x2-x0)
517
518
519    def compute_minmod_gradients(self):
520        """Compute gradients of piecewise linear function defined by centroids of
521        neighbouring volumes.
522        """
523
524        #print 'compute_minmod_gradients'
525       
526        from Numeric import Float,sign
527        from numpy import array, zeros
528       
529        def xmin(a,b):
530            return 0.5*(sign(a)+sign(b))*min(abs(a),abs(b))
531
532        def xmic(t,a,b):
533            return xmin(t*xmin(a,b), 0.50*(a+b) )
534
535
536
537        N = self.centroid_values.shape[0]
538
539
540        G = self.gradients
541        Q = self.centroid_values
542        X = self.domain.centroids
543
544        for k in range(N):
545
546            # first and last elements have boundaries
547
548            if k == 0:
549
550                #Get data
551                k0 = k
552                k1 = k+1
553                k2 = k+2
554
555                q0 = Q[k0]
556                q1 = Q[k1]
557                q2 = Q[k2]
558
559                x0 = X[k0] #V0 centroid
560                x1 = X[k1] #V1 centroid
561                x2 = X[k2]
562
563                #Gradient
564                #G[k] = (q1 - q0)/(x1 - x0)
565               
566                G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0)
567                G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1)
568
569            elif k == N-1:
570
571                #Get data
572                k0 = k
573                k1 = k-1
574                k2 = k-2
575
576                q0 = Q[k0]
577                q1 = Q[k1]
578                q2 = Q[k2]
579
580                x0 = X[k0] #V0 centroid
581                x1 = X[k1] #V1 centroid
582                x2 = X[k2]
583
584                #Gradient
585                #G[k] = (q1 - q0)/(x1 - x0)
586               
587                G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0)
588                G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1)
589
590##                #Get data
591##                k0 = k
592##                k1 = k-1
593##
594##                q0 = Q[k0]
595##                q1 = Q[k1]
596##
597##                x0 = X[k0] #V0 centroid
598##                x1 = X[k1] #V1 centroid
599##
600##                #Gradient
601##                G[k] = (q1 - q0)/(x1 - x0)
602
603            elif (self.domain.wet_nodes[k,0] == 2) & (self.domain.wet_nodes[k,1] == 2):
604                G[k] = 0.0
605
606            else:
607                #Interior Volume (2 neighbours)
608
609                #Get data
610                k0 = k-1
611                k2 = k+1
612
613                q0 = Q[k0]
614                q1 = Q[k]
615                q2 = Q[k2]
616
617                x0 = X[k0] #V0 centroid
618                x1 = X[k]  #V1 centroid (Self)
619                x2 = X[k2] #V2 centroid
620
621                # assuming uniform grid
622                d1 = (q1 - q0)/(x1-x0)
623                d2 = (q2 - q1)/(x2-x1)
624
625                #Gradient
626                #G[k] = (d1+d2)*0.5
627                #G[k] = (d1*(x2-x1) - d2*(x0-x1))/(x2-x0)               
628                G[k] = xmic( self.domain.beta, d1, d2 )
629       
630
631    def extrapolate_first_order(self):
632        """Extrapolate conserved quantities from centroid to
633        vertices for each volume using
634        first order scheme.
635        """
636
637        qc = self.centroid_values
638        qv = self.vertex_values
639
640        for i in range(2):
641            qv[:,i] = qc
642
643
644    def extrapolate_second_order(self):
645        """Extrapolate conserved quantities from centroid to
646        vertices for each volume using
647        second order scheme.
648        """
649        if self.domain.limiter == "pyvolution":
650            #Z = self.gradients
651            #print "gradients 1",Z
652            self.compute_gradients()
653            #print "gradients 2",Z
654
655            #Z = self.gradients
656            #print "gradients 1",Z
657            #self.compute_minmod_gradients()
658            #print "gradients 2", Z
659
660            G = self.gradients
661            V = self.domain.vertices
662            qc = self.centroid_values
663            qv = self.vertex_values       
664
665            #Check each triangle
666            for k in range(self.domain.number_of_elements):
667                #Centroid coordinates
668                x = self.domain.centroids[k]
669
670                #vertex coordinates
671                x0, x1 = V[k,:]
672
673                #Extrapolate
674                qv[k,0] = qc[k] + G[k]*(x0-x)
675                qv[k,1] = qc[k] + G[k]*(x1-x)
676            self.limit_pyvolution()
677        elif self.domain.limiter == "minmod_steve":
678            self.limit_minmod()
679        else:
680            self.limit_range()
681       
682       
683
684    def limit_minmod(self):
685        #Z = self.gradients
686        #print "gradients 1",Z
687        self.compute_minmod_gradients()
688        #print "gradients 2", Z
689
690        G = self.gradients
691        V = self.domain.vertices
692        qc = self.centroid_values
693        qv = self.vertex_values       
694       
695        #Check each triangle
696        for k in range(self.domain.number_of_elements):
697            #Centroid coordinates
698            x = self.domain.centroids[k]
699
700            #vertex coordinates
701            x0, x1 = V[k,:]
702
703            #Extrapolate
704            qv[k,0] = qc[k] + G[k]*(x0-x)
705            qv[k,1] = qc[k] + G[k]*(x1-x)
706
707 
708    def limit_pyvolution(self):
709        """
710        Limit slopes for each volume to eliminate artificial variance
711        introduced by e.g. second order extrapolator
712
713        This is an unsophisticated limiter as it does not take into
714        account dependencies among quantities.
715
716        precondition:
717        vertex values are estimated from gradient
718        postcondition:
719        vertex values are updated
720        """
721        from Numeric import Float
722        from numpy import zeros
723
724        N = self.domain.number_of_elements
725        beta = self.domain.beta
726        #beta = 0.8
727
728        qc = self.centroid_values
729        qv = self.vertex_values
730
731        #Find min and max of this and neighbour's centroid values
732        qmax = self.qmax
733        qmin = self.qmin
734
735        for k in range(N):
736            qmax[k] = qmin[k] = qc[k]
737            for i in range(2):
738                n = self.domain.neighbours[k,i]
739                if n >= 0:
740                    qn = qc[n] #Neighbour's centroid value
741
742                    qmin[k] = min(qmin[k], qn)
743                    qmax[k] = max(qmax[k], qn)
744
745
746        #Diffences between centroids and maxima/minima
747        dqmax = qmax - qc
748        dqmin = qmin - qc
749
750        #Deltas between vertex and centroid values
751        dq = zeros(qv.shape, Float)
752        for i in range(2):
753            dq[:,i] = qv[:,i] - qc
754
755        #Phi limiter
756        for k in range(N):
757
758            #Find the gradient limiter (phi) across vertices
759            phi = 1.0
760            for i in range(2):
761                r = 1.0
762                if (dq[k,i] > 0): r = dqmax[k]/dq[k,i]
763                if (dq[k,i] < 0): r = dqmin[k]/dq[k,i]
764
765                phi = min( min(r*beta, 1), phi )
766
767            #Then update using phi limiter
768            for i in range(2):
769                qv[k,i] = qc[k] + phi*dq[k,i]
770
771    def limit_range(self):
772        import sys
773        from Numeric import Float
774        from numpy import zeros
775        from util import minmod, minmod_kurganov, maxmod, vanleer, vanalbada
776        limiter = self.domain.limiter
777        #print limiter
778       
779        #print 'limit_range'
780        N = self.domain.number_of_elements
781        qc = self.centroid_values
782        qv = self.vertex_values
783        C = self.domain.centroids
784        X = self.domain.vertices
785        beta_p = zeros(N,Float)
786        beta_m = zeros(N,Float)
787        beta_x = zeros(N,Float)
788       
789        for k in range(N):
790            #Note: neighbours=[-1, {0, 1, 2, ..., N-1}, -2]
791            n0 = self.domain.neighbours[k,0]
792            n1 = self.domain.neighbours[k,1]           
793            if ( n0 >= 0) & (n1 >= 0):
794                #SLOPE DERIVATIVE LIMIT
795                beta_p[k] = (qc[k]-qc[k-1])/(C[k]-C[k-1])
796                beta_m[k] = (qc[k+1]-qc[k])/(C[k+1]-C[k])
797                beta_x[k] = (qc[k+1]-qc[k-1])/(C[k+1]-C[k-1])
798       
799        dq = zeros(qv.shape, Float)
800        for i in range(2):
801            dq[:,i] =self.domain.vertices[:,i]-self.domain.centroids
802           
803        #Phi limiter
804        for k in range(N):
805            #Note: neighbours=[-1, {0, 1, 2, ..., N-1}, -2]           
806            n0 = self.domain.neighbours[k,0]
807            n1 = self.domain.neighbours[k,1]
808            if n0 < 0:
809                phi = (qc[k+1] - qc[k])/(C[k+1] - C[k])
810            elif n1 < 0:
811                phi = (qc[k] - qc[k-1])/(C[k] - C[k-1])
812            else:
813                if limiter == "minmod":
814                    phi = minmod(beta_p[k],beta_m[k])
815                elif limiter == "minmod_kurganov":
816                    # Also known as monotonized central difference limiter if theta = 2.0
817                    theta = 2.0 
818                    phi = minmod_kurganov(theta*beta_p[k],theta*beta_m[k],beta_x[k])
819                elif limiter == "superbee":
820                    slope1 = minmod(beta_m[k],2.0*beta_p[k])
821                    slope2 = minmod(2.0*beta_m[k],beta_p[k])
822                    phi = maxmod(slope1,slope2)
823                elif limiter == "vanleer":
824                    phi = vanleer(beta_p[k],beta_m[k])
825                elif limiter == "vanalbada":
826                    phi = vanalbada(beta_m[k],beta_p[k])
827            for i in range(2):
828                qv[k,i] = qc[k] + phi*dq[k,i]
829
830    def limit_steve_slope(self):   
831
832        import sys
833        from Numeric import Float
834        from numpy import zeros
835        from util import minmod, minmod_kurganov, maxmod, vanleer
836
837        N = self.domain.number_of_elements
838        limiter = self.domain.limiter
839        limiter_type = self.domain.limiter_type
840           
841        qc = self.centroid_values
842        qv = self.vertex_values
843
844        #Find min and max of this and neighbour's centroid values
845        beta_p = zeros(N,Float)
846        beta_m = zeros(N,Float)
847        beta_x = zeros(N,Float)
848        C = self.domain.centroids
849        X = self.domain.vertices
850
851        for k in range(N):
852       
853            n0 = self.domain.neighbours[k,0]
854            n1 = self.domain.neighbours[k,1]
855           
856            if (n0 >= 0) & (n1 >= 0):
857                # Check denominator not zero
858                if (qc[k+1]-qc[k]) == 0.0:
859                    beta_p[k] = float(sys.maxint)
860                    beta_m[k] = float(sys.maxint)
861                else:
862                    #STEVE LIMIT
863                    beta_p[k] = (qc[k]-qc[k-1])/(qc[k+1]-qc[k])
864                    beta_m[k] = (qc[k+2]-qc[k+1])/(qc[k+1]-qc[k])
865
866        #Deltas between vertex and centroid values
867        dq = zeros(qv.shape, Float)
868        for i in range(2):
869            dq[:,i] =self.domain.vertices[:,i]-self.domain.centroids
870           
871        #Phi limiter
872        for k in range(N):
873               
874            phi = 0.0
875            if limiter == "flux_minmod":
876                #FLUX MINMOD
877                phi = minmod_kurganov(1.0,beta_m[k],beta_p[k])
878            elif limiter == "flux_superbee":
879                #FLUX SUPERBEE
880                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
881            elif limiter == "flux_muscl":
882                #FLUX MUSCL
883                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])))
884            elif limiter == "flux_vanleer":
885                #FLUX VAN LEER
886                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
887               
888                #Then update using phi limiter
889                n = self.domain.neighbours[k,1]
890                if n>=0:
891                    #qv[k,0] = qc[k] - 0.5*phi*(qc[k+1]-qc[k])
892                    #qv[k,1] = qc[k] + 0.5*phi*(qc[k+1]-qc[k])
893                    qv[k,0] = qc[k] + 0.5*phi*(qv[k,0]-qc[k])
894                    qv[k,1] = qc[k] + 0.5*phi*(qv[k,1]-qc[k])
895                else:
896                    qv[k,i] = qc[k]
897
898    def backup_centroid_values(self):
899        # Call correct module function
900        # (either from this module or C-extension)
901        #backup_centroid_values(self)
902
903        self.centroid_backup_values[:] = (self.centroid_values).astype('f')
904
905    def saxpy_centroid_values(self,a,b):
906        # Call correct module function
907        # (either from this module or C-extension)
908        self.centroid_values[:] = (a*self.centroid_values + b*self.centroid_backup_values).astype('f')
909       
910class Conserved_quantity(Quantity):
911    """Class conserved quantity adds to Quantity:
912
913    storage and method for updating, and
914    methods for extrapolation from centropid to vertices inluding
915    gradients and limiters
916    """
917
918    def __init__(self, domain, vertex_values=None):
919        Quantity.__init__(self, domain, vertex_values)
920
921        print "Use Quantity instead of Conserved_quantity"
922
923
924##
925##def newLinePlot(title='Simple Plot'):
926##    import Gnuplot
927##    g = Gnuplot.Gnuplot()
928##    g.title(title)
929##    g('set data style linespoints')
930##    g.xlabel('x')
931##    g.ylabel('y')
932##    return g
933##
934##def linePlot(g,x,y):
935##    import Gnuplot
936##    g.plot(Gnuplot.PlotItems.Data(x.flat,y.flat))
937
938def newLinePlot(title='Simple Plot'):
939    import pylab as g
940    g.ion()
941    g.hold(False)
942    g.title(title)
943    g.xlabel('x')
944    g.ylabel('y')
945   
946
947def linePlot(x,y):
948    import pylab as g
949    g.plot(x.flat,y.flat)
950
951
952def closePlots():
953    import pylab as g
954    g.close('all')
955   
956if __name__ == "__main__":
957    #from domain import Domain
958    from domain_johns import Domain   
959    from numpy import arange
960   
961    points1 = [0.0, 1.0, 2.0, 3.0]
962    vertex_values = [[1.0,2.0],[4.0,5.0],[-1.0,2.0]]
963
964    D1 = Domain(points1)
965
966    Q1 = Quantity(D1, vertex_values)
967
968    print Q1.vertex_values
969    print Q1.centroid_values
970
971    new_vertex_values = [[2.0,1.0],[3.0,4.0],[-2.0,4.0]]
972
973    Q1.set_values(new_vertex_values)
974
975    print Q1.vertex_values
976    print Q1.centroid_values
977
978    new_centroid_values = [20,30,40]
979    Q1.set_values(new_centroid_values,'centroids')
980
981    print Q1.vertex_values
982    print Q1.centroid_values
983
984    class FunClass:
985        def __init__(self,value):
986            self.value = value
987
988        def __call__(self,x):
989            return self.value*(x**2)
990
991
992    fun = FunClass(1.0)
993    Q1.set_values(fun,'vertices')
994
995    print Q1.vertex_values
996    print Q1.centroid_values
997
998    Xc = Q1.domain.vertices
999    Qc = Q1.vertex_values
1000    print Xc
1001    print Qc
1002
1003    Qc[1,0] = 3
1004
1005    Q1.extrapolate_second_order()
1006    #Q1.limit_minmod()
1007
1008    newLinePlot('plots')
1009    linePlot(Xc,Qc)
1010    raw_input('press return')
1011
1012    points2 = arange(10)
1013    D2 = Domain(points2)
1014
1015    Q2 = Quantity(D2)
1016    Q2.set_values(fun,'vertices')
1017    Xc = Q2.domain.vertices
1018    Qc = Q2.vertex_values
1019    linePlot(Xc,Qc)
1020    raw_input('press return')
1021
1022   
1023    Q2.extrapolate_second_order()
1024    #Q2.limit_minmod()
1025    Xc = Q2.domain.vertices
1026    Qc = Q2.vertex_values
1027    print Q2.centroid_values
1028    print Qc
1029    linePlot(Xc,Qc)
1030    raw_input('press return')
1031
1032
1033    for i in range(10):
1034        import pylab as g
1035        g.hold(True)
1036        fun = FunClass(i/10.0)
1037        Q2.set_values(fun,'centroids')
1038        Q2.extrapolate_second_order()
1039        #Q2.limit_minmod()
1040        Qc = Q2.vertex_values
1041        linePlot(Xc,Qc)
1042        raw_input('press return')
1043
1044    raw_input('press return to quit')
1045closePlots()
Note: See TracBrowser for help on using the repository browser.