Code reffactored to use Integrate() function. Unit tests passing. Doc…

…umentation and comments need to be added.

code/app/lagrange.py

@@ -72,8 +72,43 @@ def lobatto_weights(n, x):
 
     return gauss_lobatto_weights
 
+def gauss_nodes(N):
+    '''
+    '''
+    legendre = sp.legendre(N)
+    gauss_nodes = legendre.r
+    gauss_nodes.sort()
+
+
+    return gauss_nodes
 
-def lagrange_basis_coeffs(x):    
+def gaussian_weights(N, i):
+    '''
+    Returns the gaussian weights for :math: `N` Gaussian Nodes at index
+     :math: `i`.
+    
+    Parameters
+    ----------
+    N     : int
+            Number of Gaussian nodes for which the weight is t be calculated.
+            
+    i     : int
+            Index for which the Gaussian weight is required.
+    
+    Returns
+    -------
+    gaussian_weight : double 
+                      The gaussian weight.
+    
+    '''
+
+    gaussian_nodes = gauss_nodes(N)
+    gaussian_weight  = 2 / ((1 - (gaussian_nodes[i]) ** 2) * (np.polyder(sp.legendre(N))(gaussian_nodes[i])) ** 2)
+
+
+    return gaussian_weight
+
+def lagrange_polynomials(x):    
     '''
     A function to get the coefficients of the Lagrange basis polynomials for
     a given set of x nodes.
@@ -90,14 +125,20 @@ def lagrange_basis_coeffs(x):
 
     Returns
     -------
-    lagrange_basis_poly : numpy.ndarray
-                          A :math: `N \\times N` matrix containing the
-                          coefficients of the Lagrange basis polynomials such
-                          that :math:`i^{th}` lagrange polynomial will be the
-                          :math:`i^{th}` row of the matrix.
+    lagrange_basis_poly   : list
+                            A list of size `x.shape[0]` containing the analytical
+                            form of the Lagrange basis polynomials in numpy.poly1d
+                            form. This list is used in Integrate() function which
+                            requires the analytical form of the integrand.
+    lagrange_basis_coeffs : numpy.ndarray
+                            A :math: `N \\times N` matrix containing the
+                            coefficients of the Lagrange basis polynomials such
+                            that :math:`i^{th}` lagrange polynomial will be the
+                            :math:`i^{th}` row of the matrix.
     '''
     X = np.array(x)
-    lagrange_basis_poly = np.zeros([X.shape[0], X.shape[0]])
+    lagrange_basis_poly   = []
+    lagrange_basis_coeffs = np.zeros([X.shape[0], X.shape[0]])
 
     for j in np.arange(X.shape[0]):
         lagrange_basis_j = np.poly1d([1])
@@ -106,10 +147,10 @@ def lagrange_basis_coeffs(x):
             if m != j:
                 lagrange_basis_j *= np.poly1d([1, -X[m]]) \
                                     / (X[j] - X[m])
-        lagrange_basis_poly[j] = lagrange_basis_j.c
+        lagrange_basis_poly.append(lagrange_basis_j)
+        lagrange_basis_coeffs[j] = lagrange_basis_j.c
 
-    return lagrange_basis_poly
-
+    return lagrange_basis_poly, lagrange_basis_coeffs
 
 def lagrange_basis(i, x):
     '''
@@ -187,10 +228,47 @@ def dLp_xi_LGL(lagrange_coeff_array):
 
     return dLp_xi
 
+def product_lagrange_poly(x):
+    '''
+    Calculates the product of Lagrange basis polynomials in 'np.poly1d' in a 
+    2D array. This analytical form of the product of the Lagrange basis is used
+    in the calculation of A matrix using the Integrate() function.
+    '''
+    poly1d_list = lagrange_polynomials(x)[0]
+
+    poly1d_product_list = []
+
+    for i in range (x.shape[0]):
+        for j in range(x.shape[0]):
+            poly1d_product_list.append(poly1d_list[i] * poly1d_list[j])
+
+
+    return poly1d_product_list
+
 
 
 def Integrate(integrand, N_quad, scheme):
     '''
+    Integrates an analytical form of the integrand and integrates it using either
+    gaussian or lobatto quadrature.
+    
+    Parameters
+    ----------
+    integrand : numpy.poly1d
+                The analytical form of the integrand in numpy.poly1d form
+
+    N_quad    : int
+                The number of quadrature points to be used for Integration using
+                either Gaussian or Lobatto quadrature.    
+    scheme    : string
+                Specifies the method of integration to be used. Can take values
+                'gauss_quadrature' and 'lobatto_quadrature'.
+
+    Returns
+    -------
+    Integral : numpy.float64
+               The value o the definite integration performed using the specified
+               quadrature method.
     '''
     if (scheme == 'gauss_quadrature'):
         gaussian_nodes = gauss_nodes(N_quad)
@@ -203,9 +281,25 @@ def Integrate(integrand, N_quad, scheme):
 
     if (scheme == 'lobatto_quadrature'):
         lobatto_nodes = LGL_points(N_quad)
-        lobatto_weights  = lobatto_weight_function(N_quad)
+        weights  = af.np_to_af_array(lobatto_weights(N_quad, lobatto_nodes))
         value_at_lobatto_nodes = af.np_to_af_array(integrand(lobatto_nodes))
-        Integral = af.sum(value_at_lobatto_nodes * lobatto_weights)
+        Integral = af.sum(value_at_lobatto_nodes * weights)
 
 
     return Integral
+
+
+def wave_equation_lagrange_basis(u, element_no):
+    '''
+    Calculates the wave equation for a single element using the amplitude of the wave at
+    the mapped LGL points.
+    [TODO] - Explain math.
+    '''
+    amplitude_at_element_LGL = u[:, element_no]
+    lagrange_basis_polynomials = params.lagrange_poly1d_list
+
+    wave_equation_element = np.poly1d([0])
+    for i in range(0, params.N_LGL):
+        wave_equation_element += af.sum(amplitude_at_element_LGL[i]) * lagrange_basis_polynomials[i]
+
+    return wave_equation_element

code/app/params.py

@@ -26,7 +26,7 @@
 c          = 1
 
 # The total time for which the wave is to be evolved by the simulation. 
-total_time = 1
+total_time = 10
 
 # The c_lax to be used in the Lax-Friedrichs flux.
 c_lax      = 1
@@ -38,7 +38,7 @@
 lobatto_weights = af.np_to_af_array(lagrange.lobatto_weights(N_LGL, xi_LGL))
 
 # An array containing the coefficients of the lagrange basis polynomials.
-lagrange_coeffs = af.np_to_af_array(lagrange.lagrange_basis_coeffs(xi_LGL))
+lagrange_coeffs = af.np_to_af_array(lagrange.lagrange_polynomials(xi_LGL)[1])
 
 # Refer corresponding functions.
 dLp_xi               = lagrange.dLp_xi_LGL(lagrange_coeffs)
@@ -76,3 +76,26 @@
 u          = af.constant(0, N_LGL, N_Elements, time.shape[0],\
                                  dtype = af.Dtype.f64)
 u[:, :, 0] = u_init
+
+
+
+
+
+lagrange_poly1d_list = lagrange.lagrange_polynomials(xi_LGL)[0]
+
+poly1d_product_list = lagrange.product_lagrange_poly(xi_LGL)
+
+
+
+N_Gauss = 10
+
+gauss_points  = lagrange.gauss_nodes(N_Gauss)
+gauss_weights = af.np_to_af_array(np.zeros([N_Gauss]))
+
+for i in range(0, N_Gauss):
+    gauss_weights[i] = lagrange.gaussian_weights(N_Gauss, i)
+
+differential_lagrange_polynomial = []
+for i in range (0, N_LGL):
+    test_diff = np.poly1d.deriv(lagrange_poly1d_list[i])
+    differential_lagrange_polynomial.append(test_diff)

code/app/wave_equation.py

@@ -139,16 +139,26 @@ def A_matrix():
                obtained by LGL points, using Gauss-Lobatto quadrature method
                using :math: `N_LGL` points.
     '''
-    dx_dxi          = dx_dxi_numerical((params.element_mesh_nodes[0 : 2]),
-                                                                   params.xi_LGL)
-    dx_dxi_tile     = af.tile(dx_dxi, 1, params.N_LGL)
-    identity_matrix = af.identity(params.N_LGL, params.N_LGL, dtype = af.Dtype.f64)
-    A_matrix        = af.broadcast(utils.multiply, params.lobatto_weights,
-                                                 identity_matrix) * dx_dxi_tile
+    A_matrix = np.zeros([params.N_LGL, params.N_LGL])
+
+    for i in range (params.N_LGL):
+        for j in range(params.N_LGL):
+            A_matrix[i][j] = lagrange.Integrate(\
+                    params.poly1d_product_list[params.N_LGL * i + j],\
+                    9, 'lobatto_quadrature')
+
+    dx_dxi = af.mean(dx_dxi_numerical((params.element_mesh_nodes[0 : 2]),\
+                                                        params.xi_LGL))
+
+    A_matrix *= dx_dxi
+
+    A_matrix = af.np_to_af_array(A_matrix)
+
 
     return A_matrix
 
 
+
 def flux_x(u):
     '''
     A function which returns the value of flux for a given wave function u.
@@ -168,8 +178,19 @@ def flux_x(u):
 
     return flux
 
+def wave_equation_lagrange_polynomials(u):
+    '''
+    '''
+    wave_equation_in_lagrange_basis = []
+
+    for i in range(0, params.N_Elements):
+        element_wave_equation = lagrange.wave_equation_lagrange_basis(u, i)
+        wave_equation_in_lagrange_basis.append(element_wave_equation)  
 
-def volume_integral_flux(u):
+    return wave_equation_in_lagrange_basis
+
+
+def volume_integral_flux(u_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`
@@ -189,14 +210,17 @@ def volume_integral_flux(u):
                     A 1-D array of the value of the flux integral calculated
                     for various lagrange basis functions.
     '''
-
-    dLp_xi        = params.dLp_xi
-    weight_tile   = af.transpose(af.tile(params.lobatto_weights, 1, params.N_LGL))
-    dLp_xi       *= weight_tile
-    flux          = flux_x(u)
-    weight_flux   = flux
-    flux_integral = af.blas.matmul(dLp_xi, weight_flux)
-
+    wave_equation_in_lagrange_polynomials = wave_equation_lagrange_polynomials(u_n)
+    differential_lagrange_poly = params.differential_lagrange_polynomial
+    flux_integral = np.zeros([params.N_LGL, params.N_Elements])
+
+    for i in range(params.N_LGL):
+        for j in range(params.N_Elements):
+            integrand = wave_equation_in_lagrange_polynomials[j] * differential_lagrange_poly[i]
+            flux_integral[i][j] = lagrange.Integrate(integrand, 9, 'gauss_quadrature')
+
+    flux_integral = af.np_to_af_array(flux_integral)
+
     return flux_integral
 
 
@@ -209,8 +233,8 @@ def lax_friedrichs_flux(u_n):
     Parameters
     ----------
     u_n : arrayfire.Array [N_LGL N_Elements 1 1]
-            A 2D array consisting of the amplitude of the wave at the LGL nodes
-            at each element.
+          A 2D array consisting of the amplitude of the wave at the LGL nodes
+          at each element.
     '''
 
     u_iplus1_0    = af.shift(u_n[0, :], 0, -1)

code/main.py

@@ -11,5 +11,4 @@
 
 
 if __name__ == '__main__':
-
-    wave_equation.time_evolution()
+   wave_equation.time_evolution()

code/unit_test/test_waveEqn.py

@@ -189,6 +189,52 @@ def test_dLp_xi():
 
     assert af.max(reference_d_Lp_xi - params.dLp_xi) < threshold
 
+def test_A_matrix():
+    '''
+    Test function to check the A matrix obtained from wave_equation module with
+    one obtained by numerical integral solvers.
+    '''
+    threshold = 1e-8
+
+    reference_A_matrix = 0.1 * af.interop.np_to_af_array(np.array([\
+    [0.03333333333332194, 0.005783175201965206, -0.007358427761753982, \
+    0.008091331778355441, -0.008091331778233877, 0.007358427761705623, \
+    -0.00578317520224949, 0.002380952380963754], \
+
+    [0.005783175201965206, 0.19665727866729804, 0.017873132323192046,\
+    -0.01965330750343234, 0.019653307503020866, -0.017873132322725152,\
+    0.014046948476303067, -0.005783175202249493], \
+
+    [-0.007358427761753982, 0.017873132323192046, 0.31838117965137114, \
+    0.025006581762566073, -0.025006581761945083, 0.022741512832051156,\
+    -0.017873132322725152, 0.007358427761705624], \
+
+    [0.008091331778355441, -0.01965330750343234, 0.025006581762566073, \
+    0.3849615416814164, 0.027497252976343693, -0.025006581761945083, \
+    0.019653307503020863, -0.008091331778233875],  
+
+    [-0.008091331778233877, 0.019653307503020866, -0.025006581761945083, \
+    0.027497252976343693, 0.3849615416814164, 0.025006581762566073, \
+    -0.019653307503432346, 0.008091331778355443], \
+
+    [0.007358427761705623, -0.017873132322725152, 0.022741512832051156, \
+    -0.025006581761945083, 0.025006581762566073, 0.31838117965137114, \
+    0.017873132323192046, -0.0073584277617539835], \
+
+    [-0.005783175202249493, 0.014046948476303067, -0.017873132322725152, \
+    0.019653307503020863, -0.019653307503432346, 0.017873132323192046, \
+    0.19665727866729804, 0.0057831752019652065], \
+
+    [0.002380952380963754, -0.005783175202249493, 0.007358427761705624, \
+    -0.008091331778233875, 0.008091331778355443, -0.0073584277617539835, \
+    0.0057831752019652065, 0.033333333333321946]
+    ]))
+
+    test_A_matrix = wave_equation.A_matrix()
+    print(test_A_matrix.shape, reference_A_matrix.shape)
+    error_array = af.abs(reference_A_matrix - test_A_matrix)
+
+    assert af.max(error_array) < threshold
 
 def test_volume_integral_flux():
     '''