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

Last change on this file since 5565 was 5536, checked in by steve, 15 years ago

Added unit test files for quantity and shallow water

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