Merge pull request #74 from amanabt/advection_1d_multiple_u

Update 1D advection equation to work for multiple `u` and `flux`

.gitignore

@@ -197,4 +197,6 @@ Thumbs.db
 !examples/read_and_plot_mesh/mesh/*
 !dg_maxwell/tests/mesh/*
 !dg_maxwell/tests/msh_parser/mesh/*
+!examples/maxwells_equations/*
+!examples/lagrange_interpolation_general/*
 !docs/quazar_logo.png

.travis.yml

@@ -54,4 +54,6 @@ install:
 
 script:
     - python -c "import dg_maxwell"
-    - pytest --verbose -r P --color=yes
+    - pytest --verbose -r P --color=yes --cov dg_maxwell
+after_success:
+  - coveralls

README.md

@@ -3,7 +3,7 @@
 [![Build Status](https://travis-ci.org/QuazarTech/DG_Maxwell.svg?branch=master)](https://travis-ci.org/QuazarTech/DG_Maxwell)
 [![Codacy Badge](https://api.codacy.com/project/badge/Grade/e8733cdbf1454af0ac35ae5b2d017d9f)](https://www.codacy.com/app/aman2official/DG_Maxwell_2?utm_source=github.com&utm_medium=referral&utm_content=QuazarTech/DG_Maxwell&utm_campaign=Badge_Grade)
 [![Documentation Status](http://readthedocs.org/projects/dg-maxwell/badge/?version=latest)](http://dg-maxwell.readthedocs.io/en/latest/?badge=latest)
-
+[![Coverage Status](https://coveralls.io/repos/github/QuazarTech/DG_Maxwell/badge.svg?branch=master)](https://coveralls.io/github/QuazarTech/DG_Maxwell?branch=master)
 ```
 _____   _____   __  __                          _ _ 
 |  __ \ / ____| |  \/  |                        | | |
@@ -67,6 +67,14 @@ solving the wave equation in a rectangular domain in the
 - [arrayfire](http://arrayfire.org)
 - [gmshtranslate](https://github.com/amanabt/gmshtranslator)
 - [texlive](https://www.tug.org/texlive/)
+- h5py
+- tqdm
+- scipy
+- pytest
+- pytest-cov
+- coveralls
+- python-coveralls
+
 
 ## Documentation
 The documenation for the project can be generated using `sphinx`
@@ -93,6 +101,6 @@ $ pytest
   The movie will be created automatically and stored in `results` folder.
 
 ## Maintainers
-- Aman Abhishek Tiwari - ![aman@quazartech.com](aman@quazartech.com)
-- Balavarun P - ![f2013462@pilani.bits-pilani.ac.in](f2013462@pilani.bits-pilani.ac.in)
-- Manichandra Morumpudi - ![mani@quazartech.com](mani@quazartech.com)
+- Aman Abhishek Tiwari - [aman@quazartech.com](aman@quazartech.com)
+- Balavarun P - [f2013462@pilani.bits-pilani.ac.in](f2013462@pilani.bits-pilani.ac.in)
+- Manichandra Morumpudi - [mani@quazartech.com](mani@quazartech.com)

dg_maxwell/__init__.py

@@ -3,5 +3,6 @@
 
 import arrayfire as af
 af.set_backend('cpu')
+af.set_device(0)
 
 from dg_maxwell import params

dg_maxwell/isoparam.py

@@ -2,7 +2,10 @@
 # -*- coding: utf-8 -*-
 
 import arrayfire as af
+
 af.set_backend('cpu')
+af.set_device(0)
+
 from dg_maxwell import utils
 
 def isoparam_1D(x_nodes, xi):

dg_maxwell/lagrange.py

@@ -4,7 +4,9 @@
 import numpy as np
 from scipy import special as sp
 import arrayfire as af
+
 af.set_backend('cpu')
+af.set_device(0)
 
 from dg_maxwell import utils
 from dg_maxwell import params
@@ -376,3 +378,28 @@ def L1_norm(u):
     L1_norm = af.sum(integrate(interpolated_coeffs))
 
     return L1_norm
+
+
+def lagrange_interpolation(fn_i):
+    '''
+    Finds the general interpolation of a function.
+    
+    Parameters
+    ----------
+    fn_i : af.Array [N N_LGL 1 1]
+           Value of :math:`N` functions at the LGL points.
+    
+    Returns
+    -------
+    lagrange_interpolation : af.Array [N N_LGL 1 1]
+                             :math:`N` interpolated polynomials for
+                             :math:`N` functions.
+    '''
+
+    fn_i = af.transpose(af.reorder(fn_i, d0 = 2, d1 = 1, d2 = 0))
+    lagrange_interpolation = af.broadcast(utils.multiply,
+                                          params.lagrange_coeffs, fn_i)
+    lagrange_interpolation = af.reorder(af.sum(lagrange_interpolation, dim = 0),
+                                        d0 = 2, d1 = 1, d2 = 0)
+
+    return lagrange_interpolation

dg_maxwell/msh_parser.py

@@ -5,7 +5,9 @@
 import matplotlib.lines as lines
 import gmshtranslator.gmshtranslator as gmsh
 import arrayfire as af
+
 af.set_backend('cpu')
+af.set_device(0)
 
 from dg_maxwell import msh_parser
 from dg_maxwell import isoparam

dg_maxwell/params.py

@@ -3,7 +3,9 @@
 
 import numpy as np
 import arrayfire as af
+
 af.set_backend('cpu')
+af.set_device(0)
 
 from dg_maxwell import lagrange
 from dg_maxwell import utils
@@ -31,7 +33,7 @@
 N_quad     = 8
 
 # Wave speed.
-c          = 1.
+c          = 1.0
 
 # The total time for which the wave is to be evolved by the simulation. 
 total_time = 2.01
@@ -83,13 +85,14 @@
                   af.sum(x_nodes[1] - element_size), N_Elements)
 
 np_element_array   = np.concatenate((af.transpose(elements),
-                               af.transpose(elements + element_size)))
+                                     af.transpose(elements + element_size)))
 element_mesh_nodes = utils.linspace(af.sum(x_nodes[0]),
-                                    af.sum(x_nodes[1]), N_Elements + 1)
+                                    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.
 delta_x = af.min((element_LGL - af.shift(element_LGL, 1, 0))[1:, :])
@@ -104,7 +107,7 @@
 
 # Array of timesteps seperated by delta_t.
 time    = utils.linspace(0, int(total_time / delta_t) * delta_t,
-                                                    int(total_time / delta_t))
+                         int(total_time / delta_t))
 
 
 # The wave to be advected is either a sin or a Gaussian wave.

dg_maxwell/tests/convergence_tests.py

@@ -1,7 +1,10 @@
 from dg_maxwell import params
 from dg_maxwell import wave_equation
 
+import arrayfire as af
 
+af.set_backend('cpu')
+af.set_device(0)
 
 def L1_norm(u):
     '''

dg_maxwell/tests/isoparam/test_isoparam.py

@@ -6,7 +6,10 @@
 sys.path.insert(0, os.path.abspath('./'))
 
 import arrayfire as af
+
 af.set_backend('cpu')
+af.set_device(0)
+
 import numpy as np
 
 from dg_maxwell import isoparam

dg_maxwell/tests/msh_parser/test_msh_parser.py

@@ -6,7 +6,9 @@
 sys.path.insert(0, os.path.abspath('./'))
 
 import arrayfire as af
+
 af.set_backend('cpu')
+af.set_device(0)
 
 import numpy as np
 

dg_maxwell/tests/test_waveEqn.py

@@ -7,7 +7,9 @@
 
 import numpy as np
 import arrayfire as af
+
 af.set_backend('cpu')
+af.set_device(0)
 
 from dg_maxwell import params
 from dg_maxwell import lagrange
@@ -370,9 +372,12 @@ def test_volume_integral_flux():
 
          ])))
 
-    numerical_flux = wave_equation.volume_integral_flux(params.u[:, :, 0])
+    numerical_flux = wave_equation.volume_integral_flux(
+        params.u[:, :, 0])
     assert (af.mean(af.abs(numerical_flux - referenceFluxIntegral)) < threshold)
 
+
+
 def test_lax_friedrichs_flux():
     '''
     A test function to test the lax_friedrichs_flux function in wave_equation
@@ -387,6 +392,7 @@ def test_lax_friedrichs_flux():
     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
@@ -415,6 +421,7 @@ def test_surface_term():
     return analytical_surface_term
 
 
+
 def test_b_vector():
     '''
     A test function to check the b vector obtained analytically and compare it
@@ -427,14 +434,17 @@ def test_b_vector():
 
     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])
+    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)) * params.delta_t
     b_vector_array       = wave_equation.b_vector(params.u[:, :, 0])
 
     assert (b_vector_analytical - b_vector_array) < threshold
 
+
+
 def test_integrate():
     '''
     Testing the integrate() function by passing coefficients

dg_maxwell/tests/utils/test_utils.py

@@ -0,0 +1,73 @@
+#! /usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+import numpy as np
+import arrayfire as af
+
+from dg_maxwell import params
+from dg_maxwell import utils
+from dg_maxwell import lagrange
+
+af.set_backend('cpu')
+af.set_device(0)
+
+def test_matmul_3D():
+    '''
+    '''
+    M = 3
+    N = 2
+    P = 4
+    Q = 2
+
+    a = af.range(M * N * Q, dtype = af.Dtype.u32)
+    b = af.range(N * P * Q, dtype = af.Dtype.u32)
+
+    a = af.moddims(a, d0 = M, d1 = N, d2 = Q)
+    b = af.moddims(b, d0 = N, d1 = P, d2 = Q)
+
+    a_init = a
+    b_init = b
+
+    ref_a_0 = np.matmul(np.array(a_init[:, :, 0]),
+                        np.array(b_init[:, :, 0]))
+
+    ref_a_1 = np.matmul(np.array(a_init[:, :, 1]),
+                        np.array(b_init[:, :, 1]))
+
+    test_matmul = np.array(utils.matmul_3D(a, b))
+
+    diff_mat_0 = np.abs(test_matmul[:, :, 0] - ref_a_0)
+    diff_mat_1 = np.abs(test_matmul[:, :, 1] - ref_a_1)
+
+    assert np.all(diff_mat_0 == 0) and np.all(diff_mat_1 == 0)
+
+
+
+def test_poly1d_prod():
+    '''
+    Checks the product of the polynomials of different degrees using the
+    poly1d_product function and compares it to the analytically calculated
+    product coefficients.
+    '''
+
+    N      = 3
+
+    N_a    = 3
+    poly_a = af.range(N * N_a, dtype = af.Dtype.u32)
+    poly_a = af.moddims(poly_a, d0 = N, d1 = N_a)
+
+    N_b    = 2
+    poly_b = af.range(N * N_b, dtype = af.Dtype.u32)
+    poly_b = af.moddims(poly_b, d0 = N, d1 = N_b)
+
+    ref_poly = af.np_to_af_array(np.array([[0., 0., 9., 18.],
+                                           [1., 8., 23., 28.],
+                                           [4., 20., 41., 40.]]))
+
+    test_poly1d_prod = utils.poly1d_product(poly_a, poly_b)
+    test_poly1d_prod_commutative = utils.poly1d_product(poly_b, poly_a)
+
+    diff     = af.abs(test_poly1d_prod - ref_poly)
+    diff_commutative = af.abs(test_poly1d_prod_commutative - ref_poly)
+
+    assert af.all_true(diff == 0.) and af.all_true(diff_commutative == 0.)