# source:anuga_work/development/shallow_water_1d_old/quantity_domain.py@5533

Last change on this file since 5533 was 4959, checked in by steve, 16 years ago

Copied John Jakeman's and Sudi Mungkasi's 1d shallow water code

File size: 15.3 KB
Line
1"""Class Quantity - Implements values at each triangular element
2
3To create:
4
5   Quantity(domain, vertex_values)
6
7   domain: Associated domain structure. Required.
8
9   vertex_values: N x 3 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_t2 import Domain
23        from Numeric import array, zeros, Float
24
25        if vertex_values is None:
26            N = domain.number_of_elements
27            self.vertex_values = zeros((N, 2), Float)
28        else:
29            self.vertex_values = array(vertex_values).astype(Float)
30
31            N, V = self.vertex_values.shape
32            assert V == 2,\
33                   'Two vertex values per element must be specified'
34
35
36            msg = 'Number of vertex values (%d) must be consistent with'\
37                  %N
38            msg += 'number of elements in specified domain (%d).'\
39                   %domain.number_of_elements
40
41            assert N == domain.number_of_elements, msg
42
43        self.domain = domain
44
45        #Allocate space for other quantities
46        self.centroid_values = zeros(N, Float)
47
48        #Intialise centroid values
49        self.interpolate()
50
51
52
54    def __len__(self):
55        return self.centroid_values.shape[0]
56
57
58    def __neg__(self):
59        """Negate all values in this quantity giving meaning to the
60        expression -Q where Q is an instance of class Quantity
61        """
62
63        Q = Quantity(self.domain)
64        Q.set_values(-self.vertex_values)
65        return Q
66
67
69        """Add to self anything that could populate a quantity
70
71        E.g other can be a constant, an array, a function, another quantity
72        (except for a filename or points, attributes (for now))
73        - see set_values for details
74        """
75
76        Q = Quantity(self.domain)
77        Q.set_values(other)
78
79        result = Quantity(self.domain)
80        result.set_values(self.vertex_values + Q.vertex_values)
81        return result
82
84        """Handle cases like 7+Q, where Q is an instance of class Quantity
85        """
86        return self + other
87
88
89    def __sub__(self, other):
90        return self + -other  #Invoke __neg__
91
92    def __mul__(self, other):
93        """Multiply self with anything that could populate a quantity
94
95        E.g other can be a constant, an array, a function, another quantity
96        (except for a filename or points, attributes (for now))
97        - see set_values for details
98
99        Note that if two quantitites q1 and q2 are multiplied,
100        vertex values are multiplied entry by entry
101        while centroid and edge values are re-interpolated.
102        Hence they won't be the product of centroid or edge values
103        from q1 and q2.
104        """
105
106        Q = Quantity(self.domain)
107        Q.set_values(other)
108
109        result = Quantity(self.domain)
110        result.set_values(self.vertex_values * Q.vertex_values)
111        return result
112
113    def __rmul__(self, other):
114        """Handle cases like 3*Q, where Q is an instance of class Quantity
115        """
116        return self * other
117
118    def __pow__(self, other):
119        """Raise quantity to (numerical) power
120
121        As with __mul__ vertex values are processed entry by entry
122        while centroid and edge values are re-interpolated.
123
124        Example using __pow__:
125          Q = (Q1**2 + Q2**2)**0.5
126
127        """
128
129        result = Quantity(self.domain)
130        result.set_values(self.vertex_values**other)
131        return result
132
133
134
135    def interpolate(self):
136        """Compute interpolated values at edges and centroid
137        Pre-condition: vertex_values have been set
138        """
139
140        N = self.vertex_values.shape[0]
141        for i in range(N):
142            v0 = self.vertex_values[i, 0]
143            v1 = self.vertex_values[i, 1]
144
145            self.centroid_values[i] = (v0 + v1)/2.0
146
147
148    def set_values(self, X, location='vertices'):
149        """Set values for quantity
150
151        X: Compatible list, Numeric array (see below), constant or function
152        location: Where values are to be stored.
153                  Permissible options are: vertices, centroid
154                  Default is "vertices"
155
156        In case of location == 'centroid' the dimension values must
157        be a list of a Numerical array of length N, N being the number
158        of elements in the mesh. Otherwise it must be of dimension Nx3
159
160        The values will be stored in elements following their
161        internal ordering.
162
163        If values are described a function, it will be evaluated at specified points
164
165        If selected location is vertices, values for centroid and edges
166        will be assigned interpolated values.
167        In any other case, only values for the specified locations
168        will be assigned and the others will be left undefined.
169        """
170
171        if location not in ['vertices', 'centroids']:
172            msg = 'Invalid location: %s, (possible choices vertices, centroids)' %location
173            raise msg
174
175        if X is None:
176            msg = 'Given values are None'
177            raise msg
178
179        import types
180
181        if callable(X):
182            #Use function specific method
183            self.set_function_values(X, location)
184        elif type(X) in [types.FloatType, types.IntType, types.LongType]:
185            if location == 'centroids':
186                self.centroid_values[:] = X
187            else:
188                self.vertex_values[:] = X
189
190        else:
191            #Use array specific method
192            self.set_array_values(X, location)
193
194        if location == 'vertices':
195            #Intialise centroid values
196            self.interpolate()
197
198
199
200
201    def set_function_values(self, f, location='vertices'):
202        """Set values for quantity using specified function
203
204        f: x -> z Function where x and z are arrays
205        location: Where values are to be stored.
206                  Permissible options are: vertices, centroid
207                  Default is "vertices"
208        """
209
210        if location == 'centroids':
211
212            P = self.domain.centroids
213            self.set_values(f(P), location)
214        else:
215            #Vertices
216            P = self.domain.get_vertices()
217
218            for i in range(2):
219                self.vertex_values[:,i] = f(P[:,i])
220
221    def set_array_values(self, values, location='vertices'):
222        """Set values for quantity
223
224        values: Numeric array
225        location: Where values are to be stored.
226                  Permissible options are: vertices, centroid, edges
227                  Default is "vertices"
228
229        In case of location == 'centroid' the dimension values must
230        be a list of a Numerical array of length N, N being the number
231        of elements in the mesh. Otherwise it must be of dimension Nx2
232
233        The values will be stored in elements following their
234        internal ordering.
235
236        If selected location is vertices, values for centroid
237        will be assigned interpolated values.
238        In any other case, only values for the specified locations
239        will be assigned and the others will be left undefined.
240        """
241
242        from Numeric import array, Float
243
244        values = array(values).astype(Float)
245
246        N = self.centroid_values.shape[0]
247
248        msg = 'Number of values must match number of elements'
249        assert values.shape[0] == N, msg
250
251        if location == 'centroids':
252            assert len(values.shape) == 1, 'Values array must be 1d'
253            self.centroid_values = values
254        else:
255            assert len(values.shape) == 2, 'Values array must be 2d'
256            msg = 'Array must be N x 2'
257            assert values.shape[1] == 2, msg
258
259            self.vertex_values = values
260
261
262    def extrapolate_first_order(self):
263        """Extrapolate conserved quantities from centroid to
264        vertices for each volume using
265        first order scheme.
266        """
267
268        qc = self.centroid_values
269        qv = self.vertex_values
270
271        for i in range(2):
272            qv[:,i] = qc
273
274
275    def get_integral(self):
276        """Compute the integral of quantity across entire domain
277        """
278        integral = 0
279        for k in range(self.domain.number_of_elements):
280            area = self.domain.areas[k]
281            qc = self.centroid_values[k]
282            integral += qc*area
283
284        return integral
285
286
287
288
289class Conserved_quantity(Quantity):
290    """Class conserved quantity adds to Quantity:
291
292    boundary values, storage and method for updating, and
293    methods for (second order) extrapolation from centroid to vertices inluding
295    """
296
297    def __init__(self, domain, vertex_values=None):
298        Quantity.__init__(self, domain, vertex_values)
299
300        from Numeric import zeros, Float
301
302        #Allocate space for boundary values
303        L = len(domain.boundary)
304        self.boundary_values = zeros(L, Float)
305
306        #Allocate space for updates of conserved quantities by
307        #flux calculations and forcing functions
308
309        N = domain.number_of_elements
311        self.explicit_update = zeros(N, Float )
312        self.semi_implicit_update = zeros(N, Float )
313
314
315    def update(self, timestep):
316        #Call correct module function
317        #(either from this module or C-extension)
318        return update(self, timestep)
319
320
322        #Call correct module function
323        #(either from this module or C-extension)
325
326
327    def limit(self):
328        #Call correct module function
329        #(either from this module or C-extension)
330        limit(self)
331
332
333    def extrapolate_second_order(self):
334        #Call correct module function
335        #(either from this module or C-extension)
336        extrapolate_second_order(self)
337
338
339def update(quantity, timestep):
340    """Update centroid values based on values stored in
341    explicit_update and semi_implicit_update as well as given timestep
342
343    Function implementing forcing terms must take on argument
344    which is the domain and they must update either explicit
346
347    def gravity(domain):
348        ....
349        domain.quantities['xmomentum'].explicit_update = ...
350        domain.quantities['ymomentum'].explicit_update = ...
351
352
353
354    Explicit terms must have the form
355
356        G(q, t)
357
358    and explicit scheme is
359
360       q^{(n+1}) = q^{(n)} + delta_t G(q^{n}, n delta_t)
361
362
363    Semi implicit forcing terms are assumed to have the form
364
365       G(q, t) = H(q, t) q
366
367    and the semi implicit scheme will then be
368
369      q^{(n+1}) = q^{(n)} + delta_t H(q^{n}, n delta_t) q^{(n+1})
370
371
372    """
373
374    from Numeric import sum, equal, ones, exp, Float
375
376    N = quantity.centroid_values.shape[0]
377
378
379    #Divide H by conserved quantity to obtain G (see docstring above)
380
381
382    for k in range(N):
383        x = quantity.centroid_values[k]
384        if x == 0.0:
385            #FIXME: Is this right
386            quantity.semi_implicit_update[k] = 0.0
387        else:
388            quantity.semi_implicit_update[k] /= x
389
390
392    denominator = ones(N, Float)-timestep*quantity.semi_implicit_update
393
394    if sum(equal(denominator, 0.0)) > 0.0:
395        msg = 'Zero division in semi implicit update. Call Stephen :-)'
396        raise msg
397    else:
398        #Update conserved_quantities from semi implicit updates
399        quantity.centroid_values /= denominator
400
401#    quantity.centroid_values = exp(timestep*quantity.semi_implicit_update)*quantity.centroid_values
402
403
405    quantity.centroid_values += timestep*quantity.explicit_update
406
407def extrapolate_second_order(quantity):
408    """Extrapolate conserved quantities from centroid to
409    vertices for each volume using
410    second order scheme.
411    """
412
415
416    X = quantity.domain.vertices
417    qc = quantity.centroid_values
418    qv = quantity.vertex_values
419
420    #Check each triangle
421    for k in range(quantity.domain.number_of_elements):
422        #Centroid coordinates
423        x = quantity.domain.centroids[k]
424
425        #vertex coordinates
426        x0, x1 = X[k,:]
427
428        #Extrapolate
429        qv[k,0] = qc[k] + a[k]*(x0-x)
430        qv[k,1] = qc[k] + a[k]*(x1-x)
431
432    quantity.limit()
433
435        """Compute gradients of piecewise linear function defined by centroids of
436        neighbouring volumes.
437        """
438
439
440        from Numeric import array, zeros, Float
441
442        N = self.centroid_values.shape[0]
443
444
446        Q = self.centroid_values
447        X = self.domain.centroids
448
449        for k in range(N):
450
451            # first and last elements have boundaries
452
453            if k == 0:
454
455                #Get data
456                k0 = k
457                k1 = k+1
458
459                q0 = Q[k0]
460                q1 = Q[k1]
461
462                x0 = X[k0] #V0 centroid
463                x1 = X[k1] #V1 centroid
464
466                G[k] = (q1 - q0)/(x1 - x0)
467
468            elif k == N-1:
469
470                #Get data
471                k0 = k
472                k1 = k-1
473
474                q0 = Q[k0]
475                q1 = Q[k1]
476
477                x0 = X[k0] #V0 centroid
478                x1 = X[k1] #V1 centroid
479
481                G[k] = (q1 - q0)/(x1 - x0)
482
483            else:
484                #Interior Volume (2 neighbours)
485
486                #Get data
487                k0 = k-1
488                k2 = k+1
489
490                q0 = Q[k0]
491                q1 = Q[k]
492                q2 = Q[k2]
493
494                x0 = X[k0] #V0 centroid
495                x1 = X[k]  #V1 centroid (Self)
496                x2 = X[k2] #V2 centroid
497
499                #G[k] = (q2-q0)/(x2-x0)
500                G[k] = ((q0-q1)/(x0-x1)*(x2-x1) - (q2-q1)/(x2-x1)*(x0-x1))/(x2-x0)
501
502def limit(quantity):
503    """Limit slopes for each volume to eliminate artificial variance
504    introduced by e.g. second order extrapolator
505
506    This is an unsophisticated limiter as it does not take into
507    account dependencies among quantities.
508
509    precondition:
510    vertex values are estimated from gradient
511    postcondition:
512    vertex values are updated
513    """
514
515    from Numeric import zeros, Float
516
517    N = quantity.domain.number_of_elements
518
519    beta = quantity.domain.beta
520
521    qc = quantity.centroid_values
522    qv = quantity.vertex_values
523
524    #Find min and max of this and neighbour's centroid values
525    qmax = zeros(qc.shape, Float)
526    qmin = zeros(qc.shape, Float)
527
528    for k in range(N):
529        qmax[k] = qmin[k] = qc[k]
530        for i in range(2):
531            n = quantity.domain.neighbours[k,i]
532            if n >= 0:
533                qn = qc[n] #Neighbour's centroid value
534
535                qmin[k] = min(qmin[k], qn)
536                qmax[k] = max(qmax[k], qn)
537
538    #Diffences between centroids and maxima/minima
539    dqmax = qmax - qc
540    dqmin = qmin - qc
541
542    #Deltas between vertex and centroid values
543    dq = zeros(qv.shape, Float)
544    for i in range(2):
545        dq[:,i] = qv[:,i] - qc
546
547    #Phi limiter
548    for k in range(N):
549
550        #Find the gradient limiter (phi) across vertices
551        phi = 1.0
552        for i in range(2):
553            r = 1.0
554            if (dq[k,i] > 0): r = dqmax[k]/dq[k,i]
555            if (dq[k,i] < 0): r = dqmin[k]/dq[k,i]
556
557            phi = min( min(r*beta, 1), phi )
558
559        #Then update using phi limiter
560        for i in range(2):
561            qv[k,i] = qc[k] + phi*dq[k,i]
Note: See TracBrowser for help on using the repository browser.