Skip to content

Commit

Permalink
A matrix function written.
Browse files Browse the repository at this point in the history 
- Documented.
- Algo not properly documented
- Unit test written
- Example written
Issues
- Unit test for test_wave_equation_2d.py::test_A_matrix and
  test_waveEqn.py::test_volume_integral_flux FAILING when using opencl as backend.
  Works fine on cpu and cuda backend.
  • Loading branch information
@amanabt
amanabt committed Oct 22, 2017
1 parent 5794c88 commit 90e4f16
Show file tree
Hide file tree
Showing 8 changed files with 1,832 additions and 10 deletions.
2 changes: 1 addition & 1 deletion dg_maxwell/params.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@
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
N_LGL = 4

# Number of elements the domain is to be divided into.
N_Elements = 10
Expand Down
16 changes: 16 additions & 0 deletions dg_maxwell/tests/wave_equation_2d/files/A_matrix_ref.csv
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,16 @@
2.040837009684099909e-02,7.607021305186759813e-03,-7.605211827216239694e-03,3.401274391609880078e-03,7.607021305186719914e-03,2.835443147246799814e-03,-2.834768681946190153e-03,1.267791922578240032e-03,-7.605211827216239694e-03,-2.834768681946230051e-03,2.834094377080259922e-03,-1.267490353611489939e-03,3.401274391609880078e-03,1.267791922578240032e-03,-1.267490353611489939e-03,5.668589618928889613e-04
7.607021305186759813e-03,1.020408165439520010e-01,1.700378677733640001e-02,-7.605211827216239694e-03,2.835443147246799814e-03,3.803472113476859956e-02,6.337996011943730114e-03,-2.834768681946190153e-03,-2.834768681946230051e-03,-3.802567383658889827e-02,-6.336488396321220129e-03,2.834094377080249947e-03,1.267791922578240032e-03,1.700619964078519883e-02,2.833863961292090033e-03,-1.267490353611489939e-03
-7.605211827216239694e-03,1.700378677733640001e-02,1.020408165439520010e-01,7.607021305186730323e-03,-2.834768681946190153e-03,6.337996011943730114e-03,3.803472113476859956e-02,2.835443147246820197e-03,2.834094377080259922e-03,-6.336488396321220129e-03,-3.802567383658889827e-02,-2.834768681946200127e-03,-1.267490353611489939e-03,2.833863961292090033e-03,1.700619964078519883e-02,1.267791922578240032e-03
3.401274391609880078e-03,-7.605211827216239694e-03,7.607021305186730323e-03,2.040837009684099909e-02,1.267791922578240032e-03,-2.834768681946190153e-03,2.835443147246820197e-03,7.607021305186719914e-03,-1.267490353611489939e-03,2.834094377080249947e-03,-2.834768681946200127e-03,-7.605211827216239694e-03,5.668589618928889613e-04,-1.267490353611489939e-03,1.267791922578240032e-03,3.401274391609880078e-03
7.607021305186719914e-03,2.835443147246799814e-03,-2.834768681946190153e-03,1.267791922578240032e-03,1.020408165439520010e-01,3.803472113476869670e-02,-3.802567383658889827e-02,1.700619964078510168e-02,1.700378677733640001e-02,6.337996011943739655e-03,-6.336488396321209721e-03,2.833863961292090033e-03,-7.605211827216239694e-03,-2.834768681946230051e-03,2.834094377080259922e-03,-1.267490353611489939e-03
2.835443147246799814e-03,3.803472113476859956e-02,6.337996011943730114e-03,-2.834768681946190153e-03,3.803472113476869670e-02,5.101989130708830533e-01,8.501807243133240044e-02,-3.802567383658889827e-02,6.337996011943739655e-03,8.501807243133249758e-02,1.416716589307020005e-02,-6.336488396321209721e-03,-2.834768681946230051e-03,-3.802567383658889827e-02,-6.336488396321220129e-03,2.834094377080249947e-03
-2.834768681946190153e-03,6.337996011943730114e-03,3.803472113476859956e-02,2.835443147246820197e-03,-3.802567383658889827e-02,8.501807243133240044e-02,5.101989130708830533e-01,3.803472113476859956e-02,-6.336488396321209721e-03,1.416716589307020005e-02,8.501807243133249758e-02,6.337996011943739655e-03,2.834094377080259922e-03,-6.336488396321220129e-03,-3.802567383658889827e-02,-2.834768681946200127e-03
1.267791922578240032e-03,-2.834768681946190153e-03,2.835443147246820197e-03,7.607021305186719914e-03,1.700619964078510168e-02,-3.802567383658889827e-02,3.803472113476859956e-02,1.020408165439520010e-01,2.833863961292090033e-03,-6.336488396321209721e-03,6.337996011943739655e-03,1.700378677733640001e-02,-1.267490353611489939e-03,2.834094377080249947e-03,-2.834768681946200127e-03,-7.605211827216239694e-03
-7.605211827216239694e-03,-2.834768681946230051e-03,2.834094377080259922e-03,-1.267490353611489939e-03,1.700378677733640001e-02,6.337996011943739655e-03,-6.336488396321209721e-03,2.833863961292090033e-03,1.020408165439520010e-01,3.803472113476869670e-02,-3.802567383658889827e-02,1.700619964078510168e-02,7.607021305186719914e-03,2.835443147246799814e-03,-2.834768681946190153e-03,1.267791922578240032e-03
-2.834768681946230051e-03,-3.802567383658889827e-02,-6.336488396321220129e-03,2.834094377080249947e-03,6.337996011943739655e-03,8.501807243133249758e-02,1.416716589307020005e-02,-6.336488396321209721e-03,3.803472113476869670e-02,5.101989130708830533e-01,8.501807243133240044e-02,-3.802567383658889827e-02,2.835443147246799814e-03,3.803472113476859956e-02,6.337996011943730114e-03,-2.834768681946190153e-03
2.834094377080259922e-03,-6.336488396321220129e-03,-3.802567383658889827e-02,-2.834768681946200127e-03,-6.336488396321209721e-03,1.416716589307020005e-02,8.501807243133249758e-02,6.337996011943739655e-03,-3.802567383658889827e-02,8.501807243133240044e-02,5.101989130708830533e-01,3.803472113476859956e-02,-2.834768681946190153e-03,6.337996011943730114e-03,3.803472113476859956e-02,2.835443147246820197e-03
-1.267490353611489939e-03,2.834094377080249947e-03,-2.834768681946200127e-03,-7.605211827216239694e-03,2.833863961292090033e-03,-6.336488396321209721e-03,6.337996011943739655e-03,1.700378677733640001e-02,1.700619964078510168e-02,-3.802567383658889827e-02,3.803472113476859956e-02,1.020408165439520010e-01,1.267791922578240032e-03,-2.834768681946190153e-03,2.835443147246820197e-03,7.607021305186719914e-03
3.401274391609880078e-03,1.267791922578240032e-03,-1.267490353611489939e-03,5.668589618928889613e-04,-7.605211827216239694e-03,-2.834768681946230051e-03,2.834094377080259922e-03,-1.267490353611489939e-03,7.607021305186719914e-03,2.835443147246799814e-03,-2.834768681946190153e-03,1.267791922578240032e-03,2.040837009684099909e-02,7.607021305186759813e-03,-7.605211827216239694e-03,3.401274391609880078e-03
1.267791922578240032e-03,1.700619964078519883e-02,2.833863961292090033e-03,-1.267490353611489939e-03,-2.834768681946230051e-03,-3.802567383658889827e-02,-6.336488396321220129e-03,2.834094377080249947e-03,2.835443147246799814e-03,3.803472113476859956e-02,6.337996011943730114e-03,-2.834768681946190153e-03,7.607021305186759813e-03,1.020408165439520010e-01,1.700378677733640001e-02,-7.605211827216239694e-03
-1.267490353611489939e-03,2.833863961292090033e-03,1.700619964078519883e-02,1.267791922578240032e-03,2.834094377080259922e-03,-6.336488396321220129e-03,-3.802567383658889827e-02,-2.834768681946200127e-03,-2.834768681946190153e-03,6.337996011943730114e-03,3.803472113476859956e-02,2.835443147246820197e-03,-7.605211827216239694e-03,1.700378677733640001e-02,1.020408165439520010e-01,7.607021305186730323e-03
5.668589618928889613e-04,-1.267490353611489939e-03,1.267791922578240032e-03,3.401274391609880078e-03,-1.267490353611489939e-03,2.834094377080249947e-03,-2.834768681946200127e-03,-7.605211827216239694e-03,1.267791922578240032e-03,-2.834768681946190153e-03,2.835443147246820197e-03,7.607021305186719914e-03,3.401274391609880078e-03,-7.605211827216239694e-03,7.607021305186730323e-03,2.040837009684099909e-02
25 changes: 24 additions & 1 deletion dg_maxwell/tests/wave_equation_2d/test_wave_equation_2d.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,13 +8,15 @@

import numpy as np
import arrayfire as af
af.set_backend('cpu')

from dg_maxwell import params
from dg_maxwell import utils
from dg_maxwell import msh_parser
from dg_maxwell import lagrange
from dg_maxwell import wave_equation_2d

af.set_backend(params.backend)

def test_dx_dxi():
'''
This test checks the derivative :math:`\\frac{\\partial x}{\\partial \\xi}`
Expand Down Expand Up @@ -171,3 +173,24 @@ def test_jacobian():
check = af.abs(jacobian - jacobian_reference) < threshold

assert af.all_true(check)


def test_A_matrix():
'''
Compares the tensor product calculated using the ``A_matrix`` function
with an analytic value of the tensor product for :math:`N_{LGL} = 4`.
The analytic value of the tensor product is calculated in this
`document`_
.. _document: https://goo.gl/QNWxXp
'''
threshold = 1e-4

A_matrix_ref = af.np_to_af_array(utils.csv_to_numpy(
'dg_maxwell/tests/wave_equation_2d/files/A_matrix_ref.csv'))

params.N_LGL = 4

A_matrix_test = wave_equation_2d.A_matrix()

assert af.all_true(af.abs(A_matrix_test - A_matrix_ref) < threshold)
81 changes: 77 additions & 4 deletions dg_maxwell/wave_equation_2d.py
View file
Original file line number Diff line number Diff line change
@@ -1,8 +1,14 @@
#! /usr/bin/env python3
# -*- coding: utf-8 -*-

import arrayfire as af

from dg_maxwell import params
from dg_maxwell import isoparam
from dg_maxwell import lagrange
from dg_maxwell import utils

af.set_backend(params.backend)

def dx_dxi(x_nodes, xi, eta):
'''
Expand Down Expand Up @@ -208,10 +214,77 @@ def jacobian(x_nodes, y_nodes, xi, eta):

def A_matrix():
'''
Calculates the tensor product for the given ``params.N_LGL``.
A tensor product element is given by:
.. math:: [A^{pq}_{ij}] = \\iint L_p(\\xi) L_q(\\eta) \\
L_i(\\xi) L_j(\\eta) d\\xi d\\eta
This function finds :math:`L_p(\\xi) L_i(\\xi)` and
:math:`L_q(\\eta) L_j(\\eta)` and passes it to the ``integrate_2d``
function.
Returns
-------
A : af.Array [N_LGL^2 N_LGL^2 1 1]
The tensor product.
'''
N_LGL_quad = 9
xi_LGL_quad = lagrange.gauss_nodes(N_LGL_quad)
eta_LGL_quad = lagrange.gauss_nodes(N_LGL_quad)

N_LGL = params.N_LGL
xi_LGL = lagrange.LGL_points(N_LGL)
eta_LGL = lagrange.LGL_points(N_LGL)

_, Lp = lagrange.lagrange_polynomials(xi_LGL)
Lp = af.np_to_af_array(Lp)
Li = Lp.copy()

_, Lq = lagrange.lagrange_polynomials(eta_LGL)
Lq = af.np_to_af_array(Lq)
Lj = Lq.copy()

Li = af.reorder(Li, d0 = 0, d1 = 2, d2 = 1)
Li = af.transpose(af.tile(Li, d0 = 1, d1 = N_LGL))
Li = af.moddims(Li, d0 = N_LGL * N_LGL, d1 = 1, d2 = N_LGL)
Li = af.transpose(af.tile(Li, d0 = 1, d1 = N_LGL * N_LGL))
Li = af.transpose(af.moddims(af.transpose(Li),
d0 = N_LGL * N_LGL * N_LGL * N_LGL,
d1 = 1, d2 = N_LGL))
Li = af.reorder(Li, d0 = 2, d1 = 1, d2 = 0)

Lp = af.reorder(Lp, d0 = 0, d1 = 2, d2 = 1)
Lp = af.transpose(af.tile(Lp, d0 = 1, d1 = N_LGL))
Lp = af.moddims(Lp, d0 = N_LGL * N_LGL, d1 = 1, d2 = N_LGL)
Lp = af.tile(Lp, d0 = 1, d1 = N_LGL * N_LGL)
Lp = af.moddims(af.transpose(Lp),
d0 = N_LGL * N_LGL * N_LGL * N_LGL,
d1 = 1, d2 = N_LGL)
Lp = af.reorder(af.transpose(Lp), d0 = 2, d1 = 1, d2 = 0)

Lp_Li = af.transpose(af.convolve1(Li, Lp, conv_mode = af.CONV_MODE.EXPAND))

Lj = af.reorder(Lj, d0 = 0, d1 = 2, d2 = 1)
Lj = af.tile(Lj, d0 = 1, d1 = N_LGL)
Lj = af.moddims(Lj, d0 = N_LGL * N_LGL, d1 = 1, d2 = N_LGL)
Lj = af.transpose(af.tile(Lj, d0 = 1, d1 = N_LGL * N_LGL))
Lj = af.moddims(af.transpose(Lj),
d0 = N_LGL * N_LGL * N_LGL * N_LGL,
d1 = 1, d2 = N_LGL)
Lj = af.reorder(Lj, d0 = 2, d1 = 0, d2 = 1)

Lq = af.reorder(Lq, d0 = 0, d1 = 2, d2 = 1)
Lq = af.tile(Lq, d0 = 1, d1 = N_LGL)
Lq = af.moddims(Lq, d0 = N_LGL * N_LGL, d1 = 1, d2 = N_LGL)
Lq = af.tile(Lq, d0 = 1, d1 = N_LGL * N_LGL)
Lq = af.moddims(af.transpose(Lq),
d0 = N_LGL * N_LGL * N_LGL * N_LGL,
d1 = 1, d2 = N_LGL)
Lq = af.reorder(af.transpose(Lq), d0 = 2, d1 = 1, d2 = 0)

Lq_Lj = af.transpose(af.convolve1(Lj, Lq, conv_mode = af.CONV_MODE.EXPAND))

A = af.moddims(utils.integrate_2d(Lp_Li, Lq_Lj,
order = 9,
scheme = 'gauss'),
d0 = N_LGL * N_LGL, d1 = N_LGL * N_LGL)

return
return A
617 changes: 617 additions & 0 deletions examples/A_matrix_2d_from_integrate.ipynb
View file

Large diffs are not rendered by default.

11 changes: 10 additions & 1 deletion examples/A_matrix.ipynb → examples/A_matrix_2d_old_derivation.ipynb
View file
Original file line number Diff line number Diff line change
@@ -1,8 +1,17 @@
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Calculates A matrix given by\n",
"\n",
"$$[A^{pq}_{ij}] = \\iint L_p(\\xi) L_q(\\eta) L_i(\\xi) L_j(\\eta) J \\Big(\\frac{x, y}{\\xi, \\eta} \\Big) d\\xi d\\eta$$"
]
},
{
"cell_type": "code",
"execution_count": 15,
"execution_count": 2,
"metadata": {},
"outputs": [
{
Expand Down
Loading

0 comments on commit 90e4f16

Please sign in to comment.