This ipython notebook reproduces the bug discussed in github issue Quazartech/DG_Maxwell Issue #6
While evolving the wave equation using forward Euler method, The amplitude of the wave function (which should be a conserved quantity) inceasases with time and eventually blows up exponentially.
<function>?
In [1]:
#Required libraries are imported and plot parameters are set in this cell
import subprocess
import os
from matplotlib import pyplot as plt
import numpy as np
from scipy import integrate
from scipy import special as sp
from tqdm import trange
import arrayfire as af
%matplotlib inline
af.set_backend('cuda')
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['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'
In [2]:
#Utility functions used in the code, used with af.broadcast
def multiply(a, b):
'''
'''
return a * b
def power(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
The 1D wave equation is given by the equation
\begin{align} \frac{\partial u(x, t)}{\partial t} + \frac{\partial F(u)}{\partial x} = 0 \end{align}The solutions for the wave equation will be found for each of the elements separately. For each element, the $x$ space is mapped to $\xi$ space to and the calculations are done in the $\xi$ space.
$\xi \epsilon [-1, 1]$.
This $\xi$ space is divided into n Legendre-Gauss-Lobatto (LGL) points which are the roots of the equation
\begin{align} (1 - \xi^2)P_{n - 1}'(\xi) = 0 \end{align}Where $P_n(\xi)$ is the Legendre polynomial with index n.
The function to obtain the LGL points along with a test function is given below
In [3]:
#This cell contains the function which returns LGL points
#It is tested with values obtained from http://mathworld.wolfram.com/LobattoQuadrature.html
def LGL_points(N):
"""
Calculates and returns the LGL points for a given N, These are roots of
the formula.
(1 - xi ** 2) P_{n - 1}'(xi) = 0
Legendre polynomials satisfy the recurrence relation
(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.
Examples
--------
LGL_points(4)
Returns 4 LGL points in the space (-1, 1)
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 test_LGL_points():
'''
Comparing the LGL nodes obtained by LGL_points with
the reference nodes for N = 6
'''
reference_nodes = \
af.np_to_af_array(np.array([-1., -0.7650553239294647,\
-0.28523151648064504, 0.28523151648064504,\
0.7650553239294647, 1. \
] \
) \
)
calculated_nodes = (LGL_points(6))
assert(af.max(af.abs(reference_nodes - calculated_nodes)) <= 1e-14)
In [4]:
#The domain of the function.
x_boundary_nodes = af.np_to_af_array(np.array([-1., 1.]))
#Number of elements the domain is to be divided into.
N_Elements = 10 # Rename this to N_elements
#The number of LGL points into which an element is split.
N_LGL = 8
#Array containing the LGL points in xi space
xi_LGL = LGL_points(N_LGL)
An element is defined using two points. One of the points, say $x_i$ is mapped to $-1$ in the $\xi$ space and the second point, say $x_{i + 1}$ is mapped to $1$ in the $\xi$ space. Using this information we can map all the points in the $x$ space to $\xi$ space with the formula:
\begin{align} x = \frac{(1 - \xi)}{2} x_i + \frac{(1 + \xi)}{2} x_{i + 1} \end{align}The function to map a given $xi$ to $x$ space along with a test function in the cell below.
In [5]:
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`
Parameters
----------
x_nodes : arrayfire.Array [2 1 1 1]
Element nodes.
xi : arrayfire.Array [N 1 1 1]
Value of :math: `\\xi`coordinate for which the corresponding
:math: `x` coordinate is to be found.
Returns
-------
x : arrayfire.Array
:math: `x` value in the element corresponding to :math:`\\xi`.
Examples
--------
mapping_xi_to_x(x_nodes, 0)
Using x_nodes as [x1, x2], gives a value (x1 + x2) / 2
'''
N_0 = (1 - xi) / 2
N_1 = (1 + xi) / 2
N0_x0 = af.broadcast(multiply, N_0, x_nodes[0])
N1_x1 = af.broadcast(multiply, N_1, x_nodes[1])
x = N0_x0 + N1_x1
return x
def test_mapping_xi_to_x():
'''
A test function to check the mapping_xi_to_x function,
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
test_element_nodes = af.np_to_af_array(np.array([7, 11]))
test_xi = 0
analytical_x_value = 9
numerical_x_value = mapping_xi_to_x(test_element_nodes, test_xi)
assert af.abs(analytical_x_value - numerical_x_value) <= threshold
In [6]:
#Size of equally divided elements
element_size = af.sum((x_boundary_nodes[1] - x_boundary_nodes[0]) / N_Elements)
boundaries = af.constant(0, 2, N_Elements, dtype=af.Dtype.f64)
boundaries[0, :] = af.range(N_Elements)
boundaries[1, :] = af.range(N_Elements) + 1
#Array consisting of the element boundaries
boundaries = boundaries * element_size + af.sum(x_boundary_nodes[0])
#Array consisting of the LGL points mapped to each element
element_LGL = mapping_xi_to_x((boundaries), xi_LGL) # N_LGL x N_Elements
In [7]:
print(boundaries, element_LGL)
In [ ]:
Lagrange Basis polynomials calculated over the LGL points is used in the code. The Lagrange basis polynomial $L_p(\xi)$ is taken as
\begin{align} L_p(\xi) = \prod_{m = 0, m \neq p}^{m = N - 1} \frac{\xi - \xi_m}{\xi_m - \xi_p} \end{align}The function to calculate the Lagrange basis polynomials calculated using N LGL points is shown in the cell below
lagrange_basis_coeffs returns a matrix
\begin{align} \begin{bmatrix} a_0^{N - 1} & a_0^{N - 2} & \cdots & a_0^0\\ a_1^{N - 1} & a_1^{N - 2} & \cdots & a_1^0\\ \vdots & \vdots & \ddots & \vdots \\ a_{N - 1}^{N - 1} & a_{N - 1}^{N - 2} & \cdots & a_{N - 1}^0 \end{bmatrix} \end{align}Where the Lagrange basis function is given by \begin{align} L_i'(\xi_k) = a_i^{N - 1} \xi_k^{N - 1} + a_i^{N - 2} \xi_k^{N - 2} + \cdots + a_i^0 \xi_k^0 \end{align}
In [8]:
def lagrange_polynomials(x):
'''
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::
`L_i = \\prod_{m = 0, m \\notin i}^{N - 1}\\frac{(x - x_m)}{(x_i - x_m)}`
The analytical form of the Lagrange polynomials is required to compute
A_matrix and the volume integral flux (defined later in the code.)
Parameters
----------
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.
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.
Examples
--------
lagrange_polynomials(4)[0] gives the lagrange polynomials obtained using
4 LGL points in poly1d form
lagrange_polynomials(4)[0][2] is :math: `L_2(\\xi)`
lagrange_polynomials(4)[1] gives the coefficients of the above mentioned
lagrange basis polynomials in a 2D array.
lagrange_polynomials(4)[1][2] gives the coefficients of :math:`L_2(\\xi)`
in the form [a^2_3, a^2_2, a^2_1, a^2_0]
'''
X = np.array(x)
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])
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.append(lagrange_basis_j)
lagrange_basis_coeffs[j] = lagrange_basis_j.c
return lagrange_basis_poly, lagrange_basis_coeffs
def test_lagrange_polynomials():
'''
Function to test the lagrange_basis_coeffs function in global_variables
module by passing 8 LGL points and comparing the numerically obtained
basis function coefficients to analytically calculated ones. This array
of the coefficients is required later on in the program.
Reference
---------
The link to the sage worksheet where the calculations were carried out.
https://goo.gl/wMALRG
'''
threshold = 1e-12
reference_basis_coeffs = np.zeros([8, 8])
reference_basis_coeffs = \
np.array([
[-3.351562500008004, 3.351562500008006, \
3.867187500010295, -3.867187500010297, \
-1.054687500002225, 1.054687500002225, \
0.03906249999993106, -0.03906249999993102 \
], \
[ 8.140722718246403, -7.096594831382852, \
-11.34747768400062, 9.89205188146461, \
3.331608712119162, -2.904297073479968, \
-0.1248537463649464, 0.1088400233982081 \
], \
[-10.35813682892759, 6.128911440984293, \
18.68335515838398, -11.05494463699297, \
-8.670037141196786, 5.130062549476987, \
0.3448188117404021, -0.2040293534683072 \
], \
[ 11.38981374849497, -2.383879109609436, \
-24.03296250200938, 5.030080255538657, \
15.67350804691132, -3.28045297599924, \
-3.030359293396907, 0.6342518300700298 \
], \
[-11.38981374849497, -2.383879109609437, \
24.03296250200939, 5.030080255538648, \
-15.67350804691132, -3.28045297599924, \
3.030359293396907, 0.6342518300700299 \
], \
[ 10.35813682892759, 6.128911440984293, \
-18.68335515838398, -11.05494463699297, \
8.670037141196786, 5.130062549476987, \
-0.3448188117404021, -0.2040293534683072 \
], \
[-8.140722718246403, -7.096594831382852, \
11.34747768400062, 9.89205188146461, \
-3.331608712119162, -2.904297073479968, \
0.1248537463649464, 0.1088400233982081 \
],
[ 3.351562500008004, 3.351562500008005, \
-3.867187500010295, -3.867187500010298, \
1.054687500002225, 1.054687500002224, \
-0.039062499999931, -0.03906249999993102 \
] \
] \
)
reference_basis_coeffs = af.np_to_af_array(reference_basis_coeffs)
calculated_basis_coeffs = lagrange_basis_coeffs(LGL_points(8))
assert af.max(af.abs(reference_basis_coeffs - calculated_basis_coeffs)) < 1e-10
In [9]:
# A list of the Lagrange polynomials in poly1d form.
lagrange_poly1d_list = lagrange_polynomials(xi_LGL)[0]
# An array containing the coefficients of the lagrange basis polynomials.
lagrange_coeffs = af.np_to_af_array(lagrange_polynomials(xi_LGL)[1])
In [ ]:
In [ ]:
In the code, to evaluate the value of polynomials $f_k(x)$ at certain points $x_i$, the following algorithm is used.
Let the order of the function $f_k(x)$ be $n$ and the number of points $x_i$ be $m$
The function $f_k(x)$ can be represented as
\begin{align} f_k(x) = a_{n-1} x^{n-1} + a_{n-2} x^{n-2} + \cdots + a_0 x^0 \end{align}The coefficients of this polynomial can be arranged in a $1 \times n$ array
\begin{align} \begin{bmatrix} a_{n-1} & a_{n-2} & \cdots & a_1 & a_0 \end{bmatrix}_{1 \times n} \end{align}Now the points $x_i$ can be arranged in a $n \times m$ array.
\begin{align} \begin{bmatrix} x_0 & x_1 & \cdots & x_{m-2} & x_{m-1} \\ x_0 & x_1 & \cdots & x_{m-2} & x_{m-1} \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ x_0 & x_1 & \cdots & x_{m-2} & x_{m-1} \\ \end{bmatrix}_{n \times m} \end{align}An index array given below is to be obtained
\begin{align} \begin{bmatrix} n-1 \\ n-2 \\ \vdots \\ 1 \\ 0 \\ \end{bmatrix}_{n \times 1} \end{align}Using af.broadcast and the power() function declared earlier in the code, The array shown below can be obtained.
\begin{align} \begin{bmatrix} x_0^{n-1}& x_1^{n-1} & \cdots & x_{m-2}^{n-1} & x_{m-1}^{n-1} \\ x_0^{n-2}& x_1^{n-2} & \cdots & x_{m-2}^{n-2} & x_{m-1}^{n-2} \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ x_0^0 & x_1^0 & \cdots & x_{m-2}^0 & x_{m-1}^0 \\ \end{bmatrix}_{n \times m} \end{align}multiplying this $n \times m$ array to the coefficient array obtained earlier gives the value of the function at the required $x_i$
\begin{align} \begin{bmatrix} a_{n-1} & a_{n-2} & \cdots & a_1 & a_0 \end{bmatrix}_{1 \times n} \times \begin{bmatrix} x_0^{n-1}& x_1^{n-1} & \cdots & x_{m-2}^{n-1} & x_{m-1}^{n-1} \\ x_0^{n-2}& x_1^{n-2} & \cdots & x_{m-2}^{n-2} & x_{m-1}^{n-2} \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ x_0^0 & x_1^0 & \cdots & x_{m-2}^0 & x_{m-1}^0 \\ \end{bmatrix}_{n \times m} = \begin{bmatrix} f(x_0) & f(x_1)& \cdots & f(x_{m-1}) \end{bmatrix} \end{align}
In [10]:
def lagrange_function_value(lagrange_coeff_array):
'''
Funtion to calculate the value of lagrange basis functions over LGL
nodes.
Parameters
----------
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.
Also the value lagrange basis functions at LGL points has the property,
L_i(xi_k) = 0 for i != k
= 1 for i = k
It follows then that lagrange_function_value returns an identity matrix.
'''
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(lagrange_coeff_array[index], xi_pow)
return L_i
def test_lagrange_function_value():
'''
The Lagrange basis polynomials L_i(\\xi_k) will be zero when
'i' is not equal to 'k' and is 1 in all other cases. An identity matrix
is expected
'''
threshold = 1e-13
calculated_lagrange_basis_value = lagrange_function_value(8)
analytical_lagrange_basis_value = af.identity(8, 8, dtype=af.Dtype.f64)
error = af.max(af.abs(calculated_lagrange_basis_value - analytical_lagrange_basis_value))
assert error <= threshold
In [11]:
# Refer corresponding functions.
lagrange_basis_value = lagrange_function_value(lagrange_coeffs)
Integration is performed using the two quadrature methods, The Gauss-Lobatto quadrature and the Gauss-Legendre quadrature
Gauss Lobatto $\href{https://en.wikipedia.org/wiki/Gaussian_quadrature#Gauss.E2.80.93Lobatto_rules}{\mathrm{quadrature}}$. is a method to find integral of the type,
\begin{align} \int_{-1}^{1}f(\xi) d\xi \end{align}According to Lobatto quadrature, the integral is given by
\begin{align} \int_{-1}^{1}f(\xi) d\xi = \frac{2}{N (N - 1)}[f(1) + f(-1)] + \sum_{i = 1}^{N - 1}w_i f(x_i) \end{align}where, $x_i$ are the solution of $L'_{N - 1}$, the LGL points. $w_i$ are the weights calculated at $x_i$, it is given by,
\begin{align} w_i = \frac{2}{N(N - 1)[P'_{n - 1}(x_i)]^2} \end{align}
In [12]:
# The number quadrature points to be used for integration.
N_quad = 8
In [13]:
def lobatto_weight(n):
'''
Calculates and returns the weight function for an index n
and points x.
Parameters
----------
n : int
Lobatto weights for n quadrature points.
Returns
-------
Lobatto_weights : arrayfire.Array
An array of lobatto weight functions for
the given x points and index.
Reference
---------
Gauss-Lobatto weights Wikipedia link-
https://en.wikipedia.org/wiki/
Gaussian_quadrature#Gauss.E2.80.93Lobatto_rules
Examples
--------
lobatto_weight_function(4) returns the Gauss-Lobatto weights
which are to be used with the Lobatto nodes 'LGL_points(4)'
to integrate using Lobatto quadrature.
'''
xi_LGL = LGL_points(n)
P = sp.legendre(n - 1)
Lobatto_weights = (2 / (n * (n - 1)) / (P(xi_LGL))**2)
Lobatto_weights = af.np_to_af_array(Lobatto_weights)
return Lobatto_weights
def test_lobatto_weight_function():
'''
Compares the Lobatto weights obtained by the function and reference values.
'''
threshold = 1e-8
N_LGL = 6
lobatto_weights = lobatto_weight_function(N_LGL)
analytical_weights = af.Array([0.06666666666666667, 0.378474956297847,\
0.5548583770354863, 0.5548583770354863,\
0.378474956297847, 0.06666666666666667])
error = af.max(af.abs(lobatto_weights - analytical_weights))
print(error)
assert error <= threshold
The Gauss Legendre quadrature is a method to find integral of the type,
\begin{align} \int_{-1}^{1}f(\xi) d\xi \end{align}According to Gauss-Legendre quadrature, the integral is given by
\begin{align} \int_{-1}^{1}f(\xi) d\xi = \sum_{i=0}^{N - 1}w_i f(x_i) \end{align}where, $x_i$, the Gauss quadrature points are the roots of $P_n(x) = 1$.
$w_i$ are the weights calculated at $x_i$, it is given by,
\begin{align} w_i = \frac{2}{(1 - x_i^2)[P'_{n}(x_i)]^2} \end{align}
In [14]:
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)
gauss_nodes = legendre.r
gauss_nodes.sort()
return gauss_nodes
def test_gauss_nodes():
'''
The Gauss points obtained by the function above is compared to
analytical values.
Reference
---------
`https://goo.gl/9gqLpe`
'''
threshold = 1e-10
analytical_gauss_nodes = np.array([-0.906179845938664, -0.5384693101056831,\
0, 0.5384693101056831, 0.906179845938664])
calculated_gauss_nodes = gauss_nodes(5)
assert np.max(abs(analytical_gauss_nodes - calculated_gauss_nodes)) <= threshold
In [15]:
def gaussian_weights(N):
'''
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.
Returns
-------
gaussian_weight : arrayfire.Array [N_quad 1 1 1]
The gaussian weights.
'''
index = np.arange(N) # Index `i` in `w_i`, varies from 0 to N_quad - 1
gaussian_nodes = gauss_nodes(N)
gaussian_weight = 2 / ((1 - (gaussian_nodes[index]) ** 2) *\
(np.polyder(sp.legendre(N))(gaussian_nodes[index])) ** 2)
gaussian_weight = af.np_to_af_array(gaussian_weight)
return gaussian_weight
def test_gaussian_weights():
'''
Test to check the gaussian weights calculated.
'''
threshold = 2e-8
analytical_gauss_weights = af.Array([0.23692688505618908, 0.47862867049936647,\
0.5688888888888889, 0.47862867049936647, \
0.23692688505618908
]
)
calculated_gauss_weights = lagrange.gaussian_weights(5)
In [16]:
# N_Gauss number of Gauss nodes.
gauss_points = af.np_to_af_array(gauss_nodes(N_quad))
# The Gaussian weights.
gauss_weights = gaussian_weights(N_quad)
# The lobatto weights to be used for integration.
lobatto_weights = lobatto_weight(N_quad)
# The lobatto nodes to be used for integration.
lobatto_quadrature_nodes = LGL_points(N_quad)
Integration in the code is performed using an Integrate() function.
This function takes in the coefficients of the integrand and depending on the scheme, either lobatto_quadrature or gauss_quadrature, The integration is carried out.
To calculate the value of the integrand at the quadrature points is calculated using the algorithm discussed earlier.
In [17]:
def Integrate(integrand_coeffs):
'''
Performs integration according to the given quadrature method
by taking in the coefficients of the polynomial and the number of
quadrature points.
The number of quadrature points and the quadrature scheme are set
earlier in the code.
Parameters
----------
integrand_coeffs : arrayfire.Array [M N 1 1]
The coefficients of M number of polynomials of order N
arranged in a 2D array.
Returns
-------
Integral : arrayfire.Array [M 1 1 1]
The value of the definite integration performed using the
specified quadrature method for M polynomials.
'''
if (scheme == 'gauss_quadrature'):
integrand = integrand_coeffs
gaussian_nodes = gauss_points
Gauss_weights = gauss_weights
nodes_tile = af.transpose(af.tile(gaussian_nodes, 1, integrand.shape[1]))
power = af.flip(af.range(integrand.shape[1]))
power_tile = af.tile(power, 1, N_quad)
nodes_power = nodes_tile ** power_tile
weights_tile = af.transpose(af.tile(Gauss_weights, 1, integrand.shape[1]))
nodes_weight = nodes_power * weights_tile
value_at_gauss_nodes = af.matmul(integrand, nodes_weight)
Integral = af.sum(value_at_gauss_nodes, 1)
if (scheme == 'lobatto_quadrature'):
integrand = integrand_coeffs
lobatto_nodes = lobatto_quadrature_nodes
Lobatto_weights = lobatto_weights
nodes_tile = af.transpose(af.tile(lobatto_nodes, 1, integrand.shape[1]))
power = af.flip(af.range(integrand.shape[1]))
power_tile = af.tile(power, 1, N_quad)
nodes_power = nodes_tile ** power_tile
weights_tile = af.transpose(af.tile(Lobatto_weights, 1, integrand.shape[1]))
nodes_weight = nodes_power * weights_tile
value_at_lobatto_nodes = af.matmul(integrand, nodes_weight)
Integral = af.sum(value_at_lobatto_nodes, 1)
return Integral
def test_Integrate():
'''
Testing the Integrate() function by passing coefficients of a polynomial and comparing
it to the analytical result.
'''
threshold = 1e-14
test_coeffs = af.np_to_af_array(np.array([7., 6, 4, 2, 1, 3, 9, 2]))
# The coefficients of a test polynomial
# `7x^7 + 6x^6 + 4x^5 + 2x^4 + x^3 + 3x^2 + 9x + 2`
# Using Integrate() function. Integrates using N_quad quadrature points
calculated_integral = lagrange.Integrate(af.transpose(test_coeffs))
analytical_integral = 8.514285714285714
assert (calculated_integral - analytical_integral) <= threshold
In [18]:
# Set the scheme for integration by quadrature method.
scheme = 'lobatto_quadrature'
$\frac{dx}{d\xi}$ needs to be used throughout the code. A function to calculate the same for an element is in the cell below.
In [19]:
def dx_dxi(x_nodes, xi):
'''
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
xi : float
Value of :math: `\\xi`
Returns
-------
dx_dxi : arrayfire.Array
: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)
dx_dxi = (x2 - x1) / (2 * dxi)
return dx_dxi
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-9.
'''
threshold = 1e-9
nodes = np.array([7, 10], dtype = np.float64)
test_nodes = af.interop.np_to_af_array(nodes)
analytical_dx_dxi = 1.5 #calculating (10 - 7) / 2
check_dx_dxi = abs(af.mean(dx_dxi
(test_nodes, xi_LGL)) - analytical_dx_dxi) <= threshold
assert check_dx_dxi
In [20]:
print(dx_dxi(boundaries, 0))
The strong form (differential form) of the wave equation is converted to the weak form (integral form). To do this, A test function $v(x)$ is multiplied with the wave equation and the resultant expresion is integrated w.r.t. $x$.
\begin{align} v(x)\frac{\partial u(x, t)}{\partial t} &+ v(x)\frac{\partial F(u)}{\partial x} = 0\\ \int v(x)\frac{\partial u(x, t)}{\partial t} dx &+ \int v(x)\frac{\partial F(u)}{\partial x} dx = 0 \\ \int v(x)\frac{\partial u(x, t)}{\partial t} dx &+ v(x) F(u) - \int F(u)\frac{\partial v(x)}{\partial x} dx = 0 \end{align}Taking $u^{n}(x)$ to be the wave function at the $n^{th}$ time step. The aim is to find $u^{n + 1}(x)$ (i.e. wave function at $(n + 1)^{th}$ time step).
Taking an increment of $\Delta t$ in time between the consecutive timesteps $n$ and $n + 1$
\begin{align} \frac{\int v(x)u^{n + 1}(x) dx - \int v(x)u^{n}(x) dx} {\Delta t} + v(x) F(u^n) - \int F(u^n)\frac{\partial v(x)}{\partial x} dx = 0 \end{align}\begin{align} \int v(x)u^{n + 1}(x) dx &- \int v(x)u^{n}(x) dx + \Delta t \{v(x) F(u^n) - \int F(u^n) \frac{\partial v(x)}{\partial x} dx\} = 0 \\ \int v(x)u^{n + 1}(x) dx &= \int v(x)u^{n}(x) dx - \Delta t \{v(x) F(u^n) - \int F(u^n) \frac{\partial v(x)}{\partial x} dx\} \end{align}Mapping $x$ into $\xi$ space (whose domain is [-1, 1]), for a single element,
\begin{align} \int^1_{-1} v(\xi)u^{n + 1}(\xi) \frac{dx}{d\xi} d\xi &= \int^1_{-1} v(\xi)u^{n}(\xi) \frac{dx}{d\xi}d\xi - \Delta t \{v(\xi) F(u^n) - \int^1_{-1} F(u^n) \frac{\partial v(\xi)}{\partial \xi} d\xi\} \end{align}Choosing $v(\xi)$ as the Lagrange basis $L_p$ polynomials created using the $N$ Lagrange-Gauss-Lobatto(LGL) points
\begin{align} L_p(\xi) = \prod_{m = 0, m \neq p}^{m = N - 1} \frac{\xi - \xi_m}{\xi_m - \xi_p} \end{align}For each element, the wave function can be written in $\xi$ space as a linear combination of the Lagrange basis polynomials created using the LGL points in that space, i.e.,
\begin{align} u^n(\xi) = \sum_{i = 0}^{N - 1} u^n_i L_i(\xi) \end{align}In the above equation, $u^n_i$ are the coefficients for $L_i$ and $N$ is the total number LGL points inside the $\xi$ space.
Substituting equation $u^n(\xi) = \sum_{i = 0}^{N - 1} u^n_i L_i(\xi)$ and $v(\xi) = L_p(\xi)$ into the modified weak form of the wave equation,
\begin{align} \sum_{i = 0}^{N - 1} u^{n + 1}_i \int L_p(\xi)L_i(\xi) \frac{dx}{d\xi} d\xi &= \sum_{i = 0}^{N - 1} u^n_i \int L_p(\xi)L_i(\xi) \frac{dx}{d\xi}d\xi - \Delta t \{L_p(\xi) F(u^n) - \int F(u^n) \frac{\partial L_p(\xi)}{\partial \xi} d\xi\} \end{align}On varying $p$ from $0$ to $N - 1$, $N$ such linear equations for $N$ variable $u^{n + 1}_i$, $i \epsilon {0, 1, \cdots, N - 1}$. Writing all the system of linear equations in matrix representation.
\begin{align} \begin{bmatrix} A_{00} & A_{01} & \cdots & A_{0{N-1}} \\ A_{10} & A_{11} & \cdots & A_{1{N-1}} \\ \vdots & \vdots & \ddots & \vdots \\ A_{(N - 1)0} & A_{(N - 1)1} & \cdots & A_{(N - 1){(N-1)}} \end{bmatrix} \begin{bmatrix} u^{n + 1}_0 \\ u^{n + 1}_1 \\ \vdots \\ u^{n + 1}_{N - 1} \end{bmatrix} = \begin{bmatrix} b_0 \\ b_1 \\ \vdots \\ b_{N - 1} \\ \end{bmatrix} \end{align}or
\begin{align} Au^{n + 1} = b \\ u^{n + 1} = A^{-1} b \end{align}As seen above, $u^{n + 1}$ can be obtained if $u^n$ is known. This system of linear equations is for a single element, arranging the wave functions at LGL points mapped for all elements in the form
\begin{align} U^n = \begin{bmatrix} u^{n,0}_0 & u^{n,1}_0 & \cdots & u^{n,{N_{Elements} - 1}}_0\\ u^{n,0}_1 & u^{n,1}_1 & \cdots & u^{n,{N_{Elements} - 1}}_1\\ u^{n,0}_{N - 1} & u^{n,1}_{N - 1} & \cdots & u^{n,{N_{Elements} - 1}}_{N - 1}\\ \end{bmatrix} \end{align}An element of the array $u^{n,i}_j$ is the amplitude at the $i^{th}$ element
So the equation which describes the time evolution for all elements can be written as
\begin{align} AU^{n + 1} = b \\ U^{n + 1} = A^{-1} b \end{align}Where $u^{n, i}$ is the wave function at the mapped LGL points at timestep $n$ for an element $i$
To find the solution of the wave equation in the next time step, the $A$ matrix has to be calculated, which has elements $A_{pi}$ described as. \begin{align} A_{p i} = \int_{-1}^{1} L_p(\xi)L_i(\xi) \frac{dx}{d\xi} d\xi \end{align}
This integral can be evaluated using the Gauss-Lobatto Quadrature.
\begin{align} \sum_{k = 0}^{N - 1} w_k L_p(\xi_k)L_i(\xi_k) \frac{dx}{d\xi}|_{\xi_k} \end{align}
In [21]:
def product_lagrange_poly(x):
'''
Used to obtain the coefficients of the product of Lagrange polynomials.
A matrix involves integrals of the product of the Lagrange polynomials.
The Integrate() function requires the coefficients of the integrand to
compute the integral.
This function takes the poly1d form of the Lagrange basis polynomials,
multiplies them and stores the coefficients in a 2D array.
Parameters
----------
x : arrayfire.Array[N_LGL 1 1 1]
Contains N_LGL Gauss-Lobatto nodes.
Returns
-------
lagrange_product_coeffs : arrayfire.Array [N_LGL**2 N_LGL*2-1 1 1]
Contains the coefficients of the product of the
Lagrange polynomials.
Examples
--------
product_lagrange_poly(xi_LGL)[0] gives the coefficients of the product
`L_0(\\xi) * L_0(\\xi)`.
product_lagrange_poly(xi_LGL)[1] gives the coefficients of the product
`L_0(\\xi) * L_1(\\xi)`.
'''
poly1d_list = lagrange_polynomials(xi_LGL)[0]
lagrange_product_coeffs = np.zeros([N_LGL ** 2, N_LGL * 2 - 1])
for i in range (N_LGL):
for j in range (N_LGL):
lagrange_product_coeffs[N_LGL * i + j] = ((poly1d_list[i] * poly1d_list[j]).c)
lagrange_product_coeffs = af.np_to_af_array(lagrange_product_coeffs)
return lagrange_product_coeffs
In [22]:
lagrange_product_coefficients = product_lagrange_poly(xi_LGL)
A function to calculate the A matrix is given below
In [23]:
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}`
The integrals are computed using the Integrate() function.
Since elements are taken to be of equal size, :math: `\\frac {dx}{dxi}
is same everywhere
Returns
-------
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 the Integrate() function
'''
int_Li_Lp = Integrate(lagrange_product_coefficients)
dx_dxi_element = af.mean(dx_dxi((boundaries[0 : 2]),\
xi_LGL))
A_matrix_flat = dx_dxi_element * int_Li_Lp
A_matrix = af.moddims(A_matrix_flat, N_LGL, N_LGL)
return A_matrix
def test_A_matrix():
'''
The A matrix is expected to be a diagonal matrix with weights multiplied with
:math: `\frac{dx}{d\\xi}` as the diagonal elements.
'''
threshold = 1e-10
weights = lobatto_weights
identity_matrix = af.identity(N_LGL, N_LGL,\
dtype = af.Dtype.f64)
dx_dxi_element = af.mean(dx_dxi((boundaries[0 : 2, 0]),\
xi_LGL))
analytical_A_matrix = af.broadcast(multiply,\
identity_matrix*dx_dxi_element,\
weights)
calculated_A_matrix = A_matrix()
error = af.max(af.abs(analytical_A_matrix - calculated_A_matrix))
print(error)
assert error <= threshold
In [24]:
scheme = 'gauss_quadrature'
print(A_matrix())
To obtain the b vector, Its elements $b_p$ described below need to be calculated
\begin{align} b_p &= \Sigma_{i = 0}^{N - 1} U^{n}_i \int^1_{-1} L_p(\xi)L_i(\xi) \frac{dx}{d\xi}d\xi - \Delta t \{L_p(\xi) F(U^n)|^1_{-1} - \int^1_{-1} F(U^n) \frac{\partial L_p(\xi)}{\partial \xi} d\xi\} \end{align}The final b vector for all elements would be,
\begin{align} b = \begin{bmatrix} b_0^0 & b_0^1 & \cdots & b_0^{N_{Elements} - 1} \\ b_1^0 & b_1^1 & \cdots & b_1^{N_{Elements} - 1} \\ \vdots & \vdots & \ddots & \vdots \\ b_{N - 1}^0 & b_{N - 1}^1 & \cdots & b_{N - 1}^{N_{Elements} - 1} \\ \end{bmatrix}_{N \times 1} \end{align}$b_p$ consists of three different terms, $\Sigma_{i = 0}^{N - 1} u^n_i \int^1_{-1} L_p(\xi)L_i(\xi)\frac{dx}{d\xi}d\xi$, $\Delta t L_p(\xi) F(u^n)|^1_{-1}$ and $\Delta t \int^1_{-1} F(u^n) \frac{\partial L_p(\xi)}{\partial \xi}d\xi$.
The first term of $b_p$, $\Sigma_{i = 0}^{N - 1} U^n_i \int^1_{-1} L_p(\xi)L_i(\xi)\frac{dx}{d\xi}d\xi$ with p varying from $0$ to $N - 1$ is $A U^n$
By varying p and i from $0$ to $N - 1$ in $\int L_p(\xi)L_i(\xi)\frac{dx} {d\xi}d\xi$, The A matrix is obtained.
$U^n$ (The displacement of the mapped LGL points across the domain) at timestep n is the obtained from the variable 'U' declared earlier.
In [25]:
def first_term(t_n, A):
'''
'''
A_u_n = af.matmul(A, U[:, :, t_n])
return A_u_n
In [26]:
#The speed of the wave
c = 1
#The minimum distance between 2 mapped LGL points
delta_x = af.min((element_LGL -
af.shift(element_LGL, 1, 0))[1:, :])
delta_t = delta_x / (20 * c)
#The total time for which the simulation is to be carried out.
total_time = 10
no_of_steps = int(total_time / delta_t)
end_time = delta_t * no_of_steps
#Time array consisting of timesteps
time = linspace(0, end_time, no_of_steps)
#c_lax, used in calculation of Lax-Friedrichs flux is taken to be the wave speed.
c_lax = c
# Storing the amplitude of the wave at the mapped LGL points as time evolves
# in a 3D array where timesteps vary along dimension 3.
U = af.constant(0, N_LGL, N_Elements, time.shape[0], dtype = af.Dtype.f64)
#Initial wave amplitude
U_init = np.e ** (-(element_LGL)**2 / 0.4**2)
U[:, :, 0] = U_init
The third term of the b vector is $\int^1_{-1} F(u^n) \frac{\partial L_p(\xi)} {\partial \xi} d\xi$. This integral is to be evaluated using Gauss-Lobatto rules.
\begin{align} \int^1_{-1} F(u^n) \frac{\partial L_p(\xi)} {\partial \xi} d\xi = \Sigma_{k = 0}^{N - 1} w_k \frac{\partial L_p(\xi)}{\partial \xi}|_{\xi_k} F(u^n) \end{align}The integrand needs to be calculated and the coefficients need to be passed to the Integrate() function to obtain the volume integral.
To get the analytical form of $u$ can be calculated for each element using the amplitude of the wave at the LGL points.
To do this, first the function is represented as a linear combination of the lagrange basis polynomials.
$f(x) = a_i L_i(\xi)$
Where the coefficients a_i are the value of the function at the LGL points.
The coefficients of the integrand are then passed into Integrate() function to obtain the volume integral flux
The analytical form of $\frac{\partial L_p(\xi)} {\partial \xi} d\xi$ is calculated below.
In [27]:
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, N_LGL):
test_diff = np.poly1d.deriv(lagrange_poly1d_list[i])
diff_lagrange_poly1d.append(test_diff)
return diff_lagrange_poly1d
Functions to obtain the analytical form of the wave equation are declared below.
In [28]:
def wave_equation_lagrange_basis_single_element(u, element_no):
'''
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]
lagrange_basis_polynomials = lagrange_poly1d_list
wave_equation_element = np.poly1d([0])
for i in range(0, N_LGL):
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, 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 flux_x(u):
'''
A function which returns the value of flux for a given wave function u.
:math:`f(u) = c u^k`
Parameters
----------
u : list [N_Elements]
The analytical form of the wave equation for each element arranged in
a list of numpy.poly1d polynomials.
Returns
-------
flux : list [N_Elements]
The analytical value of the flux for each element arranged in a list
of numpy.poly1d polynomials.
'''
flux = c * u
return flux
Now to calculate the volume integral flux using the analytical form of the integrand, obtaining it's coefficients and using the Integrate() function
In [29]:
def volume_integral_flux(u_n, N_quad = 8):
'''
Calculates the volume integral of flux in the wave equation.
:math:`\\int_{-1}^1 f(u) \\frac{d L_p}{d\\xi} d\\xi`
This will give N values of flux integral as p varies from 0 to N - 1.
This integral is carried out using the analytical form of the integrand
This integrand is the used in the Integrate() function.
Parameters
----------
u : arrayfire.Array [N_LGL N_Elements 1 1]
Amplitude of the wave at the mapped LGL nodes of each element.
Returns
-------
flux_integral : arrayfire.Array [N_LGL N_Elements 1 1]
Value of the volume integral flux. It contains the integral
of all N_LGL * N_Element integrands.
'''
analytical_form_flux = flux_x(wave_equation_lagrange(u_n))
differential_lagrange_poly = differential_lagrange_polynomial
integrand = np.zeros(([N_LGL * N_Elements, 2 * N_LGL - 2]))
for i in range(N_LGL):
for j in range(N_Elements):
integrand[i + N_LGL * j] = (analytical_form_flux[j] *\
differential_lagrange_poly[i]).c
integrand = af.np_to_af_array(integrand)
flux_integral = Integrate(integrand)
flux_integral = af.moddims(flux_integral, N_LGL, N_Elements)
return flux_integral
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://goo.gl/5Mub8M`
'''
threshold = 8e-9
params.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.volume_integral_flux(params.u[:, :, 0])
assert (af.max(af.abs(numerical_flux - referenceFluxIntegral)) < threshold)
In [30]:
# list containing the poly1d forms of the differential of Lagrange
# basis polynomials.
differential_lagrange_polynomial = differential_lagrange_poly1d()
In [31]:
print(volume_integral_flux(U[:, :, 0]))
By varying p from $0$ to $N - 1$ in the second term (i.e. $F(U^n) L_p(\xi)|^{1}_{-1}$), Surface term is obtained.
The surface term is used because of the discontinuities in element boundaries. Since the value of the flux would be different for adjacent elements, There are two values for $F(u)$ at the element boundaries. To resolve this ambiguity of flux, Lax-Friedrichs flux is introduced.
The Lax-Friedrichs flux $f_M$ for an element boundary between element $M$ and $M + 1$ at a timestep $n$ is given by,
\begin{equation} f_M = \frac{F(u^{n, M+1}_{0}) + F(u^{n, M}_{N - 1})}{2} - \frac{\Delta x}{2\Delta t} (u^{n, M + 1}_{0} - u^{n, M}_{N - 1}) \end{equation}In the above equation, $u^{n, M + 1}_{0}$ is the amplitude at the first (leftmost) LGL point of the element with index $M + 1$, Similarly $u^{n, M}_{N - 1}$ is the $(N - 1)^{th}$ (rightmost) LGL point of the element with index M, $\Delta x$ is the size of the element M. This is the flux for element boundaries.
In [32]:
def lax_friedrichs_flux(u_n):
'''
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]
Amplitude of the wave at the mapped LGL nodes of each element.
Returns
-------
boundary_flux : arrayfire.Array [1 N_Elements 1 1]
Contains the value of the flux at the boundary elements.
Periodic boundary conditions are used.
'''
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 \
- c_lax * (u_iplus1_0 - u_i_N_LGL) / 2
return boundary_flux
The Surface term is obtained by varying $p$ in $F(u^n) L_p(\xi)|^{1}_{-1}$ from $0$ to $N - 1$ for an element with an index $M$
\begin{equation} F(u^n) L_p(\xi)|^{1}_{-1} = L_p (1) f_M - L_p (-1) f_{M - 1} \end{equation}The equation above has $L_p (1)$ and $L_p (-1)$. These are elements of 1D arrays of shape $N \times 1$. Since the index p varies from $0$ to $N -1$ multiplying them with the flux at the element boundaries (a number), A vector of shape $N \times 1$ is obtained for a single element.
This when calculated for $N_{Elements}$ elements, An array of shape $N \times N_{Elements}$ is obtained.
In [33]:
def surface_term(u_n):
'''
Calculates the surface term,
:math:`L_p(1) f_i - L_p(-1) f_{i - 1}`
using the lax_friedrichs_flux function and lagrange_basis_value
from params module.
Parameters
----------
u_n : arrayfire.Array [N_LGL N_Elements 1 1]
Amplitude of the wave at the mapped LGL nodes of each element.
Returns
-------
surface_term : arrayfire.Array [N_LGL N_Elements 1 1]
The surface term represented in the form of an array,
:math:`L_p (1) f_i - L_p (-1) f_{i - 1}`, where p varies
from zero to :math:`N_{LGL}` and i from zero to
:math:`N_{Elements}`. p varies along the rows and i along
columns.
Reference
---------
Link to PDF describing the algorithm to obtain the surface term.
`https://goo.gl/Nhhgzx`
'''
L_p_minus1 = lagrange_basis_value[:, 0]
L_p_1 = lagrange_basis_value[:, -1]
f_i = lax_friedrichs_flux(u_n)
f_iminus1 = af.shift(f_i, 0, 1)
surface_term = af.blas.matmul(L_p_1, f_i) - af.blas.matmul(L_p_minus1,\
f_iminus1)
return surface_term
def test_surface_term():
'''
A test function to test the surface_term function using
analytical Lax-Friedrichs flux.
'''
threshold = 1e-13
analytical_f_i = (U[-1, :, 0])
analytical_f_i_minus1 = (af.shift(U[-1, :, 0], 0, 1))
L_p_1 = af.constant(0, N_LGL, dtype = af.Dtype.f64)
L_p_1[N_LGL - 1] = 1
L_p_minus1 = af.constant(0, 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 = (surface_term(U[:, :, 0]))
assert af.max(af.abs(analytical_surface_term - numerical_surface_term)) \
< threshold
return analytical_surface_term
In [34]:
print(surface_term(U[:, :, 0]))
In [35]:
def b_vector(u_n):
'''
Calculates the b vector for N_Elements number of elements.
Parameters
----------
u_n : arrayfire.Array [N_LGL N_Elements 1 1]
Amplitude of the wave at the mapped LGL nodes of each element.
Returns
-------
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)
Surface_term = surface_term(u_n)
b_vector_array = delta_t * (volume_integral - Surface_term)
return b_vector_array
In [36]:
print(b_vector(U[:, :, 0]) * 1e9)
In [37]:
def time_evolution():
'''
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. The movie created
by stitching these images can be found in the code folder.
'''
A_inverse = af.inverse(A_matrix())
for t_n in trange(0, time.shape[0] - 1):
U[:, :, t_n + 1] = af.matmul(A_inverse, b_vector(U[:, :, t_n]))
print('u calculated!')
subprocess.run(['mkdir', 'results/1D_Wave_images'])
for t_n in trange(0, time.shape[0] - 1):
if t_n % 100 == 0:
fig = plt.figure()
axes = plt.gca()
x = element_LGL
y = U[:, :, t_n]
plt.plot(x, y)
axes.set_ylim(-1, 2)
plt.xlabel('x')
plt.ylabel('Amplitude')
plt.title('Time = %f' % (t_n * delta_t))
fig.savefig('results/1D_Wave_images/%04d' %(t_n / 100) + '.png')
plt.close('all')
# Creating a movie with the images created.
os.system("cd results/1D_Wave_images && ffmpeg -f image2 -i %04d.png -vcodec mpeg4\
-mbd rd -trellis 2 -cmp 2 -g 300 -pass 1 -r 25 -b 18000000 movie.mp4\
&& rm -f ffmpeg2pass-0.log && mv movie.mp4 /home/f2013305/DG_Maxwell/code && rm *")
return
In [38]:
time_evolution()
In [ ]: