From ea3a7b531241b6241f8daa924fbb54e9e20b3728 Mon Sep 17 00:00:00 2001 From: Balavarun5 Date: Tue, 1 Aug 2017 11:05:20 +0530 Subject: [PATCH 01/19] Push changes into lax_friedrichs_flux --- code/app/global_variables.py | 5 ++++- code/app/wave_equation.py | 23 +++++++++++------------ code/main.py | 2 +- 3 files changed, 16 insertions(+), 14 deletions(-) diff --git a/code/app/global_variables.py b/code/app/global_variables.py index 7914b27..0e2c694 100644 --- a/code/app/global_variables.py +++ b/code/app/global_variables.py @@ -57,6 +57,7 @@ time = None c = None dLp_xi = None +c_lax = None def populateGlobalVariables(Number_of_LGL_pts = 8, Number_of_elements = 10): ''' @@ -99,7 +100,6 @@ def populateGlobalVariables(Number_of_LGL_pts = 8, Number_of_elements = 10): element_nodes = af.transpose(wave_equation.mappingXiToX(\ af.transpose(element_array), xi_LGL)) - global u global time u_init = np.e ** (-(element_nodes) ** 2 / 0.4 ** 2) @@ -112,7 +112,10 @@ def populateGlobalVariables(Number_of_LGL_pts = 8, Number_of_elements = 10): dLp_xi = dLp_xi_LGL() global c + global c_lax c = 1.0 + #delta_x = af.min(element_nodes, 1) + print(element_nodes) return diff --git a/code/app/wave_equation.py b/code/app/wave_equation.py index d061d40..40ee758 100644 --- a/code/app/wave_equation.py +++ b/code/app/wave_equation.py @@ -213,18 +213,17 @@ def elementFluxIntegral(n): return volumeIntegralFlux(element_n_x_nodes, gvar.u[n, :, 0]) -def lax_friedrichs_flux(u): - ''' - ''' - - u_n_0 = u[1:, 0] - u_nminus1_N_LGL = u[:gvar.N_Elements, -1] - flux_n_0 = flux_x(u_n_0) - flux_nminus1_N_LGL = flux_x(u[:gvar.N_Elements, -1]) - c_lax = gvar.c_lax - - lax_friedrichs_flux = (flux_n_0 + flux_nminus1_N_LGL) / 2 \ - - c_lax * (u_n_0 - u_nminus1_N_LGL) +def lax_friedrichs_flux(u, t_n): + ''' + ''' + + u_n_0 = u[1:, 0, t_n] + u_nminus1_N_LGL = u[:gvar.N_Elements, -1, t_n] + flux_n_0 = flux_x(u_n_0) + flux_nminus1_N_LGL = flux_x(u_nminus1_N_LGL) + + lax_friedrichs_flux = (flux_n_0 + flux_nminus1_N_LGL) / 2 \ + - gvar.c_lax * (u_n_0 - u_nminus1_N_LGL) return lax_friedrichs_flux diff --git a/code/main.py b/code/main.py index 2ac03c3..da84582 100644 --- a/code/main.py +++ b/code/main.py @@ -10,7 +10,7 @@ if __name__ == '__main__': gvar.populateGlobalVariables(8) - print(af.transpose(wave_equation.elementFluxIntegral\ + print((wave_equation.elementFluxIntegral\ (af.range(gvar.N_Elements)))) print(af.timeit(wave_equation.elementFluxIntegral,\ af.range(gvar.N_Elements))) From b8526a2747fe1b489d283f0fa1cc917cd7d55167 Mon Sep 17 00:00:00 2001 From: Balavarun5 Date: Wed, 2 Aug 2017 17:19:37 +0530 Subject: [PATCH 02/19] Added functions to obtain the lax-friedrichs-flux and surface term. --- code/app/global_variables.py | 5 ++- code/app/wave_equation.py | 75 ++++++++++++++++++++++++++++------ code/main.py | 3 +- code/unit_test/test_waveEqn.py | 14 +++++++ 4 files changed, 81 insertions(+), 16 deletions(-) diff --git a/code/app/global_variables.py b/code/app/global_variables.py index 0e2c694..3c6b61c 100644 --- a/code/app/global_variables.py +++ b/code/app/global_variables.py @@ -104,7 +104,8 @@ def populateGlobalVariables(Number_of_LGL_pts = 8, Number_of_elements = 10): global time u_init = np.e ** (-(element_nodes) ** 2 / 0.4 ** 2) time = utils.linspace(0, 2, 200) - u = af.constant(0, (N_Elements), N_LGL, time.shape[0]) + u = af.constant(0, N_Elements, N_LGL, time.shape[0],\ + dtype = af.Dtype.f64) u[:, :, 0] = u_init @@ -114,8 +115,8 @@ def populateGlobalVariables(Number_of_LGL_pts = 8, Number_of_elements = 10): global c global c_lax c = 1.0 + c_lax = 0.1 #delta_x = af.min(element_nodes, 1) - print(element_nodes) return diff --git a/code/app/wave_equation.py b/code/app/wave_equation.py index 40ee758..caf7b2f 100644 --- a/code/app/wave_equation.py +++ b/code/app/wave_equation.py @@ -11,8 +11,8 @@ def Li_Lp_xi(L_i_xi, L_p_xi): Parameters ---------- L_i_xi : arrayfire.Array [1 N N 1] - A 2D array :math:`L_i` obtained at LGL points calculated at the - LGL nodes :math:`N_LGL`. + A 2D array :math:`L_i` obtained at LGL points calculated at the + LGL nodes :math:`N_LGL`. L_p_xi : arrayfire.Array [N 1 N 1] A 2D array :math:`L_p` obtained at LGL points calculated at the @@ -213,19 +213,70 @@ def elementFluxIntegral(n): return volumeIntegralFlux(element_n_x_nodes, gvar.u[n, :, 0]) -def lax_friedrichs_flux(u, t_n): +def lax_friedrichs_flux(t_n): ''' + A function which calculates the lax-friedrichs_flux :math:`f_i` using. + :math:`f_i = \\frac{F(u^{i + 1}_0) + F(u^i_{N_{LGL} - 1})}{2} - \frac + {\Delta x}{2\Delta t} (u^{i + 1}_0 - u^i_{N_{LGL} - 1})` + + Parameters + ---------- + u : arrayfire.Array[N_Elements N_LGL 1 1] + A 2D array consisting of the amplitude of the wave at the LGL nodes + at each element. + + t_n : float + The timestep at which the lax-friedrichs-flux is to be calculated. ''' - - u_n_0 = u[1:, 0, t_n] - u_nminus1_N_LGL = u[:gvar.N_Elements, -1, t_n] - flux_n_0 = flux_x(u_n_0) - flux_nminus1_N_LGL = flux_x(u_nminus1_N_LGL) - lax_friedrichs_flux = (flux_n_0 + flux_nminus1_N_LGL) / 2 \ - - gvar.c_lax * (u_n_0 - u_nminus1_N_LGL) - - return lax_friedrichs_flux + u_iplus1_0 = af.shift(gvar.u[:, 0, t_n], -1) + u_i_N_LGL = gvar.u[:, -1, t_n] + flux_iplus1_0 = flux_x(u_iplus1_0) + flux_i_N_LGL = flux_x(u_i_N_LGL) + + lax_friedrichs_flux = (flux_iplus1_0 + flux_i_N_LGL) / 2 \ + - gvar.c_lax * (u_iplus1_0 - u_i_N_LGL) + + return lax_friedrichs_flux + + +def surface_term(t_n): + ''' + A function which is used to calculate the surface term, + :math:`L_p (1) f_i - L_p (-1) f_{i - 1}` + using the lax_friedrichs_flux function and the dLp_xi_LGL function in gvar + module. + + Parameters + ---------- + t_n : float + The timestep at which the surface term is to be calculated. + + Returns + ------- + surface_term : 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. + + Reference + --------- + Link to PDF describing the algorithm to obtain the surface term. + + https://cocalc.com/projects/1b7f404c-87ba-40d0-816c-2eba17466aa8/files\ + /PM\_2\_5/wave\_equation/documents/surface\_term/surface\_term.pdf + + ''' + L_p_minus1 = gvar.lagrange_basis_function()[:, 0] + L_p_1 = gvar.lagrange_basis_function()[:, -1] + f_i = af.transpose(lax_friedrichs_flux(t_n)) + f_iminus1 = af.transpose(af.shift(lax_friedrichs_flux(t_n), 1)) + surface_term = af.blas.matmul(L_p_1, f_i) - af.blas.matmul(L_p_minus1,\ + f_iminus1) + + return surface_term + def b_vector(u_n): ''' diff --git a/code/main.py b/code/main.py index da84582..9b685d2 100644 --- a/code/main.py +++ b/code/main.py @@ -10,7 +10,6 @@ if __name__ == '__main__': gvar.populateGlobalVariables(8) - print((wave_equation.elementFluxIntegral\ - (af.range(gvar.N_Elements)))) print(af.timeit(wave_equation.elementFluxIntegral,\ af.range(gvar.N_Elements))) + print(wave_equation.surface_term(0)) diff --git a/code/unit_test/test_waveEqn.py b/code/unit_test/test_waveEqn.py index dc78b1e..5690ace 100644 --- a/code/unit_test/test_waveEqn.py +++ b/code/unit_test/test_waveEqn.py @@ -340,3 +340,17 @@ def test_volume_integral_flux(): assert (af.max(af.abs(numerical_flux_sum - referenceFluxIntegral_sum))\ < threshold) + +def test_lax_friedrichs_flux(): + ''' + A test function to test the lax_friedrichs_flux function in wave_equation + module. [TODO] : confirm. + ''' + threshold = 1e-14 + gvar.populateGlobalVariables(8, 10) + + f_i = wave_equation.lax_friedrichs_flux(0) + #The lax friedrichs flux at timestep 0 should just be a list of the element + #boundaries of the elements at LGL nodes. + analytical_lax_friedrichs_flux = gvar.u[:, -1, 0] + assert af.max(analytical_lax_friedrichs_flux - f_i) < threshold From 3f92f4abe40d268d7aba93616e8cfcc790856832 Mon Sep 17 00:00:00 2001 From: Balavarun5 Date: Fri, 4 Aug 2017 18:05:01 +0530 Subject: [PATCH 03/19] Time evolution of wave equation obtained. --- code/app/global_variables.py | 21 ++++--- code/app/wave_equation.py | 105 ++++++++++++++++++++++++++++----- code/main.py | 7 ++- code/unit_test/test_waveEqn.py | 65 ++++++++++++++++---- 4 files changed, 162 insertions(+), 36 deletions(-) diff --git a/code/app/global_variables.py b/code/app/global_variables.py index 3c6b61c..d9e1cba 100644 --- a/code/app/global_variables.py +++ b/code/app/global_variables.py @@ -100,23 +100,26 @@ def populateGlobalVariables(Number_of_LGL_pts = 8, Number_of_elements = 10): element_nodes = af.transpose(wave_equation.mappingXiToX(\ af.transpose(element_array), xi_LGL)) + global c + global c_lax + global delta_t + c = 1.0 + delta_x = af.min((element_nodes - af.shift(element_nodes, 0, 1))[:, 1:]) + delta_t = delta_x / (10 * c) + c_lax = delta_x / (2 * delta_t) + 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_Elements, N_LGL, time.shape[0],\ - dtype = af.Dtype.f64) + time = utils.linspace(0, 2, int(2 / delta_t)) + u_init = np.e ** (-(element_nodes) ** 2 / 0.4 ** 2) + u = af.constant(0, N_Elements, N_LGL, time.shape[0],\ + dtype = af.Dtype.f64) u[:, :, 0] = u_init global dLp_xi dLp_xi = dLp_xi_LGL() - global c - global c_lax - c = 1.0 - c_lax = 0.1 - #delta_x = af.min(element_nodes, 1) return diff --git a/code/app/wave_equation.py b/code/app/wave_equation.py index caf7b2f..3730473 100644 --- a/code/app/wave_equation.py +++ b/code/app/wave_equation.py @@ -5,6 +5,36 @@ from app import lagrange from utils import utils from app import global_variables as gvar +from matplotlib import pyplot as plt +import pylab as pl + +plt.rcParams['figure.figsize'] = 9.6, 6. +plt.rcParams['figure.dpi'] = 300 +plt.rcParams['image.cmap'] = 'jet' +plt.rcParams['lines.linewidth'] = 1.5 +plt.rcParams['font.family'] = 'serif' +plt.rcParams['font.weight'] = 'bold' +plt.rcParams['font.size'] = 20 +plt.rcParams['font.sans-serif'] = 'serif' +plt.rcParams['text.usetex'] = True +plt.rcParams['axes.linewidth'] = 1.5 +plt.rcParams['axes.titlesize'] = 'medium' +plt.rcParams['axes.labelsize'] = 'medium' +plt.rcParams['xtick.major.size'] = 8 +plt.rcParams['xtick.minor.size'] = 4 +plt.rcParams['xtick.major.pad'] = 8 +plt.rcParams['xtick.minor.pad'] = 8 +plt.rcParams['xtick.color'] = 'k' +plt.rcParams['xtick.labelsize'] = 'medium' +plt.rcParams['xtick.direction'] = 'in' +plt.rcParams['ytick.major.size'] = 8 +plt.rcParams['ytick.minor.size'] = 4 +plt.rcParams['ytick.major.pad'] = 8 +plt.rcParams['ytick.minor.pad'] = 8 +plt.rcParams['ytick.color'] = 'k' +plt.rcParams['ytick.labelsize'] = 'medium' +plt.rcParams['ytick.direction'] = 'in' + def Li_Lp_xi(L_i_xi, L_p_xi): ''' @@ -189,7 +219,7 @@ def volumeIntegralFlux(element_nodes, u): return af.reorder(flux_integral, 2, 1, 0) -def elementFluxIntegral(n): +def elementFluxIntegral(n, t_n): ''' Function which reorders the element numbers which can then be passed into volumeIntegralFlux. @@ -210,7 +240,7 @@ def elementFluxIntegral(n): ''' element_n_x_nodes = af.reorder(gvar.element_nodes[n], 1, 0, 2) - return volumeIntegralFlux(element_n_x_nodes, gvar.u[n, :, 0]) + return volumeIntegralFlux(element_n_x_nodes, gvar.u[n, :, t_n]) def lax_friedrichs_flux(t_n): @@ -221,7 +251,7 @@ def lax_friedrichs_flux(t_n): Parameters ---------- - u : arrayfire.Array[N_Elements N_LGL 1 1] + u : arrayfire.Array [N_Elements N_LGL 1 1] A 2D array consisting of the amplitude of the wave at the LGL nodes at each element. @@ -249,12 +279,13 @@ def surface_term(t_n): Parameters ---------- - t_n : float + t_n : double The timestep at which the surface term is to be calculated. Returns ------- - surface_term : The surface term represented in the form of an array, + 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 @@ -272,18 +303,64 @@ def surface_term(t_n): L_p_1 = gvar.lagrange_basis_function()[:, -1] f_i = af.transpose(lax_friedrichs_flux(t_n)) f_iminus1 = af.transpose(af.shift(lax_friedrichs_flux(t_n), 1)) - surface_term = af.blas.matmul(L_p_1, f_i) - af.blas.matmul(L_p_minus1,\ - f_iminus1) + surface_term = af.blas.matmul(L_p_1, f_i) - af.blas.matmul(L_p_minus1, f_iminus1) + + return af.transpose(surface_term) + + +def b_vector(t_n): + ''' + A function which returns the b vector for N_Elements number of elements. + + Parameters + ---------- + t_n : double + + Returns + ------- + b_vector_array : arrayfire.array + ''' + u_n_A_matrix = af.blas.matmul(gvar.u[:, :, t_n], A_matrix()) + volume_integral = elementFluxIntegral(af.range(gvar.N_Elements), t_n) + surfaceTerm = surface_term(t_n) + b_vector_array = u_n_A_matrix + gvar.delta_t * (volume_integral -\ + surfaceTerm) + - return surface_term + return b_vector_array -def b_vector(u_n): +def time_evolution(): ''' - NOTE - Incomplete. + Function which solves the wave equation + :math: `u^{t_n + 1} = b(t_n) \\times A` + iterated over time steps t_n and then plots :math: `x` against the amplitude + of the wave. The images are then stored in Wave folder. ''' - int_u_ni_Lp_Li = af.blas.matmul(A_matrix(), af.transpose(u_n)) - int_flux_dLp_dxi = volume_integral_flux(u_n) - return + A_inverse = af.lapack.inverse(A_matrix()) + element_nodes = gvar.element_nodes + + for t_n in range(0, gvar.time.shape[0] - 1): + + gvar.u[:, :, t_n + 1] = af.blas.matmul(b_vector(t_n), A_inverse) + + print('u calculated!') + + #for t_n in range(0, gvar.time.shape[0] - 1): + + #if (t_n % 6) == 0: + + #fig = plt.figure() + #x = (af.transpose(gvar.element_nodes)) + #y = (af.transpose(gvar.u[:, :, t_n])) + + #plt.plot(x, y) + #plt.xlabel('x') + #plt.ylabel('Amplitude') + #fig.savefig('Wave/%04d' %(t_n / 6) + '.png') + #plt.close('all') + + #print(t_n) + + return diff --git a/code/main.py b/code/main.py index 9b685d2..a7e5a45 100644 --- a/code/main.py +++ b/code/main.py @@ -2,14 +2,15 @@ import arrayfire as af af.set_backend('opencl') +af.set_device(1) import numpy as np from app import global_variables as gvar from unit_test import test_waveEqn from app import wave_equation from app import lagrange + + if __name__ == '__main__': gvar.populateGlobalVariables(8) - print(af.timeit(wave_equation.elementFluxIntegral,\ - af.range(gvar.N_Elements))) - print(wave_equation.surface_term(0)) + print(gvar.time.shape) diff --git a/code/unit_test/test_waveEqn.py b/code/unit_test/test_waveEqn.py index 5690ace..1d3caf9 100644 --- a/code/unit_test/test_waveEqn.py +++ b/code/unit_test/test_waveEqn.py @@ -290,14 +290,14 @@ def test_volume_integral_flux(): -0.003445512010153811, 0.0176615086879261])) numerical_flux0 = af.transpose(wave_equation.elementFluxIntegral\ - (af.range(gvar.N_Elements))[0]) + (af.range(gvar.N_Elements), 0)[0]) referenceFluxIntegral1 = af.interop.np_to_af_array(np.array([-0.018969769374 ,-0.00431252844519, -0.00882630935977, -0.0144355176966, -0.019612124119 ,-0.0209837936827, -0.0154359890788, 0.102576031756])) numerical_flux1 = af.transpose(wave_equation.elementFluxIntegral\ - (af.range(gvar.N_Elements))[1]) + (af.range(gvar.N_Elements), 0)[1]) referenceFluxIntegral2 = af.interop.np_to_af_array(np.array([-0.108222418798 @@ -305,7 +305,7 @@ def test_volume_integral_flux(): -0.0557970236273, -0.0374764132459, 0.361310165819])) numerical_flux2 = af.transpose(wave_equation.elementFluxIntegral\ - (af.range(gvar.N_Elements))[2]) + (af.range(gvar.N_Elements), 0)[2]) referenceFluxIntegral3 = af.interop.np_to_af_array(np.array([-0.374448714304 @@ -313,14 +313,14 @@ def test_volume_integral_flux(): , -0.0714664112839, -0.0422339853622, 0.771847201979])) numerical_flux3 = af.transpose(wave_equation.elementFluxIntegral\ - (af.range(gvar.N_Elements))[3]) + (af.range(gvar.N_Elements), 0)[3]) referenceFluxIntegral4 = af.interop.np_to_af_array(np.array([-0.785754362849 , -0.0396035640187, -0.0579313769517, -0.0569022801117, -0.0392041960688, -0.0172295769141, -0.00337464521455, 1.00000000213])) numerical_flux4 = af.transpose(wave_equation.elementFluxIntegral\ - (af.range(gvar.N_Elements))[4]) + (af.range(gvar.N_Elements), 0)[4]) referenceFluxIntegral5 = af.interop.np_to_af_array(np.array([-1.00000000213, @@ -328,7 +328,7 @@ def test_volume_integral_flux(): 0.0579313769517, 0.0396035640187, 0.785754362849])) numerical_flux5 = af.transpose(wave_equation.elementFluxIntegral\ - (af.range(gvar.N_Elements))[5]) + (af.range(gvar.N_Elements), 0)[5]) numerical_flux_sum = numerical_flux0 + numerical_flux1 + numerical_flux2 \ @@ -344,13 +344,58 @@ def test_volume_integral_flux(): def test_lax_friedrichs_flux(): ''' A test function to test the lax_friedrichs_flux function in wave_equation - module. [TODO] : confirm. + module. ''' threshold = 1e-14 gvar.populateGlobalVariables(8, 10) f_i = wave_equation.lax_friedrichs_flux(0) - #The lax friedrichs flux at timestep 0 should just be a list of the element - #boundaries of the elements at LGL nodes. + #The lax friedrichs flux at timestep 0 should just be a list of the + #amplitude at element boundaries. analytical_lax_friedrichs_flux = gvar.u[:, -1, 0] - assert af.max(analytical_lax_friedrichs_flux - f_i) < threshold + assert af.max(af.abs(analytical_lax_friedrichs_flux - f_i)) < threshold + + +def test_surface_term(): + ''' + A test function to test the surface_term function in the wave_equation + module using analytical Lax-Friedrichs flux. + ''' + threshold = 1e-13 + gvar.populateGlobalVariables(8, 10) + analytical_f_i = af.transpose(gvar.u[:, -1, 0]) + analytical_f_i_minus1 = af.transpose(af.shift(gvar.u[:, -1, 0], 1)) + + L_p_1 = af.constant(0, gvar.N_LGL, dtype = af.Dtype.f64) + L_p_1[gvar.N_LGL - 1] = 1 + + L_p_minus1 = af.constant(0, gvar.N_LGL, dtype = af.Dtype.f64) + L_p_minus1[0] = 1 + + analytical_surface_term = af.blas.matmul(L_p_1, analytical_f_i)\ + - af.blas.matmul(L_p_minus1, analytical_f_i_minus1) + + numerical_surface_term = af.transpose(wave_equation.surface_term(0)) + assert af.max(af.abs(analytical_surface_term - numerical_surface_term)) \ + < threshold + return analytical_surface_term + + +def test_b_vector(): + ''' + A test function to check the b vector obtained analytically and compare it + with the one returned by b_vector function in wave_equation module + ''' + threshold = 1e-13 + gvar.populateGlobalVariables(8, 10) + + u_n_A_matrix = af.blas.matmul(gvar.u[:, :, 0],\ + wave_equation.A_matrix()) + volume_integral_flux = wave_equation.elementFluxIntegral\ + (af.range(gvar.N_Elements), 0) + surface_term = test_surface_term() + b_vector_analytical = u_n_A_matrix + (volume_integral_flux -\ + af.transpose(surface_term)) * gvar.delta_t + b_vector_array = wave_equation.b_vector(0) + + assert (b_vector_analytical - b_vector_array) < threshold From 1014ef472bc3746dd3e0acdb46327708629d6f8f Mon Sep 17 00:00:00 2001 From: Balavarun5 Date: Sat, 5 Aug 2017 14:28:02 +0530 Subject: [PATCH 04/19] Simulating for higher wave speeds. --- code/app/global_variables.py | 4 ++-- code/app/wave_equation.py | 26 ++++++++++++++------------ code/main.py | 1 + 3 files changed, 17 insertions(+), 14 deletions(-) diff --git a/code/app/global_variables.py b/code/app/global_variables.py index d9e1cba..75886b2 100644 --- a/code/app/global_variables.py +++ b/code/app/global_variables.py @@ -103,14 +103,14 @@ def populateGlobalVariables(Number_of_LGL_pts = 8, Number_of_elements = 10): global c global c_lax global delta_t - c = 1.0 + c = 5.0 delta_x = af.min((element_nodes - af.shift(element_nodes, 0, 1))[:, 1:]) delta_t = delta_x / (10 * c) c_lax = delta_x / (2 * delta_t) global u global time - time = utils.linspace(0, 2, int(2 / delta_t)) + time = utils.linspace(0, 8, int(8 / delta_t)) u_init = np.e ** (-(element_nodes) ** 2 / 0.4 ** 2) u = af.constant(0, N_Elements, N_LGL, time.shape[0],\ dtype = af.Dtype.f64) diff --git a/code/app/wave_equation.py b/code/app/wave_equation.py index 3730473..15655e5 100644 --- a/code/app/wave_equation.py +++ b/code/app/wave_equation.py @@ -9,7 +9,7 @@ import pylab as pl plt.rcParams['figure.figsize'] = 9.6, 6. -plt.rcParams['figure.dpi'] = 300 +plt.rcParams['figure.dpi'] = 100 plt.rcParams['image.cmap'] = 'jet' plt.rcParams['lines.linewidth'] = 1.5 plt.rcParams['font.family'] = 'serif' @@ -340,6 +340,7 @@ def time_evolution(): A_inverse = af.lapack.inverse(A_matrix()) element_nodes = gvar.element_nodes + delta_t = gvar.delta_t for t_n in range(0, gvar.time.shape[0] - 1): @@ -347,20 +348,21 @@ def time_evolution(): print('u calculated!') - #for t_n in range(0, gvar.time.shape[0] - 1): + for t_n in range(0, gvar.time.shape[0] - 1): - #if (t_n % 6) == 0: + if (t_n % 500) == 0: - #fig = plt.figure() - #x = (af.transpose(gvar.element_nodes)) - #y = (af.transpose(gvar.u[:, :, t_n])) + fig = plt.figure() + x = (af.transpose(gvar.element_nodes)) + y = (af.transpose(gvar.u[:, :, t_n])) - #plt.plot(x, y) - #plt.xlabel('x') - #plt.ylabel('Amplitude') - #fig.savefig('Wave/%04d' %(t_n / 6) + '.png') - #plt.close('all') + plt.plot(x, y) + plt.xlabel('x') + plt.ylabel('Amplitude') + plt.title('Time = %f' %(t_n * delta_t)) + fig.savefig('Wave/%04d' %(t_n / 500) + '.png') + plt.close('all') - #print(t_n) + print(t_n) return diff --git a/code/main.py b/code/main.py index a7e5a45..59bf871 100644 --- a/code/main.py +++ b/code/main.py @@ -14,3 +14,4 @@ gvar.populateGlobalVariables(8) print(gvar.time.shape) + wave_equation.time_evolution() From c5d0a2f05185cdc34763d773a036a33ed546225e Mon Sep 17 00:00:00 2001 From: Balavarun5 Date: Sat, 5 Aug 2017 17:31:56 +0530 Subject: [PATCH 05/19] The images would be created in 1D_Wave_images folder, stitch them to form a movie. --- .gitignore | 23 ++++++++++++++++++++++- code/app/global_variables.py | 4 ++-- code/app/wave_equation.py | 7 ++++--- 3 files changed, 28 insertions(+), 6 deletions(-) diff --git a/.gitignore b/.gitignore index f5bc6e3..fe8afbf 100644 --- a/.gitignore +++ b/.gitignore @@ -2,7 +2,28 @@ __pycache__/ *.py[cod] *$py.class - +*.mp4 +*.m4p +*.m4v +*.mpg +*.mp2 +*.mpeg +*.mpe +*.mpv +*.mpg +*.mpeg +*.m2v +*.m4v +*.3gp +*.flv +*.f4v +*.f4p +*.f4a +*.f4b +*.jpeg +*.jpg +*.png +*.svg # C extensions *.so diff --git a/code/app/global_variables.py b/code/app/global_variables.py index 75886b2..e088b52 100644 --- a/code/app/global_variables.py +++ b/code/app/global_variables.py @@ -103,14 +103,14 @@ def populateGlobalVariables(Number_of_LGL_pts = 8, Number_of_elements = 10): global c global c_lax global delta_t - c = 5.0 + c = 1.0 delta_x = af.min((element_nodes - af.shift(element_nodes, 0, 1))[:, 1:]) delta_t = delta_x / (10 * c) c_lax = delta_x / (2 * delta_t) global u global time - time = utils.linspace(0, 8, int(8 / delta_t)) + time = utils.linspace(0, 30, int(30 / delta_t)) u_init = np.e ** (-(element_nodes) ** 2 / 0.4 ** 2) u = af.constant(0, N_Elements, N_LGL, time.shape[0],\ dtype = af.Dtype.f64) diff --git a/code/app/wave_equation.py b/code/app/wave_equation.py index 15655e5..a761e49 100644 --- a/code/app/wave_equation.py +++ b/code/app/wave_equation.py @@ -303,7 +303,8 @@ def surface_term(t_n): L_p_1 = gvar.lagrange_basis_function()[:, -1] f_i = af.transpose(lax_friedrichs_flux(t_n)) f_iminus1 = af.transpose(af.shift(lax_friedrichs_flux(t_n), 1)) - surface_term = af.blas.matmul(L_p_1, f_i) - af.blas.matmul(L_p_minus1, f_iminus1) + surface_term = af.blas.matmul(L_p_1, f_i) - af.blas.matmul(L_p_minus1,\ + f_iminus1) return af.transpose(surface_term) @@ -350,7 +351,7 @@ def time_evolution(): for t_n in range(0, gvar.time.shape[0] - 1): - if (t_n % 500) == 0: + if (t_n % 80) == 0: fig = plt.figure() x = (af.transpose(gvar.element_nodes)) @@ -360,7 +361,7 @@ def time_evolution(): plt.xlabel('x') plt.ylabel('Amplitude') plt.title('Time = %f' %(t_n * delta_t)) - fig.savefig('Wave/%04d' %(t_n / 500) + '.png') + fig.savefig('1D_Wave_images/%04d' %(t_n / 80) + '.png') plt.close('all') print(t_n) From bc14e39f6fb1629f54edd4215c807534cf5dc646 Mon Sep 17 00:00:00 2001 From: Balavarun5 Date: Mon, 7 Aug 2017 21:54:09 +0530 Subject: [PATCH 06/19] Changed amplitude from being represented as an N_Elements x N_LGL x time to N_LGL x N_Elements x time --- code/app/global_variables.py | 10 ++-- code/app/wave_equation.py | 76 ++++++++++----------------- code/main.py | 5 +- code/unit_test/test_waveEqn.py | 95 +++++++++++----------------------- 4 files changed, 67 insertions(+), 119 deletions(-) diff --git a/code/app/global_variables.py b/code/app/global_variables.py index e088b52..afe8907 100644 --- a/code/app/global_variables.py +++ b/code/app/global_variables.py @@ -97,22 +97,22 @@ def populateGlobalVariables(Number_of_LGL_pts = 8, Number_of_elements = 10): af.transpose(elements + element_size))) element_array = (af.transpose(af.interop.np_to_af_array(np_element_array))) - element_nodes = af.transpose(wave_equation.mappingXiToX(\ + element_nodes = (wave_equation.mappingXiToX(\ af.transpose(element_array), xi_LGL)) global c global c_lax global delta_t - c = 1.0 - delta_x = af.min((element_nodes - af.shift(element_nodes, 0, 1))[:, 1:]) + c = 4.0 + delta_x = af.min((element_nodes - af.shift(element_nodes, 1, 0))[1:, :]) delta_t = delta_x / (10 * c) c_lax = delta_x / (2 * delta_t) global u global time - time = utils.linspace(0, 30, int(30 / delta_t)) + time = utils.linspace(0, 20, int(20 / delta_t)) u_init = np.e ** (-(element_nodes) ** 2 / 0.4 ** 2) - u = af.constant(0, N_Elements, N_LGL, time.shape[0],\ + u = af.constant(0, N_LGL, N_Elements, time.shape[0],\ dtype = af.Dtype.f64) u[:, :, 0] = u_init diff --git a/code/app/wave_equation.py b/code/app/wave_equation.py index a761e49..60e7716 100644 --- a/code/app/wave_equation.py +++ b/code/app/wave_equation.py @@ -174,7 +174,7 @@ def flux_x(u): Parameters ---------- - u : arrayfire.Array [N 1 1 1] + u : arrayfire.Array A 1-D array which contains the value of wave function. Returns @@ -210,37 +210,13 @@ def volumeIntegralFlux(element_nodes, u): for various lagrange basis functions. ''' - dLp_xi = af.tile(af.transpose(gvar.dLp_xi), 1, 1, gvar.N_Elements) - weight_tile = af.tile(gvar.lobatto_weights, 1, gvar.N_LGL,\ - gvar.N_Elements) - flux = af.reorder(flux_x(u), 1, 2, 0) - flux_u_tile = af.tile(flux, 1, gvar.N_LGL, 1) - flux_integral = af.sum((weight_tile * dLp_xi * flux_u_tile), 0) + dLp_xi = gvar.dLp_xi + weight_tile = af.tile(gvar.lobatto_weights, 1, gvar.N_Elements) + flux = flux_x(u) + weight_flux = weight_tile * flux + flux_integral = af.blas.matmul(dLp_xi, weight_flux) - return af.reorder(flux_integral, 2, 1, 0) - -def elementFluxIntegral(n, t_n): - ''' - Function which reorders the element numbers which can then be passed - into volumeIntegralFlux. - - Parameters - ---------- - n : Element numbers for which the flux integral is to be calculated, - Passing an array of 0 to N_Elements - 1 would give the flux integral - for all elements at all :math:`p` - - Returns - ------- - volumeIntegralFlux(element_n_x_nodes,\ - gvar.u[n, :, 0]) : arrayfire.Array - An array of :math:`\\int_{-1}^1 f(u) - \\frac{d L_p}{d\\xi} d\\xi` for all - elements - ''' - element_n_x_nodes = af.reorder(gvar.element_nodes[n], 1, 0, 2) - - return volumeIntegralFlux(element_n_x_nodes, gvar.u[n, :, t_n]) + return flux_integral def lax_friedrichs_flux(t_n): @@ -251,7 +227,7 @@ def lax_friedrichs_flux(t_n): Parameters ---------- - u : arrayfire.Array [N_Elements N_LGL 1 1] + u : 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. @@ -259,8 +235,8 @@ def lax_friedrichs_flux(t_n): The timestep at which the lax-friedrichs-flux is to be calculated. ''' - u_iplus1_0 = af.shift(gvar.u[:, 0, t_n], -1) - u_i_N_LGL = gvar.u[:, -1, t_n] + u_iplus1_0 = af.shift(gvar.u[0, :, t_n], 0, -1) + u_i_N_LGL = gvar.u[-1, :, t_n] flux_iplus1_0 = flux_x(u_iplus1_0) flux_i_N_LGL = flux_x(u_i_N_LGL) @@ -301,12 +277,12 @@ def surface_term(t_n): ''' L_p_minus1 = gvar.lagrange_basis_function()[:, 0] L_p_1 = gvar.lagrange_basis_function()[:, -1] - f_i = af.transpose(lax_friedrichs_flux(t_n)) - f_iminus1 = af.transpose(af.shift(lax_friedrichs_flux(t_n), 1)) + f_i = lax_friedrichs_flux(t_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 af.transpose(surface_term) + return surface_term def b_vector(t_n): @@ -321,11 +297,9 @@ def b_vector(t_n): ------- b_vector_array : arrayfire.array ''' - u_n_A_matrix = af.blas.matmul(gvar.u[:, :, t_n], A_matrix()) - volume_integral = elementFluxIntegral(af.range(gvar.N_Elements), t_n) + volume_integral = volumeIntegralFlux(gvar.element_nodes, gvar.u[:, :, t_n]) surfaceTerm = surface_term(t_n) - b_vector_array = u_n_A_matrix + gvar.delta_t * (volume_integral -\ - surfaceTerm) + b_vector_array = gvar.delta_t * (volume_integral - surfaceTerm) return b_vector_array @@ -339,31 +313,37 @@ def time_evolution(): of the wave. The images are then stored in Wave folder. ''' - A_inverse = af.lapack.inverse(A_matrix()) + A_inverse = af.transpose(af.lapack.inverse(A_matrix())) element_nodes = gvar.element_nodes delta_t = gvar.delta_t for t_n in range(0, gvar.time.shape[0] - 1): - gvar.u[:, :, t_n + 1] = af.blas.matmul(b_vector(t_n), A_inverse) - + gvar.u[:, :, t_n + 1] = gvar.u[:, :, t_n] + af.blas.matmul(A_inverse, b_vector(t_n)) + if (t_n % 100 == 0): + print('timestep', t_n) print('u calculated!') + for t_n in range(0, gvar.time.shape[0] - 1): - if (t_n % 80) == 0: + if t_n % 100 == 0: fig = plt.figure() - x = (af.transpose(gvar.element_nodes)) - y = (af.transpose(gvar.u[:, :, t_n])) + x = gvar.element_nodes + y = gvar.u[:, :, t_n] plt.plot(x, y) plt.xlabel('x') plt.ylabel('Amplitude') plt.title('Time = %f' %(t_n * delta_t)) - fig.savefig('1D_Wave_images/%04d' %(t_n / 80) + '.png') + fig.savefig('1D_Wave_images/%04d' %(t_n / 100) + '.png') plt.close('all') print(t_n) + + + + return diff --git a/code/main.py b/code/main.py index 59bf871..9678c4e 100644 --- a/code/main.py +++ b/code/main.py @@ -13,5 +13,6 @@ if __name__ == '__main__': gvar.populateGlobalVariables(8) - print(gvar.time.shape) - wave_equation.time_evolution() + print(gvar.time.shape[0]) + print(wave_equation.volumeIntegralFlux(gvar.element_nodes, gvar.u[:, :, 0])) + print(wave_equation.time_evolution()) diff --git a/code/unit_test/test_waveEqn.py b/code/unit_test/test_waveEqn.py index 1d3caf9..e2b0ae6 100644 --- a/code/unit_test/test_waveEqn.py +++ b/code/unit_test/test_waveEqn.py @@ -284,62 +284,31 @@ def test_volume_integral_flux(): threshold = 4 * 1e-8 gvar.populateGlobalVariables(8) - referenceFluxIntegral0 = af.interop.np_to_af_array(np.array( + referenceFluxIntegral = af.transpose(af.interop.np_to_af_array(np.array([ [-0.002016634876668093, -0.000588597708116113, -0.0013016773719126333,\ - -0.002368387579324652, -0.003620502047659841, -0.004320197094090966, - -0.003445512010153811, 0.0176615086879261])) - - numerical_flux0 = af.transpose(wave_equation.elementFluxIntegral\ - (af.range(gvar.N_Elements), 0)[0]) - - referenceFluxIntegral1 = af.interop.np_to_af_array(np.array([-0.018969769374 - ,-0.00431252844519, -0.00882630935977, -0.0144355176966, -0.019612124119 - ,-0.0209837936827, -0.0154359890788, 0.102576031756])) - - numerical_flux1 = af.transpose(wave_equation.elementFluxIntegral\ - (af.range(gvar.N_Elements), 0)[1]) - - - referenceFluxIntegral2 = af.interop.np_to_af_array(np.array([-0.108222418798 - ,-0.0179274222595, -0.0337807018822, -0.0492589052599, -0.0588472807471, - -0.0557970236273, -0.0374764132459, 0.361310165819])) - - numerical_flux2 = af.transpose(wave_equation.elementFluxIntegral\ - (af.range(gvar.N_Elements), 0)[2]) - - - referenceFluxIntegral3 = af.interop.np_to_af_array(np.array([-0.374448714304 - , -0.0399576371245, -0.0683852285846, -0.0869229749357, -0.0884322503841 - , -0.0714664112839, -0.0422339853622, 0.771847201979])) - - numerical_flux3 = af.transpose(wave_equation.elementFluxIntegral\ - (af.range(gvar.N_Elements), 0)[3]) - - referenceFluxIntegral4 = af.interop.np_to_af_array(np.array([-0.785754362849 - , -0.0396035640187, -0.0579313769517, -0.0569022801117, - -0.0392041960688, -0.0172295769141, -0.00337464521455, 1.00000000213])) - - numerical_flux4 = af.transpose(wave_equation.elementFluxIntegral\ - (af.range(gvar.N_Elements), 0)[4]) - - - referenceFluxIntegral5 = af.interop.np_to_af_array(np.array([-1.00000000213, - 0.00337464521455, 0.0172295769141, 0.0392041960688, 0.0569022801117, \ - 0.0579313769517, 0.0396035640187, 0.785754362849])) - - numerical_flux5 = af.transpose(wave_equation.elementFluxIntegral\ - (af.range(gvar.N_Elements), 0)[5]) - - - numerical_flux_sum = numerical_flux0 + numerical_flux1 + numerical_flux2 \ - + numerical_flux3 + numerical_flux4 + numerical_flux5 - - referenceFluxIntegral_sum = referenceFluxIntegral0 + referenceFluxIntegral1\ - + referenceFluxIntegral2 + referenceFluxIntegral3 \ - + referenceFluxIntegral4 + referenceFluxIntegral5 - - assert (af.max(af.abs(numerical_flux_sum - referenceFluxIntegral_sum))\ - < threshold) + -0.002368387579324652, -0.003620502047659841, -0.004320197094090966, + -0.003445512010153811, 0.0176615086879261],\ + [-0.018969769374, -0.00431252844519,-0.00882630935977,-0.0144355176966,\ + -0.019612124119, -0.0209837936827, -0.0154359890788, 0.102576031756], \ + [-0.108222418798, -0.0179274222595, -0.0337807018822, -0.0492589052599,\ + -0.0588472807471, -0.0557970236273, -0.0374764132459, 0.361310165819],\ + [-0.374448714304, -0.0399576371245, -0.0683852285846, -0.0869229749357,\ + -0.0884322503841, -0.0714664112839, -0.0422339853622, 0.771847201979], \ + [-0.785754362849, -0.0396035640187, -0.0579313769517, -0.0569022801117,\ + -0.0392041960688, -0.0172295769141, -0.00337464521455, 1.00000000213],\ + [-1.00000000213, 0.00337464521455, 0.0172295769141, 0.0392041960688,\ + 0.0569022801117, 0.0579313769517, 0.0396035640187, 0.785754362849],\ + [-0.771847201979, 0.0422339853622, 0.0714664112839, 0.0884322503841, \ + 0.0869229749357, 0.0683852285846, 0.0399576371245, 0.374448714304],\ + [-0.361310165819, 0.0374764132459, 0.0557970236273, 0.0588472807471,\ + 0.0492589052599, 0.0337807018822, 0.0179274222595, 0.108222418798], \ + [-0.102576031756, 0.0154359890788, 0.0209837936827, 0.019612124119, \ + 0.0144355176966, 0.00882630935977, 0.00431252844519, 0.018969769374],\ + [-0.0176615086879, 0.00344551201015 ,0.00432019709409, 0.00362050204766,\ + 0.00236838757932, 0.00130167737191, 0.000588597708116, 0.00201663487667]]))) + + numerical_flux = wave_equation.volumeIntegralFlux(gvar.element_nodes, gvar.u[:, :, 0]) + assert (af.max(af.abs(numerical_flux - referenceFluxIntegral)) < threshold) def test_lax_friedrichs_flux(): ''' @@ -352,7 +321,7 @@ def test_lax_friedrichs_flux(): f_i = wave_equation.lax_friedrichs_flux(0) #The lax friedrichs flux at timestep 0 should just be a list of the #amplitude at element boundaries. - analytical_lax_friedrichs_flux = gvar.u[:, -1, 0] + analytical_lax_friedrichs_flux = gvar.u[-1, :, 0] assert af.max(af.abs(analytical_lax_friedrichs_flux - f_i)) < threshold @@ -363,8 +332,8 @@ def test_surface_term(): ''' threshold = 1e-13 gvar.populateGlobalVariables(8, 10) - analytical_f_i = af.transpose(gvar.u[:, -1, 0]) - analytical_f_i_minus1 = af.transpose(af.shift(gvar.u[:, -1, 0], 1)) + analytical_f_i = (gvar.u[-1, :, 0]) + analytical_f_i_minus1 = (af.shift(gvar.u[-1, :, 0], 0, 1)) L_p_1 = af.constant(0, gvar.N_LGL, dtype = af.Dtype.f64) L_p_1[gvar.N_LGL - 1] = 1 @@ -375,7 +344,7 @@ def test_surface_term(): analytical_surface_term = af.blas.matmul(L_p_1, analytical_f_i)\ - af.blas.matmul(L_p_minus1, analytical_f_i_minus1) - numerical_surface_term = af.transpose(wave_equation.surface_term(0)) + numerical_surface_term = (wave_equation.surface_term(0)) assert af.max(af.abs(analytical_surface_term - numerical_surface_term)) \ < threshold return analytical_surface_term @@ -389,13 +358,11 @@ def test_b_vector(): threshold = 1e-13 gvar.populateGlobalVariables(8, 10) - u_n_A_matrix = af.blas.matmul(gvar.u[:, :, 0],\ - wave_equation.A_matrix()) - volume_integral_flux = wave_equation.elementFluxIntegral\ - (af.range(gvar.N_Elements), 0) + u_n_A_matrix = af.blas.matmul(wave_equation.A_matrix(), gvar.u[:, :, 0]) + volume_integral_flux = wave_equation.volumeIntegralFlux(gvar.element_nodes, gvar.u[:, :, 0]) surface_term = test_surface_term() b_vector_analytical = u_n_A_matrix + (volume_integral_flux -\ - af.transpose(surface_term)) * gvar.delta_t + (surface_term)) * gvar.delta_t b_vector_array = wave_equation.b_vector(0) assert (b_vector_analytical - b_vector_array) < threshold From 0011aaf4f411560c378de7105c1d97bbfd126e4e Mon Sep 17 00:00:00 2001 From: Balavarun5 Date: Tue, 8 Aug 2017 23:02:39 +0530 Subject: [PATCH 07/19] Resolving wave speed abnormality. --- code/app/global_variables.py | 35 ++++++++++++++++---------- code/app/wave_equation.py | 48 +++++++++++++++++------------------- code/main.py | 4 +-- 3 files changed, 46 insertions(+), 41 deletions(-) diff --git a/code/app/global_variables.py b/code/app/global_variables.py index afe8907..2556532 100644 --- a/code/app/global_variables.py +++ b/code/app/global_variables.py @@ -46,7 +46,7 @@ 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.])) +x_nodes = af.interop.np_to_af_array(np.array([-2., 2.])) N_LGL = 16 xi_LGL = None lBasisArray = None @@ -86,33 +86,42 @@ def populateGlobalVariables(Number_of_LGL_pts = 8, Number_of_elements = 10): lobatto_weight_function(N_LGL, xi_LGL)) global N_Elements - global element_nodes + global element_nodes #[TODO] Change nodes notation. + global elementMeshNodes + 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) + 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)) + elementMeshNodes = utils.linspace(af.sum(x_nodes[0]), + af.sum(x_nodes[1]), + N_Elements + 1) + + + 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 c global c_lax global delta_t - c = 4.0 + c = 1.0 delta_x = af.min((element_nodes - af.shift(element_nodes, 1, 0))[1:, :]) - delta_t = delta_x / (10 * c) - c_lax = delta_x / (2 * delta_t) + delta_t = delta_x / (50 * c) + c_lax = 0 global u global time - time = utils.linspace(0, 20, int(20 / delta_t)) - u_init = np.e ** (-(element_nodes) ** 2 / 0.4 ** 2) - u = af.constant(0, N_LGL, N_Elements, time.shape[0],\ + total_time = 3 + time = utils.linspace(0, total_time, int(total_time / delta_t)) + u_init = np.e ** (-(element_nodes) ** 2 / 0.4 ** 2) + u = af.constant(0, N_LGL, N_Elements, time.shape[0],\ dtype = af.Dtype.f64) u[:, :, 0] = u_init diff --git a/code/app/wave_equation.py b/code/app/wave_equation.py index 60e7716..7892c1e 100644 --- a/code/app/wave_equation.py +++ b/code/app/wave_equation.py @@ -7,6 +7,7 @@ from app import global_variables as gvar from matplotlib import pyplot as plt import pylab as pl +from tqdm import trange plt.rcParams['figure.figsize'] = 9.6, 6. plt.rcParams['figure.dpi'] = 100 @@ -144,11 +145,11 @@ def A_matrix(): using :math: `N_LGL` points. ''' - x_tile = af.transpose(af.tile(gvar.xi_LGL, 1, gvar.N_LGL)) - power = af.flip(af.range(gvar.N_LGL)) - power_tile = af.tile(power, 1, gvar.N_LGL) - x_pow = af.arith.pow(x_tile, power_tile) - lobatto_weights = gvar.lobatto_weights + x_tile = af.transpose(af.tile(gvar.xi_LGL, 1, gvar.N_LGL)) + power = af.flip(af.range(gvar.N_LGL)) + power_tile = af.tile(power, 1, gvar.N_LGL) + x_pow = af.arith.pow(x_tile, power_tile) + lobatto_weights = gvar.lobatto_weights lobatto_weights_tile = af.tile(af.reorder(lobatto_weights, 1, 2, 0),\ gvar.N_LGL, gvar.N_LGL) @@ -158,11 +159,11 @@ def A_matrix(): L_p = af.reorder(L_i, 0, 2, 1) L_i = af.reorder(L_i, 2, 0, 1) - dx_dxi = dx_dxi_numerical((gvar.x_nodes),gvar.xi_LGL) + dx_dxi = dx_dxi_numerical((gvar.elementMeshNodes[0 : 2]), gvar.xi_LGL) dx_dxi_tile = af.tile(dx_dxi, 1, gvar.N_LGL, gvar.N_LGL) - Li_Lp_array = Li_Lp_xi(L_p, L_i) - L_element = (Li_Lp_array * lobatto_weights_tile * dx_dxi_tile) - A_matrix = af.sum(L_element, dim = 2) + Li_Lp_array = Li_Lp_xi(L_p, L_i) + L_element = (Li_Lp_array * lobatto_weights_tile * dx_dxi_tile) + A_matrix = af.sum(L_element, dim = 2) return A_matrix @@ -241,7 +242,7 @@ def lax_friedrichs_flux(t_n): flux_i_N_LGL = flux_x(u_i_N_LGL) lax_friedrichs_flux = (flux_iplus1_0 + flux_i_N_LGL) / 2 \ - - gvar.c_lax * (u_iplus1_0 - u_i_N_LGL) + - gvar.c_lax * (u_iplus1_0 - u_i_N_LGL) / 2 return lax_friedrichs_flux @@ -307,7 +308,7 @@ def b_vector(t_n): def time_evolution(): ''' - Function which solves the wave equation + Function which solves the wave equation :math: `u^{t_n + 1} = b(t_n) \\times A` iterated over time steps t_n and then plots :math: `x` against the amplitude of the wave. The images are then stored in Wave folder. @@ -317,15 +318,16 @@ def time_evolution(): element_nodes = gvar.element_nodes delta_t = gvar.delta_t - for t_n in range(0, gvar.time.shape[0] - 1): + for t_n in trange(0, gvar.time.shape[0] - 1): + + gvar.u[:, :, t_n + 1] = gvar.u[:, :, t_n] + af.blas.matmul(A_inverse,\ + b_vector(t_n)) + - gvar.u[:, :, t_n + 1] = gvar.u[:, :, t_n] + af.blas.matmul(A_inverse, b_vector(t_n)) - if (t_n % 100 == 0): - print('timestep', t_n) - print('u calculated!') + print('u calculated!') - for t_n in range(0, gvar.time.shape[0] - 1): + for t_n in trange(0, gvar.time.shape[0]): if t_n % 100 == 0: @@ -336,14 +338,8 @@ def time_evolution(): plt.plot(x, y) plt.xlabel('x') plt.ylabel('Amplitude') - plt.title('Time = %f' %(t_n * delta_t)) + plt.title('Time = %f' % (t_n * delta_t)) fig.savefig('1D_Wave_images/%04d' %(t_n / 100) + '.png') plt.close('all') - - print(t_n) - - - - - - return + + return diff --git a/code/main.py b/code/main.py index 9678c4e..79211ac 100644 --- a/code/main.py +++ b/code/main.py @@ -12,7 +12,7 @@ if __name__ == '__main__': - gvar.populateGlobalVariables(8) + gvar.populateGlobalVariables(16) print(gvar.time.shape[0]) - print(wave_equation.volumeIntegralFlux(gvar.element_nodes, gvar.u[:, :, 0])) print(wave_equation.time_evolution()) + #print(gvar.element_nodes) From 64f414c78d483e55870c998de733815599680677 Mon Sep 17 00:00:00 2001 From: Balavarun5 Date: Fri, 11 Aug 2017 16:19:23 +0530 Subject: [PATCH 08/19] -Wave is moving at the expected speed by taking c_lax to be the wave speed. -testwaveEqn has issues accessing modules. --- code/app/__init__.py | 0 code/app/global_variables.py | 27 ++++++++++-------- code/app/wave_equation.py | 51 +++++++++++++++++++++------------- code/main.py | 6 ++-- code/unit_test/__init__.py | 0 code/unit_test/test_waveEqn.py | 47 +++++++++++++++++++++---------- code/utils/utils.py | 2 ++ 7 files changed, 83 insertions(+), 50 deletions(-) delete mode 100644 code/app/__init__.py delete mode 100644 code/unit_test/__init__.py diff --git a/code/app/__init__.py b/code/app/__init__.py deleted file mode 100644 index e69de29..0000000 diff --git a/code/app/global_variables.py b/code/app/global_variables.py index 2556532..d1e49a6 100644 --- a/code/app/global_variables.py +++ b/code/app/global_variables.py @@ -1,5 +1,7 @@ import numpy as np import arrayfire as af +af.set_backend('opencl') +af.set_device(1) from scipy import special as sp from app import lagrange from app import wave_equation @@ -46,13 +48,13 @@ 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([-2., 2.])) +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 +element_LGL = None u = None time = None c = None @@ -86,7 +88,7 @@ def populateGlobalVariables(Number_of_LGL_pts = 8, Number_of_elements = 10): lobatto_weight_function(N_LGL, xi_LGL)) global N_Elements - global element_nodes #[TODO] Change nodes notation. + global element_LGL global elementMeshNodes N_Elements = Number_of_elements @@ -105,27 +107,30 @@ def populateGlobalVariables(Number_of_LGL_pts = 8, Number_of_elements = 10): element_array = af.transpose(af.interop.np_to_af_array(np_element_array)) - element_nodes = wave_equation.mappingXiToX(\ + element_LGL = wave_equation.mappingXiToX(\ af.transpose(element_array), xi_LGL) global c global c_lax global delta_t - c = 1.0 - delta_x = af.min((element_nodes - af.shift(element_nodes, 1, 0))[1:, :]) - delta_t = delta_x / (50 * c) - c_lax = 0 + c = 4 + delta_x = af.min((element_LGL - af.shift(element_LGL, 1, 0))[1:, :]) + delta_t = delta_x / (80 * c) + c_lax = c / 2 # Was previously taken to be 0.1 even works if it's + # taken to be c. + global u global time - total_time = 3 + total_time = 1.05 time = utils.linspace(0, total_time, int(total_time / delta_t)) - u_init = np.e ** (-(element_nodes) ** 2 / 0.4 ** 2) + u_init = np.e ** (-(element_LGL) ** 2 / 0.4 ** 2) 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() @@ -139,7 +144,7 @@ def lobatto_weight_function(n, x): and points :math: `x`. :math:: - w_{n} = \\frac{2 P(x)^2}{n (n - 1)}, + `w_{n} = \\frac{2 P(x)^2}{n (n - 1)}`, Where P(x) is $ (n - 1)^th $ index. Parameters diff --git a/code/app/wave_equation.py b/code/app/wave_equation.py index 7892c1e..c80ed1e 100644 --- a/code/app/wave_equation.py +++ b/code/app/wave_equation.py @@ -1,6 +1,8 @@ #! /usr/bin/env python3 import arrayfire as af +af.set_backend('opencl') +af.set_device(1) import numpy as np from app import lagrange from utils import utils @@ -143,6 +145,11 @@ def A_matrix(): The value of integral of product of lagrange basis functions obtained by LGL points, using Gauss-Lobatto quadrature method using :math: `N_LGL` points. + + [NOTE]: + + The A matrix will vary for each element. The one calculatedis for the case + of 1D elements which are of equal size. ''' x_tile = af.transpose(af.tile(gvar.xi_LGL, 1, gvar.N_LGL)) @@ -186,7 +193,7 @@ def flux_x(u): return gvar.c * u -def volumeIntegralFlux(element_nodes, u): +def volumeIntegralFlux(element_LGL, u): ''' A function to calculate the volume integral of flux in the wave equation. :math:`\\int_{-1}^1 f(u) \\frac{d L_p}{d\\xi} d\\xi` @@ -196,11 +203,11 @@ def volumeIntegralFlux(element_nodes, u): Parameters ---------- - element_nodes : arrayfire.Array [N 1 1 1] - A 1-D array consisting of the LGL nodes mapped onto the + element_LGL : arrayfire.Array [N_LGL N_Elements 1 1] + A 2-D array consisting of the LGL nodes mapped onto the element's domain. - u : arrayfire.Array [1 N 1 1] + u : arrayfire.Array [N_LGL N_Elements 1 1] A 1-D array containing the value of the wave function at the mapped LGL nodes in the element. @@ -220,7 +227,7 @@ def volumeIntegralFlux(element_nodes, u): return flux_integral -def lax_friedrichs_flux(t_n): +def laxFriedrichsFlux(t_n): ''' A function which calculates the lax-friedrichs_flux :math:`f_i` using. :math:`f_i = \\frac{F(u^{i + 1}_0) + F(u^i_{N_{LGL} - 1})}{2} - \frac @@ -241,17 +248,17 @@ def lax_friedrichs_flux(t_n): flux_iplus1_0 = flux_x(u_iplus1_0) flux_i_N_LGL = flux_x(u_i_N_LGL) - lax_friedrichs_flux = (flux_iplus1_0 + flux_i_N_LGL) / 2 \ + laxFriedrichsFlux = (flux_iplus1_0 + flux_i_N_LGL) / 2 \ - gvar.c_lax * (u_iplus1_0 - u_i_N_LGL) / 2 - return lax_friedrichs_flux + return laxFriedrichsFlux def surface_term(t_n): ''' A function which is used to calculate the surface term, :math:`L_p (1) f_i - L_p (-1) f_{i - 1}` - using the lax_friedrichs_flux function and the dLp_xi_LGL function in gvar + using the laxFriedrichsFlux function and the dLp_xi_LGL function in gvar module. Parameters @@ -274,11 +281,10 @@ def surface_term(t_n): https://cocalc.com/projects/1b7f404c-87ba-40d0-816c-2eba17466aa8/files\ /PM\_2\_5/wave\_equation/documents/surface\_term/surface\_term.pdf - ''' L_p_minus1 = gvar.lagrange_basis_function()[:, 0] L_p_1 = gvar.lagrange_basis_function()[:, -1] - f_i = lax_friedrichs_flux(t_n) + f_i = laxFriedrichsFlux(t_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) @@ -296,9 +302,9 @@ def b_vector(t_n): Returns ------- - b_vector_array : arrayfire.array + b_vector_array : arrayfire.Array ''' - volume_integral = volumeIntegralFlux(gvar.element_nodes, gvar.u[:, :, t_n]) + volume_integral = volumeIntegralFlux(gvar.element_LGL, gvar.u[:, :, t_n]) surfaceTerm = surface_term(t_n) b_vector_array = gvar.delta_t * (volume_integral - surfaceTerm) @@ -314,25 +320,30 @@ def time_evolution(): of the wave. The images are then stored in Wave folder. ''' - A_inverse = af.transpose(af.lapack.inverse(A_matrix())) - element_nodes = gvar.element_nodes - delta_t = gvar.delta_t + A_inverse = af.lapack.inverse(A_matrix()) + element_LGL = gvar.element_LGL + delta_t = gvar.delta_t + for t_n in trange(0, gvar.time.shape[0] - 1): gvar.u[:, :, t_n + 1] = gvar.u[:, :, t_n] + af.blas.matmul(A_inverse,\ - b_vector(t_n)) - - + b_vector(t_n)) print('u calculated!') - for t_n in trange(0, gvar.time.shape[0]): + approximate_1_s = (int(1 / gvar.delta_t) * gvar.delta_t) + analytical_u_after_1s = np.e ** (-(gvar.element_LGL - gvar.c * (1 - approximate_1_s)) ** 2 / 0.4 ** 2) + + af.display(analytical_u_after_1s, 10) + af.display(gvar.u[:, :, int(1 / gvar.delta_t)], 10) + af.display(gvar.u[:, :, 0], 10) + for t_n in trange(0, gvar.time.shape[0] - 1): if t_n % 100 == 0: fig = plt.figure() - x = gvar.element_nodes + x = gvar.element_LGL y = gvar.u[:, :, t_n] plt.plot(x, y) diff --git a/code/main.py b/code/main.py index 79211ac..a7e52e6 100644 --- a/code/main.py +++ b/code/main.py @@ -5,14 +5,12 @@ af.set_device(1) import numpy as np from app import global_variables as gvar -from unit_test import test_waveEqn +import test_waveEqn from app import wave_equation from app import lagrange if __name__ == '__main__': - gvar.populateGlobalVariables(16) - print(gvar.time.shape[0]) + gvar.populateGlobalVariables(8) print(wave_equation.time_evolution()) - #print(gvar.element_nodes) diff --git a/code/unit_test/__init__.py b/code/unit_test/__init__.py deleted file mode 100644 index e69de29..0000000 diff --git a/code/unit_test/test_waveEqn.py b/code/unit_test/test_waveEqn.py index e2b0ae6..c4f1362 100644 --- a/code/unit_test/test_waveEqn.py +++ b/code/unit_test/test_waveEqn.py @@ -1,12 +1,13 @@ import numpy as np import arrayfire as af -af.set_backend('opencl') -from app import lagrange -from app import global_variables as gvar -from app import wave_equation -from matplotlib import pyplot as plt -from utils import utils import math +from matplotlib import pyplot as plt +from code.app.lagrange import lagrange +from code.app.global_variables import global_variables as gvar +from code.app.wave_equation import wave_equation +from utils.utils import utils +af.set_backend('opencl') +af.set_device(1) def test_mappingXiToX(): @@ -186,7 +187,7 @@ def test_Integral_Li_Lp(): Obtaining the A_matrix from the function and setting the value of all elements above a certain threshold to be 1 and plotting it. ''' - threshold = 1e-5 + threshold = 1e-5 gvar.populateGlobalVariables(8, 8) A_matrix_structure = np.zeros([gvar.N_LGL, gvar.N_LGL]) non_zero_indices = np.where(np.array(wave_equation.A_matrix()) > threshold) @@ -214,14 +215,18 @@ def test_A_matrix(): ''' threshold = 1e-8 - gvar.populateGlobalVariables(8, 8) + gvar.populateGlobalVariables(8) reference_A_matrix = af.tile(gvar.lobatto_weights, 1, gvar.N_LGL)\ - * af.identity(gvar.N_LGL, gvar.N_LGL, dtype = af.Dtype.f64) + * af.identity(gvar.N_LGL, gvar.N_LGL, dtype = af.Dtype.f64)\ + * af.mean(wave_equation.dx_dxi_numerical((gvar.elementMeshNodes[0 : 2])\ + ,gvar.xi_LGL)) test_A_matrix = wave_equation.A_matrix() error_array = af.abs(reference_A_matrix - test_A_matrix) + print(test_A_matrix, reference_A_matrix) + assert af.algorithm.max(error_array) < threshold @@ -237,7 +242,10 @@ def test_dLp_xi(): /PM_2_5/wave_equation/worksheets/dLp_xi.sagews ''' threshold = 1e-13 + gvar.x_nodes = af.interop.np_to_af_array(np.array([-1., 1.])) gvar.populateGlobalVariables(8) + gvar.c = 1 + reference_d_Lp_xi = af.interop.np_to_af_array(np.array([\ [-14.0000000000226,-3.20991570302344,0.792476681323880,-0.372150435728984,\ @@ -283,6 +291,7 @@ def test_volume_integral_flux(): ''' threshold = 4 * 1e-8 gvar.populateGlobalVariables(8) + gvar.c = 1 referenceFluxIntegral = af.transpose(af.interop.np_to_af_array(np.array([ [-0.002016634876668093, -0.000588597708116113, -0.0013016773719126333,\ @@ -307,18 +316,19 @@ def test_volume_integral_flux(): [-0.0176615086879, 0.00344551201015 ,0.00432019709409, 0.00362050204766,\ 0.00236838757932, 0.00130167737191, 0.000588597708116, 0.00201663487667]]))) - numerical_flux = wave_equation.volumeIntegralFlux(gvar.element_nodes, gvar.u[:, :, 0]) + numerical_flux = wave_equation.volumeIntegralFlux(gvar.element_LGL, gvar.u[:, :, 0]) assert (af.max(af.abs(numerical_flux - referenceFluxIntegral)) < threshold) def test_lax_friedrichs_flux(): ''' - A test function to test the lax_friedrichs_flux function in wave_equation + A test function to test the laxFriedrichsFlux function in wave_equation module. ''' threshold = 1e-14 gvar.populateGlobalVariables(8, 10) + gvar.c = 1 - f_i = wave_equation.lax_friedrichs_flux(0) + f_i = wave_equation.laxFriedrichsFlux(0) #The lax friedrichs flux at timestep 0 should just be a list of the #amplitude at element boundaries. analytical_lax_friedrichs_flux = gvar.u[-1, :, 0] @@ -332,6 +342,9 @@ def test_surface_term(): ''' threshold = 1e-13 gvar.populateGlobalVariables(8, 10) + gvar.c = 1 + + analytical_f_i = (gvar.u[-1, :, 0]) analytical_f_i_minus1 = (af.shift(gvar.u[-1, :, 0], 0, 1)) @@ -353,16 +366,20 @@ def test_surface_term(): def test_b_vector(): ''' A test function to check the b vector obtained analytically and compare it - with the one returned by b_vector function in wave_equation module + with the one returned by b_vector function in wave_equation module ''' threshold = 1e-13 - gvar.populateGlobalVariables(8, 10) + gvar.populateGlobalVariables(8) + gvar.c = 1 u_n_A_matrix = af.blas.matmul(wave_equation.A_matrix(), gvar.u[:, :, 0]) - volume_integral_flux = wave_equation.volumeIntegralFlux(gvar.element_nodes, gvar.u[:, :, 0]) + volume_integral_flux = wave_equation.volumeIntegralFlux(gvar.element_LGL, gvar.u[:, :, 0]) surface_term = test_surface_term() b_vector_analytical = u_n_A_matrix + (volume_integral_flux -\ (surface_term)) * gvar.delta_t b_vector_array = wave_equation.b_vector(0) assert (b_vector_analytical - b_vector_array) < threshold + + +#def test_time_evo diff --git a/code/utils/utils.py b/code/utils/utils.py index 72dc08b..af31bc8 100644 --- a/code/utils/utils.py +++ b/code/utils/utils.py @@ -1,4 +1,6 @@ import arrayfire as af +af.set_backend('opencl') +af.set_device(1) def add(a, b): ''' From 7507b3d694934c08efcde25121a721b421fb064e Mon Sep 17 00:00:00 2001 From: Balavarun5 Date: Fri, 11 Aug 2017 16:27:56 +0530 Subject: [PATCH 09/19] Unit tests working. --- code/unit_test/test_waveEqn.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/code/unit_test/test_waveEqn.py b/code/unit_test/test_waveEqn.py index c4f1362..7298b0b 100644 --- a/code/unit_test/test_waveEqn.py +++ b/code/unit_test/test_waveEqn.py @@ -2,10 +2,10 @@ import arrayfire as af import math from matplotlib import pyplot as plt -from code.app.lagrange import lagrange -from code.app.global_variables import global_variables as gvar -from code.app.wave_equation import wave_equation -from utils.utils import utils +from app import lagrange +from app import global_variables as gvar +from app import wave_equation +from utils import utils af.set_backend('opencl') af.set_device(1) From 756c8f6eeb4916c6d0f408c0b8793c7e6f684d24 Mon Sep 17 00:00:00 2001 From: Balavarun5 Date: Fri, 11 Aug 2017 16:29:13 +0530 Subject: [PATCH 10/19] Unit tests working. --- code/main.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/code/main.py b/code/main.py index a7e52e6..588ffd4 100644 --- a/code/main.py +++ b/code/main.py @@ -5,7 +5,7 @@ af.set_device(1) import numpy as np from app import global_variables as gvar -import test_waveEqn +from unit_test import test_waveEqn from app import wave_equation from app import lagrange From bdffc6cec57abb0ebe510111d1cdc0628b78dbe0 Mon Sep 17 00:00:00 2001 From: Balavarun5 Date: Sat, 26 Aug 2017 16:58:35 +0530 Subject: [PATCH 11/19] Edits in the documentation. --- code/app/global_variables.py | 19 +++++++++--------- code/app/wave_equation.py | 35 ++++++++++++++++------------------ code/main.py | 5 +++-- code/unit_test/test_waveEqn.py | 22 +++++++++------------ 4 files changed, 37 insertions(+), 44 deletions(-) diff --git a/code/app/global_variables.py b/code/app/global_variables.py index d1e49a6..52cb363 100644 --- a/code/app/global_variables.py +++ b/code/app/global_variables.py @@ -6,7 +6,7 @@ from app import lagrange from app import wave_equation from utils import utils - +import math LGL_list = [ \ [-1.0,1.0], \ @@ -84,7 +84,8 @@ def populateGlobalVariables(Number_of_LGL_pts = 8, Number_of_elements = 10): lBasisArray = af.interop.np_to_af_array( \ lagrange.lagrange_basis_coeffs(xi_LGL)) - lobatto_weights = af.interop.np_to_af_array(\ + lobatto_weights = af.interop.\ + np_to_af_array(\ lobatto_weight_function(N_LGL, xi_LGL)) global N_Elements @@ -113,24 +114,22 @@ def populateGlobalVariables(Number_of_LGL_pts = 8, Number_of_elements = 10): global c global c_lax global delta_t - c = 4 + c = 1 delta_x = af.min((element_LGL - af.shift(element_LGL, 1, 0))[1:, :]) - delta_t = delta_x / (80 * c) - c_lax = c / 2 # Was previously taken to be 0.1 even works if it's - # taken to be c. - + delta_t = delta_x / (20 * c) + c_lax = c global u global time - total_time = 1.05 + total_time = 5 time = utils.linspace(0, total_time, int(total_time / delta_t)) u_init = np.e ** (-(element_LGL) ** 2 / 0.4 ** 2) 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() @@ -193,7 +192,7 @@ def lagrange_basis_function(): 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 + as a 2D array of :math: `L_i' (\\xi_{LGL})`. Where i varies along rows and the nodes vary along the columns. Returns diff --git a/code/app/wave_equation.py b/code/app/wave_equation.py index c80ed1e..adcea6d 100644 --- a/code/app/wave_equation.py +++ b/code/app/wave_equation.py @@ -41,6 +41,9 @@ def Li_Lp_xi(L_i_xi, L_p_xi): ''' + A function which broadcasts and multiplies two 2D arrays of + shapes [1 N N 1] and [N 1 N 1] + Parameters ---------- L_i_xi : arrayfire.Array [1 N N 1] @@ -126,12 +129,14 @@ def dx_dxi_analytical(x_nodes, xi): Returns ------- - (x_nodes[1] - x_nodes[0]) / 2 : arrayfire.Array - The analytical solution of - \\frac{dx}{d\\xi} for an element. + dx_xi : arrayfire.Array + The analytical solution of + \\frac{dx}{d\\xi} for an element. ''' - return (x_nodes[1] - x_nodes[0]) / 2 + dx_dxi = (x_nodes[1] - x_nodes[0]) / 2 + + return dx_dxi def A_matrix(): @@ -148,7 +153,7 @@ def A_matrix(): [NOTE]: - The A matrix will vary for each element. The one calculatedis for the case + The A matrix will vary for each element. The one calculated is for the case of 1D elements which are of equal size. ''' @@ -172,6 +177,7 @@ def A_matrix(): L_element = (Li_Lp_array * lobatto_weights_tile * dx_dxi_tile) A_matrix = af.sum(L_element, dim = 2) + return A_matrix @@ -193,7 +199,7 @@ def flux_x(u): return gvar.c * u -def volumeIntegralFlux(element_LGL, u): +def volumeIntegralFlux(u): ''' A function to calculate the volume integral of flux in the wave equation. :math:`\\int_{-1}^1 f(u) \\frac{d L_p}{d\\xi} d\\xi` @@ -203,17 +209,13 @@ def volumeIntegralFlux(element_LGL, u): Parameters ---------- - element_LGL : arrayfire.Array [N_LGL N_Elements 1 1] - A 2-D array consisting of the LGL nodes mapped onto the - element's domain. - u : arrayfire.Array [N_LGL N_Elements 1 1] A 1-D array containing the value of the wave function at the mapped LGL nodes in the element. Returns ------- - flux_integral : arrayfire.Array [1 N 1 1] + flux_integral : arrayfire.Array [N_LGL N_Elements 1 1] A 1-D array of the value of the flux integral calculated for various lagrange basis functions. ''' @@ -251,6 +253,7 @@ def laxFriedrichsFlux(t_n): laxFriedrichsFlux = (flux_iplus1_0 + flux_i_N_LGL) / 2 \ - gvar.c_lax * (u_iplus1_0 - u_i_N_LGL) / 2 + return laxFriedrichsFlux @@ -304,7 +307,7 @@ def b_vector(t_n): ------- b_vector_array : arrayfire.Array ''' - volume_integral = volumeIntegralFlux(gvar.element_LGL, gvar.u[:, :, t_n]) + volume_integral = volumeIntegralFlux(gvar.u[:, :, t_n]) surfaceTerm = surface_term(t_n) b_vector_array = gvar.delta_t * (volume_integral - surfaceTerm) @@ -320,7 +323,7 @@ def time_evolution(): of the wave. The images are then stored in Wave folder. ''' - A_inverse = af.lapack.inverse(A_matrix()) + A_inverse = af.inverse(A_matrix()) element_LGL = gvar.element_LGL delta_t = gvar.delta_t @@ -332,12 +335,6 @@ def time_evolution(): print('u calculated!') - approximate_1_s = (int(1 / gvar.delta_t) * gvar.delta_t) - analytical_u_after_1s = np.e ** (-(gvar.element_LGL - gvar.c * (1 - approximate_1_s)) ** 2 / 0.4 ** 2) - - af.display(analytical_u_after_1s, 10) - af.display(gvar.u[:, :, int(1 / gvar.delta_t)], 10) - af.display(gvar.u[:, :, 0], 10) for t_n in trange(0, gvar.time.shape[0] - 1): if t_n % 100 == 0: diff --git a/code/main.py b/code/main.py index 588ffd4..d076b55 100644 --- a/code/main.py +++ b/code/main.py @@ -12,5 +12,6 @@ if __name__ == '__main__': - gvar.populateGlobalVariables(8) - print(wave_equation.time_evolution()) + gvar.populateGlobalVariables(8,10) + wave_equation.time_evolution() + print(wave_equation.laxFriedrichsFlux(0)) diff --git a/code/unit_test/test_waveEqn.py b/code/unit_test/test_waveEqn.py index 7298b0b..ac8dfa6 100644 --- a/code/unit_test/test_waveEqn.py +++ b/code/unit_test/test_waveEqn.py @@ -33,7 +33,7 @@ def test_Li_Lp_xi(): ''' A Test function to check the Li_Lp_xi function in wave_equation module by passing two test arrays and comparing the analytical product with the - numerically calculated one with a tolerance of 1e-14. + numerically calculated one with to a tolerance of 1e-14. ''' gvar.populateGlobalVariables(3) @@ -66,7 +66,7 @@ def test_dx_dxi(): passing nodes of an element and using the LGL points. Analytically, the differential would be a constant. The check has a tolerance 1e-7. ''' - threshold = 1e-7 + threshold = 1e-14 gvar.populateGlobalVariables(8) nodes = np.array([7, 10], dtype = np.float64) test_nodes = af.interop.np_to_af_array(nodes) @@ -214,7 +214,7 @@ def test_A_matrix(): The A matrix calculated analytically gives a different matrix. ''' - threshold = 1e-8 + threshold = 2e-10 gvar.populateGlobalVariables(8) reference_A_matrix = af.tile(gvar.lobatto_weights, 1, gvar.N_LGL)\ @@ -225,7 +225,7 @@ def test_A_matrix(): test_A_matrix = wave_equation.A_matrix() error_array = af.abs(reference_A_matrix - test_A_matrix) - print(test_A_matrix, reference_A_matrix) + print(af.max(error_array)) assert af.algorithm.max(error_array) < threshold @@ -289,7 +289,7 @@ def test_volume_integral_flux(): /PM_2_5/wave_equation/worksheets/volume_integral_flux.sagews ''' - threshold = 4 * 1e-8 + threshold = 4e-8 gvar.populateGlobalVariables(8) gvar.c = 1 @@ -314,9 +314,10 @@ def test_volume_integral_flux(): [-0.102576031756, 0.0154359890788, 0.0209837936827, 0.019612124119, \ 0.0144355176966, 0.00882630935977, 0.00431252844519, 0.018969769374],\ [-0.0176615086879, 0.00344551201015 ,0.00432019709409, 0.00362050204766,\ - 0.00236838757932, 0.00130167737191, 0.000588597708116, 0.00201663487667]]))) + 0.00236838757932, 0.00130167737191, 0.000588597708116, 0.00201663487667\ + ]]))) - numerical_flux = wave_equation.volumeIntegralFlux(gvar.element_LGL, gvar.u[:, :, 0]) + numerical_flux = wave_equation.volumeIntegralFlux(gvar.u[:, :, 0]) assert (af.max(af.abs(numerical_flux - referenceFluxIntegral)) < threshold) def test_lax_friedrichs_flux(): @@ -329,8 +330,6 @@ def test_lax_friedrichs_flux(): gvar.c = 1 f_i = wave_equation.laxFriedrichsFlux(0) - #The lax friedrichs flux at timestep 0 should just be a list of the - #amplitude at element boundaries. analytical_lax_friedrichs_flux = gvar.u[-1, :, 0] assert af.max(af.abs(analytical_lax_friedrichs_flux - f_i)) < threshold @@ -373,13 +372,10 @@ def test_b_vector(): gvar.c = 1 u_n_A_matrix = af.blas.matmul(wave_equation.A_matrix(), gvar.u[:, :, 0]) - volume_integral_flux = wave_equation.volumeIntegralFlux(gvar.element_LGL, gvar.u[:, :, 0]) + volume_integral_flux = wave_equation.volumeIntegralFlux(gvar.u[:, :, 0]) surface_term = test_surface_term() b_vector_analytical = u_n_A_matrix + (volume_integral_flux -\ (surface_term)) * gvar.delta_t b_vector_array = wave_equation.b_vector(0) assert (b_vector_analytical - b_vector_array) < threshold - - -#def test_time_evo From 88b04d0d65885428b579b59a5693688d64123aa9 Mon Sep 17 00:00:00 2001 From: Balavarun5 Date: Sat, 26 Aug 2017 17:05:41 +0530 Subject: [PATCH 12/19] Minor edits --- code/app/wave_equation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/code/app/wave_equation.py b/code/app/wave_equation.py index adcea6d..22cfe02 100644 --- a/code/app/wave_equation.py +++ b/code/app/wave_equation.py @@ -233,7 +233,7 @@ def laxFriedrichsFlux(t_n): ''' A function which calculates the lax-friedrichs_flux :math:`f_i` using. :math:`f_i = \\frac{F(u^{i + 1}_0) + F(u^i_{N_{LGL} - 1})}{2} - \frac - {\Delta x}{2\Delta t} (u^{i + 1}_0 - u^i_{N_{LGL} - 1})` + {\Delta x}{2\Delta t} (u^{i + 1}_0 - u^i_{N_{LGL} - 1})` Parameters ---------- From 6c5c3d040905536aa119bb7e48a96b76261721fc Mon Sep 17 00:00:00 2001 From: Balavarun5 Date: Thu, 31 Aug 2017 11:08:30 +0530 Subject: [PATCH 13/19] Documentations edited. Tabs converted to spaces. --- code/app/global_variables.py | 346 ++++++++++------- code/app/lagrange.py | 172 ++++----- code/app/wave_equation.py | 535 +++++++++++++------------- code/main.py | 8 +- code/unit_test/test_waveEqn.py | 666 ++++++++++++++++----------------- code/utils/utils.py | 36 +- 6 files changed, 909 insertions(+), 854 deletions(-) diff --git a/code/app/global_variables.py b/code/app/global_variables.py index 52cb363..eb70cd8 100644 --- a/code/app/global_variables.py +++ b/code/app/global_variables.py @@ -45,8 +45,8 @@ for idx in np.arange(len(LGL_list)): - LGL_list[idx] = np.array(LGL_list[idx], dtype = np.float64) - LGL_list[idx] = af.interop.np_to_af_array(LGL_list[idx]) + 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 @@ -62,153 +62,211 @@ c_lax = None def populateGlobalVariables(Number_of_LGL_pts = 8, Number_of_elements = 10): - ''' - Initialize the global variables. - 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_LGL - global elementMeshNodes - - 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))) - - elementMeshNodes = utils.linspace(af.sum(x_nodes[0]), - af.sum(x_nodes[1]), - N_Elements + 1) - - - element_array = af.transpose(af.interop.np_to_af_array(np_element_array)) - element_LGL = wave_equation.mappingXiToX(\ - af.transpose(element_array), xi_LGL) - - global c - global c_lax - global delta_t - c = 1 - delta_x = af.min((element_LGL - af.shift(element_LGL, 1, 0))[1:, :]) - delta_t = delta_x / (20 * c) - c_lax = c - - global u - global time - total_time = 5 - time = utils.linspace(0, total_time, int(total_time / delta_t)) - u_init = np.e ** (-(element_LGL) ** 2 / 0.4 ** 2) - 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() - - - return + ''' + Time evolution of the wave function requires the use of 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_LGL + global elementMeshNodes + + 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))) + + elementMeshNodes = utils.linspace(af.sum(x_nodes[0]), + af.sum(x_nodes[1]), + N_Elements + 1) + + + element_array = af.transpose(af.interop.np_to_af_array(np_element_array)) + element_LGL = wave_equation.mappingXiToX(\ + af.transpose(element_array), xi_LGL) + + global c + global c_lax + global delta_t + c = 4 + delta_x = af.min((element_LGL - af.shift(element_LGL, 1, 0))[1:, :]) + delta_t = delta_x / (20 * c) + c_lax = 1 + + global u + global time + total_time = 5 + time = utils.linspace(0, total_time, int(total_time / delta_t)) + u_init = np.e ** (-(element_LGL) ** 2 / 0.4 ** 2) + 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() + + + 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. - - .. lobatto weight function - - https://en.wikipedia.org/wiki/Gaussian_quadrature#Gauss.E2.80.93Lobatto_rules - - Returns - ------- - (2 / (n * (n - 1)) / (P(x))**2) : arrayfire.Array - An array of lobatto weight functions for - the given :math: `x` points and index. - - ''' - P = sp.legendre(n - 1) - - return (2 / (n * (n - 1)) / (P(x))**2) + ''' + 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 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 + ''' + 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 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 + ''' + 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 diff --git a/code/app/lagrange.py b/code/app/lagrange.py index 91d41f2..14ef566 100644 --- a/code/app/lagrange.py +++ b/code/app/lagrange.py @@ -5,98 +5,98 @@ from app import global_variables as gvar def LGL_points(N): - ''' - Returns the :math: `N` Legendre-Gauss-Laguere points which are - the roots of the equation - :math:: - (1 - x^2)L'_N = 0 - - Parameters - ---------- - N : int - Number of LGL-points to be generated. - 2 < N < 16 - - Returns - ------- - lgl : arrayfire.Array - An array of :math: `N` LGL points. - ''' - if N > 16 or N < 2: - print('Skipping! This function can only return from ', - '2 to 16 LGL points.') - - n = N - 2 + ''' + Returns the :math: `N` Legendre-Gauss-Laguere points which are + the roots of the equation + :math:: + (1 - x^2)L'_N = 0 + + Parameters + ---------- + N : int + Number of LGL-points to be generated. + 2 < N < 16 + + Returns + ------- + lgl : arrayfire.Array + An array of :math: `N` LGL points. + ''' + if N > 16 or N < 2: + print('Skipping! This function can only return from ', + '2 to 16 LGL points.') + + n = N - 2 - lgl = af.Array(gvar.LGL_list[n]) - return lgl + lgl = af.Array(gvar.LGL_list[n]) + return lgl def lagrange_basis_coeffs(x): - ''' - A function to get the coefficients of the Lagrange basis polynomials for - a given set of x nodes. + ''' + A function to get the coefficients of the Lagrange basis polynomials for + a given set of x nodes. + + This function calculates the Lagrange basis + polynomials by this formula: + :math:: + `L_i = \\prod_{m = 0, m \\notin i}^{N - 1}\\frac{(x - x_m)}{(x_i - x_m)}` - This function calculates the Lagrange basis - polynomials by this formula: - :math:: - `L_i = \\prod_{m = 0, m \\notin i}^{N - 1}\\frac{(x - x_m)}{(x_i - x_m)}` + Parameters + ---------- + x : numpy.array + A numpy array consisting of the :math: `x` nodes using which the + lagrange basis functions need to be evaluated. - Parameters - ---------- - x : numpy.array - A numpy array consisting of the :math: `x` nodes using which the - lagrange basis functions need to be evaluated. - - 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. - ''' - X = np.array(x) - lagrange_basis_poly = np.zeros([X.shape[0], X.shape[0]]) - - for j in np.arange(X.shape[0]): - lagrange_basis_j = np.poly1d([1]) - - for m in np.arange(X.shape[0]): - if m != j: - lagrange_basis_j *= np.poly1d([1, -X[m]]) \ - / (X[j] - X[m]) - lagrange_basis_poly[j] = lagrange_basis_j.c - - return lagrange_basis_poly + 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. + ''' + X = np.array(x) + lagrange_basis_poly = np.zeros([X.shape[0], X.shape[0]]) + + for j in np.arange(X.shape[0]): + lagrange_basis_j = np.poly1d([1]) + + for m in np.arange(X.shape[0]): + if m != j: + lagrange_basis_j *= np.poly1d([1, -X[m]]) \ + / (X[j] - X[m]) + lagrange_basis_poly[j] = lagrange_basis_j.c + + return lagrange_basis_poly def lagrange_basis(i, x): - ''' - Calculates the value of the :math: `i^{th}` Lagrange basis (calculated - using the :math:`\\xi_{LGL}` points) at the :math: `x` coordinates. - - Parameters - ---------- - - i : int - The Lagrange basis which is to be calculated. - - x : arrayfire.Array - The coordinates at which the `i^{th}` Lagrange polynomial is to be - calculated. - - Returns - ------- - - l_xi_j : arrayfire.Array - Array of values of the :math: `i^{th}` Lagrange basis - calculated at the given :math: `x` coordinates. - ''' - x_tile = af.transpose(af.tile(x, 1, gvar.N_LGL)) - power = utils.linspace((gvar.N_LGL-1), 0, gvar.N_LGL) - power_tile = af.tile(power, 1, x.shape[0]) - x_pow = af.arith.pow(x_tile, power_tile) - l_xi_j = af.blas.matmul(gvar.lBasisArray[i], x_pow) - - return l_xi_j + ''' + Calculates the value of the :math: `i^{th}` Lagrange basis (calculated + using the :math:`\\xi_{LGL}` points) at the :math: `x` coordinates. + + Parameters + ---------- + + i : int + The Lagrange basis which is to be calculated. + + x : arrayfire.Array + The coordinates at which the `i^{th}` Lagrange polynomial is to be + calculated. + + Returns + ------- + + l_xi_j : arrayfire.Array + Array of values of the :math: `i^{th}` Lagrange basis + calculated at the given :math: `x` coordinates. + ''' + x_tile = af.transpose(af.tile(x, 1, gvar.N_LGL)) + power = utils.linspace((gvar.N_LGL-1), 0, gvar.N_LGL) + power_tile = af.tile(power, 1, x.shape[0]) + x_pow = af.arith.pow(x_tile, power_tile) + l_xi_j = af.blas.matmul(gvar.lBasisArray[i], x_pow) + + return l_xi_j diff --git a/code/app/wave_equation.py b/code/app/wave_equation.py index 22cfe02..8e55035 100644 --- a/code/app/wave_equation.py +++ b/code/app/wave_equation.py @@ -40,145 +40,141 @@ def Li_Lp_xi(L_i_xi, L_p_xi): - ''' - A function which broadcasts and multiplies two 2D arrays of - shapes [1 N N 1] and [N 1 N 1] - - Parameters - ---------- - L_i_xi : arrayfire.Array [1 N N 1] - A 2D array :math:`L_i` obtained at LGL points calculated at the - LGL nodes :math:`N_LGL`. - - L_p_xi : arrayfire.Array [N 1 N 1] - A 2D array :math:`L_p` obtained at LGL points calculated at the - LGL nodes :math:`N_LGL`. - - Returns - ------- - Li_Lp_xi : arrayfire.Array [N N N 1] - Matrix of :math:`L_p (N_LGL) L_i (N_LGL)`. - - ''' - Li_Lp_xi = af.bcast.broadcast(utils.multiply, L_i_xi, L_p_xi) - + ''' + Broadcasts and multiplies two arrays of different ordering. Used in the + calculation of A matrix. + + Parameters + ---------- + L_i_xi : arrayfire.Array [1 N N 1] + A 2D array :math:`L_i` obtained at LGL points calculated at + the LGL nodes :math:`N_LGL`. + + L_p_xi : arrayfire.Array [N 1 N 1] + A 2D array :math:`L_p` obtained at LGL points calculated at the + LGL nodes :math:`N_LGL`. + + Returns + ------- + Li_Lp_xi : arrayfire.Array [N N N 1] + Matrix of :math:`L_p (N_LGL) L_i (N_LGL)`. + + ''' + Li_Lp_xi = af.bcast.broadcast(utils.multiply, L_i_xi, L_p_xi) + return Li_Lp_xi def mappingXiToX(x_nodes, xi): - ''' - Parameters - ---------- - - x_nodes : arrayfire.Array - Element nodes. - - xi : numpy.float64 - Value of :math: `\\xi`coordinate for which the corresponding - :math: `x` coordinate is to be found. + ''' + Maps points in :math: `\\xi` space to :math:`x` space using the formula + :math: `x = \\frac{1 - \\xi}{2} x_0 + \\frac{1 + \\xi}{2} x_1` + + Parameters + ---------- + + x_nodes : arrayfire.Array + Element nodes. + + xi : numpy.float64 + Value of :math: `\\xi`coordinate for which the corresponding + :math: `x` coordinate is to be found. . - - Returns - ------- - N0_x0 + N1_x1 : arrayfire.Array - :math: `x` value in the element corresponding to - :math:`\\xi`. - ''' - N_0 = (1 - xi) / 2 - N_1 = (1 + xi) / 2 - - N0_x0 = af.bcast.broadcast(utils.multiply, N_0, x_nodes[0]) - N1_x1 = af.bcast.broadcast(utils.multiply, N_1, x_nodes[1]) - - return N0_x0 + N1_x1 + + Returns + ------- + x : arrayfire.Array + :math: `x` value in the element corresponding to :math:`\\xi`. + ''' + N_0 = (1 - xi) / 2 + N_1 = (1 + xi) / 2 + + N0_x0 = af.bcast.broadcast(utils.multiply, N_0, x_nodes[0]) + N1_x1 = af.bcast.broadcast(utils.multiply, N_1, x_nodes[1]) + + x = N0_x0 + N1_x1 + + return x + def dx_dxi_numerical(x_nodes, xi): - ''' - Differential calculated by central differential method about :math: `\\xi` - using the mappingXiToX function. - - Parameters - ---------- - - x_nodes : arrayfire.Array - Contains the nodes of elements - - xi : float - Value of :math: `\\xi` - - Returns - ------- - (x2 - x1) / (2 * dxi) : arrayfire.Array - :math:`\\frac{dx}{d \\xi}`. - ''' - dxi = 1e-7 - x2 = mappingXiToX(x_nodes, xi + dxi) - x1 = mappingXiToX(x_nodes, xi - dxi) - - return (x2 - x1) / (2 * dxi) + ''' + Differential :math: `\\frac{dx}{d \\xi}` calculated by central differential + method about xi using the mappingXiToX function. + + Parameters + ---------- + + x_nodes : arrayfire.Array + Contains the nodes of elements + + xi : float + Value of :math: `\\xi` + + Returns + ------- + dx_dxi : arrayfire.Array + :math:`\\frac{dx}{d \\xi}`. + ''' + dxi = 1e-7 + x2 = mappingXiToX(x_nodes, xi + dxi) + x1 = mappingXiToX(x_nodes, xi - dxi) + + dx_dxi = (x2 - x1) / (2 * dxi) + + return dx_dxi def dx_dxi_analytical(x_nodes, xi): - ''' - - Parameters - ---------- - x_nodes : arrayfire.Array - An array containing the nodes of an element. - - Returns - ------- - dx_xi : arrayfire.Array - The analytical solution of - \\frac{dx}{d\\xi} for an element. - - ''' - dx_dxi = (x_nodes[1] - x_nodes[0]) / 2 - - return dx_dxi + ''' + The analytical result for :math:`\\frac{dx}{d \\xi}` for a 1D element is + :math: `\\frac{x_1 - x_0}{2} + Parameters + ---------- + x_nodes : arrayfire.Array + An array containing the nodes of an element. + + Returns + ------- + analytical_dx_dxi : arrayfire.Array + The analytical solution of + \\frac{dx}{d\\xi} for an element. + + ''' + analytical_dx_dxi = (x_nodes[1] - x_nodes[0]) / 2 + + return analytical_dx_dxi def A_matrix(): + ''' + Calculates A matrix whose elements :math:`A_{p i}` are given by + :math: `A_{p i} &= \\int^1_{-1} L_p(\\xi)L_i(\\xi) \\frac{dx}{d\\xi}` + These elements are to be arranged in an :math:`N \times N` array with p + varying from 0 to N - 1 along the rows and i along the columns. + The integration is carried out using Gauss-Lobatto quadrature. + + Full description of the algorithm can be found here- + https://cocalc.com/projects/1b7f404c-87ba-40d0-816c-2eba17466aa8/files/ + PM_2_5/wave_equation/documents/wave_equation_report + /A_matrix.pdf?session=default + + Returns + ------- + A_matrix : arrayfire.Array + The value of integral of product of lagrange basis functions + obtained by LGL points, using Gauss-Lobatto quadrature method + using :math: `N_LGL` points. ''' - Calculates the value of lagrange basis functions obtained for :math: `N_LGL` - points at the LGL nodes. - - Returns - ------- - A_matrix : arrayfire.Array - The value of integral of product of lagrange basis functions - obtained by LGL points, using Gauss-Lobatto quadrature method - using :math: `N_LGL` points. - - [NOTE]: - - The A matrix will vary for each element. The one calculated is for the case - of 1D elements which are of equal size. - ''' - - x_tile = af.transpose(af.tile(gvar.xi_LGL, 1, gvar.N_LGL)) - power = af.flip(af.range(gvar.N_LGL)) - power_tile = af.tile(power, 1, gvar.N_LGL) - x_pow = af.arith.pow(x_tile, power_tile) - lobatto_weights = gvar.lobatto_weights - - lobatto_weights_tile = af.tile(af.reorder(lobatto_weights, 1, 2, 0),\ - gvar.N_LGL, gvar.N_LGL) - - index = af.range(gvar.N_LGL) - L_i = af.blas.matmul(gvar.lBasisArray[index], x_pow) - L_p = af.reorder(L_i, 0, 2, 1) - L_i = af.reorder(L_i, 2, 0, 1) - - dx_dxi = dx_dxi_numerical((gvar.elementMeshNodes[0 : 2]), gvar.xi_LGL) - dx_dxi_tile = af.tile(dx_dxi, 1, gvar.N_LGL, gvar.N_LGL) - Li_Lp_array = Li_Lp_xi(L_p, L_i) - L_element = (Li_Lp_array * lobatto_weights_tile * dx_dxi_tile) - A_matrix = af.sum(L_element, dim = 2) - - - return A_matrix + dx_dxi = dx_dxi_numerical((gvar.elementMeshNodes[0 : 2]), + gvar.xi_LGL) + dx_dxi_tile = af.tile(dx_dxi, 1, gvar.N_LGL) + identity_matrix = af.identity(gvar.N_LGL, gvar.N_LGL, dtype = af.Dtype.f64) + A_matrix = af.bcast.broadcast(utils.multiply, gvar.lobatto_weights, + identity_matrix) * dx_dxi_tile + + return A_matrix def flux_x(u): @@ -189,165 +185,166 @@ def flux_x(u): Parameters ---------- u : arrayfire.Array - A 1-D array which contains the value of wave function. - - Returns - ------- - c * u : arrayfire.Array - The value of the flux for given u. + A 1-D array which contains the value of wave function. + + Returns + ------- + c * u : arrayfire.Array + The value of the flux for given u. ''' return gvar.c * u def volumeIntegralFlux(u): - ''' - A function to calculate 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 over an element with LGL nodes mapped onto it. - - Parameters - ---------- - u : arrayfire.Array [N_LGL N_Elements 1 1] - A 1-D array containing the value of the wave function at the - mapped LGL nodes in the element. - - Returns - ------- - flux_integral : arrayfire.Array [N_LGL N_Elements 1 1] - A 1-D array of the value of the flux integral calculated - for various lagrange basis functions. - ''' - - dLp_xi = gvar.dLp_xi - weight_tile = af.tile(gvar.lobatto_weights, 1, gvar.N_Elements) - flux = flux_x(u) - weight_flux = weight_tile * flux - flux_integral = af.blas.matmul(dLp_xi, weight_flux) - - return flux_integral + ''' + A function to calculate 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 over an element with LGL nodes mapped onto it. + + Parameters + ---------- + u : arrayfire.Array [N_LGL N_Elements 1 1] + A 1-D array containing the value of the wave function at the + mapped LGL nodes in the element. + + Returns + ------- + flux_integral : arrayfire.Array [N_LGL N_Elements 1 1] + A 1-D array of the value of the flux integral calculated + for various lagrange basis functions. + ''' + + dLp_xi = gvar.dLp_xi + weight_tile = af.transpose(af.tile(gvar.lobatto_weights, 1, gvar.N_LGL)) + dLp_xi *= weight_tile + flux = flux_x(u) + weight_flux = flux + flux_integral = af.blas.matmul(dLp_xi, weight_flux) + + return flux_integral def laxFriedrichsFlux(t_n): - ''' - A function which calculates the lax-friedrichs_flux :math:`f_i` using. - :math:`f_i = \\frac{F(u^{i + 1}_0) + F(u^i_{N_{LGL} - 1})}{2} - \frac - {\Delta x}{2\Delta t} (u^{i + 1}_0 - u^i_{N_{LGL} - 1})` - - Parameters - ---------- - u : 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. - - t_n : float - The timestep at which the lax-friedrichs-flux is to be calculated. - ''' - - u_iplus1_0 = af.shift(gvar.u[0, :, t_n], 0, -1) - u_i_N_LGL = gvar.u[-1, :, t_n] - flux_iplus1_0 = flux_x(u_iplus1_0) - flux_i_N_LGL = flux_x(u_i_N_LGL) - - laxFriedrichsFlux = (flux_iplus1_0 + flux_i_N_LGL) / 2 \ - - gvar.c_lax * (u_iplus1_0 - u_i_N_LGL) / 2 - - - return laxFriedrichsFlux + ''' + A function which calculates the lax-friedrichs_flux :math:`f_i` using. + :math:`f_i = \\frac{F(u^{i + 1}_0) + F(u^i_{N_{LGL} - 1})}{2} - \frac + {\Delta x}{2\Delta t} (u^{i + 1}_0 - u^i_{N_{LGL} - 1})` + + Parameters + ---------- + u : 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. + + t_n : float + The timestep at which the lax-friedrichs-flux is to be calculated. + ''' + + u_iplus1_0 = af.shift(gvar.u[0, :, t_n], 0, -1) + u_i_N_LGL = gvar.u[-1, :, t_n] + flux_iplus1_0 = flux_x(u_iplus1_0) + flux_i_N_LGL = flux_x(u_i_N_LGL) + + laxFriedrichsFlux = (flux_iplus1_0 + flux_i_N_LGL) / 2 \ + - gvar.c_lax * (u_iplus1_0 - u_i_N_LGL) / 2 + + + return laxFriedrichsFlux def surface_term(t_n): - ''' - A function which is used to calculate the surface term, - :math:`L_p (1) f_i - L_p (-1) f_{i - 1}` - using the laxFriedrichsFlux function and the dLp_xi_LGL function in gvar - module. - - Parameters - ---------- - t_n : double - The timestep at which the surface term is to be calculated. - - 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. - - Reference - --------- - Link to PDF describing the algorithm to obtain the surface term. - - https://cocalc.com/projects/1b7f404c-87ba-40d0-816c-2eba17466aa8/files\ - /PM\_2\_5/wave\_equation/documents/surface\_term/surface\_term.pdf - ''' - L_p_minus1 = gvar.lagrange_basis_function()[:, 0] - L_p_1 = gvar.lagrange_basis_function()[:, -1] - f_i = laxFriedrichsFlux(t_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 + ''' + A function which is used to calculate the surface term, + :math:`L_p (1) f_i - L_p (-1) f_{i - 1}` + using the laxFriedrichsFlux function and the dLp_xi_LGL function in gvar + module. + + Parameters + ---------- + t_n : double + The timestep at which the surface term is to be calculated. + + 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. + + Reference + --------- + Link to PDF describing the algorithm to obtain the surface term. + + https://cocalc.com/projects/1b7f404c-87ba-40d0-816c-2eba17466aa8/files\ + /PM\_2\_5/wave\_equation/documents/surface\_term/surface\_term.pdf + ''' + L_p_minus1 = gvar.lagrange_basis_function()[:, 0] + L_p_1 = gvar.lagrange_basis_function()[:, -1] + f_i = laxFriedrichsFlux(t_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 def b_vector(t_n): - ''' - A function which returns the b vector for N_Elements number of elements. - - Parameters - ---------- - t_n : double - - Returns - ------- - b_vector_array : arrayfire.Array - ''' - volume_integral = volumeIntegralFlux(gvar.u[:, :, t_n]) - surfaceTerm = surface_term(t_n) - b_vector_array = gvar.delta_t * (volume_integral - surfaceTerm) - - - return b_vector_array + ''' + A function which returns the b vector for N_Elements number of elements. + + Parameters + ---------- + t_n : double + + Returns + ------- + b_vector_array : arrayfire.Array + ''' + volume_integral = volumeIntegralFlux(gvar.u[:, :, t_n]) + surfaceTerm = surface_term(t_n) + b_vector_array = gvar.delta_t * (volume_integral - surfaceTerm) + + + return b_vector_array def time_evolution(): - ''' - Function which solves the wave equation - :math: `u^{t_n + 1} = b(t_n) \\times A` - iterated over time steps t_n and then plots :math: `x` against the amplitude - of the wave. The images are then stored in Wave folder. - ''' - - A_inverse = af.inverse(A_matrix()) - element_LGL = gvar.element_LGL - delta_t = gvar.delta_t - - - for t_n in trange(0, gvar.time.shape[0] - 1): - - gvar.u[:, :, t_n + 1] = gvar.u[:, :, t_n] + af.blas.matmul(A_inverse,\ - b_vector(t_n)) - - print('u calculated!') - - for t_n in trange(0, gvar.time.shape[0] - 1): - - if t_n % 100 == 0: - - fig = plt.figure() - x = gvar.element_LGL - y = gvar.u[:, :, t_n] - - plt.plot(x, y) - plt.xlabel('x') - plt.ylabel('Amplitude') - plt.title('Time = %f' % (t_n * delta_t)) - fig.savefig('1D_Wave_images/%04d' %(t_n / 100) + '.png') - plt.close('all') - - return + ''' + Function which solves the wave equation + :math: `u^{t_n + 1} = b(t_n) \\times A` + iterated over time steps t_n and then plots :math: `x` against the amplitude + of the wave. The images are then stored in Wave folder. + ''' + + A_inverse = af.inverse(A_matrix()) + element_LGL = gvar.element_LGL + delta_t = gvar.delta_t + + + for t_n in trange(0, gvar.time.shape[0] - 1): + + gvar.u[:, :, t_n + 1] = gvar.u[:, :, t_n] + af.blas.matmul(A_inverse,\ + b_vector(t_n)) + + print('u calculated!') + + for t_n in trange(0, gvar.time.shape[0] - 1): + + if t_n % 100 == 0: + + fig = plt.figure() + x = gvar.element_LGL + y = gvar.u[:, :, t_n] + + plt.plot(x, y) + plt.xlabel('x') + plt.ylabel('Amplitude') + plt.title('Time = %f' % (t_n * delta_t)) + fig.savefig('1D_Wave_images/%04d' %(t_n / 100) + '.png') + plt.close('all') + + return diff --git a/code/main.py b/code/main.py index d076b55..ede5969 100644 --- a/code/main.py +++ b/code/main.py @@ -11,7 +11,7 @@ if __name__ == '__main__': - - gvar.populateGlobalVariables(8,10) - wave_equation.time_evolution() - print(wave_equation.laxFriedrichsFlux(0)) + + gvar.populateGlobalVariables(8,10) + wave_equation.time_evolution() + #print(gvar.lobatto_weights) diff --git a/code/unit_test/test_waveEqn.py b/code/unit_test/test_waveEqn.py index ac8dfa6..7e62fa0 100644 --- a/code/unit_test/test_waveEqn.py +++ b/code/unit_test/test_waveEqn.py @@ -11,371 +11,371 @@ def test_mappingXiToX(): - ''' - A test function to check the mappingXiToX function in wave_equation module, - The test involves passing trial element nodes and :math: `\\xi` and - comparing it with the x obatined by passing the trial parameters to - mappingXiToX function. - ''' - threshold = 1e-14 - gvar.populateGlobalVariables() - - test_element_nodes = af.interop.np_to_af_array(np.array([7, 11])) - test_xi = 0 - analytical_x_value = 9 - numerical_x_value = wave_equation.mappingXiToX(test_element_nodes, test_xi) - - assert af.abs(analytical_x_value - numerical_x_value) <= threshold + ''' + A test function to check the mappingXiToX function in wave_equation module, + The test involves passing trial element nodes and :math: `\\xi` and + comparing it with the x obatined by passing the trial parameters to + mappingXiToX function. + ''' + threshold = 1e-14 + gvar.populateGlobalVariables() + + test_element_nodes = af.interop.np_to_af_array(np.array([7, 11])) + test_xi = 0 + analytical_x_value = 9 + numerical_x_value = wave_equation.mappingXiToX(test_element_nodes, test_xi) + + assert af.abs(analytical_x_value - numerical_x_value) <= threshold def test_Li_Lp_xi(): - ''' - A Test function to check the Li_Lp_xi function in wave_equation module - by passing two test arrays and comparing the analytical product with the - numerically calculated one with to a tolerance of 1e-14. - ''' - - gvar.populateGlobalVariables(3) - - threshold = 1e-14 - - test_array = af.interop.np_to_af_array(np.array([[1, 2, 3], \ - [4, 5, 6], \ - [7, 8, 9]], \ - dtype = np.float64)) - - test_array1 = af.reorder(test_array, 2, 0, 1) - - product = np.zeros([3,3,3]) - product[0] = np.array([[1,2,3],[8,10,12],[21,24,27]]) - product[1] = np.array([[4,8,12],[20,25,30],[42,48,54]]) - product[2] = np.array([[7,14,21],[32,40,48],[63,72,81]]) - - analytical_product = af.interop.np_to_af_array(product) - - check_order3 = af.sum(af.abs(wave_equation.Li_Lp_xi\ - (test_array, test_array1) - analytical_product)) <= threshold - - assert check_order3 + ''' + A Test function to check the Li_Lp_xi function in wave_equation module + by passing two test arrays and comparing the analytical product with the + numerically calculated one with to a tolerance of 1e-14. + ''' + + gvar.populateGlobalVariables(3) + + threshold = 1e-14 + + test_array = af.interop.np_to_af_array(np.array([[1, 2, 3], \ + [4, 5, 6], \ + [7, 8, 9]], \ + dtype = np.float64)) + + test_array1 = af.reorder(test_array, 2, 0, 1) + + product = np.zeros([3,3,3]) + product[0] = np.array([[1,2,3],[8,10,12],[21,24,27]]) + product[1] = np.array([[4,8,12],[20,25,30],[42,48,54]]) + product[2] = np.array([[7,14,21],[32,40,48],[63,72,81]]) + + analytical_product = af.interop.np_to_af_array(product) + + check_order3 = af.sum(af.abs(wave_equation.Li_Lp_xi\ + (test_array, test_array1) - analytical_product)) <= threshold + + assert check_order3 def test_dx_dxi(): - ''' - A Test function to check the dx_xi function in wave_equation module by - passing nodes of an element and using the LGL points. Analytically, the - differential would be a constant. The check has a tolerance 1e-7. - ''' - threshold = 1e-14 - gvar.populateGlobalVariables(8) - nodes = np.array([7, 10], dtype = np.float64) - test_nodes = af.interop.np_to_af_array(nodes) - analytical_dx_dxi = 1.5 - - check_dx_dxi = (af.statistics.mean(wave_equation.dx_dxi_numerical - (test_nodes,gvar.xi_LGL)) - analytical_dx_dxi) <= threshold - - assert check_dx_dxi + ''' + A Test function to check the dx_xi function in wave_equation module by + passing nodes of an element and using the LGL points. Analytically, the + differential would be a constant. The check has a tolerance 1e-7. + ''' + threshold = 1e-14 + gvar.populateGlobalVariables(8) + nodes = np.array([7, 10], dtype = np.float64) + test_nodes = af.interop.np_to_af_array(nodes) + analytical_dx_dxi = 1.5 + + check_dx_dxi = (af.statistics.mean(wave_equation.dx_dxi_numerical + (test_nodes,gvar.xi_LGL)) - analytical_dx_dxi) <= threshold + + assert check_dx_dxi def test_dx_dxi_analytical(): - ''' - Test to check the dx_dxi_analytical in wave equation module for an element - and compare it with an analytical value - ''' - threshold = 1e-14 - nodes = af.Array([2,6]) - check_analytical_dx_dxi = af.sum(af.abs(wave_equation.dx_dxi_analytical - (nodes, 0) - 2)) <= threshold - assert check_analytical_dx_dxi + ''' + Test to check the dx_dxi_analytical in wave equation module for an element + and compare it with an analytical value + ''' + threshold = 1e-14 + nodes = af.Array([2,6]) + check_analytical_dx_dxi = af.sum(af.abs(wave_equation.dx_dxi_analytical + (nodes, 0) - 2)) <= threshold + assert check_analytical_dx_dxi def test_lBasisArray(): - ''' - Function to test the lBasisArray function in global_variables module by - passing 8 LGL points and comparing the numerically obtained basis function - coefficients to analytically calculated ones. - - Reference - --------- - The link to the sage worksheet where the calculations were carried out. - - https://cocalc.com/projects/1b7f404c-87ba-40d0-816c-2eba17466aa8/files\ - /PM_2_5/wave_equation/worksheets/l_basis_array.sagews - ''' - threshold = 1e-12 - gvar.populateGlobalVariables(8) - basis_array_analytical = np.zeros([8, 8]) - - basis_array_analytical[0] = np.array([-3.351562500008004,\ - 3.351562500008006, \ - 3.867187500010295,\ - -3.867187500010297,\ - - 1.054687500002225, \ - 1.054687500002225, \ - 0.03906249999993106,\ - - 0.03906249999993102]) - basis_array_analytical[1] = np.array([8.140722718246403,\ - - 7.096594831382852,\ - - 11.34747768400062,\ - 9.89205188146461, \ - 3.331608712119162, \ - - 2.904297073479968,\ - - 0.1248537463649464,\ - 0.1088400233982081]) - basis_array_analytical[2] = np.array([-10.35813682892759,\ - 6.128911440984293,\ - 18.68335515838398,\ - - 11.05494463699297,\ - - 8.670037141196786,\ - 5.130062549476987,\ - 0.3448188117404021,\ - - 0.2040293534683072]) + ''' + Function to test the lBasisArray function in global_variables module by + passing 8 LGL points and comparing the numerically obtained basis function + coefficients to analytically calculated ones. + + Reference + --------- + The link to the sage worksheet where the calculations were carried out. + + https://cocalc.com/projects/1b7f404c-87ba-40d0-816c-2eba17466aa8/files\ + /PM_2_5/wave_equation/worksheets/l_basis_array.sagews + ''' + threshold = 1e-12 + gvar.populateGlobalVariables(8) + basis_array_analytical = np.zeros([8, 8]) + + basis_array_analytical[0] = np.array([-3.351562500008004,\ + 3.351562500008006, \ + 3.867187500010295,\ + -3.867187500010297,\ + - 1.054687500002225, \ + 1.054687500002225, \ + 0.03906249999993106,\ + - 0.03906249999993102]) + basis_array_analytical[1] = np.array([8.140722718246403,\ + - 7.096594831382852,\ + - 11.34747768400062,\ + 9.89205188146461, \ + 3.331608712119162, \ + - 2.904297073479968,\ + - 0.1248537463649464,\ + 0.1088400233982081]) + basis_array_analytical[2] = np.array([-10.35813682892759,\ + 6.128911440984293,\ + 18.68335515838398,\ + - 11.05494463699297,\ + - 8.670037141196786,\ + 5.130062549476987,\ + 0.3448188117404021,\ + - 0.2040293534683072]) - basis_array_analytical[3] = np.array([11.38981374849497,\ - - 2.383879109609436,\ - - 24.03296250200938,\ - 5.030080255538657,\ - 15.67350804691132,\ - - 3.28045297599924,\ - - 3.030359293396907,\ - 0.6342518300700298]) + basis_array_analytical[3] = np.array([11.38981374849497,\ + - 2.383879109609436,\ + - 24.03296250200938,\ + 5.030080255538657,\ + 15.67350804691132,\ + - 3.28045297599924,\ + - 3.030359293396907,\ + 0.6342518300700298]) - basis_array_analytical[4] = np.array([-11.38981374849497,\ - - 2.383879109609437,\ - 24.03296250200939,\ - 5.030080255538648,\ - - 15.67350804691132,\ - - 3.28045297599924,\ - 3.030359293396907,\ - 0.6342518300700299]) + basis_array_analytical[4] = np.array([-11.38981374849497,\ + - 2.383879109609437,\ + 24.03296250200939,\ + 5.030080255538648,\ + - 15.67350804691132,\ + - 3.28045297599924,\ + 3.030359293396907,\ + 0.6342518300700299]) - basis_array_analytical[5] = np.array([10.35813682892759,\ - 6.128911440984293,\ - -18.68335515838398,\ - - 11.05494463699297,\ - 8.670037141196786,\ - 5.130062549476987,\ - - 0.3448188117404021,\ - - 0.2040293534683072]) - basis_array_analytical[6] = np.array([-8.140722718246403,\ - - 7.096594831382852,\ - 11.34747768400062,\ - 9.89205188146461, \ - -3.331608712119162, \ - - 2.904297073479968,\ - 0.1248537463649464,\ - 0.1088400233982081]) - basis_array_analytical[7] = np.array([3.351562500008004,\ - 3.351562500008005, \ - - 3.867187500010295,\ - - 3.867187500010298,\ - 1.054687500002225, \ - 1.054687500002224, \ - - 0.039062499999931,\ - - 0.03906249999993102]) - - basis_array_analytical = af.interop.np_to_af_array(basis_array_analytical) - - assert af.sum(af.abs(basis_array_analytical - gvar.lBasisArray)) < threshold + basis_array_analytical[5] = np.array([10.35813682892759,\ + 6.128911440984293,\ + -18.68335515838398,\ + - 11.05494463699297,\ + 8.670037141196786,\ + 5.130062549476987,\ + - 0.3448188117404021,\ + - 0.2040293534683072]) + basis_array_analytical[6] = np.array([-8.140722718246403,\ + - 7.096594831382852,\ + 11.34747768400062,\ + 9.89205188146461, \ + -3.331608712119162, \ + - 2.904297073479968,\ + 0.1248537463649464,\ + 0.1088400233982081]) + basis_array_analytical[7] = np.array([3.351562500008004,\ + 3.351562500008005, \ + - 3.867187500010295,\ + - 3.867187500010298,\ + 1.054687500002225, \ + 1.054687500002224, \ + - 0.039062499999931,\ + - 0.03906249999993102]) + + basis_array_analytical = af.interop.np_to_af_array(basis_array_analytical) + + assert af.sum(af.abs(basis_array_analytical - gvar.lBasisArray)) < threshold def test_Integral_Li_Lp(): - ''' - Test function to check the A_matrix function in wave_equation module. - - Obtaining the A_matrix from the function and setting the value of - all elements above a certain threshold to be 1 and plotting it. - ''' - threshold = 1e-5 - gvar.populateGlobalVariables(8, 8) - A_matrix_structure = np.zeros([gvar.N_LGL, gvar.N_LGL]) - non_zero_indices = np.where(np.array(wave_equation.A_matrix()) > threshold) - - A_matrix_structure[non_zero_indices] = 1. - - plt.gca().invert_yaxis() - plt.contourf(A_matrix_structure, 100, cmap = 'Blues') - plt.axes().set_aspect('equal') - plt.colorbar() - plt.show() - - return + ''' + Test function to check the A_matrix function in wave_equation module. + + Obtaining the A_matrix from the function and setting the value of + all elements above a certain threshold to be 1 and plotting it. + ''' + threshold = 1e-5 + gvar.populateGlobalVariables(8, 8) + A_matrix_structure = np.zeros([gvar.N_LGL, gvar.N_LGL]) + non_zero_indices = np.where(np.array(wave_equation.A_matrix()) > threshold) + + A_matrix_structure[non_zero_indices] = 1. + + plt.gca().invert_yaxis() + plt.contourf(A_matrix_structure, 100, cmap = 'Blues') + plt.axes().set_aspect('equal') + plt.colorbar() + plt.show() + + return def test_A_matrix(): - ''' - Test function to check the A matrix obtained from wave_equation module with - one obtained by numerical integral solvers. - - NOTE - ---- - A diagonal A-matrix is used as a reference in this unit test. - The A matrix calculated analytically gives a different matrix. - - ''' - threshold = 2e-10 - gvar.populateGlobalVariables(8) - - reference_A_matrix = af.tile(gvar.lobatto_weights, 1, gvar.N_LGL)\ - * af.identity(gvar.N_LGL, gvar.N_LGL, dtype = af.Dtype.f64)\ - * af.mean(wave_equation.dx_dxi_numerical((gvar.elementMeshNodes[0 : 2])\ - ,gvar.xi_LGL)) - - test_A_matrix = wave_equation.A_matrix() - error_array = af.abs(reference_A_matrix - test_A_matrix) - - print(af.max(error_array)) - - assert af.algorithm.max(error_array) < threshold + ''' + Test function to check the A matrix obtained from wave_equation module with + one obtained by numerical integral solvers. + + NOTE + ---- + A diagonal A-matrix is used as a reference in this unit test. + The A matrix calculated analytically gives a different matrix. + + ''' + threshold = 2e-10 + gvar.populateGlobalVariables(8) + + reference_A_matrix = af.tile(gvar.lobatto_weights, 1, gvar.N_LGL)\ + * af.identity(gvar.N_LGL, gvar.N_LGL, dtype = af.Dtype.f64)\ + * af.mean(wave_equation.dx_dxi_numerical((gvar.elementMeshNodes[0 : 2])\ + ,gvar.xi_LGL)) + + test_A_matrix = wave_equation.A_matrix() + error_array = af.abs(reference_A_matrix - test_A_matrix) + + print(af.max(error_array)) + + assert af.algorithm.max(error_array) < threshold def test_dLp_xi(): - ''' - Test function to check the dLp_xi calculated in gvar mdule with a - numerically obtained one. - - Refrence - -------- - The link to the sage worksheet where the calculations were carried out. - https://cocalc.com/projects/1b7f404c-87ba-40d0-816c-2eba17466aa8/files\ - /PM_2_5/wave_equation/worksheets/dLp_xi.sagews - ''' - threshold = 1e-13 - gvar.x_nodes = af.interop.np_to_af_array(np.array([-1., 1.])) - gvar.populateGlobalVariables(8) - gvar.c = 1 - - - reference_d_Lp_xi = af.interop.np_to_af_array(np.array([\ - [-14.0000000000226,-3.20991570302344,0.792476681323880,-0.372150435728984,\ - 0.243330712724289,-0.203284568901545,0.219957514771985, -0.500000000000000], - - [18.9375986071129, 3.31499272476776e-11, -2.80647579473469,1.07894468878725\ - ,-0.661157350899271,0.537039586158262, -0.573565414940005,1.29768738831567], - - [-7.56928981931106, 4.54358506455201, -6.49524878326702e-12, \ - -2.37818723350641, 1.13535801687865, -0.845022556506714, 0.869448098330221,\ - -1.94165942553537], - - [4.29790816425547,-2.11206121431525,2.87551740597844,-1.18896004153157e-11,\ - -2.38892435916370, 1.37278583181113, -1.29423205091574, 2.81018898925442], - - [-2.81018898925442, 1.29423205091574, -1.37278583181113, 2.38892435916370, \ - 1.18892673484083e-11,-2.87551740597844, 2.11206121431525,-4.29790816425547], - - [1.94165942553537, -0.869448098330221, 0.845022556506714,-1.13535801687865,\ - 2.37818723350641, 6.49524878326702e-12,-4.54358506455201,7.56928981931106],\ - - [-1.29768738831567, 0.573565414940005,-0.537039586158262,0.661157350899271,\ - -1.07894468878725,2.80647579473469,-3.31498162253752e-11,-18.9375986071129], - - [0.500000000000000,-0.219957514771985,0.203284568901545,-0.243330712724289,\ - 0.372150435728984, -0.792476681323880, 3.20991570302344, 14.0000000000226] - ])) - - assert af.max(reference_d_Lp_xi - gvar.dLp_xi) < threshold + ''' + Test function to check the dLp_xi calculated in gvar mdule with a + numerically obtained one. + + Refrence + -------- + The link to the sage worksheet where the calculations were carried out. + https://cocalc.com/projects/1b7f404c-87ba-40d0-816c-2eba17466aa8/files\ + /PM_2_5/wave_equation/worksheets/dLp_xi.sagews + ''' + threshold = 1e-13 + gvar.x_nodes = af.interop.np_to_af_array(np.array([-1., 1.])) + gvar.populateGlobalVariables(8) + gvar.c = 1 + + + reference_d_Lp_xi = af.interop.np_to_af_array(np.array([\ + [-14.0000000000226,-3.20991570302344,0.792476681323880,-0.372150435728984,\ + 0.243330712724289,-0.203284568901545,0.219957514771985, -0.500000000000000], + + [18.9375986071129, 3.31499272476776e-11, -2.80647579473469,1.07894468878725\ + ,-0.661157350899271,0.537039586158262, -0.573565414940005,1.29768738831567], + + [-7.56928981931106, 4.54358506455201, -6.49524878326702e-12, \ + -2.37818723350641, 1.13535801687865, -0.845022556506714, 0.869448098330221,\ + -1.94165942553537], + + [4.29790816425547,-2.11206121431525,2.87551740597844,-1.18896004153157e-11,\ + -2.38892435916370, 1.37278583181113, -1.29423205091574, 2.81018898925442], + + [-2.81018898925442, 1.29423205091574, -1.37278583181113, 2.38892435916370, \ + 1.18892673484083e-11,-2.87551740597844, 2.11206121431525,-4.29790816425547], + + [1.94165942553537, -0.869448098330221, 0.845022556506714,-1.13535801687865,\ + 2.37818723350641, 6.49524878326702e-12,-4.54358506455201,7.56928981931106],\ + + [-1.29768738831567, 0.573565414940005,-0.537039586158262,0.661157350899271,\ + -1.07894468878725,2.80647579473469,-3.31498162253752e-11,-18.9375986071129], + + [0.500000000000000,-0.219957514771985,0.203284568901545,-0.243330712724289,\ + 0.372150435728984, -0.792476681323880, 3.20991570302344, 14.0000000000226] + ])) + + assert af.max(reference_d_Lp_xi - gvar.dLp_xi) < threshold def test_volume_integral_flux(): - ''' - A test function to check the volume_integral_flux function in wave_equation - module by analytically calculated Gauss-Lobatto quadrature. - - Reference - --------- - The link to the sage worksheet where the calculations were caried out is - given below. - https://cocalc.com/projects/1b7f404c-87ba-40d0-816c-2eba17466aa8/files\ - /PM_2_5/wave_equation/worksheets/volume_integral_flux.sagews - - ''' - threshold = 4e-8 - gvar.populateGlobalVariables(8) - gvar.c = 1 - - referenceFluxIntegral = af.transpose(af.interop.np_to_af_array(np.array([ - [-0.002016634876668093, -0.000588597708116113, -0.0013016773719126333,\ - -0.002368387579324652, -0.003620502047659841, -0.004320197094090966, - -0.003445512010153811, 0.0176615086879261],\ - [-0.018969769374, -0.00431252844519,-0.00882630935977,-0.0144355176966,\ - -0.019612124119, -0.0209837936827, -0.0154359890788, 0.102576031756], \ - [-0.108222418798, -0.0179274222595, -0.0337807018822, -0.0492589052599,\ - -0.0588472807471, -0.0557970236273, -0.0374764132459, 0.361310165819],\ - [-0.374448714304, -0.0399576371245, -0.0683852285846, -0.0869229749357,\ - -0.0884322503841, -0.0714664112839, -0.0422339853622, 0.771847201979], \ - [-0.785754362849, -0.0396035640187, -0.0579313769517, -0.0569022801117,\ - -0.0392041960688, -0.0172295769141, -0.00337464521455, 1.00000000213],\ - [-1.00000000213, 0.00337464521455, 0.0172295769141, 0.0392041960688,\ - 0.0569022801117, 0.0579313769517, 0.0396035640187, 0.785754362849],\ - [-0.771847201979, 0.0422339853622, 0.0714664112839, 0.0884322503841, \ - 0.0869229749357, 0.0683852285846, 0.0399576371245, 0.374448714304],\ - [-0.361310165819, 0.0374764132459, 0.0557970236273, 0.0588472807471,\ - 0.0492589052599, 0.0337807018822, 0.0179274222595, 0.108222418798], \ - [-0.102576031756, 0.0154359890788, 0.0209837936827, 0.019612124119, \ - 0.0144355176966, 0.00882630935977, 0.00431252844519, 0.018969769374],\ - [-0.0176615086879, 0.00344551201015 ,0.00432019709409, 0.00362050204766,\ - 0.00236838757932, 0.00130167737191, 0.000588597708116, 0.00201663487667\ - ]]))) - - numerical_flux = wave_equation.volumeIntegralFlux(gvar.u[:, :, 0]) - assert (af.max(af.abs(numerical_flux - referenceFluxIntegral)) < threshold) + ''' + A test function to check the volume_integral_flux function in wave_equation + module by analytically calculated Gauss-Lobatto quadrature. + + Reference + --------- + The link to the sage worksheet where the calculations were caried out is + given below. + https://cocalc.com/projects/1b7f404c-87ba-40d0-816c-2eba17466aa8/files\ + /PM_2_5/wave_equation/worksheets/volume_integral_flux.sagews + + ''' + threshold = 4e-8 + gvar.populateGlobalVariables(8) + gvar.c = 1 + + referenceFluxIntegral = af.transpose(af.interop.np_to_af_array(np.array([ + [-0.002016634876668093, -0.000588597708116113, -0.0013016773719126333,\ + -0.002368387579324652, -0.003620502047659841, -0.004320197094090966, + -0.003445512010153811, 0.0176615086879261],\ + [-0.018969769374, -0.00431252844519,-0.00882630935977,-0.0144355176966,\ + -0.019612124119, -0.0209837936827, -0.0154359890788, 0.102576031756], \ + [-0.108222418798, -0.0179274222595, -0.0337807018822, -0.0492589052599,\ + -0.0588472807471, -0.0557970236273, -0.0374764132459, 0.361310165819],\ + [-0.374448714304, -0.0399576371245, -0.0683852285846, -0.0869229749357,\ + -0.0884322503841, -0.0714664112839, -0.0422339853622, 0.771847201979], \ + [-0.785754362849, -0.0396035640187, -0.0579313769517, -0.0569022801117,\ + -0.0392041960688, -0.0172295769141, -0.00337464521455, 1.00000000213],\ + [-1.00000000213, 0.00337464521455, 0.0172295769141, 0.0392041960688,\ + 0.0569022801117, 0.0579313769517, 0.0396035640187, 0.785754362849],\ + [-0.771847201979, 0.0422339853622, 0.0714664112839, 0.0884322503841, \ + 0.0869229749357, 0.0683852285846, 0.0399576371245, 0.374448714304],\ + [-0.361310165819, 0.0374764132459, 0.0557970236273, 0.0588472807471,\ + 0.0492589052599, 0.0337807018822, 0.0179274222595, 0.108222418798], \ + [-0.102576031756, 0.0154359890788, 0.0209837936827, 0.019612124119, \ + 0.0144355176966, 0.00882630935977, 0.00431252844519, 0.018969769374],\ + [-0.0176615086879, 0.00344551201015 ,0.00432019709409, 0.00362050204766,\ + 0.00236838757932, 0.00130167737191, 0.000588597708116, 0.00201663487667\ + ]]))) + + numerical_flux = wave_equation.volumeIntegralFlux(gvar.u[:, :, 0]) + assert (af.max(af.abs(numerical_flux - referenceFluxIntegral)) < threshold) def test_lax_friedrichs_flux(): - ''' - A test function to test the laxFriedrichsFlux function in wave_equation - module. - ''' - threshold = 1e-14 - gvar.populateGlobalVariables(8, 10) - gvar.c = 1 - - f_i = wave_equation.laxFriedrichsFlux(0) - analytical_lax_friedrichs_flux = gvar.u[-1, :, 0] - assert af.max(af.abs(analytical_lax_friedrichs_flux - f_i)) < threshold + ''' + A test function to test the laxFriedrichsFlux function in wave_equation + module. + ''' + threshold = 1e-14 + gvar.populateGlobalVariables(8, 10) + gvar.c = 1 + + f_i = wave_equation.laxFriedrichsFlux(0) + analytical_lax_friedrichs_flux = gvar.u[-1, :, 0] + assert af.max(af.abs(analytical_lax_friedrichs_flux - f_i)) < threshold def test_surface_term(): - ''' - A test function to test the surface_term function in the wave_equation - module using analytical Lax-Friedrichs flux. - ''' - threshold = 1e-13 - gvar.populateGlobalVariables(8, 10) - gvar.c = 1 - - - analytical_f_i = (gvar.u[-1, :, 0]) - analytical_f_i_minus1 = (af.shift(gvar.u[-1, :, 0], 0, 1)) - - L_p_1 = af.constant(0, gvar.N_LGL, dtype = af.Dtype.f64) - L_p_1[gvar.N_LGL - 1] = 1 - - L_p_minus1 = af.constant(0, gvar.N_LGL, dtype = af.Dtype.f64) - L_p_minus1[0] = 1 - - analytical_surface_term = af.blas.matmul(L_p_1, analytical_f_i)\ - - af.blas.matmul(L_p_minus1, analytical_f_i_minus1) - - numerical_surface_term = (wave_equation.surface_term(0)) - assert af.max(af.abs(analytical_surface_term - numerical_surface_term)) \ - < threshold - return analytical_surface_term + ''' + A test function to test the surface_term function in the wave_equation + module using analytical Lax-Friedrichs flux. + ''' + threshold = 1e-13 + gvar.populateGlobalVariables(8, 10) + gvar.c = 1 + + + analytical_f_i = (gvar.u[-1, :, 0]) + analytical_f_i_minus1 = (af.shift(gvar.u[-1, :, 0], 0, 1)) + + L_p_1 = af.constant(0, gvar.N_LGL, dtype = af.Dtype.f64) + L_p_1[gvar.N_LGL - 1] = 1 + + L_p_minus1 = af.constant(0, gvar.N_LGL, dtype = af.Dtype.f64) + L_p_minus1[0] = 1 + + analytical_surface_term = af.blas.matmul(L_p_1, analytical_f_i)\ + - af.blas.matmul(L_p_minus1, analytical_f_i_minus1) + + numerical_surface_term = (wave_equation.surface_term(0)) + assert af.max(af.abs(analytical_surface_term - numerical_surface_term)) \ + < threshold + return analytical_surface_term def test_b_vector(): - ''' - A test function to check the b vector obtained analytically and compare it - with the one returned by b_vector function in wave_equation module - ''' - threshold = 1e-13 - gvar.populateGlobalVariables(8) - gvar.c = 1 - - u_n_A_matrix = af.blas.matmul(wave_equation.A_matrix(), gvar.u[:, :, 0]) - volume_integral_flux = wave_equation.volumeIntegralFlux(gvar.u[:, :, 0]) - surface_term = test_surface_term() - b_vector_analytical = u_n_A_matrix + (volume_integral_flux -\ - (surface_term)) * gvar.delta_t - b_vector_array = wave_equation.b_vector(0) - - assert (b_vector_analytical - b_vector_array) < threshold + ''' + A test function to check the b vector obtained analytically and compare it + with the one returned by b_vector function in wave_equation module + ''' + threshold = 1e-13 + gvar.populateGlobalVariables(8) + gvar.c = 1 + + u_n_A_matrix = af.blas.matmul(wave_equation.A_matrix(), gvar.u[:, :, 0]) + volume_integral_flux = wave_equation.volumeIntegralFlux(gvar.u[:, :, 0]) + surface_term = test_surface_term() + b_vector_analytical = u_n_A_matrix + (volume_integral_flux -\ + (surface_term)) * gvar.delta_t + b_vector_array = wave_equation.b_vector(0) + + assert (b_vector_analytical - b_vector_array) < threshold diff --git a/code/utils/utils.py b/code/utils/utils.py index af31bc8..ef06f2f 100644 --- a/code/utils/utils.py +++ b/code/utils/utils.py @@ -3,30 +3,30 @@ af.set_device(1) def add(a, b): - ''' - ''' - return a + b + ''' + ''' + return a + b def divide(a, b): - ''' - ''' - return a / b + ''' + ''' + return a / b def multiply(a, b): - ''' - ''' - return a * b + ''' + ''' + return a * b def linspace(start, end, number_of_points): - ''' - Linspace implementation using arrayfire. - ''' - X = af.range(number_of_points, dtype = af.Dtype.f64) - d = (end - start) / (number_of_points - 1) - X = X * d - X = X + start - - return X + ''' + Linspace implementation using arrayfire. + ''' + X = af.range(number_of_points, dtype = af.Dtype.f64) + d = (end - start) / (number_of_points - 1) + X = X * d + X = X + start + + return X From 4b5b4a7d523ec939d2924c538cdb414c7cbc5e35 Mon Sep 17 00:00:00 2001 From: Balavarun5 Date: Fri, 8 Sep 2017 13:19:07 +0530 Subject: [PATCH 14/19] Cleanup of the code, Unit tests working. Iteration speed unchanged (decreases with larger array size). --- code/app/global_variables.py | 314 ++++++++++----------------------- code/app/lagrange.py | 75 +++++--- code/app/wave_equation.py | 136 ++++++-------- code/main.py | 3 - code/unit_test/test_waveEqn.py | 101 +++-------- 5 files changed, 229 insertions(+), 400 deletions(-) diff --git a/code/app/global_variables.py b/code/app/global_variables.py index eb70cd8..bf3f424 100644 --- a/code/app/global_variables.py +++ b/code/app/global_variables.py @@ -1,194 +1,47 @@ import numpy as np import arrayfire as af af.set_backend('opencl') -af.set_device(1) from scipy import special as sp from app import lagrange from app import wave_equation from utils import utils import math -LGL_list = [ \ -[-1.0,1.0], \ -[-1.0,0.0,1.0], \ -[-1.0,-0.4472135955,0.4472135955,1.0], \ -[-1.0,-0.654653670708,0.0,0.654653670708,1.0], \ -[-1.0,-0.765055323929,-0.285231516481,0.285231516481,0.765055323929,1.0], \ -[-1.0,-0.830223896279,-0.468848793471,0.0,0.468848793471,0.830223896279, \ -1.0], \ -[-1.0,-0.87174014851,-0.591700181433,-0.209299217902,0.209299217902, \ -0.591700181433,0.87174014851,1.0], \ -[-1.0,-0.899757995411,-0.677186279511,-0.363117463826,0.0,0.363117463826, \ -0.677186279511,0.899757995411,1.0], \ -[-1.0,-0.919533908167,-0.738773865105,-0.47792494981,-0.165278957666, \ -0.165278957666,0.47792494981,0.738773865106,0.919533908166,1.0], \ -[-1.0,-0.934001430408,-0.784483473663,-0.565235326996,-0.295758135587, \ -0.0,0.295758135587,0.565235326996,0.784483473663,0.934001430408,1.0], \ -[-1.0,-0.944899272223,-0.819279321644,-0.632876153032,-0.399530940965, \ --0.136552932855,0.136552932855,0.399530940965,0.632876153032, \ -0.819279321644,0.944899272223,1.0], \ -[-1.0,-0.953309846642,-0.846347564652,-0.686188469082,-0.482909821091, \ --0.249286930106,0.0,0.249286930106,0.482909821091,0.686188469082, \ -0.846347564652,0.953309846642,1.0], \ -[-0.999999999996,-0.959935045274,-0.867801053826,-0.728868599093, \ --0.550639402928,-0.342724013343,-0.116331868884,0.116331868884, \ -0.342724013343,0.550639402929,0.728868599091,0.86780105383, \ -0.959935045267,1.0], \ -[-0.999999999996,-0.965245926511,-0.885082044219,-0.763519689953, \ --0.60625320547,-0.420638054714,-0.215353955364,0.0,0.215353955364, \ -0.420638054714,0.60625320547,0.763519689952,0.885082044223, \ -0.965245926503,1.0], \ -[-0.999999999984,-0.9695680463,-0.899200533072,-0.792008291871, \ --0.65238870288,-0.486059421887,-0.299830468901,-0.101326273522, \ -0.101326273522,0.299830468901,0.486059421887,0.652388702882, \ -0.792008291863,0.899200533092,0.969568046272,0.999999999999]] +# This module is used to change the parameters of the simulation. -for idx in np.arange(len(LGL_list)): - LGL_list[idx] = np.array(LGL_list[idx], dtype = np.float64) - LGL_list[idx] = af.interop.np_to_af_array(LGL_list[idx]) +# The functions to calculate Legendre-Gauss-Lobatto points as well +# as the Lobatto weights used for integration are given below. -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_LGL = None -u = None -time = None -c = None -dLp_xi = None -c_lax = None - -def populateGlobalVariables(Number_of_LGL_pts = 8, Number_of_elements = 10): +def LGL_points(N): ''' - Time evolution of the wave function requires the use of 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`. - - + Calculates and returns the LGL points for a given N, These are roots of + the formula. + :math: `(1 - xi ** 2) P_{n - 1}'(xi) = 0` + Legendre polynomials satisfy the recurrence relation + :math:`(1 - x ** 2) P_n' (x) = -n x P_n(x) + n P_{n - 1} (x)` + 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_LGL - global elementMeshNodes - - 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))) - - elementMeshNodes = utils.linspace(af.sum(x_nodes[0]), - af.sum(x_nodes[1]), - N_Elements + 1) - - - element_array = af.transpose(af.interop.np_to_af_array(np_element_array)) - element_LGL = wave_equation.mappingXiToX(\ - af.transpose(element_array), xi_LGL) - - global c - global c_lax - global delta_t - c = 4 - delta_x = af.min((element_LGL - af.shift(element_LGL, 1, 0))[1:, :]) - delta_t = delta_x / (20 * c) - c_lax = 1 - - global u - global time - total_time = 5 - time = utils.linspace(0, total_time, int(total_time / delta_t)) - u_init = np.e ** (-(element_LGL) ** 2 / 0.4 ** 2) - u = af.constant(0, N_LGL, N_Elements, time.shape[0],\ - dtype = af.Dtype.f64) - u[:, :, 0] = u_init - + N : Int + Number of LGL nodes required - - global dLp_xi - dLp_xi = dLp_xi_LGL() - - - return + Returns + ------- + lgl : arrayfire.Array [N 1 1 1] + An arrayfire array consisting of the Lagrange-Gauss-Lobatto Nodes. + + Reference + --------- + http://mathworld.wolfram.com/LobattoQuadrature.html + ''' + xi = np.poly1d([1, 0]) + legendre_N_minus_1 = N * (xi * sp.legendre(N - 1) - sp.legendre(N)) + lgl_points = legendre_N_minus_1.r + lgl_points.sort() + lgl_points = af.np_to_af_array(lgl_points) + + return lgl_points def lobatto_weight_function(n, x): @@ -197,7 +50,7 @@ def lobatto_weight_function(n, x): and points :math: `x`. :math:: - w_{n} = \\frac{2 P(x)^2}{n (n - 1)}, + `w_{n} = \\frac{2 P(x)^2}{n (n - 1)}`, Where P(x) is $ (n - 1)^th $ index. Parameters @@ -205,13 +58,13 @@ def lobatto_weight_function(n, x): n : int Index for which lobatto weight function - x : arrayfire.Array + x : arrayfire.Array [N 1 1 1] 1D array of points where weight function is to be calculated. Returns ------- - gaussLobattoWeights : arrayfire.Array + gauss_lobatto_weights : arrayfire.Array An array of lobatto weight functions for the given :math: `x` points and index. Reference @@ -223,50 +76,75 @@ def lobatto_weight_function(n, x): ''' P = sp.legendre(n - 1) - gaussLobattoWeights = (2 / (n * (n - 1)) / (P(x))**2) + gauss_lobatto_weights = (2 / (n * (n - 1)) / (P(x))**2) - return gaussLobattoWeights + return gauss_lobatto_weights -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 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] + +# The domain of the function which is of interest. +x_nodes = af.np_to_af_array(np.array([-1., 1.])) + +# The number of LGL points into which an element is split +N_LGL = 8 + +# Number of elements the domain is to be divided into +N_Elements = 10 + +# The speed of the wave. +c = 1 + +# The total time for which the simulation is to be carried out. +total_time = 1 + +# The c_lax to be used in the Lax-Friedrichs flux. +c_lax = 1 + +# Array containing the LGL points in xi space. +xi_LGL = LGL_points(N_LGL) + + +# Array of the lobatto weights used during integration. +lobatto_weights = af.np_to_af_array(lobatto_weight_function(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)) + +# Refer corresponding functions. +dLp_xi = lagrange.dLp_xi_LGL(lagrange_coeffs) +lagrange_basis_value = lagrange.lagrange_basis_function(lagrange_coeffs) + +# Obtaining an array consisting of the LGL points mapped onto the 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) - 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) +np_element_array = np.concatenate((af.transpose(elements), + af.transpose(elements + element_size))) - dLp_xi = af.blas.matmul(differentiation_coeffs, nodes_power_tile) +element_mesh_nodes = utils.linspace(af.sum(x_nodes[0]), + af.sum(x_nodes[1]), + N_Elements + 1) - return dLp_xi +element_array = af.transpose(af.interop.np_to_af_array(np_element_array)) +element_LGL = wave_equation.mapping_xi_to_x(af.transpose(element_array), xi_LGL) + + +# The minimum distance between 2 mapped LGL points. +delta_x = af.min((element_LGL - af.shift(element_LGL, 1, 0))[1:, :]) + +# The value of time-step. +delta_t = delta_x / (20 * c) + +# Array of timesteps seperated by delta_t. +time = utils.linspace(0, int(total_time / delta_t) * delta_t, + int(total_time / delta_t)) + +# Initializing the amplitudes. Change u_init to required initial conditions. +u_init = np.e ** (-(element_LGL) ** 2 / 0.4 ** 2) +u = af.constant(0, N_LGL, N_Elements, time.shape[0],\ + dtype = af.Dtype.f64) +u[:, :, 0] = u_init diff --git a/code/app/lagrange.py b/code/app/lagrange.py index 14ef566..2301bb4 100644 --- a/code/app/lagrange.py +++ b/code/app/lagrange.py @@ -4,33 +4,6 @@ from utils import utils from app import global_variables as gvar -def LGL_points(N): - ''' - Returns the :math: `N` Legendre-Gauss-Laguere points which are - the roots of the equation - :math:: - (1 - x^2)L'_N = 0 - - Parameters - ---------- - N : int - Number of LGL-points to be generated. - 2 < N < 16 - - Returns - ------- - lgl : arrayfire.Array - An array of :math: `N` LGL points. - ''' - if N > 16 or N < 2: - print('Skipping! This function can only return from ', - '2 to 16 LGL points.') - - n = N - 2 - - lgl = af.Array(gvar.LGL_list[n]) - return lgl - def lagrange_basis_coeffs(x): ''' @@ -100,3 +73,51 @@ def lagrange_basis(i, x): l_xi_j = af.blas.matmul(gvar.lBasisArray[i], x_pow) return l_xi_j + + +def lagrange_basis_function(lagrange_coeff_array): + ''' + 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(gvar.xi_LGL, 1, gvar.N_LGL)) + power = af.flip(af.range(gvar.N_LGL)) + power_tile = af.tile(power, 1, gvar.N_LGL) + xi_pow = af.arith.pow(xi_tile, power_tile) + index = af.range(gvar.N_LGL) + L_i = af.blas.matmul(lagrange_coeff_array[index], xi_pow) + + return L_i + + +def dLp_xi_LGL(lagrange_coeff_array): + ''' + 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(gvar.N_LGL), 1, gvar.N_LGL))) + * lagrange_coeff_array)[:, :-1] + + nodes_tile = af.transpose(af.tile(gvar.xi_LGL, 1, gvar.N_LGL - 1)) + power_tile = af.flip(af.tile(af.range(gvar.N_LGL - 1), 1, gvar.N_LGL)) + nodes_power_tile = af.pow(nodes_tile, power_tile) + + dLp_xi = af.blas.matmul(differentiation_coeffs, nodes_power_tile) + + return dLp_xi + + diff --git a/code/app/wave_equation.py b/code/app/wave_equation.py index 8e55035..a28919a 100644 --- a/code/app/wave_equation.py +++ b/code/app/wave_equation.py @@ -2,11 +2,10 @@ import arrayfire as af af.set_backend('opencl') -af.set_device(1) import numpy as np +from app import global_variables as gvar from app import lagrange from utils import utils -from app import global_variables as gvar from matplotlib import pyplot as plt import pylab as pl from tqdm import trange @@ -39,33 +38,8 @@ plt.rcParams['ytick.direction'] = 'in' -def Li_Lp_xi(L_i_xi, L_p_xi): - ''' - Broadcasts and multiplies two arrays of different ordering. Used in the - calculation of A matrix. - - Parameters - ---------- - L_i_xi : arrayfire.Array [1 N N 1] - A 2D array :math:`L_i` obtained at LGL points calculated at - the LGL nodes :math:`N_LGL`. - - L_p_xi : arrayfire.Array [N 1 N 1] - A 2D array :math:`L_p` obtained at LGL points calculated at the - LGL nodes :math:`N_LGL`. - - Returns - ------- - Li_Lp_xi : arrayfire.Array [N N N 1] - Matrix of :math:`L_p (N_LGL) L_i (N_LGL)`. - - ''' - Li_Lp_xi = af.bcast.broadcast(utils.multiply, L_i_xi, L_p_xi) - - return Li_Lp_xi - -def mappingXiToX(x_nodes, xi): +def mapping_xi_to_x(x_nodes, xi): ''' Maps points in :math: `\\xi` space to :math:`x` space using the formula :math: `x = \\frac{1 - \\xi}{2} x_0 + \\frac{1 + \\xi}{2} x_1` @@ -89,19 +63,19 @@ def mappingXiToX(x_nodes, xi): N_0 = (1 - xi) / 2 N_1 = (1 + xi) / 2 - N0_x0 = af.bcast.broadcast(utils.multiply, N_0, x_nodes[0]) - N1_x1 = af.bcast.broadcast(utils.multiply, N_1, x_nodes[1]) + N0_x0 = af.broadcast(utils.multiply, N_0, x_nodes[0]) + N1_x1 = af.broadcast(utils.multiply, N_1, x_nodes[1]) x = N0_x0 + N1_x1 - return x + return x def dx_dxi_numerical(x_nodes, xi): ''' Differential :math: `\\frac{dx}{d \\xi}` calculated by central differential - method about xi using the mappingXiToX function. + method about xi using the mapping_xi_to_x function. Parameters ---------- @@ -118,12 +92,12 @@ def dx_dxi_numerical(x_nodes, xi): :math:`\\frac{dx}{d \\xi}`. ''' dxi = 1e-7 - x2 = mappingXiToX(x_nodes, xi + dxi) - x1 = mappingXiToX(x_nodes, xi - dxi) + x2 = mapping_xi_to_x(x_nodes, xi + dxi) + x1 = mapping_xi_to_x(x_nodes, xi - dxi) dx_dxi = (x2 - x1) / (2 * dxi) - return dx_dxi + return dx_dxi def dx_dxi_analytical(x_nodes, xi): @@ -144,7 +118,7 @@ def dx_dxi_analytical(x_nodes, xi): ''' analytical_dx_dxi = (x_nodes[1] - x_nodes[0]) / 2 - return analytical_dx_dxi + return analytical_dx_dxi def A_matrix(): @@ -166,12 +140,12 @@ def A_matrix(): The value of integral of product of lagrange basis functions obtained by LGL points, using Gauss-Lobatto quadrature method using :math: `N_LGL` points. - ''' - dx_dxi = dx_dxi_numerical((gvar.elementMeshNodes[0 : 2]), - gvar.xi_LGL) + ''' + dx_dxi = dx_dxi_numerical((gvar.element_mesh_nodes[0 : 2]), + gvar.xi_LGL) dx_dxi_tile = af.tile(dx_dxi, 1, gvar.N_LGL) identity_matrix = af.identity(gvar.N_LGL, gvar.N_LGL, dtype = af.Dtype.f64) - A_matrix = af.bcast.broadcast(utils.multiply, gvar.lobatto_weights, + A_matrix = af.broadcast(utils.multiply, gvar.lobatto_weights, identity_matrix) * dx_dxi_tile return A_matrix @@ -184,18 +158,20 @@ def flux_x(u): Parameters ---------- - u : arrayfire.Array - A 1-D array which contains the value of wave function. + u : arrayfire.Array + A 1-D array which contains the value of wave function. Returns ------- - c * u : arrayfire.Array - The value of the flux for given u. + flux : arrayfire.Array + The value of the flux for given u. ''' - return gvar.c * u + flux = gvar.c * u + return flux -def volumeIntegralFlux(u): + +def volume_integral_flux(u): ''' A function to calculate the volume integral of flux in the wave equation. :math:`\\int_{-1}^1 f(u) \\frac{d L_p}{d\\xi} d\\xi` @@ -206,8 +182,8 @@ def volumeIntegralFlux(u): Parameters ---------- u : arrayfire.Array [N_LGL N_Elements 1 1] - A 1-D array containing the value of the wave function at the - mapped LGL nodes in the element. + A 1-D array containing the value of the wave function at + the mapped LGL nodes in the element. Returns ------- @@ -226,7 +202,7 @@ def volumeIntegralFlux(u): return flux_integral -def laxFriedrichsFlux(t_n): +def lax_friedrichs_flux(u_t_n): ''' A function which calculates the lax-friedrichs_flux :math:`f_i` using. :math:`f_i = \\frac{F(u^{i + 1}_0) + F(u^i_{N_{LGL} - 1})}{2} - \frac @@ -234,36 +210,33 @@ def laxFriedrichsFlux(t_n): Parameters ---------- - u : 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. - - t_n : float - The timestep at which the lax-friedrichs-flux is to be calculated. + u_t_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. ''' - u_iplus1_0 = af.shift(gvar.u[0, :, t_n], 0, -1) - u_i_N_LGL = gvar.u[-1, :, t_n] + u_iplus1_0 = af.shift(u_t_n[0, :], 0, -1) + u_i_N_LGL = u_t_n[-1, :] flux_iplus1_0 = flux_x(u_iplus1_0) flux_i_N_LGL = flux_x(u_i_N_LGL) - laxFriedrichsFlux = (flux_iplus1_0 + flux_i_N_LGL) / 2 \ + boundary_flux = (flux_iplus1_0 + flux_i_N_LGL) / 2 \ - gvar.c_lax * (u_iplus1_0 - u_i_N_LGL) / 2 - return laxFriedrichsFlux + return boundary_flux -def surface_term(t_n): +def surface_term(u_t_n): ''' A function which is used to calculate the surface term, :math:`L_p (1) f_i - L_p (-1) f_{i - 1}` - using the laxFriedrichsFlux function and the dLp_xi_LGL function in gvar - module. + using the lax_friedrichs_flux function and lagrange_basis_value + from gvar module. Parameters ---------- - t_n : double + u_t_n : arrayfire.Array [N M 1 1] The timestep at which the surface term is to be calculated. Returns @@ -279,12 +252,13 @@ def surface_term(t_n): --------- Link to PDF describing the algorithm to obtain the surface term. - https://cocalc.com/projects/1b7f404c-87ba-40d0-816c-2eba17466aa8/files\ + https://cocalc.com/projects/1b7f404c-87ba-40d0-816c-2eba17466aa8/files /PM\_2\_5/wave\_equation/documents/surface\_term/surface\_term.pdf ''' - L_p_minus1 = gvar.lagrange_basis_function()[:, 0] - L_p_1 = gvar.lagrange_basis_function()[:, -1] - f_i = laxFriedrichsFlux(t_n) + + L_p_minus1 = gvar.lagrange_basis_value[:, 0] + L_p_1 = gvar.lagrange_basis_value[:, -1] + f_i = lax_friedrichs_flux(u_t_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) @@ -292,22 +266,22 @@ def surface_term(t_n): return surface_term -def b_vector(t_n): +def b_vector(u_t_n): ''' A function which returns the b vector for N_Elements number of elements. Parameters ---------- - t_n : double + t_n : double Returns ------- b_vector_array : arrayfire.Array ''' - volume_integral = volumeIntegralFlux(gvar.u[:, :, t_n]) - surfaceTerm = surface_term(t_n) - b_vector_array = gvar.delta_t * (volume_integral - surfaceTerm) - + + volume_integral = volume_integral_flux(u_t_n) + Surface_term = surface_term(u_t_n) + b_vector_array = gvar.delta_t * (volume_integral - Surface_term) return b_vector_array @@ -323,28 +297,30 @@ def time_evolution(): A_inverse = af.inverse(A_matrix()) element_LGL = gvar.element_LGL delta_t = gvar.delta_t + amplitude = gvar.u + time = gvar.time - - for t_n in trange(0, gvar.time.shape[0] - 1): + for t_n in trange(0, time.shape[0] - 1): - gvar.u[:, :, t_n + 1] = gvar.u[:, :, t_n] + af.blas.matmul(A_inverse,\ - b_vector(t_n)) + amplitude[:, :, t_n + 1] = amplitude[:, :, t_n] + af.blas.matmul(A_inverse,\ + b_vector(amplitude[:, :, t_n])) print('u calculated!') - for t_n in trange(0, gvar.time.shape[0] - 1): + for t_n in trange(0, time.shape[0] - 1): if t_n % 100 == 0: fig = plt.figure() x = gvar.element_LGL - y = gvar.u[:, :, t_n] + y = amplitude[:, :, t_n] plt.plot(x, y) plt.xlabel('x') plt.ylabel('Amplitude') plt.title('Time = %f' % (t_n * delta_t)) - fig.savefig('1D_Wave_images/%04d' %(t_n / 100) + '.png') + fig.savefig('results/1D_Wave_images/%04d' %(t_n / 100) + '.png') plt.close('all') return + diff --git a/code/main.py b/code/main.py index ede5969..13a9f6d 100644 --- a/code/main.py +++ b/code/main.py @@ -2,7 +2,6 @@ import arrayfire as af af.set_backend('opencl') -af.set_device(1) import numpy as np from app import global_variables as gvar from unit_test import test_waveEqn @@ -12,6 +11,4 @@ if __name__ == '__main__': - gvar.populateGlobalVariables(8,10) wave_equation.time_evolution() - #print(gvar.lobatto_weights) diff --git a/code/unit_test/test_waveEqn.py b/code/unit_test/test_waveEqn.py index 7e62fa0..da0704e 100644 --- a/code/unit_test/test_waveEqn.py +++ b/code/unit_test/test_waveEqn.py @@ -2,64 +2,30 @@ import arrayfire as af import math from matplotlib import pyplot as plt -from app import lagrange from app import global_variables as gvar +from app import lagrange from app import wave_equation from utils import utils af.set_backend('opencl') -af.set_device(1) -def test_mappingXiToX(): +def test_mapping_xi_to_x(): ''' - A test function to check the mappingXiToX function in wave_equation module, + A test function to check the mapping_xi_to_x function in wave_equation module, The test involves passing trial element nodes and :math: `\\xi` and comparing it with the x obatined by passing the trial parameters to - mappingXiToX function. + mapping_xi_to_x function. ''' threshold = 1e-14 - gvar.populateGlobalVariables() test_element_nodes = af.interop.np_to_af_array(np.array([7, 11])) test_xi = 0 analytical_x_value = 9 - numerical_x_value = wave_equation.mappingXiToX(test_element_nodes, test_xi) + numerical_x_value = wave_equation.mapping_xi_to_x(test_element_nodes, test_xi) assert af.abs(analytical_x_value - numerical_x_value) <= threshold - -def test_Li_Lp_xi(): - ''' - A Test function to check the Li_Lp_xi function in wave_equation module - by passing two test arrays and comparing the analytical product with the - numerically calculated one with to a tolerance of 1e-14. - ''' - - gvar.populateGlobalVariables(3) - - threshold = 1e-14 - - test_array = af.interop.np_to_af_array(np.array([[1, 2, 3], \ - [4, 5, 6], \ - [7, 8, 9]], \ - dtype = np.float64)) - - test_array1 = af.reorder(test_array, 2, 0, 1) - - product = np.zeros([3,3,3]) - product[0] = np.array([[1,2,3],[8,10,12],[21,24,27]]) - product[1] = np.array([[4,8,12],[20,25,30],[42,48,54]]) - product[2] = np.array([[7,14,21],[32,40,48],[63,72,81]]) - - analytical_product = af.interop.np_to_af_array(product) - - check_order3 = af.sum(af.abs(wave_equation.Li_Lp_xi\ - (test_array, test_array1) - analytical_product)) <= threshold - - assert check_order3 - - def test_dx_dxi(): ''' A Test function to check the dx_xi function in wave_equation module by @@ -67,7 +33,6 @@ def test_dx_dxi(): differential would be a constant. The check has a tolerance 1e-7. ''' threshold = 1e-14 - gvar.populateGlobalVariables(8) nodes = np.array([7, 10], dtype = np.float64) test_nodes = af.interop.np_to_af_array(nodes) analytical_dx_dxi = 1.5 @@ -84,6 +49,7 @@ def test_dx_dxi_analytical(): and compare it with an analytical value ''' threshold = 1e-14 + nodes = af.Array([2,6]) check_analytical_dx_dxi = af.sum(af.abs(wave_equation.dx_dxi_analytical (nodes, 0) - 2)) <= threshold @@ -103,8 +69,7 @@ def test_lBasisArray(): https://cocalc.com/projects/1b7f404c-87ba-40d0-816c-2eba17466aa8/files\ /PM_2_5/wave_equation/worksheets/l_basis_array.sagews ''' - threshold = 1e-12 - gvar.populateGlobalVariables(8) + threshold = 6e-10 basis_array_analytical = np.zeros([8, 8]) basis_array_analytical[0] = np.array([-3.351562500008004,\ @@ -177,7 +142,7 @@ def test_lBasisArray(): basis_array_analytical = af.interop.np_to_af_array(basis_array_analytical) - assert af.sum(af.abs(basis_array_analytical - gvar.lBasisArray)) < threshold + assert af.sum(af.abs(basis_array_analytical - gvar.lagrange_coeffs)) < threshold def test_Integral_Li_Lp(): @@ -188,7 +153,6 @@ def test_Integral_Li_Lp(): all elements above a certain threshold to be 1 and plotting it. ''' threshold = 1e-5 - gvar.populateGlobalVariables(8, 8) A_matrix_structure = np.zeros([gvar.N_LGL, gvar.N_LGL]) non_zero_indices = np.where(np.array(wave_equation.A_matrix()) > threshold) @@ -214,12 +178,11 @@ def test_A_matrix(): The A matrix calculated analytically gives a different matrix. ''' - threshold = 2e-10 - gvar.populateGlobalVariables(8) + threshold = 3e-10 reference_A_matrix = af.tile(gvar.lobatto_weights, 1, gvar.N_LGL)\ * af.identity(gvar.N_LGL, gvar.N_LGL, dtype = af.Dtype.f64)\ - * af.mean(wave_equation.dx_dxi_numerical((gvar.elementMeshNodes[0 : 2])\ + * af.mean(wave_equation.dx_dxi_numerical((gvar.element_mesh_nodes[0 : 2])\ ,gvar.xi_LGL)) test_A_matrix = wave_equation.A_matrix() @@ -234,48 +197,45 @@ def test_dLp_xi(): ''' Test function to check the dLp_xi calculated in gvar mdule with a numerically obtained one. - + Refrence -------- The link to the sage worksheet where the calculations were carried out. https://cocalc.com/projects/1b7f404c-87ba-40d0-816c-2eba17466aa8/files\ /PM_2_5/wave_equation/worksheets/dLp_xi.sagews ''' - threshold = 1e-13 - gvar.x_nodes = af.interop.np_to_af_array(np.array([-1., 1.])) - gvar.populateGlobalVariables(8) - gvar.c = 1 - + threshold = 4e-11 - reference_d_Lp_xi = af.interop.np_to_af_array(np.array([\ + reference_d_Lp_xi = af.np_to_af_array(np.array([\ [-14.0000000000226,-3.20991570302344,0.792476681323880,-0.372150435728984,\ 0.243330712724289,-0.203284568901545,0.219957514771985, -0.500000000000000], [18.9375986071129, 3.31499272476776e-11, -2.80647579473469,1.07894468878725\ ,-0.661157350899271,0.537039586158262, -0.573565414940005,1.29768738831567], - + [-7.56928981931106, 4.54358506455201, -6.49524878326702e-12, \ -2.37818723350641, 1.13535801687865, -0.845022556506714, 0.869448098330221,\ -1.94165942553537], - + [4.29790816425547,-2.11206121431525,2.87551740597844,-1.18896004153157e-11,\ -2.38892435916370, 1.37278583181113, -1.29423205091574, 2.81018898925442], - + [-2.81018898925442, 1.29423205091574, -1.37278583181113, 2.38892435916370, \ 1.18892673484083e-11,-2.87551740597844, 2.11206121431525,-4.29790816425547], - + [1.94165942553537, -0.869448098330221, 0.845022556506714,-1.13535801687865,\ - 2.37818723350641, 6.49524878326702e-12,-4.54358506455201,7.56928981931106],\ - + 2.37818723350641, 6.49524878326702e-12,-4.54358506455201,7.56928981931106],\ + [-1.29768738831567, 0.573565414940005,-0.537039586158262,0.661157350899271,\ -1.07894468878725,2.80647579473469,-3.31498162253752e-11,-18.9375986071129], - + [0.500000000000000,-0.219957514771985,0.203284568901545,-0.243330712724289,\ 0.372150435728984, -0.792476681323880, 3.20991570302344, 14.0000000000226] ])) - + assert af.max(reference_d_Lp_xi - gvar.dLp_xi) < threshold + def test_volume_integral_flux(): ''' A test function to check the volume_integral_flux function in wave_equation @@ -290,7 +250,6 @@ def test_volume_integral_flux(): ''' threshold = 4e-8 - gvar.populateGlobalVariables(8) gvar.c = 1 referenceFluxIntegral = af.transpose(af.interop.np_to_af_array(np.array([ @@ -317,19 +276,19 @@ def test_volume_integral_flux(): 0.00236838757932, 0.00130167737191, 0.000588597708116, 0.00201663487667\ ]]))) - numerical_flux = wave_equation.volumeIntegralFlux(gvar.u[:, :, 0]) + numerical_flux = wave_equation.volume_integral_flux(gvar.u[:, :, 0]) assert (af.max(af.abs(numerical_flux - referenceFluxIntegral)) < threshold) def test_lax_friedrichs_flux(): ''' - A test function to test the laxFriedrichsFlux function in wave_equation + A test function to test the lax_friedrichs_flux function in wave_equation module. ''' threshold = 1e-14 - gvar.populateGlobalVariables(8, 10) + gvar.c = 1 - f_i = wave_equation.laxFriedrichsFlux(0) + f_i = wave_equation.lax_friedrichs_flux(gvar.u[:, :, 0]) analytical_lax_friedrichs_flux = gvar.u[-1, :, 0] assert af.max(af.abs(analytical_lax_friedrichs_flux - f_i)) < threshold @@ -340,7 +299,6 @@ def test_surface_term(): module using analytical Lax-Friedrichs flux. ''' threshold = 1e-13 - gvar.populateGlobalVariables(8, 10) gvar.c = 1 @@ -356,7 +314,7 @@ def test_surface_term(): analytical_surface_term = af.blas.matmul(L_p_1, analytical_f_i)\ - af.blas.matmul(L_p_minus1, analytical_f_i_minus1) - numerical_surface_term = (wave_equation.surface_term(0)) + numerical_surface_term = (wave_equation.surface_term(gvar.u[:, :, 0])) assert af.max(af.abs(analytical_surface_term - numerical_surface_term)) \ < threshold return analytical_surface_term @@ -368,14 +326,13 @@ def test_b_vector(): with the one returned by b_vector function in wave_equation module ''' threshold = 1e-13 - gvar.populateGlobalVariables(8) gvar.c = 1 u_n_A_matrix = af.blas.matmul(wave_equation.A_matrix(), gvar.u[:, :, 0]) - volume_integral_flux = wave_equation.volumeIntegralFlux(gvar.u[:, :, 0]) + volume_integral_flux = wave_equation.volume_integral_flux(gvar.u[:, :, 0]) surface_term = test_surface_term() b_vector_analytical = u_n_A_matrix + (volume_integral_flux -\ (surface_term)) * gvar.delta_t - b_vector_array = wave_equation.b_vector(0) + b_vector_array = wave_equation.b_vector(gvar.u[:, :, 0]) assert (b_vector_analytical - b_vector_array) < threshold From 8d384577563a51c93019acb762e46f0aaba5fd5b Mon Sep 17 00:00:00 2001 From: Balavarun5 Date: Fri, 8 Sep 2017 13:38:33 +0530 Subject: [PATCH 15/19] Minor edits in documentation of functions. --- code/app/global_variables.py | 6 +++--- code/app/wave_equation.py | 11 ++++++----- code/unit_test/test_waveEqn.py | 18 +++++++++--------- 3 files changed, 18 insertions(+), 17 deletions(-) diff --git a/code/app/global_variables.py b/code/app/global_variables.py index bf3f424..645e69f 100644 --- a/code/app/global_variables.py +++ b/code/app/global_variables.py @@ -10,8 +10,8 @@ # This module is used to change the parameters of the simulation. -# The functions to calculate Legendre-Gauss-Lobatto points as well -# as the Lobatto weights used for integration are given below. +# The functions to calculate Legendre-Gauss-Lobatto points and +# the Lobatto weights used for integration are given below. def LGL_points(N): ''' @@ -140,7 +140,7 @@ def lobatto_weight_function(n, x): delta_t = delta_x / (20 * c) # Array of timesteps seperated by delta_t. -time = utils.linspace(0, int(total_time / delta_t) * delta_t, +time = utils.linspace(0, int(total_time / delta_t) * delta_t, int(total_time / delta_t)) # Initializing the amplitudes. Change u_init to required initial conditions. diff --git a/code/app/wave_equation.py b/code/app/wave_equation.py index a28919a..b50cb08 100644 --- a/code/app/wave_equation.py +++ b/code/app/wave_equation.py @@ -236,9 +236,9 @@ def surface_term(u_t_n): Parameters ---------- - u_t_n : arrayfire.Array [N M 1 1] - The timestep at which the surface term is to be calculated. - + u_t_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. Returns ------- surface_term : arrayfire.Array [N_LGL N_Elements 1 1] @@ -272,8 +272,9 @@ def b_vector(u_t_n): Parameters ---------- - t_n : double - + u_t_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. Returns ------- b_vector_array : arrayfire.Array diff --git a/code/unit_test/test_waveEqn.py b/code/unit_test/test_waveEqn.py index da0704e..5310a82 100644 --- a/code/unit_test/test_waveEqn.py +++ b/code/unit_test/test_waveEqn.py @@ -46,7 +46,7 @@ def test_dx_dxi(): def test_dx_dxi_analytical(): ''' Test to check the dx_dxi_analytical in wave equation module for an element - and compare it with an analytical value + and compare it with an analytical value. ''' threshold = 1e-14 @@ -56,9 +56,9 @@ def test_dx_dxi_analytical(): assert check_analytical_dx_dxi -def test_lBasisArray(): +def test_lagrange_coeffs(): ''' - Function to test the lBasisArray function in global_variables module by + Function to test the lagrange_coeffs in global_variables module by passing 8 LGL points and comparing the numerically obtained basis function coefficients to analytically calculated ones. @@ -66,7 +66,7 @@ def test_lBasisArray(): --------- The link to the sage worksheet where the calculations were carried out. - https://cocalc.com/projects/1b7f404c-87ba-40d0-816c-2eba17466aa8/files\ + https://cocalc.com/projects/1b7f404c-87ba-40d0-816c-2eba17466aa8/files /PM_2_5/wave_equation/worksheets/l_basis_array.sagews ''' threshold = 6e-10 @@ -145,7 +145,7 @@ def test_lBasisArray(): assert af.sum(af.abs(basis_array_analytical - gvar.lagrange_coeffs)) < threshold -def test_Integral_Li_Lp(): +def test_integral_Li_Lp(): ''' Test function to check the A_matrix function in wave_equation module. @@ -195,13 +195,13 @@ def test_A_matrix(): def test_dLp_xi(): ''' - Test function to check the dLp_xi calculated in gvar mdule with a + Test function to check the dLp_xi calculated in gvar module with a numerically obtained one. Refrence -------- The link to the sage worksheet where the calculations were carried out. - https://cocalc.com/projects/1b7f404c-87ba-40d0-816c-2eba17466aa8/files\ + https://cocalc.com/projects/1b7f404c-87ba-40d0-816c-2eba17466aa8/files /PM_2_5/wave_equation/worksheets/dLp_xi.sagews ''' threshold = 4e-11 @@ -245,7 +245,7 @@ def test_volume_integral_flux(): --------- The link to the sage worksheet where the calculations were caried out is given below. - https://cocalc.com/projects/1b7f404c-87ba-40d0-816c-2eba17466aa8/files\ + https://cocalc.com/projects/1b7f404c-87ba-40d0-816c-2eba17466aa8/files /PM_2_5/wave_equation/worksheets/volume_integral_flux.sagews ''' @@ -323,7 +323,7 @@ def test_surface_term(): def test_b_vector(): ''' A test function to check the b vector obtained analytically and compare it - with the one returned by b_vector function in wave_equation module + with the one returned by b_vector function in wave_equation module. ''' threshold = 1e-13 gvar.c = 1 From 457f42dbc31715928691bad96a609f415fcc7c16 Mon Sep 17 00:00:00 2001 From: Balavarun5 Date: Fri, 8 Sep 2017 13:41:30 +0530 Subject: [PATCH 16/19] Change in note for developers in readme.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 05d2b32..d3c0d61 100644 --- a/README.md +++ b/README.md @@ -27,4 +27,4 @@ python3 main.py * **Mani Chandra** - [GitHub Profile](https://github.com/mchandra) ## Note for developers: -* Use tab spaces for indentation. +* Use spaces for indentation. From f82df886b9d035daa9f8044afbe1340dd5c07230 Mon Sep 17 00:00:00 2001 From: Balavarun5 Date: Fri, 15 Sep 2017 17:50:01 +0530 Subject: [PATCH 17/19] Temporary commit. Not all changes have been implemented. Code is to be refactored to include `Integrate()` function. --- README.md | 4 +- code/app/global_variables.py | 150 ------------------------ code/app/lagrange.py | 122 ++++++++++++++++--- code/app/params.py | 78 ++++++++++++ code/app/wave_equation.py | 163 +++++++++++++------------- code/main.py | 8 +- code/unit_test/lagrange_polynomial.py | 7 +- code/unit_test/test_waveEqn.py | 126 ++++++-------------- code/utils/utils.py | 59 ++++++++++ 9 files changed, 368 insertions(+), 349 deletions(-) delete mode 100644 code/app/global_variables.py create mode 100644 code/app/params.py diff --git a/README.md b/README.md index d3c0d61..a8d8527 100644 --- a/README.md +++ b/README.md @@ -22,9 +22,9 @@ python3 main.py ## Authors -* **Balavarun P** - [GitHub Profile](https://github.com/Balavarun5) +* **Balavarun P** - [GitHub Profile](https://github.com/Balavarun5) * **Aman Abhishek Tiwari** - [GitHub Profile](https://github.com/amanabt) -* **Mani Chandra** - [GitHub Profile](https://github.com/mchandra) +* **Mani Chandra** - [GitHub Profile](https://github.com/mchandra) ## Note for developers: * Use spaces for indentation. diff --git a/code/app/global_variables.py b/code/app/global_variables.py deleted file mode 100644 index 6871ee6..0000000 --- a/code/app/global_variables.py +++ /dev/null @@ -1,150 +0,0 @@ -import numpy as np -import arrayfire as af -af.set_backend('opencl') -from scipy import special as sp -from app import lagrange -from app import wave_equation -from utils import utils -import math - - -# This module is used to change the parameters of the simulation. - -# The functions to calculate Legendre-Gauss-Lobatto points and -# the Lobatto weights used for integration are given below. - -def LGL_points(N): - ''' - Calculates and returns the LGL points for a given N, These are roots of - the formula. - :math: `(1 - xi ** 2) P_{n - 1}'(xi) = 0` - Legendre polynomials satisfy the recurrence relation - :math:`(1 - x ** 2) P_n' (x) = -n x P_n(x) + n P_{n - 1} (x)` - - Parameters - ---------- - N : Int - Number of LGL nodes required - - Returns - ------- - lgl : arrayfire.Array [N 1 1 1] - An arrayfire array consisting of the Lagrange-Gauss-Lobatto Nodes. - - Reference - --------- - http://mathworld.wolfram.com/LobattoQuadrature.html - ''' - xi = np.poly1d([1, 0]) - legendre_N_minus_1 = N * (xi * sp.legendre(N - 1) - sp.legendre(N)) - lgl_points = legendre_N_minus_1.r - lgl_points.sort() - lgl_points = af.np_to_af_array(lgl_points) - - return lgl_points - -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 [N 1 1 1] - 1D array of points where weight function is to be calculated. - - - Returns - ------- - - gauss_lobatto_weights : 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) - - gauss_lobatto_weights = (2 / (n * (n - 1)) / (P(x))**2) - - return gauss_lobatto_weights - - - - -# The domain of the function which is of interest. -x_nodes = af.np_to_af_array(np.array([-1., 1.])) - -# The number of LGL points into which an element is split -N_LGL = 8 - -# Number of elements the domain is to be divided into -N_Elements = 10 - -# The speed of the wave. -c = 1 - -# The total time for which the simulation is to be carried out. -total_time = 1 - -# The c_lax to be used in the Lax-Friedrichs flux. -c_lax = 1 - -# Array containing the LGL points in xi space. -xi_LGL = LGL_points(N_LGL) - - -# Array of the lobatto weights used during integration. -lobatto_weights = af.np_to_af_array(lobatto_weight_function(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)) - -# Refer corresponding functions. -dLp_xi = lagrange.dLp_xi_LGL(lagrange_coeffs) -lagrange_basis_value = lagrange.lagrange_basis_function(lagrange_coeffs) - -# Obtaining an array consisting of the LGL points mapped onto the 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_mesh_nodes = utils.linspace(af.sum(x_nodes[0]), - af.sum(x_nodes[1]), - N_Elements + 1) - -element_array = af.transpose(af.interop.np_to_af_array(np_element_array)) -element_LGL = wave_equation.mapping_xi_to_x(af.transpose(element_array), xi_LGL) - - -# The minimum distance between 2 mapped LGL points. -delta_x = af.min((element_LGL - af.shift(element_LGL, 1, 0))[1:, :]) - -# The value of time-step. -delta_t = delta_x / (20 * c) - -# Array of timesteps seperated by delta_t. -time = utils.linspace(0, int(total_time / delta_t) * delta_t, - int(total_time / delta_t)) - -# Initializing the amplitudes. Change u_init to required initial conditions. -u_init = np.e ** (-(element_LGL) ** 2 / 0.4 ** 2) -u = af.constant(0, N_LGL, N_Elements, time.shape[0],\ - dtype = af.Dtype.f64) -u[:, :, 0] = u_init diff --git a/code/app/lagrange.py b/code/app/lagrange.py index 6badfc6..e8cf6ef 100644 --- a/code/app/lagrange.py +++ b/code/app/lagrange.py @@ -1,8 +1,76 @@ #! /usr/bin/env python3 + import numpy as np import arrayfire as af +from scipy import special as sp + from utils import utils -from app import global_variables as gvar +from app import params + +def LGL_points(N): + ''' + Calculates : math: `N` Legendre-Gauss-Lobatto (LGL) points. + LGL points are the roots of the polynomial + :math: `(1 - \\xi ** 2) P_{n - 1}'(\\xi) = 0` + Where :math: `P_{n}(\\xi)` are the Legendre polynomials. + This function finds the roots of the above polynomial. + + Parameters + ---------- + N : int + Number of LGL nodes required + + Returns + ------- + lgl : arrayfire.Array [N 1 1 1] + An array consisting of the Lagrange-Gauss-Lobatto Nodes. + + Reference + --------- + `https://goo.gl/KdG2Sv` + ''' + xi = np.poly1d([1, 0]) + legendre_N_minus_1 = N * (xi * sp.legendre(N - 1) - sp.legendre(N)) + lgl_points = legendre_N_minus_1.r + lgl_points.sort() + lgl_points = af.np_to_af_array(lgl_points) + + return lgl_points + + +def lobatto_weights(n, x): + ''' + Calculates 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 [N 1 1 1] + 1D array of points where weight function is to be calculated. + + + Returns + ------- + gauss_lobatto_weights : arrayfire.Array + An array of lobatto weight functions for + the given :math: `x` points and index. + Reference + --------- + Gauss-Lobatto weights Wikipedia link- + `https://goo.gl/o7WE4K` + ''' + P = sp.legendre(n - 1) + + gauss_lobatto_weights = (2 / (n * (n - 1)) / (P(x))**2) + + return gauss_lobatto_weights def lagrange_basis_coeffs(x): @@ -10,15 +78,14 @@ def lagrange_basis_coeffs(x): A function to get the coefficients of the Lagrange basis polynomials for a given set of x nodes. - This function calculates the Lagrange basis - polynomials by this formula: + It calculates the Lagrange basis polynomials using the formula: :math:: `L_i = \\prod_{m = 0, m \\notin i}^{N - 1}\\frac{(x - x_m)}{(x_i - x_m)}` Parameters ---------- x : numpy.array - A numpy array consisting of the :math: `x` nodes using which the + An array consisting of the :math: `x` nodes using which the lagrange basis functions need to be evaluated. Returns @@ -66,18 +133,18 @@ def lagrange_basis(i, x): Array of values of the :math: `i^{th}` Lagrange basis calculated at the given :math: `x` coordinates. ''' - x_tile = af.transpose(af.tile(x, 1, gvar.N_LGL)) - power = utils.linspace((gvar.N_LGL-1), 0, gvar.N_LGL) + x_tile = af.transpose(af.tile(x, 1, params.N_LGL)) + power = utils.linspace((params.N_LGL-1), 0, params.N_LGL) power_tile = af.tile(power, 1, x.shape[0]) x_pow = af.arith.pow(x_tile, power_tile) - l_xi_j = af.blas.matmul(gvar.lBasisArray[i], x_pow) + l_xi_j = af.blas.matmul(params.lBasisArray[i], x_pow) return l_xi_j def lagrange_basis_function(lagrange_coeff_array): ''' - Funtion which calculates the value of lagrange basis functions over LGL + Funtion to calculate the value of lagrange basis functions over LGL nodes. Returns @@ -86,11 +153,11 @@ def lagrange_basis_function(lagrange_coeff_array): The value of lagrange basis functions calculated over the LGL nodes. ''' - xi_tile = af.transpose(af.tile(gvar.xi_LGL, 1, gvar.N_LGL)) - power = af.flip(af.range(gvar.N_LGL)) - power_tile = af.tile(power, 1, gvar.N_LGL) + xi_tile = af.transpose(af.tile(params.xi_LGL, 1, params.N_LGL)) + power = af.flip(af.range(params.N_LGL)) + power_tile = af.tile(power, 1, params.N_LGL) xi_pow = af.arith.pow(xi_tile, power_tile) - index = af.range(gvar.N_LGL) + index = af.range(params.N_LGL) L_i = af.blas.matmul(lagrange_coeff_array[index], xi_pow) return L_i @@ -109,13 +176,36 @@ def dLp_xi_LGL(lagrange_coeff_array): along dimension 1 and p varies along second dimension. ''' differentiation_coeffs = (af.transpose(af.flip(af.tile\ - (af.range(gvar.N_LGL), 1, gvar.N_LGL))) + (af.range(params.N_LGL), 1, params.N_LGL))) * lagrange_coeff_array)[:, :-1] - nodes_tile = af.transpose(af.tile(gvar.xi_LGL, 1, gvar.N_LGL - 1)) - power_tile = af.flip(af.tile(af.range(gvar.N_LGL - 1), 1, gvar.N_LGL)) + nodes_tile = af.transpose(af.tile(params.xi_LGL, 1, params.N_LGL - 1)) + power_tile = af.flip(af.tile(af.range(params.N_LGL - 1), 1, params.N_LGL)) nodes_power_tile = af.pow(nodes_tile, power_tile) dLp_xi = af.blas.matmul(differentiation_coeffs, nodes_power_tile) - return dLp_xi \ No newline at end of file + return dLp_xi + + + +def Integrate(integrand, N_quad, scheme): + ''' + ''' + if (scheme == 'gauss_quadrature'): + gaussian_nodes = gauss_nodes(N_quad) + gauss_weights = af.np_to_af_array(np.zeros([N_quad])) + + for i in range(0, N_quad): + gauss_weights[i] = gaussian_weights(N_quad, i) + value_at_gauss_nodes = af.np_to_af_array(integrand(gaussian_nodes)) + Integral = af.sum(value_at_gauss_nodes * gauss_weights) + + if (scheme == 'lobatto_quadrature'): + lobatto_nodes = LGL_points(N_quad) + lobatto_weights = lobatto_weight_function(N_quad) + value_at_lobatto_nodes = af.np_to_af_array(integrand(lobatto_nodes)) + Integral = af.sum(value_at_lobatto_nodes * lobatto_weights) + + + return Integral diff --git a/code/app/params.py b/code/app/params.py new file mode 100644 index 0000000..2b53c28 --- /dev/null +++ b/code/app/params.py @@ -0,0 +1,78 @@ +#! /usr/bin/env python3 + +import numpy as np +import arrayfire as af +af.set_backend('opencl') + +from app import lagrange +from utils import utils +from app import wave_equation + + +# The domain of the function. +x_nodes = af.np_to_af_array(np.array([-1., 1.])) + +# The number of LGL points into which an element is split. +N_LGL = 8 + +# Number of elements the domain is to be divided into. +N_Elements = 10 + +# The scheme to be used for integration. Values are either +# 'gauss_quadrature' or 'lobatto_quadrature' +scheme = 'lobatto_quadrature' + +# Wave speed. +c = 1 + +# The total time for which the wave is to be evolved by the simulation. +total_time = 1 + +# The c_lax to be used in the Lax-Friedrichs flux. +c_lax = 1 + +# Array containing the LGL points in xi space. +xi_LGL = lagrange.LGL_points(N_LGL) + +# Array of the lobatto weights used during integration. +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)) + +# Refer corresponding functions. +dLp_xi = lagrange.dLp_xi_LGL(lagrange_coeffs) +lagrange_basis_value = lagrange.lagrange_basis_function(lagrange_coeffs) + + +# Obtaining an array consisting of the LGL points mapped onto the 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_mesh_nodes = utils.linspace(af.sum(x_nodes[0]), + af.sum(x_nodes[1]), N_Elements + 1) + +element_array = af.transpose(af.interop.np_to_af_array(np_element_array)) +element_LGL = wave_equation.mapping_xi_to_x(af.transpose(element_array), xi_LGL) + + +# The minimum distance between 2 mapped LGL points. +delta_x = af.min((element_LGL - af.shift(element_LGL, 1, 0))[1:, :]) + +# The value of time-step. +delta_t = delta_x / (20 * c) + +# Array of timesteps seperated by delta_t. +time = utils.linspace(0, int(total_time / delta_t) * delta_t, + int(total_time / delta_t)) + +# Initializing the amplitudes. Change u_init to required initial conditions. +u_init = np.e ** (-(element_LGL) ** 2 / 0.4 ** 2) +u = af.constant(0, N_LGL, N_Elements, time.shape[0],\ + dtype = af.Dtype.f64) +u[:, :, 0] = u_init diff --git a/code/app/wave_equation.py b/code/app/wave_equation.py index 0a0da98..f4911c3 100644 --- a/code/app/wave_equation.py +++ b/code/app/wave_equation.py @@ -3,39 +3,42 @@ import arrayfire as af af.set_backend('opencl') import numpy as np -from app import global_variables as gvar -from app import lagrange -from utils import utils from matplotlib import pyplot as plt import pylab as pl from tqdm import trange -plt.rcParams['figure.figsize'] = 9.6, 6. -plt.rcParams['figure.dpi'] = 100 -plt.rcParams['image.cmap'] = 'jet' -plt.rcParams['lines.linewidth'] = 1.5 -plt.rcParams['font.family'] = 'serif' -plt.rcParams['font.weight'] = 'bold' -plt.rcParams['font.size'] = 20 -plt.rcParams['font.sans-serif'] = 'serif' -plt.rcParams['text.usetex'] = True -plt.rcParams['axes.linewidth'] = 1.5 -plt.rcParams['axes.titlesize'] = 'medium' -plt.rcParams['axes.labelsize'] = 'medium' +from app import params +from app import lagrange +from utils import utils + +plt.rcParams['figure.figsize' ] = 9.6, 6. +plt.rcParams['figure.dpi' ] = 100 +plt.rcParams['image.cmap' ] = 'jet' +plt.rcParams['lines.linewidth' ] = 1.5 +plt.rcParams['font.family' ] = 'serif' +plt.rcParams['font.weight' ] = 'bold' +plt.rcParams['font.size' ] = 20 +plt.rcParams['font.sans-serif' ] = 'serif' +plt.rcParams['text.usetex' ] = True +plt.rcParams['axes.linewidth' ] = 1.5 +plt.rcParams['axes.titlesize' ] = 'medium' +plt.rcParams['axes.labelsize' ] = 'medium' plt.rcParams['xtick.major.size'] = 8 plt.rcParams['xtick.minor.size'] = 4 -plt.rcParams['xtick.major.pad'] = 8 -plt.rcParams['xtick.minor.pad'] = 8 -plt.rcParams['xtick.color'] = 'k' -plt.rcParams['xtick.labelsize'] = 'medium' -plt.rcParams['xtick.direction'] = 'in' +plt.rcParams['xtick.major.pad' ] = 8 +plt.rcParams['xtick.minor.pad' ] = 8 +plt.rcParams['xtick.color' ] = 'k' +plt.rcParams['xtick.labelsize' ] = 'medium' +plt.rcParams['xtick.direction' ] = 'in' plt.rcParams['ytick.major.size'] = 8 plt.rcParams['ytick.minor.size'] = 4 -plt.rcParams['ytick.major.pad'] = 8 -plt.rcParams['ytick.minor.pad'] = 8 -plt.rcParams['ytick.color'] = 'k' -plt.rcParams['ytick.labelsize'] = 'medium' -plt.rcParams['ytick.direction'] = 'in' +plt.rcParams['ytick.major.pad' ] = 8 +plt.rcParams['ytick.minor.pad' ] = 8 +plt.rcParams['ytick.color' ] = 'k' +plt.rcParams['ytick.labelsize' ] = 'medium' +plt.rcParams['ytick.direction' ] = 'in' + + def mapping_xi_to_x(x_nodes, xi): @@ -52,7 +55,6 @@ def mapping_xi_to_x(x_nodes, xi): xi : numpy.float64 Value of :math: `\\xi`coordinate for which the corresponding :math: `x` coordinate is to be found. -. Returns ------- @@ -62,7 +64,6 @@ def mapping_xi_to_x(x_nodes, xi): N_0 = (1 - xi) / 2 N_1 = (1 + xi) / 2 - N0_x0 = af.broadcast(utils.multiply, N_0, x_nodes[0]) N1_x1 = af.broadcast(utils.multiply, N_1, x_nodes[1]) @@ -70,7 +71,8 @@ def mapping_xi_to_x(x_nodes, xi): return x - + + def dx_dxi_numerical(x_nodes, xi): ''' Differential :math: `\\frac{dx}{d \\xi}` calculated by central differential @@ -91,7 +93,6 @@ def dx_dxi_numerical(x_nodes, xi): :math:`\\frac{dx}{d \\xi}`. ''' dxi = 1e-7 - x2 = mapping_xi_to_x(x_nodes, xi + dxi) x1 = mapping_xi_to_x(x_nodes, xi - dxi) @@ -103,7 +104,7 @@ def dx_dxi_numerical(x_nodes, xi): def dx_dxi_analytical(x_nodes, xi): ''' The analytical result for :math:`\\frac{dx}{d \\xi}` for a 1D element is - :math: `\\frac{x_1 - x_0}{2} + :math: `\\frac{x_1 - x_0}{2}` Parameters ---------- x_nodes : arrayfire.Array @@ -112,14 +113,15 @@ def dx_dxi_analytical(x_nodes, xi): Returns ------- analytical_dx_dxi : arrayfire.Array - The analytical solution of - \\frac{dx}{d\\xi} for an element. + The analytical solution of :math: + `\\frac{dx}{d\\xi}` for an element. ''' analytical_dx_dxi = (x_nodes[1] - x_nodes[0]) / 2 return analytical_dx_dxi + def A_matrix(): ''' Calculates A matrix whose elements :math:`A_{p i}` are given by @@ -129,10 +131,7 @@ def A_matrix(): The integration is carried out using Gauss-Lobatto quadrature. Full description of the algorithm can be found here- - https://cocalc.com/projects/1b7f404c-87ba-40d0-816c-2eba17466aa8/files/ - PM_2_5/wave_equation/documents/wave_equation_report - /A_matrix.pdf?session=default - + `https://goo.gl/Cw6tnw` Returns ------- A_matrix : arrayfire.Array @@ -140,11 +139,11 @@ def A_matrix(): obtained by LGL points, using Gauss-Lobatto quadrature method using :math: `N_LGL` points. ''' - dx_dxi = dx_dxi_numerical((gvar.element_mesh_nodes[0 : 2]), - gvar.xi_LGL) - dx_dxi_tile = af.tile(dx_dxi, 1, gvar.N_LGL) - identity_matrix = af.identity(gvar.N_LGL, gvar.N_LGL, dtype = af.Dtype.f64) - A_matrix = af.broadcast(utils.multiply, gvar.lobatto_weights, + 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 return A_matrix @@ -152,27 +151,27 @@ def A_matrix(): def flux_x(u): ''' - A function which calcultes and returns the value of flux for a given - wave function u. :math:`f(u) = c u^k` + A function which returns the value of flux for a given wave function u. + :math:`f(u) = c u^k` Parameters ---------- - u : arrayfire.Array [N 1 1 1] - A 1-D array which contains the value of wave function. + u : arrayfire.Array + A 1-D array which contains the value of wave function. Returns ------- flux : arrayfire.Array The value of the flux for given u. ''' - flux = gvar.c * u + flux = params.c * u return flux def volume_integral_flux(u): ''' - A function to calculate the volume integral of flux in the wave equation. + 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. @@ -180,10 +179,9 @@ def volume_integral_flux(u): Parameters ---------- - - u : arrayfire.Array [N_LGL N_Elements 1 1] - A 1-D array containing the value of the wave function at - the mapped LGL nodes in the element. + u : arrayfire.Array [N_LGL N_Elements 1 1] + A 2D array containing the value of the wave function at + the mapped LGL nodes in all the elements. Returns ------- @@ -192,9 +190,8 @@ def volume_integral_flux(u): for various lagrange basis functions. ''' - - dLp_xi = gvar.dLp_xi - weight_tile = af.transpose(af.tile(gvar.lobatto_weights, 1, gvar.N_LGL)) + 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 @@ -203,7 +200,7 @@ def volume_integral_flux(u): return flux_integral -def lax_friedrichs_flux(u_t_n): +def lax_friedrichs_flux(u_n): ''' A function which calculates the lax-friedrichs_flux :math:`f_i` using. :math:`f_i = \\frac{F(u^{i + 1}_0) + F(u^i_{N_{LGL} - 1})}{2} - \frac @@ -211,35 +208,36 @@ def lax_friedrichs_flux(u_t_n): Parameters ---------- - u_t_n : arrayfire.Array [N_LGL N_Elements 1 1] + 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. ''' - u_iplus1_0 = af.shift(u_t_n[0, :], 0, -1) - u_i_N_LGL = u_t_n[-1, :] + u_iplus1_0 = af.shift(u_n[0, :], 0, -1) + u_i_N_LGL = u_n[-1, :] flux_iplus1_0 = flux_x(u_iplus1_0) flux_i_N_LGL = flux_x(u_i_N_LGL) boundary_flux = (flux_iplus1_0 + flux_i_N_LGL) / 2 \ - - gvar.c_lax * (u_iplus1_0 - u_i_N_LGL) / 2 + - params.c_lax * (u_iplus1_0 - u_i_N_LGL) / 2 return boundary_flux -def surface_term(u_t_n): +def surface_term(u_n): ''' A function which is used to calculate 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 gvar module. + from params module. Parameters ---------- - u_t_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. + 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. + Returns ------- surface_term : arrayfire.Array [N_LGL N_Elements 1 1] @@ -253,13 +251,12 @@ def surface_term(u_t_n): --------- Link to PDF describing the algorithm to obtain the surface term. - https://cocalc.com/projects/1b7f404c-87ba-40d0-816c-2eba17466aa8/files - /PM\_2\_5/wave\_equation/documents/surface\_term/surface\_term.pdf + `https://goo.gl/Nhhgzx` ''' - L_p_minus1 = gvar.lagrange_basis_value[:, 0] - L_p_1 = gvar.lagrange_basis_value[:, -1] - f_i = lax_friedrichs_flux(u_t_n) + 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) @@ -267,23 +264,23 @@ def surface_term(u_t_n): return surface_term -def b_vector(u_t_n): +def b_vector(u_n): ''' - A function which returns the b vector for N_Elements number of elements. + Calculates the b vector for N_Elements number of elements. Parameters ---------- - u_t_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. + 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. Returns ------- b_vector_array : arrayfire.Array ''' - volume_integral = volume_integral_flux(u_t_n) - Surface_term = surface_term(u_t_n) - b_vector_array = gvar.delta_t * (volume_integral - Surface_term) + volume_integral = volume_integral_flux(u_n) + Surface_term = surface_term(u_n) + b_vector_array = params.delta_t * (volume_integral - Surface_term) return b_vector_array @@ -297,10 +294,10 @@ def time_evolution(): ''' A_inverse = af.inverse(A_matrix()) - element_LGL = gvar.element_LGL - delta_t = gvar.delta_t - amplitude = gvar.u - time = gvar.time + element_LGL = params.element_LGL + delta_t = params.delta_t + amplitude = params.u + time = params.time for t_n in trange(0, time.shape[0] - 1): @@ -314,7 +311,7 @@ def time_evolution(): if t_n % 100 == 0: fig = plt.figure() - x = gvar.element_LGL + x = params.element_LGL y = amplitude[:, :, t_n] plt.plot(x, y) @@ -324,4 +321,4 @@ def time_evolution(): fig.savefig('results/1D_Wave_images/%04d' %(t_n / 100) + '.png') plt.close('all') - return \ No newline at end of file + return diff --git a/code/main.py b/code/main.py index fbd85e6..4b3a1ff 100644 --- a/code/main.py +++ b/code/main.py @@ -1,9 +1,10 @@ #! /usr/bin/env python3 import arrayfire as af -af.set_backend('opencl') import numpy as np -from app import global_variables as gvar +af.set_backend('opencl') + +from app import params from unit_test import test_waveEqn from app import wave_equation from app import lagrange @@ -11,5 +12,4 @@ if __name__ == '__main__': - - wave_equation.time_evolution() \ No newline at end of file + wave_equation.time_evolution() diff --git a/code/unit_test/lagrange_polynomial.py b/code/unit_test/lagrange_polynomial.py index a5fedd5..5cd42ab 100644 --- a/code/unit_test/lagrange_polynomial.py +++ b/code/unit_test/lagrange_polynomial.py @@ -1,10 +1,13 @@ -from app import lagrange -from utils import utils +#! /usr/bin/env python3 + import math + import arrayfire as af import numpy as np from matplotlib import pyplot as plt +from app import lagrange +from utils import utils from app import global_variables as gvar plt.rcParams['figure.figsize'] = 9.6, 6. diff --git a/code/unit_test/test_waveEqn.py b/code/unit_test/test_waveEqn.py index 37fbea0..65ab090 100644 --- a/code/unit_test/test_waveEqn.py +++ b/code/unit_test/test_waveEqn.py @@ -1,14 +1,17 @@ +#! /usr/bin/env python3 +import math + import numpy as np import arrayfire as af -import math from matplotlib import pyplot as plt -from app import global_variables as gvar +af.set_backend('opencl') + +from app import params from app import lagrange from app import wave_equation from utils import utils -af.set_backend('opencl') - +# This test uses the initial paramters N_LGL = 8, N_Elements = 10 def test_mapping_xi_to_x(): ''' @@ -22,26 +25,24 @@ def test_mapping_xi_to_x(): test_element_nodes = af.interop.np_to_af_array(np.array([7, 11])) test_xi = 0 analytical_x_value = 9 - numerical_x_value = wave_equation.mapping_xi_to_x(test_element_nodes, test_xi) assert af.abs(analytical_x_value - numerical_x_value) <= threshold + def test_dx_dxi(): ''' A Test function to check the dx_xi function in wave_equation module by passing nodes of an element and using the LGL points. Analytically, the differential would be a constant. The check has a tolerance 1e-7. ''' - threshold = 1e-14 - nodes = np.array([7, 10], dtype = np.float64) test_nodes = af.interop.np_to_af_array(nodes) analytical_dx_dxi = 1.5 check_dx_dxi = (af.statistics.mean(wave_equation.dx_dxi_numerical - (test_nodes,gvar.xi_LGL)) - analytical_dx_dxi) <= threshold + (test_nodes,params.xi_LGL)) - analytical_dx_dxi) <= threshold assert check_dx_dxi @@ -59,7 +60,6 @@ def test_dx_dxi_analytical(): assert check_analytical_dx_dxi - def test_lagrange_coeffs(): ''' Function to test the lagrange_coeffs in global_variables module by @@ -70,10 +70,8 @@ def test_lagrange_coeffs(): --------- The link to the sage worksheet where the calculations were carried out. - https://cocalc.com/projects/1b7f404c-87ba-40d0-816c-2eba17466aa8/files - /PM_2_5/wave_equation/worksheets/l_basis_array.sagews + `https://goo.gl/6EFX5S` ''' - threshold = 6e-10 basis_array_analytical = np.zeros([8, 8]) @@ -147,71 +145,18 @@ def test_lagrange_coeffs(): basis_array_analytical = af.interop.np_to_af_array(basis_array_analytical) - - assert af.sum(af.abs(basis_array_analytical - gvar.lagrange_coeffs)) < threshold - - -def test_integral_Li_Lp(): - ''' - Test function to check the A_matrix function in wave_equation module. - - Obtaining the A_matrix from the function and setting the value of - all elements above a certain threshold to be 1 and plotting it. - ''' - threshold = 1e-5 - - A_matrix_structure = np.zeros([gvar.N_LGL, gvar.N_LGL]) - non_zero_indices = np.where(np.array(wave_equation.A_matrix()) > threshold) - - A_matrix_structure[non_zero_indices] = 1. - - plt.gca().invert_yaxis() - plt.contourf(A_matrix_structure, 100, cmap = 'Blues') - plt.axes().set_aspect('equal') - plt.colorbar() - plt.show() - - return - - -def test_A_matrix(): - ''' - Test function to check the A matrix obtained from wave_equation module with - one obtained by numerical integral solvers. - - NOTE - ---- - A diagonal A-matrix is used as a reference in this unit test. - The A matrix calculated analytically gives a different matrix. - - ''' - - threshold = 3e-10 - - reference_A_matrix = af.tile(gvar.lobatto_weights, 1, gvar.N_LGL)\ - * af.identity(gvar.N_LGL, gvar.N_LGL, dtype = af.Dtype.f64)\ - * af.mean(wave_equation.dx_dxi_numerical((gvar.element_mesh_nodes[0 : 2])\ - ,gvar.xi_LGL)) - - test_A_matrix = wave_equation.A_matrix() - error_array = af.abs(reference_A_matrix - test_A_matrix) - - - print(af.max(error_array)) - - assert af.algorithm.max(error_array) < threshold + assert af.sum(af.abs(basis_array_analytical - params.lagrange_coeffs)) < threshold def test_dLp_xi(): ''' - Test function to check the dLp_xi calculated in gvar module with a + Test function to check the dLp_xi calculated in params module with a numerically obtained one. Refrence -------- The link to the sage worksheet where the calculations were carried out. - https://cocalc.com/projects/1b7f404c-87ba-40d0-816c-2eba17466aa8/files - /PM_2_5/wave_equation/worksheets/dLp_xi.sagews + `https://goo.gl/Xxp3PT` ''' threshold = 4e-11 @@ -242,7 +187,8 @@ def test_dLp_xi(): 0.372150435728984, -0.792476681323880, 3.20991570302344, 14.0000000000226] ])) - assert af.max(reference_d_Lp_xi - gvar.dLp_xi) < threshold + assert af.max(reference_d_Lp_xi - params.dLp_xi) < threshold + def test_volume_integral_flux(): ''' @@ -253,13 +199,10 @@ def test_volume_integral_flux(): --------- The link to the sage worksheet where the calculations were caried out is given below. - - https://cocalc.com/projects/1b7f404c-87ba-40d0-816c-2eba17466aa8/files - /PM_2_5/wave_equation/worksheets/volume_integral_flux.sagews - + `https://goo.gl/5Mub8M` ''' threshold = 4e-8 - gvar.c = 1 + params.c = 1 referenceFluxIntegral = af.transpose(af.interop.np_to_af_array(np.array([ [-0.002016634876668093, -0.000588597708116113, -0.0013016773719126333,\ @@ -282,11 +225,10 @@ def test_volume_integral_flux(): [-0.102576031756, 0.0154359890788, 0.0209837936827, 0.019612124119, \ 0.0144355176966, 0.00882630935977, 0.00431252844519, 0.018969769374],\ [-0.0176615086879, 0.00344551201015 ,0.00432019709409, 0.00362050204766,\ - 0.00236838757932, 0.00130167737191, 0.000588597708116, 0.00201663487667\ ]]))) - numerical_flux = wave_equation.volume_integral_flux(gvar.u[:, :, 0]) + numerical_flux = wave_equation.volume_integral_flux(params.u[:, :, 0]) assert (af.max(af.abs(numerical_flux - referenceFluxIntegral)) < threshold) def test_lax_friedrichs_flux(): @@ -296,10 +238,10 @@ def test_lax_friedrichs_flux(): ''' threshold = 1e-14 - gvar.c = 1 + params.c = 1 - f_i = wave_equation.lax_friedrichs_flux(gvar.u[:, :, 0]) - analytical_lax_friedrichs_flux = gvar.u[-1, :, 0] + f_i = wave_equation.lax_friedrichs_flux(params.u[:, :, 0]) + analytical_lax_friedrichs_flux = params.u[-1, :, 0] assert af.max(af.abs(analytical_lax_friedrichs_flux - f_i)) < threshold @@ -309,22 +251,22 @@ def test_surface_term(): module using analytical Lax-Friedrichs flux. ''' threshold = 1e-13 - gvar.c = 1 + params.c = 1 - analytical_f_i = (gvar.u[-1, :, 0]) - analytical_f_i_minus1 = (af.shift(gvar.u[-1, :, 0], 0, 1)) + analytical_f_i = (params.u[-1, :, 0]) + analytical_f_i_minus1 = (af.shift(params.u[-1, :, 0], 0, 1)) - L_p_1 = af.constant(0, gvar.N_LGL, dtype = af.Dtype.f64) - L_p_1[gvar.N_LGL - 1] = 1 + L_p_1 = af.constant(0, params.N_LGL, dtype = af.Dtype.f64) + L_p_1[params.N_LGL - 1] = 1 - L_p_minus1 = af.constant(0, gvar.N_LGL, dtype = af.Dtype.f64) + L_p_minus1 = af.constant(0, params.N_LGL, dtype = af.Dtype.f64) L_p_minus1[0] = 1 analytical_surface_term = af.blas.matmul(L_p_1, analytical_f_i)\ - af.blas.matmul(L_p_minus1, analytical_f_i_minus1) - numerical_surface_term = (wave_equation.surface_term(gvar.u[:, :, 0])) + numerical_surface_term = (wave_equation.surface_term(params.u[:, :, 0])) assert af.max(af.abs(analytical_surface_term - numerical_surface_term)) \ < threshold return analytical_surface_term @@ -336,13 +278,13 @@ def test_b_vector(): with the one returned by b_vector function in wave_equation module. ''' threshold = 1e-13 - gvar.c = 1 + params.c = 1 - u_n_A_matrix = af.blas.matmul(wave_equation.A_matrix(), gvar.u[:, :, 0]) - volume_integral_flux = wave_equation.volume_integral_flux(gvar.u[:, :, 0]) + u_n_A_matrix = af.blas.matmul(wave_equation.A_matrix(), params.u[:, :, 0]) + volume_integral_flux = wave_equation.volume_integral_flux(params.u[:, :, 0]) surface_term = test_surface_term() b_vector_analytical = u_n_A_matrix + (volume_integral_flux -\ - (surface_term)) * gvar.delta_t - b_vector_array = wave_equation.b_vector(gvar.u[:, :, 0]) + (surface_term)) * params.delta_t + b_vector_array = wave_equation.b_vector(params.u[:, :, 0]) - assert (b_vector_analytical - b_vector_array) < threshold \ No newline at end of file + assert (b_vector_analytical - b_vector_array) < threshold diff --git a/code/utils/utils.py b/code/utils/utils.py index aa7e00a..1827914 100644 --- a/code/utils/utils.py +++ b/code/utils/utils.py @@ -4,6 +4,22 @@ def add(a, b): ''' + For broadcasting purposes, To sum two arrays of different + shapes, A function which can sum two variables is required. + + Parameters + ---------- + a : arrayfire.Array [N M 1 1] + One of the arrays which need to be broadcasted and summed. + + b : arrayfire.Array [1 M L 1] + One of the arrays which need to be broadcasted and summed. + + Returns + ------- + sum : arrayfire.Array + returns the sum of a and b. When used along with af.broadcast + can be used to sum different size arrays. ''' sum = a + b @@ -12,6 +28,24 @@ def add(a, b): def divide(a, b): ''' + + For broadcasting purposes, To divide two arrays of different + shapes, A function which can sum two variables is required. + + Parameters + ---------- + a : arrayfire.Array [N M 1 1] + One of the arrays which need to be broadcasted and divided. + + b : arrayfire.Array [1 M L 1] + One of the arrays which need to be broadcasted and divided. + + Returns + ------- + quotient : arrayfire.Array + returns the quotient a / b. When used along with af.broadcast + can be used to give quotient of two different size arrays + by dividing elements of the broadcasted array. ''' quotient = a / b @@ -20,6 +54,24 @@ def divide(a, b): def multiply(a, b): ''' + + For broadcasting purposes, To divide two arrays of different + shapes, A function which can sum two variables is required. + + Parameters + ---------- + a : arrayfire.Array [N M 1 1] + One of the arrays which need to be broadcasted and multiplying. + + b : arrayfire.Array [1 M L 1] + One of the arrays which need to be broadcasted and multiplying. + + Returns + ------- + product : arrayfire.Array + returns the quotient a / b. When used along with af.broadcast + can be used to give quotient of two different size arrays + by multiplying elements of the broadcasted array. ''' product = a* b @@ -29,6 +81,12 @@ def multiply(a, b): def linspace(start, end, number_of_points): ''' Linspace implementation using arrayfire. + + Returns + ------- + X : arrayfire.Array + An array which contains 'number_of_points' evenly spaced points + between 'start' and 'end' ''' X = af.range(number_of_points, dtype = af.Dtype.f64) d = (end - start) / (number_of_points - 1) @@ -36,3 +94,4 @@ def linspace(start, end, number_of_points): X = X + start return X + From e6b5c71f20facfb2a351fedfe719fc7ad292f9ef Mon Sep 17 00:00:00 2001 From: Balavarun5 Date: Mon, 18 Sep 2017 00:39:16 +0530 Subject: [PATCH 18/19] Code reffactored to use Integrate() function. Unit tests passing. Documentation and comments need to be added. --- code/app/lagrange.py | 118 +++++++++++++++++++++++++++++---- code/app/params.py | 27 +++++++- code/app/wave_equation.py | 58 +++++++++++----- code/main.py | 3 +- code/unit_test/test_waveEqn.py | 46 +++++++++++++ 5 files changed, 219 insertions(+), 33 deletions(-) diff --git a/code/app/lagrange.py b/code/app/lagrange.py index e8cf6ef..9b24c53 100644 --- a/code/app/lagrange.py +++ b/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 diff --git a/code/app/params.py b/code/app/params.py index 2b53c28..39031c1 100644 --- a/code/app/params.py +++ b/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) diff --git a/code/app/wave_equation.py b/code/app/wave_equation.py index f4911c3..d26616d 100644 --- a/code/app/wave_equation.py +++ b/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) diff --git a/code/main.py b/code/main.py index 4b3a1ff..68d0905 100644 --- a/code/main.py +++ b/code/main.py @@ -11,5 +11,4 @@ if __name__ == '__main__': - - wave_equation.time_evolution() + wave_equation.time_evolution() diff --git a/code/unit_test/test_waveEqn.py b/code/unit_test/test_waveEqn.py index 65ab090..9b58c4f 100644 --- a/code/unit_test/test_waveEqn.py +++ b/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(): ''' From f79ceb91ffe9b6ca96c412bee09330ee150796e0 Mon Sep 17 00:00:00 2001 From: Balavarun5 Date: Mon, 18 Sep 2017 14:24:41 +0530 Subject: [PATCH 19/19] Edited the documentation. Examples need to be added. --- README.md | 13 +- code/app/lagrange.py | 261 +++++++++++++++++++++------------ code/app/params.py | 55 +++---- code/app/wave_equation.py | 101 +++++++------ code/main.py | 4 +- code/unit_test/test_waveEqn.py | 47 +----- 6 files changed, 267 insertions(+), 214 deletions(-) diff --git a/README.md b/README.md index a8d8527..d67649f 100644 --- a/README.md +++ b/README.md @@ -13,13 +13,23 @@ cd DG_Maxwell * Arrayfire * Numpy * Matplotlib +* tqdm +* pytest ## Usage: ``` cd DG_Maxwell/code python3 main.py ``` +* The parameters of the simulation are stored in global_variables.py in + the app folder, These can be changed accordingly. + +* The images of the wave are stored in the folder 1D_wave_images folder. +* To stitch the images in the folder and obtain a video of the simulation, + use the command in the terminal - + `ffmpeg -f image2 -i %04d.png -vcodec mpeg4 -mbd rd -trellis 2 -cmp 2 -g 300 -pass 1 -r 25 -b 18000000 movie.mp4` + ## Authors * **Balavarun P** - [GitHub Profile](https://github.com/Balavarun5) @@ -27,4 +37,5 @@ python3 main.py * **Mani Chandra** - [GitHub Profile](https://github.com/mchandra) ## Note for developers: -* Use spaces for indentation. +* Use tab spaces for indentation. + diff --git a/code/app/lagrange.py b/code/app/lagrange.py index 9b24c53..454c87c 100644 --- a/code/app/lagrange.py +++ b/code/app/lagrange.py @@ -23,7 +23,7 @@ def LGL_points(N): Returns ------- lgl : arrayfire.Array [N 1 1 1] - An array consisting of the Lagrange-Gauss-Lobatto Nodes. + The Lagrange-Gauss-Lobatto Nodes. Reference --------- @@ -53,14 +53,14 @@ def lobatto_weights(n, x): Index for which lobatto weight function x : arrayfire.Array [N 1 1 1] - 1D array of points where weight function is to be calculated. + Points where weight function is to be calculated. Returns ------- - gauss_lobatto_weights : arrayfire.Array - An array of lobatto weight functions for - the given :math: `x` points and index. + gauss_lobatto_weights : numpy.ndarray [N 1 1 1] + Lobatto weight for the given :math: `x` + points and index `i`. Reference --------- Gauss-Lobatto weights Wikipedia link- @@ -72,10 +72,32 @@ def lobatto_weights(n, x): return gauss_lobatto_weights -def gauss_nodes(N): +def gauss_nodes(n): ''' + Calculates :math: `N` Gaussian nodes used for Integration by + Gaussia quadrature. + + Gaussian node :math: `x_i` is the `i^{th}` root of + + :math: `P_n(\\xi)` + Where :math: `P_{n}(\\xi)` are the Legendre polynomials. + + Parameters + ---------- + n : int + The number of Gaussian nodes required. + + Returns + ------- + gauss_nodes : numpy.ndarray + The Gauss nodes :math: `x_i`. + + Reference + --------- + A Wikipedia article about the Gauss-Legendre quadrature + `https://goo.gl/9gqLpe` ''' - legendre = sp.legendre(N) + legendre = sp.legendre(n) gauss_nodes = legendre.r gauss_nodes.sort() @@ -84,16 +106,22 @@ def gauss_nodes(N): def gaussian_weights(N, i): ''' - Returns the gaussian weights for :math: `N` Gaussian Nodes at index - :math: `i`. + Returns the gaussian weights :math:`w_i` for :math: `N` Gaussian Nodes + at index :math: `i`. They are given by + + :math: `w_i = \\frac{2}{(1 - x_i^2) P'n(x_i)}` + + Where :math:`x_i` are the Gaussian nodes and :math: `P_{n}(\\xi)` + are the Legendre polynomials. + Parameters ---------- - N : int - Number of Gaussian nodes for which the weight is t be calculated. + N : int + Number of Gaussian nodes for which the weight is t be calculated. - i : int - Index for which the Gaussian weight is required. + i : int + Index for which the Gaussian weight is required. Returns ------- @@ -103,15 +131,16 @@ def gaussian_weights(N, i): ''' gaussian_nodes = gauss_nodes(N) - gaussian_weight = 2 / ((1 - (gaussian_nodes[i]) ** 2) * (np.polyder(sp.legendre(N))(gaussian_nodes[i])) ** 2) + 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. + A function to get the analytical form and thr coefficients of + Lagrange basis polynomials evaluated using x nodes. It calculates the Lagrange basis polynomials using the formula: :math:: @@ -119,17 +148,19 @@ def lagrange_polynomials(x): Parameters ---------- - x : numpy.array - An array consisting of the :math: `x` nodes using which the + x : numpy.array [N_LGL 1 1 1] + Contains the :math: `x` nodes using which the lagrange basis functions need to be evaluated. Returns ------- 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. + 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 @@ -152,47 +183,38 @@ def lagrange_polynomials(x): return lagrange_basis_poly, lagrange_basis_coeffs -def lagrange_basis(i, x): + +def lagrange_function_value(lagrange_coeff_array): ''' - Calculates the value of the :math: `i^{th}` Lagrange basis (calculated - using the :math:`\\xi_{LGL}` points) at the :math: `x` coordinates. - + Funtion to calculate the value of lagrange basis functions over LGL + nodes. + Parameters ---------- - - i : int - The Lagrange basis which is to be calculated. - - x : arrayfire.Array - The coordinates at which the `i^{th}` Lagrange polynomial is to be - calculated. + lagrange_coeff_array : arrayfire.Array[N_LGL N_LGL 1 1] + Contains the coefficients of the + Lagrange basis polynomials Returns ------- + L_i : arrayfire.Array [N 1 1 1] + The value of lagrange basis functions calculated over the LGL + nodes. + + Examples + -------- + lagrange_function_value(4) gives the value of the four + Lagrange basis functions evaluated over 4 LGL points + arranged in a 2D array where Lagrange polynomials + evaluated at the same LGL point are in the same column. - l_xi_j : arrayfire.Array - Array of values of the :math: `i^{th}` Lagrange basis - calculated at the given :math: `x` coordinates. - ''' - x_tile = af.transpose(af.tile(x, 1, params.N_LGL)) - power = utils.linspace((params.N_LGL-1), 0, params.N_LGL) - power_tile = af.tile(power, 1, x.shape[0]) - x_pow = af.arith.pow(x_tile, power_tile) - l_xi_j = af.blas.matmul(params.lBasisArray[i], x_pow) + Also the value lagrange basis functions at LGL points has the property, - return l_xi_j - - -def lagrange_basis_function(lagrange_coeff_array): - ''' - Funtion to calculate the value of lagrange basis functions over LGL - nodes. + L_i(xi_k) = 0 for i != k + = 1 for i = k + + It follows then that lagrange_function_value returns an identity matrix. - 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(params.xi_LGL, 1, params.N_LGL)) power = af.flip(af.range(params.N_LGL)) @@ -203,36 +225,22 @@ def lagrange_basis_function(lagrange_coeff_array): return L_i - -def dLp_xi_LGL(lagrange_coeff_array): - ''' - 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(params.N_LGL), 1, params.N_LGL))) - * lagrange_coeff_array)[:, :-1] - - nodes_tile = af.transpose(af.tile(params.xi_LGL, 1, params.N_LGL - 1)) - power_tile = af.flip(af.tile(af.range(params.N_LGL - 1), 1, params.N_LGL)) - nodes_power_tile = af.pow(nodes_tile, power_tile) - - dLp_xi = af.blas.matmul(differentiation_coeffs, nodes_power_tile) - - 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. + + Parameters + ---------- + x : arrayfire.Array[N_LGL 1 1 1] + Contains N_LGL Gauss-Lobatto nodes. + + Returns + ------- + poly1d_product_list : list [N_LGL ** 2] + Contains the poly1d form of the product of the Lagrange + basis polynomials. ''' poly1d_list = lagrange_polynomials(x)[0] @@ -249,8 +257,8 @@ def product_lagrange_poly(x): def Integrate(integrand, N_quad, scheme): ''' - Integrates an analytical form of the integrand and integrates it using either - gaussian or lobatto quadrature. + Integrates an analytical form of the integrand and integrates it using + either gaussian or lobatto quadrature. Parameters ---------- @@ -258,8 +266,9 @@ def Integrate(integrand, N_quad, scheme): 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. + 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'. @@ -267,8 +276,8 @@ def Integrate(integrand, N_quad, scheme): Returns ------- Integral : numpy.float64 - The value o the definite integration performed using the specified - quadrature method. + The value of the definite integration performed using the + specified quadrature method. ''' if (scheme == 'gauss_quadrature'): gaussian_nodes = gauss_nodes(N_quad) @@ -289,17 +298,87 @@ def Integrate(integrand, N_quad, scheme): return Integral -def wave_equation_lagrange_basis(u, element_no): +def wave_equation_lagrange_basis_single_element(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. + Calculates the function which describes the amplitude of the wave in + a particular element. + + Using the value of the amplitude at the LGL points, A function which + describes this behaviour is obtained by expressing it as a linear + combination of the Lagrange basis polynomials. + + :math: ` f(x) = '\\sigma_i a_i L_i(\\xi)` + + Where the coefficients a_i are the value of the function at the + LGL points. + + Parameters + ---------- + u : arrayfire.Array [N_LGL N_Elements 1 1] + The amplitude of the wave at the LGL points for a + single element. + + element_no : int + The element for which the analytical form of the wave equation + is required. + + Returns + ------- + wave_equation_element : numpy.poly1d + The analytical form of the function which describes + the amplitude locally. ''' - amplitude_at_element_LGL = u[:, element_no] + 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] + wave_equation_element += af.sum(amplitude_at_element_LGL[i])\ + * lagrange_basis_polynomials[i] return wave_equation_element + +def wave_equation_lagrange(u): + ''' + Calculates the local wave equation in the Lagrange basis space + for all elements using wave_equation_lagrange_basis_single_element function. + + Parameters + ---------- + u : arrayfire.Array [N_LGL N_Elements 1 1] + Contains the amplitude of the wave at the LGL points for all elements. + + Returns + ------- + wave_equation_lagrange_basis : list [N_Elements] + Contains the local approximation of the wave + function in the form of a list + ''' + wave_equation_lagrange_basis = [] + + for i in range(0, params.N_Elements): + element_wave_equation = wave_equation_lagrange_basis_single_element(u, i) + + wave_equation_lagrange_basis.append(element_wave_equation) + + return wave_equation_lagrange_basis + +def differential_lagrange_poly1d(): + ''' + Calculates the differential of the analytical form of the Lagrange basis + polynomials. + + Returns + ------- + diff_lagrange_poly1d : list [N_LGL] + Contains the differential of the Lagrange basis + polynomials in numpy.poly1d form. + ''' + diff_lagrange_poly1d = [] + + for i in range (0, params.N_LGL): + test_diff = np.poly1d.deriv(params.lagrange_poly1d_list[i]) + diff_lagrange_poly1d.append(test_diff) + + return diff_lagrange_poly1d diff --git a/code/app/params.py b/code/app/params.py index 39031c1..c434fed 100644 --- a/code/app/params.py +++ b/code/app/params.py @@ -37,12 +37,37 @@ # Array of the lobatto weights used during integration. lobatto_weights = af.np_to_af_array(lagrange.lobatto_weights(N_LGL, xi_LGL)) + +# The number of Gaussian nodes for Gauss quadrature +N_Gauss = 10 + +# N_Gauss number of Gauss nodes. +gauss_points = lagrange.gauss_nodes(N_Gauss) + +# The Gaussian weights. +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) + + +# A list of the Lagrange polynomials in poly1d form. +lagrange_poly1d_list = lagrange.lagrange_polynomials(xi_LGL)[0] + +# List of the product of the Lagrange basis polynomials. Used in +# the calculation of A matrix. +poly1d_product_list = lagrange.product_lagrange_poly(xi_LGL) + +# list containing the poly1d forms of the differential of Lagrange +# basis polynomials. +differential_lagrange_polynomial = lagrange.differential_lagrange_poly1d() + + # An array containing the coefficients of the lagrange basis polynomials. lagrange_coeffs = af.np_to_af_array(lagrange.lagrange_polynomials(xi_LGL)[1]) # Refer corresponding functions. -dLp_xi = lagrange.dLp_xi_LGL(lagrange_coeffs) -lagrange_basis_value = lagrange.lagrange_basis_function(lagrange_coeffs) +lagrange_basis_value = lagrange.lagrange_function_value(lagrange_coeffs) # Obtaining an array consisting of the LGL points mapped onto the elements. @@ -58,7 +83,8 @@ af.sum(x_nodes[1]), N_Elements + 1) element_array = af.transpose(af.interop.np_to_af_array(np_element_array)) -element_LGL = wave_equation.mapping_xi_to_x(af.transpose(element_array), xi_LGL) +element_LGL = wave_equation.mapping_xi_to_x(af.transpose(element_array),\ + xi_LGL) # The minimum distance between 2 mapped LGL points. @@ -76,26 +102,3 @@ 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) diff --git a/code/app/wave_equation.py b/code/app/wave_equation.py index d26616d..fc02e4c 100644 --- a/code/app/wave_equation.py +++ b/code/app/wave_equation.py @@ -75,21 +75,21 @@ def mapping_xi_to_x(x_nodes, xi): def dx_dxi_numerical(x_nodes, xi): ''' - Differential :math: `\\frac{dx}{d \\xi}` calculated by central differential - method about xi using the mapping_xi_to_x function. + Differential :math: `\\frac{dx}{d \\xi}` calculated by central + differential method about xi using the mapping_xi_to_x function. Parameters ---------- - x_nodes : arrayfire.Array - Contains the nodes of elements + x_nodes : arrayfire.Array [N_Elements 1 1 1] + Contains the nodes of elements. xi : float Value of :math: `\\xi` Returns ------- - dx_dxi : arrayfire.Array + dx_dxi : arrayfire.Array [N_Elements 1 1 1] :math:`\\frac{dx}{d \\xi}`. ''' dxi = 1e-7 @@ -107,12 +107,12 @@ def dx_dxi_analytical(x_nodes, xi): :math: `\\frac{x_1 - x_0}{2}` Parameters ---------- - x_nodes : arrayfire.Array - An array containing the nodes of an element. - + x_nodes : arrayfire.Array [N_Elements 1 1 1] + Contains the nodes of elements. + Returns ------- - analytical_dx_dxi : arrayfire.Array + analytical_dx_dxi : arrayfire.Array [N_Elements 1 1 1] The analytical solution of :math: `\\frac{dx}{d\\xi}` for an element. @@ -134,10 +134,9 @@ def A_matrix(): `https://goo.gl/Cw6tnw` Returns ------- - A_matrix : arrayfire.Array + A_matrix : arrayfire.Array [N_LGL N_LGL 1 1] The value of integral of product of lagrange basis functions - obtained by LGL points, using Gauss-Lobatto quadrature method - using :math: `N_LGL` points. + obtained by LGL points, using the Integrate() function ''' A_matrix = np.zeros([params.N_LGL, params.N_LGL]) @@ -145,7 +144,7 @@ def A_matrix(): 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') + params.N_LGL + 1, 'lobatto_quadrature') dx_dxi = af.mean(dx_dxi_numerical((params.element_mesh_nodes[0 : 2]),\ params.xi_LGL)) @@ -166,58 +165,48 @@ def flux_x(u): Parameters ---------- - u : arrayfire.Array - A 1-D array which contains the value of wave function. - + u : list [N_Elements] + The analytical form of the wave equation for each element arranged in + a list of numpy.poly1d polynomials. Returns ------- - flux : arrayfire.Array - The value of the flux for given u. + flux : list [N_Elements] + The analytical value of the flux for each element arranged in a list + of numpy.poly1d polynomials. ''' flux = params.c * 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) - - 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` This will give N values of flux integral as p varies from 0 to N - 1. - This integral is carried out over an element with LGL nodes mapped onto it. + This integral is carried out using the analytical form of the integrand + This integrand is the used in the Integrate() function. Parameters ---------- u : arrayfire.Array [N_LGL N_Elements 1 1] - A 2D array containing the value of the wave function at - the mapped LGL nodes in all the elements. - + Amplitude of the wave at the mapped LGL nodes of each element. + Returns ------- flux_integral : arrayfire.Array [N_LGL N_Elements 1 1] - A 1-D array of the value of the flux integral calculated - for various lagrange basis functions. + Value of the volume integral flux. It contains the integral + of all N_LGL * N_Element integrands. ''' - wave_equation_in_lagrange_polynomials = wave_equation_lagrange_polynomials(u_n) + analytical_form_flux = flux_x(lagrange.wave_equation_lagrange(u_n)) differential_lagrange_poly = params.differential_lagrange_polynomial - flux_integral = np.zeros([params.N_LGL, params.N_Elements]) + 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') + integrand = analytical_form_flux[j] * differential_lagrange_poly[i] + flux_integral[i][j] = lagrange.Integrate(integrand, params.N_LGL,\ + 'gauss_quadrature') flux_integral = af.np_to_af_array(flux_integral) @@ -229,13 +218,17 @@ def lax_friedrichs_flux(u_n): A function which calculates the lax-friedrichs_flux :math:`f_i` using. :math:`f_i = \\frac{F(u^{i + 1}_0) + F(u^i_{N_{LGL} - 1})}{2} - \frac {\Delta x}{2\Delta t} (u^{i + 1}_0 - u^i_{N_{LGL} - 1})` - + + The algorithm used is explained in the link given below + `https://goo.gl/sNsXXK` + 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. + Amplitude of the wave at the mapped LGL nodes of each element. + ''' + u_iplus1_0 = af.shift(u_n[0, :], 0, -1) u_i_N_LGL = u_n[-1, :] @@ -259,15 +252,14 @@ def surface_term(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. - + 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:`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. @@ -295,11 +287,18 @@ def b_vector(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. + Amplitude of the wave at the mapped LGL nodes of each element. + Returns ------- - b_vector_array : arrayfire.Array + b_vector_array : arrayfire.Array [N_LGL N_Elements 1 1] + Contains the b vector(of shape [N_LGL 1 1 1]) + for each element. + + Reference + --------- + A report for the b-vector can be found here + `https://goo.gl/sNsXXK` ''' volume_integral = volume_integral_flux(u_n) diff --git a/code/main.py b/code/main.py index 68d0905..889ef44 100644 --- a/code/main.py +++ b/code/main.py @@ -11,4 +11,6 @@ if __name__ == '__main__': - wave_equation.time_evolution() + #print(lagrange.product_lagrange_poly(params.xi_LGL)) + wave_equation.time_evolution() + print(params.lagrange_basis_value) diff --git a/code/unit_test/test_waveEqn.py b/code/unit_test/test_waveEqn.py index 9b58c4f..e3962b7 100644 --- a/code/unit_test/test_waveEqn.py +++ b/code/unit_test/test_waveEqn.py @@ -1,4 +1,5 @@ #! /usr/bin/env python3 + import math import numpy as np @@ -11,7 +12,7 @@ from app import wave_equation from utils import utils -# This test uses the initial paramters N_LGL = 8, N_Elements = 10 +# This test uses the initial paramters N_LGL = 8, N_Elements = 10 and c = 1. def test_mapping_xi_to_x(): ''' @@ -147,48 +148,6 @@ def test_lagrange_coeffs(): assert af.sum(af.abs(basis_array_analytical - params.lagrange_coeffs)) < threshold - -def test_dLp_xi(): - ''' - Test function to check the dLp_xi calculated in params module with a - numerically obtained one. - - Refrence - -------- - The link to the sage worksheet where the calculations were carried out. - `https://goo.gl/Xxp3PT` - ''' - threshold = 4e-11 - - reference_d_Lp_xi = af.np_to_af_array(np.array([\ - [-14.0000000000226,-3.20991570302344,0.792476681323880,-0.372150435728984,\ - 0.243330712724289,-0.203284568901545,0.219957514771985, -0.500000000000000], - - [18.9375986071129, 3.31499272476776e-11, -2.80647579473469,1.07894468878725\ - ,-0.661157350899271,0.537039586158262, -0.573565414940005,1.29768738831567], - - [-7.56928981931106, 4.54358506455201, -6.49524878326702e-12, \ - -2.37818723350641, 1.13535801687865, -0.845022556506714, 0.869448098330221,\ - -1.94165942553537], - - [4.29790816425547,-2.11206121431525,2.87551740597844,-1.18896004153157e-11,\ - -2.38892435916370, 1.37278583181113, -1.29423205091574, 2.81018898925442], - - [-2.81018898925442, 1.29423205091574, -1.37278583181113, 2.38892435916370, \ - 1.18892673484083e-11,-2.87551740597844, 2.11206121431525,-4.29790816425547], - - [1.94165942553537, -0.869448098330221, 0.845022556506714,-1.13535801687865,\ - 2.37818723350641, 6.49524878326702e-12,-4.54358506455201,7.56928981931106],\ - - [-1.29768738831567, 0.573565414940005,-0.537039586158262,0.661157350899271,\ - -1.07894468878725,2.80647579473469,-3.31498162253752e-11,-18.9375986071129], - - [0.500000000000000,-0.219957514771985,0.203284568901545,-0.243330712724289,\ - 0.372150435728984, -0.792476681323880, 3.20991570302344, 14.0000000000226] - ])) - - 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 @@ -247,7 +206,7 @@ def test_volume_integral_flux(): given below. `https://goo.gl/5Mub8M` ''' - threshold = 4e-8 + threshold = 8e-9 params.c = 1 referenceFluxIntegral = af.transpose(af.interop.np_to_af_array(np.array([