Making code faster w/o using global variables

[amanabt]

Commit: 89296e1
PR #47

Issue:
Calculation of value of a constant variable for every time loop decreases execution speed.
One instance is the function wave_equation.surface_term. Function is shown below

def surface_term(u_n):
    '''

    Calculates the surface term,
    :math:`L_p(1) f_i - L_p(-1) f_{i - 1}`
    using the lax_friedrichs_flux function and lagrange_basis_value
    from params module.
    
    Parameters
    ----------
    u_n : arrayfire.Array [N_LGL N_Elements 1 1]
          Amplitude of the wave at the mapped LGL nodes of each element.
          
    Returns
    -------
    surface_term : arrayfire.Array [N_LGL N_Elements 1 1]
                   The surface term represented in the form of an array,
                   :math:`L_p (1) f_i - L_p (-1) f_{i - 1}`, where p varies
                   from zero to :math:`N_{LGL}` and i from zero to
                   :math:`N_{Elements}`. p varies along the rows and i along
                   columns.
    
    **See:** `PDF`_ describing the algorithm to obtain the surface term.
    
    .. _PDF: https://goo.gl/Nhhgzx

    '''

    L_p_minus1   = params.lagrange_basis_value[:, 0]
    L_p_1        = params.lagrange_basis_value[:, -1]
    f_i          = lax_friedrichs_flux(u_n)
    f_iminus1    = af.shift(f_i, 0, 1)
    surface_term = af.blas.matmul(L_p_1, f_i) - af.blas.matmul(L_p_minus1,\
                                                                    f_iminus1)
    
    return surface_term

On replacing the params.lagrange_basis_value with lagrange_basis_value_local, where

lagrange_coeffs_local = af.np_to_af_array(\
                                lagrange.lagrange_polynomials( \
                                    lagrange.LGL_points(N_LGL))[1])

lagrange_basis_value_local = lagrange.lagrange_function_value(lagrange_coeffs_local)

This change reduces the time iteration rate from 280iter/s to 50iter/s.

Similarly, instead of using params.volume_integrand_N_LGL in the wave_equation.volume_integral_flux, if we calculate it's value in the time loop, it reduces the speed from 280iter/s to 160iter/s.
Below is the code for wave_equation.volume_integral_flux

def volume_integral_flux(u_n, t_n):
    '''

    Calculates the volume integral of flux in the wave equation.

    :math:`\\int_{-1}^1 f(u) \\frac{d L_p}{d\\xi} d\\xi`

    This will give N values of flux integral as p varies from 0 to N - 1.
    
    This integral is carried out using the analytical form of the integrand
    obtained as a linear combination of Lagrange basis polynomials.

    This integrand is the used in the Integrate() function.

    Calculation of volume integral flux using N_LGL Lobatto quadrature points
    can be vectorized and is much faster.
    
    Parameters
    ----------

    u : arrayfire.Array [N_LGL N_Elements 1 1]
        Amplitude of the wave at the mapped LGL nodes of each element.
            
    Returns
    -------

    flux_integral : arrayfire.Array [N_LGL N_Elements 1 1]
                    Value of the volume integral flux. It contains the integral
                    of all N_LGL * N_Element integrands.

    '''
    # The coefficients of dLp / d\xi
    diff_lag_coeff  = params.volume_integrand_N_LGL

    lobatto_nodes   = params.lobatto_quadrature_nodes
    Lobatto_weights = params.lobatto_weights_quadrature

    nodes_tile   = af.transpose(af.tile(lobatto_nodes, 1, diff_lag_coeff.shape[1]))
    power        = af.flip(af.range(diff_lag_coeff.shape[1]))
    power_tile   = af.tile(power, 1, params.N_quad)
    nodes_power  = nodes_tile ** power_tile
    weights_tile = af.transpose(af.tile(Lobatto_weights, 1, diff_lag_coeff.shape[1]))
    nodes_weight = nodes_power * weights_tile

    dLp_dxi      = af.matmul(diff_lag_coeff, nodes_weight)

Reduction in iteration speed for other global variables in the time loop was not noticeable.

Possible Implementation:
Write a container class for storing constant variables. The class object will be initialized before the time loop begins, and passed to the functions that need it. Below is the code for the class

class constant_container
{
    def __init__(self):
        '''
        '''
        self.lagrange_coeffs = some_value
        self.lagrange_basis_value = some_value
        #.......
        # Other variables which remains constant throughout the program
        # ......
        return
}
# Now the functions which need these constant variables are redefined to accept
# an object of the class `constant_containor` as an argument.