# source:anuga_work/development/anuga_1d/quantity.py@6453

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

File size: 30.7 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 array, zeros, Float
27
28        msg = 'First argument in Quantity.__init__ '
29        msg += 'must be of class Domain (or a subclass thereof)'
30        assert isinstance(domain, Domain), msg
31
32        if vertex_values is None:
33            N = domain.number_of_elements
34            self.vertex_values = zeros((N, 2), Float)
35        else:
36            self.vertex_values = array(vertex_values, Float)
37
38            N, V = self.vertex_values.shape
39            assert V == 2,\
40                   'Two vertex values per element must be specified'
41
42
43            msg = 'Number of vertex values (%d) must be consistent with'\
44                  %N
45            msg += 'number of elements in specified domain (%d).'\
46                   %domain.number_of_elements
47
48            assert N == domain.number_of_elements, msg
49
50        self.domain = domain
51
52        #Allocate space for other quantities
53        self.centroid_values = zeros(N, Float)
54        self.centroid_backup_values = zeros(N, Float)
55        #self.edge_values = zeros((N, 2), Float)
56        #edge values are values of the ends of each interval
57
58        #Intialise centroid values
59        self.interpolate()
60
61
62        from Numeric import zeros, Float
63
64        #Allocate space for boundary values
65        #L = len(domain.boundary)
66        self.boundary_values = zeros(2, Float) #assumes no parrellism
67
68        #Allocate space for updates of conserved quantities by
69        #flux calculations and forcing functions
70
71        N = domain.number_of_elements
72        self.explicit_update = zeros(N, Float )
73        self.semi_implicit_update = zeros(N, Float )
74
76        self.qmax = zeros(self.centroid_values.shape, Float)
77        self.qmin = zeros(self.centroid_values.shape, Float)
78
79        self.beta = domain.beta
80
81
82    def __len__(self):
83        """
84        Returns number of intervals.
85        """
86        return self.centroid_values.shape[0]
87
88    def interpolate(self):
89        """
90        Compute interpolated values at centroid
91        Pre-condition: vertex_values have been set
92        """
93
94        N = self.vertex_values.shape[0]
95        for i in range(N):
96            v0 = self.vertex_values[i, 0]
97            v1 = self.vertex_values[i, 1]
98
99            self.centroid_values[i] = (v0 + v1)/2.0
100
101    def set_values(self, X, location='vertices'):
102        """Set values for quantity
103
104        X: Compatible list, Numeric array (see below), constant or function
105        location: Where values are to be stored.
106                  Permissible options are: vertices, centroid
107                  Default is "vertices"
108
109        In case of location == 'centroid' the dimension values must
110        be a list of a Numerical array of length N, N being the number
111        of elements in the mesh. Otherwise it must be of dimension Nx2
112
113        The values will be stored in elements following their
114        internal ordering.
115
116        If values are described a function, it will be evaluated at specified points
117
118        If selected location is vertices, values for centroid and edges
119        will be assigned interpolated values.
120        In any other case, only values for the specified locations
121        will be assigned and the others will be left undefined.
122        """
123
124        if location not in ['vertices', 'centroids']:
125            msg = 'Invalid location: %s, (possible choices vertices, centroids)' %location
126            raise msg
127
128        if X is None:
129            msg = 'Given values are None'
130            raise msg
131
132        import types
133
134        if callable(X):
135            #Use function specific method
136            self.set_function_values(X, location)
137
138        elif type(X) in [types.FloatType, types.IntType, types.LongType]:
139            if location == 'centroids':
140                self.centroid_values[:] = X
141            else:
142                self.vertex_values[:] = X
143
144        else:
145            #Use array specific method
146            self.set_array_values(X, location)
147
148        if location == 'vertices':
149            #Intialise centroid and edge values
150            self.interpolate()
151
152
153
154
155
156    def set_function_values(self, f, location='vertices'):
157        """Set values for quantity using specified function
158
159        f: x -> z Function where x and z are arrays
160        location: Where values are to be stored.
161                  Permissible options are: vertices, centroid
162                  Default is "vertices"
163        """
164
165        if location == 'centroids':
166
167            P = self.domain.centroids
168            self.set_values(f(P), location)
169        else:
170            #Vertices
171
172            P = self.domain.get_vertices()
173
174            for i in range(2):
175
176                self.vertex_values[:,i] = f(P[:,i])
177
178    def set_array_values(self, values, location='vertices'):
179        """Set values for quantity
180
181        values: Numeric array
182        location: Where values are to be stored.
183                  Permissible options are: vertices, centroid, edges
184                  Default is "vertices"
185
186        In case of location == 'centroid' the dimension values must
187        be a list of a Numerical array of length N, N being the number
188        of elements in the mesh. Otherwise it must be of dimension Nx2
189
190        The values will be stored in elements following their
191        internal ordering.
192
193        If selected location is vertices, values for centroid
194        will be assigned interpolated values.
195        In any other case, only values for the specified locations
196        will be assigned and the others will be left undefined.
197        """
198
199        from Numeric import array, Float
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 Numeric 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 concatenate, zeros, Float, Int, array, reshape
316
317
318        if smooth is None:
319            smooth = self.domain.smooth
320
321        if precision is None:
322            precision = Float
323
324        if reduction is None:
325            reduction = self.domain.reduction
326
327        #Create connectivity
328
329        if smooth == True:
330
331            V = self.domain.get_vertices()
332            N = len(self.domain.vertexlist)
333            #N = len(self.domain.vertices)
334            A = zeros(N, precision)
335
336            #Smoothing loop
337            for k in range(N):
338                L = self.domain.vertexlist[k]
339                #L = self.domain.vertices[k]
340
341                #Go through all triangle, vertex pairs
342                #contributing to vertex k and register vertex value
343
344                if L is None: continue #In case there are unused points
345
346                contributions = []
347                for volume_id, vertex_id in L:
348                    v = self.vertex_values[volume_id, vertex_id]
349                    contributions.append(v)
350
351                A[k] = reduction(contributions)
352
353            if x is True:
354                 #X = self.domain.coordinates[:,0].astype(precision)
355                 X = self.domain.coordinates[:].astype(precision)
356                 #Y = self.domain.coordinates[:,1].astype(precision)
357
358                 #return X, Y, A, V
359                 return X, A, V
360
361            #else:
362            return A, V
363        else:
364            #Don't smooth
365            #obj machinery moved to general_mesh
366
367            # Create a V like [[0 1 2], [3 4 5]....[3*m-2 3*m-1 3*m]]
368            # These vert_id's will relate to the verts created below
369            #m = len(self.domain)  #Number of volumes
370            #M = 3*m        #Total number of unique vertices
371            #V = reshape(array(range(M)).astype(Int), (m,3))
372
373            #V = self.domain.get_triangles(obj=True)
374            V = self.domain.get_vertices
375            #FIXME use get_vertices, when ready
376
377            A = self.vertex_values.flat
378
379            #Do vertex coordinates
380            if x is True:
381                X = self.domain.get_vertex_coordinates()
382
383                #X = C[:,0:6:2].copy()
384                #Y = C[:,1:6:2].copy()
385
386                return X.flat, A, V
387            else:
388                return A, V
389
390    def get_integral(self):
391        """Compute the integral of quantity across entire domain
392        """
393        integral = 0
394        for k in range(self.domain.number_of_elements):
395            area = self.domain.areas[k]
396            qc = self.centroid_values[k]
397            integral += qc*area
398
399        return integral
400
401
402    def update(self, timestep):
403        """Update centroid values based on values stored in
404        explicit_update and semi_implicit_update as well as given timestep
405        """
406
407        from Numeric import sum, equal, ones, Float
408
409        N = self.centroid_values.shape[0]
410
412        self.centroid_values += timestep*self.explicit_update
413
415        denominator = ones(N, Float)-timestep*self.semi_implicit_update
416
417        if sum(equal(denominator, 0.0)) > 0.0:
418            msg = 'Zero division in semi implicit update. Call Stephen :-)'
419            raise msg
420        else:
421            #Update conserved_quantities from semi implicit updates
422            self.centroid_values /= denominator
423
424
426        """Compute gradients of piecewise linear function defined by centroids of
427        neighbouring volumes.
428        """
429
431
432        from Numeric import array, zeros, Float
433
434        N = self.centroid_values.shape[0]
435
436
438        Q = self.centroid_values
439        X = self.domain.centroids
440
441        for k in range(N):
442
443            # first and last elements have boundaries
444
445            if k == 0:
446
447                #Get data
448                k0 = k
449                k1 = k+1
450                k2 = k+2
451
452                q0 = Q[k0]
453                q1 = Q[k1]
454                q2 = Q[k2]
455
456                x0 = X[k0] #V0 centroid
457                x1 = X[k1] #V1 centroid
458                x2 = X[k2]
459
461                #G[k] = (q1 - q0)/(x1 - x0)
462
463                G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0)
464                G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1)
465
466            elif k == N-1:
467
468                #Get data
469                k0 = k
470                k1 = k-1
471                k2 = k-2
472
473                q0 = Q[k0]
474                q1 = Q[k1]
475                q2 = Q[k2]
476
477                x0 = X[k0] #V0 centroid
478                x1 = X[k1] #V1 centroid
479                x2 = X[k2]
480
482                #G[k] = (q1 - q0)/(x1 - x0)
483
484                G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0)
485                G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1)
486
487##                q0 = Q[k0]
488##                q1 = Q[k1]
489##
490##                x0 = X[k0] #V0 centroid
491##                x1 = X[k1] #V1 centroid
492##
494##                G[k] = (q1 - q0)/(x1 - x0)
495
496            else:
497                #Interior Volume (2 neighbours)
498
499                #Get data
500                k0 = k-1
501                k2 = k+1
502
503                q0 = Q[k0]
504                q1 = Q[k]
505                q2 = Q[k2]
506
507                x0 = X[k0] #V0 centroid
508                x1 = X[k]  #V1 centroid (Self)
509                x2 = X[k2] #V2 centroid
510
512                #G[k] = (q2-q0)/(x2-x0)
513                G[k] = ((q0-q1)/(x0-x1)*(x2-x1) - (q2-q1)/(x2-x1)*(x0-x1))/(x2-x0)
514
515
517        """Compute gradients of piecewise linear function defined by centroids of
518        neighbouring volumes.
519        """
520
522
523        from Numeric import array, zeros, Float,sign
524
525        def xmin(a,b):
526            return 0.5*(sign(a)+sign(b))*min(abs(a),abs(b))
527
528        def xmic(t,a,b):
529            return xmin(t*xmin(a,b), 0.50*(a+b) )
530
531
532
533        N = self.centroid_values.shape[0]
534
535
537        Q = self.centroid_values
538        X = self.domain.centroids
539
540        for k in range(N):
541
542            # first and last elements have boundaries
543
544            if k == 0:
545
546                #Get data
547                k0 = k
548                k1 = k+1
549                k2 = k+2
550
551                q0 = Q[k0]
552                q1 = Q[k1]
553                q2 = Q[k2]
554
555                x0 = X[k0] #V0 centroid
556                x1 = X[k1] #V1 centroid
557                x2 = X[k2]
558
560                #G[k] = (q1 - q0)/(x1 - x0)
561
562                G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0)
563                G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1)
564
565            elif k == N-1:
566
567                #Get data
568                k0 = k
569                k1 = k-1
570                k2 = k-2
571
572                q0 = Q[k0]
573                q1 = Q[k1]
574                q2 = Q[k2]
575
576                x0 = X[k0] #V0 centroid
577                x1 = X[k1] #V1 centroid
578                x2 = X[k2]
579
581                #G[k] = (q1 - q0)/(x1 - x0)
582
583                G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0)
584                G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1)
585
586##                #Get data
587##                k0 = k
588##                k1 = k-1
589##
590##                q0 = Q[k0]
591##                q1 = Q[k1]
592##
593##                x0 = X[k0] #V0 centroid
594##                x1 = X[k1] #V1 centroid
595##
597##                G[k] = (q1 - q0)/(x1 - x0)
598
599            elif (self.domain.wet_nodes[k,0] == 2) & (self.domain.wet_nodes[k,1] == 2):
600                G[k] = 0.0
601
602            else:
603                #Interior Volume (2 neighbours)
604
605                #Get data
606                k0 = k-1
607                k2 = k+1
608
609                q0 = Q[k0]
610                q1 = Q[k]
611                q2 = Q[k2]
612
613                x0 = X[k0] #V0 centroid
614                x1 = X[k]  #V1 centroid (Self)
615                x2 = X[k2] #V2 centroid
616
617                # assuming uniform grid
618                d1 = (q1 - q0)/(x1-x0)
619                d2 = (q2 - q1)/(x2-x1)
620
622                #G[k] = (d1+d2)*0.5
623                #G[k] = (d1*(x2-x1) - d2*(x0-x1))/(x2-x0)
624                G[k] = xmic( self.domain.beta, d1, d2 )
625
626
627    def extrapolate_first_order(self):
628        """Extrapolate conserved quantities from centroid to
629        vertices for each volume using
630        first order scheme.
631        """
632
633        qc = self.centroid_values
634        qv = self.vertex_values
635
636        for i in range(2):
637            qv[:,i] = qc
638
639
640    def extrapolate_second_order(self):
641        """Extrapolate conserved quantities from centroid to
642        vertices for each volume using
643        second order scheme.
644        """
645        if self.domain.limiter == "pyvolution":
650
655
657            V = self.domain.vertices
658            qc = self.centroid_values
659            qv = self.vertex_values
660
661            #Check each triangle
662            for k in range(self.domain.number_of_elements):
663                #Centroid coordinates
664                x = self.domain.centroids[k]
665
666                #vertex coordinates
667                x0, x1 = V[k,:]
668
669                #Extrapolate
670                qv[k,0] = qc[k] + G[k]*(x0-x)
671                qv[k,1] = qc[k] + G[k]*(x1-x)
672            self.limit_pyvolution()
673        elif self.domain.limiter == "minmod_steve":
674            self.limit_minmod()
675        else:
676            self.limit_range()
677
678
679
680    def limit_minmod(self):
685
687        V = self.domain.vertices
688        qc = self.centroid_values
689        qv = self.vertex_values
690
691        #Check each triangle
692        for k in range(self.domain.number_of_elements):
693            #Centroid coordinates
694            x = self.domain.centroids[k]
695
696            #vertex coordinates
697            x0, x1 = V[k,:]
698
699            #Extrapolate
700            qv[k,0] = qc[k] + G[k]*(x0-x)
701            qv[k,1] = qc[k] + G[k]*(x1-x)
702
703
704    def limit_pyvolution(self):
705        """
706        Limit slopes for each volume to eliminate artificial variance
707        introduced by e.g. second order extrapolator
708
709        This is an unsophisticated limiter as it does not take into
710        account dependencies among quantities.
711
712        precondition:
713        vertex values are estimated from gradient
714        postcondition:
715        vertex values are updated
716        """
717        from Numeric import zeros, Float
718
719        N = self.domain.number_of_elements
720        beta = self.domain.beta
721        #beta = 0.8
722
723        qc = self.centroid_values
724        qv = self.vertex_values
725
726        #Find min and max of this and neighbour's centroid values
727        qmax = self.qmax
728        qmin = self.qmin
729
730        for k in range(N):
731            qmax[k] = qmin[k] = qc[k]
732            for i in range(2):
733                n = self.domain.neighbours[k,i]
734                if n >= 0:
735                    qn = qc[n] #Neighbour's centroid value
736
737                    qmin[k] = min(qmin[k], qn)
738                    qmax[k] = max(qmax[k], qn)
739
740
741        #Diffences between centroids and maxima/minima
742        dqmax = qmax - qc
743        dqmin = qmin - qc
744
745        #Deltas between vertex and centroid values
746        dq = zeros(qv.shape, Float)
747        for i in range(2):
748            dq[:,i] = qv[:,i] - qc
749
750        #Phi limiter
751        for k in range(N):
752
753            #Find the gradient limiter (phi) across vertices
754            phi = 1.0
755            for i in range(2):
756                r = 1.0
757                if (dq[k,i] > 0): r = dqmax[k]/dq[k,i]
758                if (dq[k,i] < 0): r = dqmin[k]/dq[k,i]
759
760                phi = min( min(r*beta, 1), phi )
761
762            #Then update using phi limiter
763            for i in range(2):
764                qv[k,i] = qc[k] + phi*dq[k,i]
765
766    def limit_range(self):
767        import sys
768        from Numeric import zeros, Float
769        from util import minmod, minmod_kurganov, maxmod, vanleer, vanalbada
770        limiter = self.domain.limiter
771        #print limiter
772
773        #print 'limit_range'
774        N = self.domain.number_of_elements
775        qc = self.centroid_values
776        qv = self.vertex_values
777        C = self.domain.centroids
778        X = self.domain.vertices
779        beta_p = zeros(N,Float)
780        beta_m = zeros(N,Float)
781        beta_x = zeros(N,Float)
782
783        for k in range(N):
784
785            n0 = self.domain.neighbours[k,0]
786            n1 = self.domain.neighbours[k,1]
787
788            if ( n0 >= 0) & (n1 >= 0):
789                #SLOPE DERIVATIVE LIMIT
790                beta_p[k] = (qc[k]-qc[k-1])/(C[k]-C[k-1])
791                beta_m[k] = (qc[k+1]-qc[k])/(C[k+1]-C[k])
792                beta_x[k] = (qc[k+1]-qc[k-1])/(C[k+1]-C[k-1])
793
794        dq = zeros(qv.shape, Float)
795        for i in range(2):
796            dq[:,i] =self.domain.vertices[:,i]-self.domain.centroids
797
798        #Phi limiter
799        for k in range(N):
800            n0 = self.domain.neighbours[k,0]
801            n1 = self.domain.neighbours[k,1]
802            if n0 < 0:
803                phi = (qc[k+1] - qc[k])/(C[k+1] - C[k])
804            elif n1 < 0:
805                phi = (qc[k] - qc[k-1])/(C[k] - C[k-1])
806            #elif (self.domain.wet_nodes[k,0] == 2) & (self.domain.wet_nodes[k,1] == 2):
807            #    phi = 0.0
808            else:
809                if limiter == "minmod":
810                    phi = minmod(beta_p[k],beta_m[k])
811
812                elif limiter == "minmod_kurganov":#Change this
813                    # Also known as monotonized central difference limiter
814                    # if theta = 2.0
815                    theta = 2.0
816                    phi = minmod_kurganov(theta*beta_p[k],theta*beta_m[k],beta_x[k])
817                elif limiter == "superbee":
818                    slope1 = minmod(beta_m[k],2.0*beta_p[k])
819                    slope2 = minmod(2.0*beta_m[k],beta_p[k])
820                    phi = maxmod(slope1,slope2)
821
822                elif limiter == "vanleer":
823                    phi = vanleer(beta_p[k],beta_m[k])
824
827
828            for i in range(2):
829                qv[k,i] = qc[k] + phi*dq[k,i]
830
831    def limit_steve_slope(self):
832
833        import sys
834        from Numeric import zeros, Float
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
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 shallow_water_domain import Domain
959    from Numeric 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