source: anuga_work/development/shallow_water_1d_old/quantity.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: 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_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.0
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        for k in range(self.domain.number_of_elements):
377            area = self.domain.areas[k]
378            qc = self.centroid_values[k]
379            integral += qc*area
380
381        return integral
382
383class Conserved_quantity(Quantity):
384    """Class conserved quantity adds to Quantity:
385
386    storage and method for updating, and
387    methods for extrapolation from centropid to vertices inluding
388    gradients and limiters
389    """
390
391    def __init__(self, domain, vertex_values=None):
392        Quantity.__init__(self, domain, vertex_values)
393
394        from Numeric import zeros, Float
395       
396        #Allocate space for boundary values
397        #L = len(domain.boundary)
398        self.boundary_values = zeros(2, Float) #assumes no parrellism
399
400        #Allocate space for updates of conserved quantities by
401        #flux calculations and forcing functions
402       
403        N = domain.number_of_elements
404        self.explicit_update = zeros(N, Float )
405        self.semi_implicit_update = zeros(N, Float )
406
407        self.gradients = zeros(N, Float)
408        self.qmax = zeros(self.centroid_values.shape, Float)
409        self.qmin = zeros(self.centroid_values.shape, Float)
410
411        self.beta = domain.beta
412
413    def update(self, timestep):
414        """Update centroid values based on values stored in
415        explicit_update and semi_implicit_update as well as given timestep
416        """
417
418        from Numeric import sum, equal, ones, Float
419
420        N = self.centroid_values.shape[0]
421
422        #print "pre update",self.centroid_values
423        #print "timestep",timestep
424        #print min(self.domain.areas)
425        #print "explicit_update", self.explicit_update
426
427        #Explicit updates
428        self.centroid_values += timestep*self.explicit_update
429
430        #print "post update",self.centroid_values
431       
432        #Semi implicit updates
433        denominator = ones(N, Float)-timestep*self.semi_implicit_update
434
435        if sum(equal(denominator, 0.0)) > 0.0:
436            msg = 'Zero division in semi implicit update. Call Stephen :-)'
437            raise msg
438        else:
439            #Update conserved_quantities from semi implicit updates
440            self.centroid_values /= denominator
441
442
443    def compute_gradients(self):
444        """Compute gradients of piecewise linear function defined by centroids of
445        neighbouring volumes.
446        """
447
448
449        from Numeric import array, zeros, Float
450
451        N = self.centroid_values.shape[0]
452
453
454        G = self.gradients
455        Q = self.centroid_values
456        X = self.domain.centroids
457
458        for k in range(N):
459
460            # first and last elements have boundaries
461
462            if k == 0:
463
464                #Get data
465                k0 = k
466                k1 = k+1
467
468                q0 = Q[k0]
469                q1 = Q[k1]
470
471                x0 = X[k0] #V0 centroid
472                x1 = X[k1] #V1 centroid
473
474                #Gradient
475                G[k] = (q1 - q0)/(x1 - x0)
476
477            elif k == N-1:
478
479                #Get data
480                k0 = k
481                k1 = k-1
482
483                q0 = Q[k0]
484                q1 = Q[k1]
485
486                x0 = X[k0] #V0 centroid
487                x1 = X[k1] #V1 centroid
488
489                #Gradient
490                G[k] = (q1 - q0)/(x1 - x0)
491
492            else:
493                #Interior Volume (2 neighbours)
494
495                #Get data
496                k0 = k-1
497                k2 = k+1
498
499                q0 = Q[k0]
500                q1 = Q[k]
501                q2 = Q[k2]
502
503                x0 = X[k0] #V0 centroid
504                x1 = X[k]  #V1 centroid (Self)
505                x2 = X[k2] #V2 centroid
506
507                #Gradient
508                #G[k] = (q2-q0)/(x2-x0)
509                G[k] = ((q0-q1)/(x0-x1)*(x2-x1) - (q2-q1)/(x2-x1)*(x0-x1))/(x2-x0)
510
511
512    def compute_minmod_gradients(self):
513        """Compute gradients of piecewise linear function defined by centroids of
514        neighbouring volumes.
515        """
516
517        from Numeric import array, zeros, Float,sign
518       
519        def xmin(a,b):
520            return 0.5*(sign(a)+sign(b))*min(abs(a),abs(b))
521
522        def xmic(t,a,b):
523            return xmin(t*xmin(a,b), 0.50*(a+b) )
524
525
526
527        N = self.centroid_values.shape[0]
528
529
530        G = self.gradients
531        Q = self.centroid_values
532        X = self.domain.centroids
533
534        for k in range(N):
535
536            # first and last elements have boundaries
537
538            if k == 0:
539
540                #Get data
541                k0 = k
542                k1 = k+1
543
544                q0 = Q[k0]
545                q1 = Q[k1]
546
547                x0 = X[k0] #V0 centroid
548                x1 = X[k1] #V1 centroid
549
550                #Gradient
551                G[k] = (q1 - q0)/(x1 - x0)
552
553            elif k == N-1:
554
555                #Get data
556                k0 = k
557                k1 = k-1
558
559                q0 = Q[k0]
560                q1 = Q[k1]
561
562                x0 = X[k0] #V0 centroid
563                x1 = X[k1] #V1 centroid
564
565                #Gradient
566                G[k] = (q1 - q0)/(x1 - x0)
567
568            elif (self.domain.wet_nodes[k,0] == 2) & (self.domain.wet_nodes[k,1] == 2):
569                G[k] = 0.0
570
571            else:
572                #Interior Volume (2 neighbours)
573
574                #Get data
575                k0 = k-1
576                k2 = k+1
577
578                q0 = Q[k0]
579                q1 = Q[k]
580                q2 = Q[k2]
581
582                x0 = X[k0] #V0 centroid
583                x1 = X[k]  #V1 centroid (Self)
584                x2 = X[k2] #V2 centroid
585
586                # assuming uniform grid
587                d1 = (q1 - q0)/(x1-x0)
588                d2 = (q2 - q1)/(x2-x1)
589
590                #Gradient
591                #G[k] = (d1+d2)*0.5
592                #G[k] = (d1*(x2-x1) - d2*(x0-x1))/(x2-x0)               
593                G[k] = xmic( self.domain.beta, d1, d2 )
594       
595
596    def extrapolate_first_order(self):
597        """Extrapolate conserved quantities from centroid to
598        vertices for each volume using
599        first order scheme.
600        """
601
602        qc = self.centroid_values
603        qv = self.vertex_values
604
605        for i in range(2):
606            qv[:,i] = qc
607
608
609    def extrapolate_second_order(self):
610        """Extrapolate conserved quantities from centroid to
611        vertices for each volume using
612        second order scheme.
613        """
614        if self.domain.limiter == "pyvolution":
615            #Z = self.gradients
616            #print "gradients 1",Z
617            self.compute_gradients()
618            #print "gradients 2",Z
619
620            #Z = self.gradients
621            #print "gradients 1",Z
622            #self.compute_minmod_gradients()
623            #print "gradients 2", Z
624
625            G = self.gradients
626            V = self.domain.vertices
627            qc = self.centroid_values
628            qv = self.vertex_values       
629
630            #Check each triangle
631            for k in range(self.domain.number_of_elements):
632                #Centroid coordinates
633                x = self.domain.centroids[k]
634
635                #vertex coordinates
636                x0, x1 = V[k,:]
637
638                #Extrapolate
639                qv[k,0] = qc[k] + G[k]*(x0-x)
640                qv[k,1] = qc[k] + G[k]*(x1-x)
641            self.limit_pyvolution()
642        elif self.domain.limiter == "steve_minmod":
643            self.limit_minmod()
644        else:
645            self.limit_range()
646       
647       
648
649    def limit_minmod(self):
650        #Z = self.gradients
651        #print "gradients 1",Z
652        self.compute_minmod_gradients()
653        #print "gradients 2", Z
654
655        G = self.gradients
656        V = self.domain.vertices
657        qc = self.centroid_values
658        qv = self.vertex_values       
659       
660        #Check each triangle
661        for k in range(self.domain.number_of_elements):
662            #Centroid coordinates
663            x = self.domain.centroids[k]
664
665            #vertex coordinates
666            x0, x1 = V[k,:]
667
668            #Extrapolate
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            else:
777                if limiter == "minmod":
778                    phi = minmod(beta_p[k],beta_m[k])
779
780                elif limiter == "minmod_kurganov":#Change this
781                    # Also known as monotonized central difference limiter
782                    # if theta = 2.0
783                    theta = 2.0 
784                    phi = minmod_kurganov(theta*beta_p[k],theta*beta_m[k],beta_x[k])
785                elif limiter == "superbee":
786                    slope1 = minmod(beta_m[k],2.0*beta_p[k])
787                    slope2 = minmod(2.0*beta_m[k],beta_p[k])
788                    phi = maxmod(slope1,slope2)
789
790                elif limiter == "vanleer":
791                    phi = vanleer(beta_p[k],beta_m[k])
792
793                elif limiter == "vanalbada":
794                    phi = vanalbada(beta_m[k],beta_p[k])
795           
796            for i in range(2):
797                qv[k,i] = qc[k] + phi*dq[k,i]
798
799    def limit_steve_slope(self):   
800
801        import sys
802        from Numeric import zeros, Float
803        from util import minmod, minmod_kurganov, maxmod, vanleer
804
805        N = self.domain.number_of_elements
806        limiter = self.domain.limiter
807        limiter_type = self.domain.limiter_type
808           
809        qc = self.centroid_values
810        qv = self.vertex_values
811
812        #Find min and max of this and neighbour's centroid values
813        beta_p = zeros(N,Float)
814        beta_m = zeros(N,Float)
815        beta_x = zeros(N,Float)
816        C = self.domain.centroids
817        X = self.domain.vertices
818
819        for k in range(N):
820       
821            n0 = self.domain.neighbours[k,0]
822            n1 = self.domain.neighbours[k,1]
823           
824            if (n0 >= 0) & (n1 >= 0):
825                # Check denominator not zero
826                if (qc[k+1]-qc[k]) == 0.0:
827                    beta_p[k] = float(sys.maxint)
828                    beta_m[k] = float(sys.maxint)
829                else:
830                    #STEVE LIMIT
831                    beta_p[k] = (qc[k]-qc[k-1])/(qc[k+1]-qc[k])
832                    beta_m[k] = (qc[k+2]-qc[k+1])/(qc[k+1]-qc[k])
833
834        #Deltas between vertex and centroid values
835        dq = zeros(qv.shape, Float)
836        for i in range(2):
837            dq[:,i] =self.domain.vertices[:,i]-self.domain.centroids
838           
839        #Phi limiter
840        for k in range(N):
841               
842            phi = 0.0
843            if limiter == "flux_minmod":
844                #FLUX MINMOD
845                phi = minmod_kurganov(1.0,beta_m[k],beta_p[k])
846            elif limiter == "flux_superbee":
847                #FLUX SUPERBEE
848                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
849            elif limiter == "flux_muscl":
850                #FLUX MUSCL
851                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])))
852            elif limiter == "flux_vanleer":
853                #FLUX VAN LEER
854                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
855               
856                #Then update using phi limiter
857                n = self.domain.neighbours[k,1]
858                if n>=0:
859                    #qv[k,0] = qc[k] - 0.5*phi*(qc[k+1]-qc[k])
860                    #qv[k,1] = qc[k] + 0.5*phi*(qc[k+1]-qc[k])
861                    qv[k,0] = qc[k] + 0.5*phi*(qv[k,0]-qc[k])
862                    qv[k,1] = qc[k] + 0.5*phi*(qv[k,1]-qc[k])
863                else:
864                    qv[k,i] = qc[k]
865                   
866
867def newLinePlot(title='Simple Plot'):
868    import Gnuplot
869    g = Gnuplot.Gnuplot()
870    g.title(title)
871    g('set data style linespoints') 
872    g.xlabel('x')
873    g.ylabel('y')
874    return g
875
876def linePlot(g,x,y):
877    import Gnuplot
878    g.plot(Gnuplot.PlotItems.Data(x.flat,y.flat))
879
880
881
882
883
884if __name__ == "__main__":
885    #from domain import Domain
886    from shallow_water_1d import Domain     
887    from Numeric import arange
888   
889    points1 = [0.0, 1.0, 2.0, 3.0]
890    vertex_values = [[1.0,2.0],[4.0,5.0],[-1.0,2.0]]
891
892    D1 = Domain(points1)
893
894    Q1 = Conserved_quantity(D1, vertex_values)
895
896    print Q1.vertex_values
897    print Q1.centroid_values
898
899    new_vertex_values = [[2.0,1.0],[3.0,4.0],[-2.0,4.0]]
900
901    Q1.set_values(new_vertex_values)
902
903    print Q1.vertex_values
904    print Q1.centroid_values
905
906    new_centroid_values = [20,30,40]
907    Q1.set_values(new_centroid_values,'centroids')
908
909    print Q1.vertex_values
910    print Q1.centroid_values
911
912    class FunClass:
913        def __init__(self,value):
914            self.value = value
915
916        def __call__(self,x):
917            return self.value*(x**2)
918
919
920    fun = FunClass(1.0)
921    Q1.set_values(fun,'vertices')
922
923    print Q1.vertex_values
924    print Q1.centroid_values
925
926    Xc = Q1.domain.vertices
927    Qc = Q1.vertex_values
928    print Xc
929    print Qc
930
931    Qc[1,0] = 3
932
933    Q1.beta = 1.0
934    Q1.extrapolate_second_order()
935    Q1.limit_minmod()
936
937    g1 = newLinePlot('plot 1')
938    linePlot(g1,Xc,Qc)
939
940    points2 = arange(10)
941    D2 = Domain(points2)
942
943    Q2 = Conserved_quantity(D2)
944    Q2.set_values(fun,'vertices')
945    Xc = Q2.domain.vertices
946    Qc = Q2.vertex_values
947
948    g2 = newLinePlot('plot 2')
949    linePlot(g2,Xc,Qc)
950
951
952   
953    Q2.extrapolate_second_order()
954    Q2.limit_minmod()
955    Xc = Q2.domain.vertices
956    Qc = Q2.vertex_values
957
958    print Q2.centroid_values
959    print Qc
960    raw_input('press_return')
961   
962    g3 = newLinePlot('plot 3')
963    linePlot(g3,Xc,Qc)
964    raw_input('press return')
965
966
967    for i in range(10):
968        fun = FunClass(i/10.0)
969        Q2.set_values(fun,'vertices')
970        Qc = Q2.vertex_values
971        linePlot(g3,Xc,Qc)
972
973    raw_input('press return')
974
Note: See TracBrowser for help on using the repository browser.