source: inundation/ga/storm_surge/pyvolution-1d/quantity.py @ 705

"""Class Quantity - Implements values at each 1d element

To create:

   Quantity(domain, vertex_values)

   domain: Associated domain structure. Required.
   
   vertex_values: N x 2 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 import Domain
        from Numeric import array, zeros, Float
        
        msg = 'First argument in Quantity.__init__ '
        msg += 'must be of class Domain (or a subclass thereof)'
        assert isinstance(domain, Domain), msg

        if vertex_values is None:
            N = domain.number_of_elements            
            self.vertex_values = zeros((N, 2), Float)
        else:    
            self.vertex_values = array(vertex_values, 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()
        
    
    def interpolate(self):
        """Compute interpolated values at 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

    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
            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
                  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




class Conserved_quantity(Quantity):
    """Class conserved quantity adds to Quantity:

    storage and method for updating, and
    methods for extrapolation from centropid 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 updates of conserved quantities by
        #flux calculations and forcing functions

        N = domain.number_of_elements
        self.explicit_update = zeros(N, Float )
        self.semi_implicit_update = zeros(N, Float )

        self.gradients = zeros(N, Float)
        self.qmax = zeros(self.centroid_values.shape, Float)
        self.qmin = zeros(self.centroid_values.shape, Float)              
                

    def update(self, timestep):
        """Update centroid values based on values stored in
        explicit_update and semi_implicit_update as well as given timestep
        """

        from Numeric import sum, equal, ones, Float
        
        N = self.centroid_values.shape[0]
        
        #Explicit updates
        self.centroid_values += timestep*self.explicit_update
            
        #Semi implicit updates
        denominator = ones(N, Float)-timestep*self.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
            self.centroid_values /= denominator


    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] = ((q0-q1)/(x0-x1)*(x2-x1) - (q2-q1)/(x2-x1)*(x0-x1))/(x2-x0)
        
        return
        
    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 extrapolate_second_order(self):
        """Extrapolate conserved quantities from centroid to
        vertices for each volume using
        second order scheme.
        """
        
        self.compute_gradients()

        G = self.gradients
        V = self.domain.vertices
        Qc = self.centroid_values
        Qv = self.vertex_values
        
        #Check each triangle
        for k in range(self.domain.number_of_elements):
            #Centroid coordinates            
            x = self.domain.centroids[k]

            #vertex coordinates
            x0, x1 = V[k,:]

            #Extrapolate
            Qv[k,0] = Qc[k] + G[k]*(x0-x)
            Qv[k,1] = Qc[k] + G[k]*(x1-x)


            

    def limit(self):        
        """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 = self.domain.number_of_elements
        beta = self.domain.beta

        qc = self.centroid_values
        qv = self.vertex_values

        #Find min and max of this and neighbour's centroid values
        qmax = self.qmax
        qmin = self.qmin

        for k in range(N):
            qmax[k] = qmin[k] = qc[k]
            
            for i in [-1,1]:
                n = k+i
                if (n >= 0) & (n <= N-1):
                    qn = qc[n] #Neighbour's centroid value

                    qmin[k] = min(qmin[k], qn)
                    qmax[k] = max(qmax[k], qn)

        #Phi limiter    
        for k in range(N):

            #Diffences between centroids and maxima/minima
            dqmax = qmax[k] - qc[k] 
            dqmin = qmin[k] - qc[k]

            #Deltas between vertex and centroid values
            dq = [0.0, 0.0]
            for i in range(2):
                dq[i] = qv[k,i] - qc[k]            

            #Find the gradient limiter (phi) across vertices
            phi = 1.0
            for i in range(2):
                r = 1.0
                if (dq[i] > 0): r = dqmax/dq[i]
                if (dq[i] < 0): r = dqmin/dq[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[i]



if __name__ == "__main__":
    from domain import Domain

    points1 = [0.0, 1.0, 2.0, 3.0]
    vertex_values = [[1.0,2.0],[4.0,5.0],[-1.0,2.0]]
    
    D1 = Domain(points1)

    Q1 = Conserved_quantity(D1, vertex_values)

    print Q1.vertex_values
    print Q1.centroid_values

    new_vertex_values = [[2.0,1.0],[3.0,4.0],[-2.0,4.0]]

    Q1.set_values(new_vertex_values)
    
    print Q1.vertex_values
    print Q1.centroid_values

    new_centroid_values = [20,30,40]
    Q1.set_values(new_centroid_values,'centroids')
            
    print Q1.vertex_values
    print Q1.centroid_values

    def fun(x):
        return x**2

    Q1.set_values(fun,'vertices')    
            
    print Q1.vertex_values
    print Q1.centroid_values

    Xc = Q1.domain.vertices
    Qc = Q1.vertex_values
    print Xc
    print Qc

    Qc[1,0] = 3

    Q1.extrapolate_second_order()
    Q1.limit()
    
    from matplotlib.matlab import *
    plot(Xc,Qc)

    show()