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

Last change on this file since 5827 was 5827, checked in by sudi, 15 years ago

Adding new versions of shallow_water_domain

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