Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

1D wave equation solver.[WIP] #16

Merged
merged 20 commits into from
Sep 18, 2017
Merged

1D wave equation solver.[WIP] #16

Show file tree
Hide file tree
Changes from 17 commits
Commits
Show all changes
20 commits
Select commit Hold shift + click to select a range
ea3a7b5
Push changes into lax_friedrichs_flux
Balavarun5 Aug 1, 2017
b8526a2
Added functions to obtain the lax-friedrichs-flux and surface term.
Balavarun5 Aug 2, 2017
3f92f4a
Time evolution of wave equation obtained.
Balavarun5 Aug 4, 2017
1014ef4
Simulating for higher wave speeds.
Balavarun5 Aug 5, 2017
c5d0a2f
The images would be created in 1D_Wave_images folder, stitch them to …
Balavarun5 Aug 5, 2017
bc14e39
Changed amplitude from being represented as an N_Elements x N_LGL x t…
Balavarun5 Aug 7, 2017
0011aaf
Resolving wave speed abnormality.
Balavarun5 Aug 8, 2017
64f414c
-Wave is moving at the expected speed by taking c_lax to be the wave …
Balavarun5 Aug 11, 2017
7507b3d
Unit tests working.
Balavarun5 Aug 11, 2017
756c8f6
Unit tests working.
Balavarun5 Aug 11, 2017
bdffc6c
Edits in the documentation.
Balavarun5 Aug 26, 2017
88b04d0
Minor edits
Balavarun5 Aug 26, 2017
6c5c3d0
Documentations edited. Tabs converted to spaces.
Balavarun5 Aug 31, 2017
4b5b4a7
Cleanup of the code, Unit tests working. Iteration speed unchanged (d…
Balavarun5 Sep 8, 2017
8d38457
Minor edits in documentation of functions.
Balavarun5 Sep 8, 2017
457f42d
Change in note for developers in readme.md
Balavarun5 Sep 8, 2017
28729c3
Merge branch 'master' into lax_friedrichs_flux
Balavarun5 Sep 8, 2017
f82df88
Temporary commit. Not all changes have been implemented. Code is to b…
Balavarun5 Sep 15, 2017
e6b5c71
Code reffactored to use Integrate() function. Unit tests passing. Doc…
Balavarun5 Sep 17, 2017
f79ceb9
Edited the documentation. Examples need to be added.
Balavarun5 Sep 18, 2017
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
23 changes: 22 additions & 1 deletion .gitignore
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
2 changes: 1 addition & 1 deletion README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,4 +27,4 @@ python3 main.py
* **Mani Chandra** - [GitHub Profile](https://github.com/mchandra)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We updated the README a few days ago. You should have the updated README here. I am attaching the
link where you can find it.
hdf5_ImplementationForStoring_u

## Note for developers:
* Use tab spaces for indentation.
* Use spaces for indentation.
Empty file removed code/app/__init__.py
View file
Open in desktop
Empty file.
315 changes: 105 additions & 210 deletions code/app/global_variables.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,201 +1,72 @@
import numpy as np
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

According to the Quazartech style guidelines, the imports should be grouped in the following order
Imports should be grouped in the following order:

  • standard library imports
  • related third party imports
  • local application/library specific imports

I would prefer if you give spacing between different groups.
I am adding this review only here but you should check for this issue in every python file.

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


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]]


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])

x_nodes = af.interop.np_to_af_array(np.array([-1., 1.]))
N_LGL = 16
xi_LGL = None
lBasisArray = None
lobatto_weights = None
N_Elements = None
element_nodes = None
u = None
time = None
c = None
dLp_xi = None

def populateGlobalVariables(Number_of_LGL_pts = 8, Number_of_elements = 10):
# 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):
'''
Time evolution of the wave function requires the use of many variables which
are time independent throughout the program. These variables have been
declared as global variables in the program. The global variables being used
in the program are listed below.

x_nodes : af.Array [2 1 1 1]
Stores the domain of the wave. Default value is set to
[-1, 1]. This means that the domain is from
:math:`x \\epsilon [-1., 1.]`. This domain will be split
into `Number_of_elements` elements.

N_Elements : int
The number of elements into which the domain is to be
divided.

N_LGL : int
The number of `LGL` nodes in an element.

xi_LGL : af.Array [N_LGL 1 1 1]
`N_LGL` LGL points.

element_LGL : af.Array [N_LGL N_Elements 1 1]
The domain is divided into `N_elements` number of elements
and each element has N_LGL points corresponding to the
each LGL point. This variable stores the :math:`x`
coordinates corresponding to each `xi_LGL` for each
element. element_LGL[n_LGL, n_element] returns the
:math:`n\\_LGL^{th}` :math:`x` coordinate for the
:math:`n_element^{th}` element.

lBasisArray : af.Array [N_LGL N_LGL 1 1]
Contains the coefficients of the Lagrange basis functions
created using LGL nodes stored in xi_LGL.
lBasisArray[i, j] is the coefficient corresponding to the
:math:`(N_LGL - j - 1)^{th}` power of :math:`x` in the
:mat:`i^th` Lagrange polynomial.

dLp_xi : af.Array [N_LGL N_LGL 1 1]
Stores the value of the derivative
:math:`\\frac{dL_p}{d\\xi}` at all the LGL points.

lobatto_weights : af.Array [N_LGL 1 1 1]
Lobatto weights to integrate using Gauss-Lobatto
quadrature.

c : float
Number denoting the wave speed. Used in flux calculation.

time : af.Array
Array containing the time for which the :math:`u` is to be
calculated.

u : af.array [N_LGL N_Elements time.shape[0] 1]
Stores :math:`u` calculated at the :math:`x` cordinates
stored in the variable `element_LGL` at time :math:`t`
stored in stored in variable `time`.


Calculates and returns the LGL points for a given N, These are roots of
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You can write,

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

the formula.
:math: `(1 - xi ** 2) P_{n - 1}'(xi) = 0`
Legendre polynomials satisfy the recurrence relation
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What is the point adding this recurrence relation?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The variable legendre_N_minus_1 uses this recurrence relation.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It would have made sense to give the recurrence relation if you had explained the algorithm. But since you have not explained the algorithm, the reader won't know why you have written 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_nodes
N_Elements = Number_of_elements
element_size = af.sum((x_nodes[1] - x_nodes[0]) / N_Elements)
elements_xi_LGL = af.constant(0, N_Elements, N_LGL)
elements = utils.linspace(af.sum(x_nodes[0]), \
af.sum(x_nodes[1] - element_size), N_Elements)
N : Int
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Int is not a keyword
int is the keyword.
Type Int() with upper case I and tell me the output.

Number of LGL nodes required

np_element_array = np.concatenate((af.transpose(elements),
af.transpose(elements + element_size)))

element_array = af.transpose(af.interop.np_to_af_array(np_element_array))
element_nodes = wave_equation.mappingXiToX(
af.transpose(element_array), xi_LGL)

global u
global time
u_init = np.e ** (-(element_nodes) ** 2 / 0.4 ** 2)
time = utils.linspace(0, 2, 200)
u = af.constant(0, N_LGL, N_Elements, time.shape[0],\
dtype = af.Dtype.f64)

u[:, :, 0] = u_init

global dLp_xi
dLp_xi = dLp_xi_LGL()

global c
c = 1.0

return
Returns
-------
lgl : arrayfire.Array [N 1 1 1]
An arrayfire array consisting of the Lagrange-Gauss-Lobatto Nodes.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Avoid redundant information. You have already mentioned that it is an arrayfire.Array.
I don't think that you should state that again.


Reference
---------
http://mathworld.wolfram.com/LobattoQuadrature.html
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Use markdown to state a link.

'''
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`
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

  • What is a weight function? It's not returning a weight function, it's returning weights.
  • You should also mention what weights is it returning, Gaussian weights, Lobatto weight etc.
  • And you do not need to write Returns, if the reader wants to know what it returns, then
    he/she could look at the Returns section of the doc-string.

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
----------
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
An array of lobatto weight functions for
the given :math: `x` points and index.

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-
Expand All @@ -205,51 +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]

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

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

why are you still using global variables?
You mentioned in Issue #17 that this code is not using global variables. What's going on?

# The domain of the function which is of interest.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

  • What function?
  • No need to write which is of interest.
  • You may also write what each of the element of the array stores.

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.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Be minimalistic. Just say Wave speed

c = 1

# The total time for which the simulation is to be carried out.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Total time for which the wave is 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 = 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
Loading