Merge pull request #2 from Balavarun5/gaussian_quadrature_experiment

Gaussian quadrature experiment

.gitignore

@@ -99,5 +99,7 @@ ENV/
 
 # mypy
 .mypy_cache/
+
+# Kdev
 *kdev*
 *.ipynb

code/app/global_variables.py

@@ -1,7 +1,9 @@
 import numpy as np
 import arrayfire as af
-
+from scipy import special as sp
 from app import lagrange
+from app import wave_equation
+from utils import utils
 
 
 LGL_list = [ \
@@ -41,31 +43,213 @@
 
 
 for idx in np.arange(len(LGL_list)):
-	LGL_list[idx] = af.arith.cast(af.Array(LGL_list[idx]), af.Dtype.f64)
-	pass
+    LGL_list[idx] = np.array(LGL_list[idx], dtype = np.float64)
+    LGL_list[idx] = af.interop.np_to_af_array(LGL_list[idx])
+
+x_nodes         = af.interop.np_to_af_array(np.array([-1., 1.]))
+N_LGL           = 16
+xi_LGL          = None
+lBasisArray     = None
+lobatto_weights = None
+N_Elements      = None
+element_nodes   = None
+u               = None
+time            = None
+c               = None
+dLp_xi          = None
+
+def populateGlobalVariables(Number_of_LGL_pts = 8, Number_of_elements = 10):
+    '''
+    Time evolution of the wave function requires the use of many variables which
+    are time independent throughout the program. These variables have been
+    declared as global variables in the program. The global variables being used
+    in the program are listed below.
+    
+    x_nodes         : af.Array [2 1 1 1]
+                      Stores the domain of the wave. Default value is set to
+                      [-1, 1]. This means that the domain is from
+                      :math:`x \\epsilon [-1., 1.]`. This domain will be split
+                      into `Number_of_elements` elements.
+    
+    N_Elements      : int
+                      The number of elements into which the domain is to be
+                      divided.
+                 
+    N_LGL           : int
+                      The number of `LGL` nodes in an element.
+              
+    xi_LGL          : af.Array [N_LGL 1 1 1]
+                      `N_LGL` LGL points.
+    
+    element_LGL     : af.Array [N_LGL N_Elements 1 1]
+                      The domain is divided into `N_elements` number of elements
+                      and each element has N_LGL points corresponding to the
+                      each LGL point. This variable stores the :math:`x`
+                      coordinates corresponding to each `xi_LGL` for each
+                      element. element_LGL[n_LGL, n_element] returns the
+                      :math:`n\\_LGL^{th}` :math:`x` coordinate for the
+                      :math:`n_element^{th}` element.
+                  
+    lBasisArray     : af.Array [N_LGL N_LGL 1 1]
+                      Contains the coefficients of the Lagrange basis functions
+                      created using LGL nodes stored in xi_LGL.
+                      lBasisArray[i, j] is the coefficient corresponding to the
+                      :math:`(N_LGL - j - 1)^{th}` power of :math:`x` in the
+                      :mat:`i^th` Lagrange polynomial.
+    
+    dLp_xi          : af.Array [N_LGL N_LGL 1 1]
+                      Stores the value of the derivative
+                      :math:`\\frac{dL_p}{d\\xi}` at all the LGL points.
+    
+    lobatto_weights : af.Array [N_LGL 1 1 1]
+                      Lobatto weights to integrate using Gauss-Lobatto
+                      quadrature.
+                      
+    c               : float
+                      Number denoting the wave speed. Used in flux calculation.
+    
+    time            : af.Array
+                      Array containing the time for which the :math:`u` is to be
+                      calculated.
+    
+    u               : af.array [N_LGL N_Elements time.shape[0] 1]
+                      Stores :math:`u` calculated at the :math:`x` cordinates
+                      stored in the variable `element_LGL` at time :math:`t`
+                      stored in stored in variable `time`.
+        
+        
+    Parameters
+    ----------
+    Number_of_LGL_pts  : int
+                         Number of LGL nodes.
+    Number_of_elements : int
+                         Number of elements into which the domain is to be
+                         split.
+
+    '''
+
+    global N_LGL
+    global xi_LGL
+    global lBasisArray
+    global lobatto_weights
+    N_LGL       = Number_of_LGL_pts
+    xi_LGL      = lagrange.LGL_points(N_LGL)
+    lBasisArray = af.interop.np_to_af_array(
+                                         lagrange.lagrange_basis_coeffs(xi_LGL))
+
+    lobatto_weights = af.interop.np_to_af_array(\
+                                         lobatto_weight_function(N_LGL, xi_LGL)) 
+
+    global N_Elements
+    global element_nodes
+    N_Elements       = Number_of_elements
+    element_size     = af.sum((x_nodes[1] - x_nodes[0]) / N_Elements)
+    elements_xi_LGL  = af.constant(0, N_Elements, N_LGL)
+    elements         = utils.linspace(af.sum(x_nodes[0]), \
+                                  af.sum(x_nodes[1] - element_size), N_Elements)
+
+    np_element_array = np.concatenate((af.transpose(elements), 
+                                         af.transpose(elements + element_size)))
+
+    element_array = af.transpose(af.interop.np_to_af_array(np_element_array))
+    element_nodes = wave_equation.mappingXiToX(
+                                            af.transpose(element_array), xi_LGL)
+
+    global u
+    global time
+    u_init     = np.e ** (-(element_nodes) ** 2 / 0.4 ** 2)
+    time       = utils.linspace(0, 2, 200)
+    u          = af.constant(0, N_LGL, N_Elements, time.shape[0],\
+                                                           dtype = af.Dtype.f64)
+
+    u[:, :, 0] = u_init
+
+    global dLp_xi
+    dLp_xi = dLp_xi_LGL()
+
+    global c
+    c = 1.0
+
+    return
+
+
+def lobatto_weight_function(n, x):
+    '''
+    Calculates and returns the weight function for an index :math:`n`
+    and points :math: `x`.
+    
+    :math::
+        w_{n} = \\frac{2 P(x)^2}{n (n - 1)},
+        Where P(x) is $ (n - 1)^th $ index.
+    
+    Parameters
+    ----------
+    n : int
+        Index for which lobatto weight function
+    
+    x : arrayfire.Array
+        1D array of points where weight function is to be calculated.
+    
+    
+    Returns
+    -------
+    gaussLobattoWeights : arrayfire.Array
+                          An array of lobatto weight functions for
+                          the given :math: `x` points and index.
+    Reference
+    ---------
+    Gauss-Lobatto weights Wikipedia link-
+    https://en.wikipedia.org/wiki/
+    Gaussian_quadrature#Gauss.E2.80.93Lobatto_rules
+    
+    '''
+    P = sp.legendre(n - 1)
+
+    gaussLobattoWeights = (2 / (n * (n - 1)) / (P(x))**2)
+
+
+    return gaussLobattoWeights
+
 
-N_LGL       = 16
-xi_LGL      = None
-lBasisArray = None
+def lagrange_basis_function():
+    '''
+    Funtion which calculates the value of lagrange basis functions over LGL
+    nodes.
+    
+    Returns
+    -------
+    L_i    : arrayfire.Array [N 1 1 1]
+             The value of lagrange basis functions calculated over the LGL
+             nodes.
+    '''
+    xi_tile    = af.transpose(af.tile(xi_LGL, 1, N_LGL))
+    power      = af.flip(af.range(N_LGL))
+    power_tile = af.tile(power, 1, N_LGL)
+    xi_pow     = af.arith.pow(xi_tile, power_tile)
+    index      = af.range(N_LGL)
+    L_i        = af.blas.matmul(lBasisArray[index], xi_pow)
+
+    return L_i
 
-def populateGlobalVariables(N = 16):
-	'''
-	Initialize the global variables
-	
-	Parameters
-	----------
-	N : int
-		Number of LGL points.
-	'''
-
-	global N_LGL
-	global xi_LGL
-	global lBasisArray
-
-	N_LGL       = N
-	xi_LGL      = lagrange.LGL_points(N_LGL)
-
-	lBasisArray = af.interop.np_to_af_array( \
-		lagrange.lagrange_basis_coeffs(xi_LGL))
-
-	return
+def dLp_xi_LGL():
+    '''
+    Function to calculate :math: `\\frac{d L_p(\\xi_{LGL})}{d\\xi}`
+    as a 2D array of :math: `L_i (\\xi_{LGL})`. Where i varies along rows
+    and the nodes vary along the columns.
+    
+    Returns
+    -------
+    dLp_xi        : arrayfire.Array [N N 1 1]
+                    A 2D array :math: `L_i (\\xi_p)`, where i varies
+                    along dimension 1 and p varies along second dimension.
+    '''
+    differentiation_coeffs = (af.transpose(af.flip(af.tile\
+        (af.range(N_LGL), 1, N_LGL))) * lBasisArray)[:, :-1]
+
+    nodes_tile       = af.transpose(af.tile(xi_LGL, 1, N_LGL - 1))
+    power_tile       = af.flip(af.tile(af.range(N_LGL - 1), 1, N_LGL))
+    nodes_power_tile = af.pow(nodes_tile, power_tile)
+
+    dLp_xi = af.blas.matmul(differentiation_coeffs, nodes_power_tile)
+
+    return dLp_xi