"""Class Quantity - Implements values at each triangular element To create: Quantity(domain, vertex_values) domain: Associated domain structure. Required. vertex_values: N x 3 array of values at each vertex for each element. Default None If vertex_values are None Create array of zeros compatible with domain. Otherwise check that it is compatible with dimenions of domain. Otherwise raise an exception """ class Quantity: def __init__(self, domain, vertex_values=None): from domain_t2 import Domain from Numeric import array, zeros, Float if vertex_values is None: N = domain.number_of_elements self.vertex_values = zeros((N, 2), Float) else: self.vertex_values = array(vertex_values).astype(Float) N, V = self.vertex_values.shape assert V == 2,\ 'Two vertex values per element must be specified' msg = 'Number of vertex values (%d) must be consistent with'\ %N msg += 'number of elements in specified domain (%d).'\ %domain.number_of_elements assert N == domain.number_of_elements, msg self.domain = domain #Allocate space for other quantities self.centroid_values = zeros(N, Float) #Intialise centroid values self.interpolate() #Methods for operator overloading def __len__(self): return self.centroid_values.shape[0] def __neg__(self): """Negate all values in this quantity giving meaning to the expression -Q where Q is an instance of class Quantity """ Q = Quantity(self.domain) Q.set_values(-self.vertex_values) return Q def __add__(self, other): """Add to self anything that could populate a quantity E.g other can be a constant, an array, a function, another quantity (except for a filename or points, attributes (for now)) - see set_values for details """ Q = Quantity(self.domain) Q.set_values(other) result = Quantity(self.domain) result.set_values(self.vertex_values + Q.vertex_values) return result def __radd__(self, other): """Handle cases like 7+Q, where Q is an instance of class Quantity """ return self + other def __sub__(self, other): return self + -other #Invoke __neg__ def __mul__(self, other): """Multiply self with anything that could populate a quantity E.g other can be a constant, an array, a function, another quantity (except for a filename or points, attributes (for now)) - see set_values for details Note that if two quantitites q1 and q2 are multiplied, vertex values are multiplied entry by entry while centroid and edge values are re-interpolated. Hence they won't be the product of centroid or edge values from q1 and q2. """ Q = Quantity(self.domain) Q.set_values(other) result = Quantity(self.domain) result.set_values(self.vertex_values * Q.vertex_values) return result def __rmul__(self, other): """Handle cases like 3*Q, where Q is an instance of class Quantity """ return self * other def __pow__(self, other): """Raise quantity to (numerical) power As with __mul__ vertex values are processed entry by entry while centroid and edge values are re-interpolated. Example using __pow__: Q = (Q1**2 + Q2**2)**0.5 """ result = Quantity(self.domain) result.set_values(self.vertex_values**other) return result def interpolate(self): """Compute interpolated values at edges and centroid Pre-condition: vertex_values have been set """ N = self.vertex_values.shape[0] for i in range(N): v0 = self.vertex_values[i, 0] v1 = self.vertex_values[i, 1] self.centroid_values[i] = (v0 + v1)/2.0 def set_values(self, X, location='vertices'): """Set values for quantity X: Compatible list, Numeric array (see below), constant or function location: Where values are to be stored. Permissible options are: vertices, centroid Default is "vertices" In case of location == 'centroid' the dimension values must be a list of a Numerical array of length N, N being the number of elements in the mesh. Otherwise it must be of dimension Nx3 The values will be stored in elements following their internal ordering. If values are described a function, it will be evaluated at specified points If selected location is vertices, values for centroid and edges will be assigned interpolated values. In any other case, only values for the specified locations will be assigned and the others will be left undefined. """ if location not in ['vertices', 'centroids']: msg = 'Invalid location: %s, (possible choices vertices, centroids)' %location raise msg if X is None: msg = 'Given values are None' raise msg import types if callable(X): #Use function specific method self.set_function_values(X, location) elif type(X) in [types.FloatType, types.IntType, types.LongType]: if location == 'centroids': self.centroid_values[:] = X else: self.vertex_values[:] = X else: #Use array specific method self.set_array_values(X, location) if location == 'vertices': #Intialise centroid values self.interpolate() def set_function_values(self, f, location='vertices'): """Set values for quantity using specified function f: x -> z Function where x and z are arrays location: Where values are to be stored. Permissible options are: vertices, centroid Default is "vertices" """ if location == 'centroids': P = self.domain.centroids self.set_values(f(P), location) else: #Vertices P = self.domain.get_vertices() for i in range(2): self.vertex_values[:,i] = f(P[:,i]) def set_array_values(self, values, location='vertices'): """Set values for quantity values: Numeric array location: Where values are to be stored. Permissible options are: vertices, centroid, edges Default is "vertices" In case of location == 'centroid' the dimension values must be a list of a Numerical array of length N, N being the number of elements in the mesh. Otherwise it must be of dimension Nx2 The values will be stored in elements following their internal ordering. If selected location is vertices, values for centroid will be assigned interpolated values. In any other case, only values for the specified locations will be assigned and the others will be left undefined. """ from Numeric import array, Float values = array(values).astype(Float) N = self.centroid_values.shape[0] msg = 'Number of values must match number of elements' assert values.shape[0] == N, msg if location == 'centroids': assert len(values.shape) == 1, 'Values array must be 1d' self.centroid_values = values else: assert len(values.shape) == 2, 'Values array must be 2d' msg = 'Array must be N x 2' assert values.shape[1] == 2, msg self.vertex_values = values def extrapolate_first_order(self): """Extrapolate conserved quantities from centroid to vertices for each volume using first order scheme. """ qc = self.centroid_values qv = self.vertex_values for i in range(2): qv[:,i] = qc def get_integral(self): """Compute the integral of quantity across entire domain """ integral = 0 for k in range(self.domain.number_of_elements): area = self.domain.areas[k] qc = self.centroid_values[k] integral += qc*area return integral class Conserved_quantity(Quantity): """Class conserved quantity adds to Quantity: boundary values, storage and method for updating, and methods for (second order) extrapolation from centroid to vertices inluding gradients and limiters """ def __init__(self, domain, vertex_values=None): Quantity.__init__(self, domain, vertex_values) from Numeric import zeros, Float #Allocate space for boundary values L = len(domain.boundary) self.boundary_values = zeros(L, Float) #Allocate space for updates of conserved quantities by #flux calculations and forcing functions N = domain.number_of_elements self.gradients = zeros(N, Float) self.explicit_update = zeros(N, Float ) self.semi_implicit_update = zeros(N, Float ) def update(self, timestep): #Call correct module function #(either from this module or C-extension) return update(self, timestep) def compute_gradients(self): #Call correct module function #(either from this module or C-extension) return compute_gradients(self) def limit(self): #Call correct module function #(either from this module or C-extension) limit(self) def extrapolate_second_order(self): #Call correct module function #(either from this module or C-extension) extrapolate_second_order(self) def update(quantity, timestep): """Update centroid values based on values stored in explicit_update and semi_implicit_update as well as given timestep Function implementing forcing terms must take on argument which is the domain and they must update either explicit or implicit updates, e,g,: def gravity(domain): .... domain.quantities['xmomentum'].explicit_update = ... domain.quantities['ymomentum'].explicit_update = ... Explicit terms must have the form G(q, t) and explicit scheme is q^{(n+1}) = q^{(n)} + delta_t G(q^{n}, n delta_t) Semi implicit forcing terms are assumed to have the form G(q, t) = H(q, t) q and the semi implicit scheme will then be q^{(n+1}) = q^{(n)} + delta_t H(q^{n}, n delta_t) q^{(n+1}) """ from Numeric import sum, equal, ones, exp, Float N = quantity.centroid_values.shape[0] #Divide H by conserved quantity to obtain G (see docstring above) for k in range(N): x = quantity.centroid_values[k] if x == 0.0: #FIXME: Is this right quantity.semi_implicit_update[k] = 0.0 else: quantity.semi_implicit_update[k] /= x #Semi implicit updates denominator = ones(N, Float)-timestep*quantity.semi_implicit_update if sum(equal(denominator, 0.0)) > 0.0: msg = 'Zero division in semi implicit update. Call Stephen :-)' raise msg else: #Update conserved_quantities from semi implicit updates quantity.centroid_values /= denominator # quantity.centroid_values = exp(timestep*quantity.semi_implicit_update)*quantity.centroid_values #Explicit updates quantity.centroid_values += timestep*quantity.explicit_update def extrapolate_second_order(quantity): """Extrapolate conserved quantities from centroid to vertices for each volume using second order scheme. """ quantity.compute_gradients() a = quantity.gradients X = quantity.domain.vertices qc = quantity.centroid_values qv = quantity.vertex_values #Check each triangle for k in range(quantity.domain.number_of_elements): #Centroid coordinates x = quantity.domain.centroids[k] #vertex coordinates x0, x1 = X[k,:] #Extrapolate qv[k,0] = qc[k] + a[k]*(x0-x) qv[k,1] = qc[k] + a[k]*(x1-x) quantity.limit() def compute_gradients(self): """Compute gradients of piecewise linear function defined by centroids of neighbouring volumes. """ from Numeric import array, zeros, Float N = self.centroid_values.shape[0] G = self.gradients Q = self.centroid_values X = self.domain.centroids for k in range(N): # first and last elements have boundaries if k == 0: #Get data k0 = k k1 = k+1 q0 = Q[k0] q1 = Q[k1] x0 = X[k0] #V0 centroid x1 = X[k1] #V1 centroid #Gradient G[k] = (q1 - q0)/(x1 - x0) elif k == N-1: #Get data k0 = k k1 = k-1 q0 = Q[k0] q1 = Q[k1] x0 = X[k0] #V0 centroid x1 = X[k1] #V1 centroid #Gradient G[k] = (q1 - q0)/(x1 - x0) else: #Interior Volume (2 neighbours) #Get data k0 = k-1 k2 = k+1 q0 = Q[k0] q1 = Q[k] q2 = Q[k2] x0 = X[k0] #V0 centroid x1 = X[k] #V1 centroid (Self) x2 = X[k2] #V2 centroid #Gradient #G[k] = (q2-q0)/(x2-x0) G[k] = ((q0-q1)/(x0-x1)*(x2-x1) - (q2-q1)/(x2-x1)*(x0-x1))/(x2-x0) def limit(quantity): """Limit slopes for each volume to eliminate artificial variance introduced by e.g. second order extrapolator This is an unsophisticated limiter as it does not take into account dependencies among quantities. precondition: vertex values are estimated from gradient postcondition: vertex values are updated """ from Numeric import zeros, Float N = quantity.domain.number_of_elements beta = quantity.domain.beta qc = quantity.centroid_values qv = quantity.vertex_values #Find min and max of this and neighbour's centroid values qmax = zeros(qc.shape, Float) qmin = zeros(qc.shape, Float) for k in range(N): qmax[k] = qmin[k] = qc[k] for i in range(2): n = quantity.domain.neighbours[k,i] if n >= 0: qn = qc[n] #Neighbour's centroid value qmin[k] = min(qmin[k], qn) qmax[k] = max(qmax[k], qn) #Diffences between centroids and maxima/minima dqmax = qmax - qc dqmin = qmin - qc #Deltas between vertex and centroid values dq = zeros(qv.shape, Float) for i in range(2): dq[:,i] = qv[:,i] - qc #Phi limiter for k in range(N): #Find the gradient limiter (phi) across vertices phi = 1.0 for i in range(2): r = 1.0 if (dq[k,i] > 0): r = dqmax[k]/dq[k,i] if (dq[k,i] < 0): r = dqmin[k]/dq[k,i] phi = min( min(r*beta, 1), phi ) #Then update using phi limiter for i in range(2): qv[k,i] = qc[k] + phi*dq[k,i]