source: anuga_work/development/sudi/sw_1d/avalanche/B_momentum/quantity.py @ 7837

Last change on this file since 7837 was 7837, checked in by mungkasi, 14 years ago

Again, adding some codes for 1d problems on debris avalanche and periodic waves.

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