source: anuga_work/development/pipeflow/old_quantity.py @ 7825

Last change on this file since 7825 was 7823, checked in by steve, 15 years ago

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