source: trunk/anuga_work/development/2010-projects/anuga_1d/avalanche-sudi/quantity.py @ 8427

Last change on this file since 8427 was 8427, checked in by davies, 13 years ago

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