source: anuga_work/development/shallow_water_1d_old/quantity_adjust.py @ 5533

Last change on this file since 5533 was 4959, checked in by steve, 16 years ago

Copied John Jakeman's and Sudi Mungkasi's 1d shallow water code

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