source: development/pyvolution-1d/quantity.py @ 2229

Last change on this file since 2229 was 2229, checked in by steve, 18 years ago

Moved directories into production and development parent directories

File size: 13.1 KB
Line 
1"""Class Quantity - Implements values at each 1d element
2
3To create:
4
5   Quantity(domain, vertex_values)
6
7   domain: Associated domain structure. Required.
8
9   vertex_values: N x 2 array of values at each vertex for each element.
10                  Default None
11
12   If vertex_values are None Create array of zeros compatible with domain.
13   Otherwise check that it is compatible with dimenions of domain.
14   Otherwise raise an exception
15"""
16
17
18class Quantity:
19
20    def __init__(self, domain, vertex_values=None):
21
22        from domain import Domain
23        from Numeric import array, zeros, Float
24
25        msg = 'First argument in Quantity.__init__ '
26        msg += 'must be of class Domain (or a subclass thereof)'
27        assert isinstance(domain, Domain), msg
28
29        if vertex_values is None:
30            N = domain.number_of_elements
31            self.vertex_values = zeros((N, 2), Float)
32        else:
33            self.vertex_values = array(vertex_values, Float)
34
35            N, V = self.vertex_values.shape
36            assert V == 2,\
37                   'Two vertex values per element must be specified'
38
39
40            msg = 'Number of vertex values (%d) must be consistent with'\
41                  %N
42            msg += 'number of elements in specified domain (%d).'\
43                   %domain.number_of_elements
44
45            assert N == domain.number_of_elements, msg
46
47        self.domain = domain
48
49        #Allocate space for other quantities
50        self.centroid_values = zeros(N, Float)
51
52        #Intialise centroid values
53        self.interpolate()
54
55
56    def interpolate(self):
57        """Compute interpolated values at centroid
58        Pre-condition: vertex_values have been set
59        """
60
61        N = self.vertex_values.shape[0]
62        for i in range(N):
63            v0 = self.vertex_values[i, 0]
64            v1 = self.vertex_values[i, 1]
65
66            self.centroid_values[i] = (v0 + v1)/2
67
68    def set_values(self, X, location='vertices'):
69        """Set values for quantity
70
71        X: Compatible list, Numeric array (see below), constant or function
72        location: Where values are to be stored.
73                  Permissible options are: vertices, centroid
74                  Default is "vertices"
75
76        In case of location == 'centroid' the dimension values must
77        be a list of a Numerical array of length N, N being the number
78        of elements in the mesh. Otherwise it must be of dimension Nx3
79
80        The values will be stored in elements following their
81        internal ordering.
82
83        If values are described a function, it will be evaluated at specified points
84
85        If selected location is vertices, values for centroid and edges
86        will be assigned interpolated values.
87        In any other case, only values for the specified locations
88        will be assigned and the others will be left undefined.
89        """
90
91        if location not in ['vertices', 'centroids']:
92            msg = 'Invalid location: %s, (possible choices vertices, centroids)' %location
93            raise msg
94
95        if X is None:
96            msg = 'Given values are None'
97            raise msg
98
99        import types
100
101        if callable(X):
102            #Use function specific method
103            self.set_function_values(X, location)
104        elif type(X) in [types.FloatType, types.IntType, types.LongType]:
105            if location == 'centroids':
106                self.centroid_values[:] = X
107            else:
108                self.vertex_values[:] = X
109
110        else:
111            #Use array specific method
112            self.set_array_values(X, location)
113
114        if location == 'vertices':
115            #Intialise centroid
116            self.interpolate()
117
118
119
120
121    def set_function_values(self, f, location='vertices'):
122        """Set values for quantity using specified function
123
124        f: x -> z Function where x and z are arrays
125        location: Where values are to be stored.
126                  Permissible options are: vertices, centroid
127                  Default is "vertices"
128        """
129
130        if location == 'centroids':
131            P = self.domain.centroids
132            self.set_values(f(P), location)
133        else:
134            #Vertices
135            P = self.domain.get_vertices()
136            for i in range(2):
137                self.vertex_values[:,i] = f(P[:,i])
138
139
140    def set_array_values(self, values, location='vertices'):
141        """Set values for quantity
142
143        values: Numeric array
144        location: Where values are to be stored.
145                  Permissible options are: vertices, centroid
146                  Default is "vertices"
147
148        In case of location == 'centroid' the dimension values must
149        be a list of a Numerical array of length N, N being the number
150        of elements in the mesh. Otherwise it must be of dimension Nx2
151
152        The values will be stored in elements following their
153        internal ordering.
154
155        If selected location is vertices, values for centroid
156        will be assigned interpolated values.
157        In any other case, only values for the specified locations
158        will be assigned and the others will be left undefined.
159        """
160
161        from Numeric import array, Float
162
163        values = array(values).astype(Float)
164
165        N = self.centroid_values.shape[0]
166
167        msg = 'Number of values must match number of elements'
168        assert values.shape[0] == N, msg
169
170        if location == 'centroids':
171            assert len(values.shape) == 1, 'Values array must be 1d'
172            self.centroid_values = values
173        else:
174            assert len(values.shape) == 2, 'Values array must be 2d'
175            msg = 'Array must be N x 2'
176            assert values.shape[1] == 2, msg
177
178            self.vertex_values = values
179
180
181
182
183class Conserved_quantity(Quantity):
184    """Class conserved quantity adds to Quantity:
185
186    storage and method for updating, and
187    methods for extrapolation from centropid to vertices inluding
188    gradients and limiters
189    """
190
191    def __init__(self, domain, vertex_values=None):
192        Quantity.__init__(self, domain, vertex_values)
193
194        from Numeric import zeros, Float
195
196        #Allocate space for updates of conserved quantities by
197        #flux calculations and forcing functions
198
199        N = domain.number_of_elements
200        self.explicit_update = zeros(N, Float )
201        self.semi_implicit_update = zeros(N, Float )
202
203        self.gradients = zeros(N, Float)
204        self.qmax = zeros(self.centroid_values.shape, Float)
205        self.qmin = zeros(self.centroid_values.shape, Float)
206
207        self.beta = domain.beta
208
209
210    def update(self, timestep):
211        """Update centroid values based on values stored in
212        explicit_update and semi_implicit_update as well as given timestep
213        """
214
215        from Numeric import sum, equal, ones, Float
216
217        N = self.centroid_values.shape[0]
218
219        #Explicit updates
220        self.centroid_values += timestep*self.explicit_update
221
222        #Semi implicit updates
223        denominator = ones(N, Float)-timestep*self.semi_implicit_update
224
225        if sum(equal(denominator, 0.0)) > 0.0:
226            msg = 'Zero division in semi implicit update. Call Stephen :-)'
227            raise msg
228        else:
229            #Update conserved_quantities from semi implicit updates
230            self.centroid_values /= denominator
231
232
233    def compute_gradients(self):
234        """Compute gradients of piecewise linear function defined by centroids of
235        neighbouring volumes.
236        """
237
238
239        from Numeric import array, zeros, Float
240
241        N = self.centroid_values.shape[0]
242
243
244        G = self.gradients
245        Q = self.centroid_values
246        X = self.domain.centroids
247
248        for k in range(N):
249
250            # first and last elements have boundaries
251
252            if k == 0:
253
254                #Get data
255                k0 = k
256                k1 = k+1
257
258                q0 = Q[k0]
259                q1 = Q[k1]
260
261                x0 = X[k0] #V0 centroid
262                x1 = X[k1] #V1 centroid
263
264                #Gradient
265                G[k] = (q1 - q0)/(x1 - x0)
266
267            elif k == N-1:
268
269                #Get data
270                k0 = k
271                k1 = k-1
272
273                q0 = Q[k0]
274                q1 = Q[k1]
275
276                x0 = X[k0] #V0 centroid
277                x1 = X[k1] #V1 centroid
278
279                #Gradient
280                G[k] = (q1 - q0)/(x1 - x0)
281
282            else:
283                #Interior Volume (2 neighbours)
284
285                #Get data
286                k0 = k-1
287                k2 = k+1
288
289                q0 = Q[k0]
290                q1 = Q[k]
291                q2 = Q[k2]
292
293                x0 = X[k0] #V0 centroid
294                x1 = X[k] #V1 centroid (Self)
295                x2 = X[k2] #V2 centroid
296
297                #Gradient
298                G[k] = ((q0-q1)/(x0-x1)*(x2-x1) - (q2-q1)/(x2-x1)*(x0-x1))/(x2-x0)
299
300        return
301
302    def extrapolate_first_order(self):
303        """Extrapolate conserved quantities from centroid to
304        vertices for each volume using
305        first order scheme.
306        """
307
308        qc = self.centroid_values
309        qv = self.vertex_values
310
311        for i in range(2):
312            qv[:,i] = qc
313
314
315    def extrapolate_second_order(self):
316        """Extrapolate conserved quantities from centroid to
317        vertices for each volume using
318        second order scheme.
319        """
320
321        self.compute_gradients()
322
323        G = self.gradients
324        V = self.domain.vertices
325        Qc = self.centroid_values
326        Qv = self.vertex_values
327
328        #Check each triangle
329        for k in range(self.domain.number_of_elements):
330            #Centroid coordinates
331            x = self.domain.centroids[k]
332
333            #vertex coordinates
334            x0, x1 = V[k,:]
335
336            #Extrapolate
337            Qv[k,0] = Qc[k] + G[k]*(x0-x)
338            Qv[k,1] = Qc[k] + G[k]*(x1-x)
339
340
341
342
343    def limit(self):
344        """Limit slopes for each volume to eliminate artificial variance
345        introduced by e.g. second order extrapolator
346
347        This is an unsophisticated limiter as it does not take into
348        account dependencies among quantities.
349
350        precondition:
351        vertex values are estimated from gradient
352        postcondition:
353        vertex values are updated
354        """
355
356        from Numeric import zeros, Float
357
358        N = self.domain.number_of_elements
359        beta = self.beta
360
361        qc = self.centroid_values
362        qv = self.vertex_values
363
364        #Find min and max of this and neighbour's centroid values
365        qmax = self.qmax
366        qmin = self.qmin
367
368        for k in range(N):
369            qmax[k] = qmin[k] = qc[k]
370
371            for i in [-1,1]:
372                n = k+i
373                if (n >= 0) & (n <= N-1):
374                    qn = qc[n] #Neighbour's centroid value
375
376                    qmin[k] = min(qmin[k], qn)
377                    qmax[k] = max(qmax[k], qn)
378
379        #Phi limiter
380        for k in range(N):
381
382            #Diffences between centroids and maxima/minima
383            dqmax = qmax[k] - qc[k]
384            dqmin = qmin[k] - qc[k]
385
386            #Deltas between vertex and centroid values
387            dq = [0.0, 0.0]
388            for i in range(2):
389                dq[i] = qv[k,i] - qc[k]
390
391            #Find the gradient limiter (phi) across vertices
392            phi = 1.0
393            for i in range(2):
394                r = 1.0
395                if (dq[i] > 0): r = dqmax/dq[i]
396                if (dq[i] < 0): r = dqmin/dq[i]
397
398                phi = min( min(r*beta, 1), phi )
399
400            #Then update using phi limiter
401            for i in range(2):
402                qv[k,i] = qc[k] + phi*dq[i]
403
404
405
406if __name__ == "__main__":
407    from domain import Domain
408
409    points1 = [0.0, 1.0, 2.0, 3.0]
410    vertex_values = [[1.0,2.0],[4.0,5.0],[-1.0,2.0]]
411
412    D1 = Domain(points1)
413
414    Q1 = Conserved_quantity(D1, vertex_values)
415
416    print Q1.vertex_values
417    print Q1.centroid_values
418
419    new_vertex_values = [[2.0,1.0],[3.0,4.0],[-2.0,4.0]]
420
421    Q1.set_values(new_vertex_values)
422
423    print Q1.vertex_values
424    print Q1.centroid_values
425
426    new_centroid_values = [20,30,40]
427    Q1.set_values(new_centroid_values,'centroids')
428
429    print Q1.vertex_values
430    print Q1.centroid_values
431
432    def fun(x):
433        return x**2
434
435    Q1.set_values(fun,'vertices')
436
437    print Q1.vertex_values
438    print Q1.centroid_values
439
440    Xc = Q1.domain.vertices
441    Qc = Q1.vertex_values
442    print Xc
443    print Qc
444
445    Qc[1,0] = 3
446
447    Q1.beta = 1.0
448    Q1.extrapolate_second_order()
449    Q1.limit()
450
451    from pylab import *
452    #plot(Xc,Qc)
453
454    #show()
455
456    points2 = arange(10)
457
458
459    D2 = Domain(points2)
460
461    Q2 = Conserved_quantity(D2)
462    Q2.set_values(fun,'vertices')
463    Xc = Q2.domain.vertices
464    Qc = Q2.vertex_values
465    plot(Xc,Qc)
466
467    show()
468    Q2.extrapolate_second_order()
469    Q2.limit()
470    Xc = Q2.domain.vertices
471    Qc = Q2.vertex_values
472
473    print Q2.centroid_values
474    print Qc
475
476    plot(Xc,Qc)
477
478    show()
479
480
481
482
483
484
485
486
487
Note: See TracBrowser for help on using the repository browser.