gauss_legendre

module for Gauss-Legendre-Lobatto and Gauss-Legendre-Radau spectral element grids

A struct for passing around elemental matrices on Gauss-Legendre points in 1D

A struct for Gauss-Legendre arrays needed for global operations in 1D, contains the struct of elemental matrices for Lobatto and Radau points, as well as some assembled 1D global matrices.

A function that takes the first derivative in each element of coord.grid, leaving the result (element-wise) in coord.scratch_2d.

Function to carry out a 1D (global) mass matrix solve.

Gauss-Legendre derivative at arbitrary x values, for boundary condition on Radau points. D0 – the vector xj – the x location where the derivative is evaluated ngrid – number of points in x x – the grid from -1, 1 Note that D0 is not scaled to the physical grid with a scaling factor.

Assign abitrary weak (nonlinear) inner product matrix Q on a 1D line with Jacobian equal to 1. matrix Q acts on two vectors x1 and x2 such that the quadratic form y = x1 * Q * x2 is also a vector. See documentation of corresponding function for linear inner product matrices.

Assign abitrary weak inner product matrix Q on a 1D line with Jacobian equal to 1 matrix Q acts on a single vector x such that y = Q * x is also a vector.

We use a projection onto Gauss-Legendre polynomials to carry out the calculation in two steps (see, e.g, S. A. Teukolsky, Short note on the mass matrix for Gauss–Lobatto grid points, J. Comput. Phys. 283 (2015) 408–413. https://doi.org/10.1016/j.jcp.2014.12.012). First, we write the desired matrix elements in terms of Legendre polynomials

\[ l_i(x) = \sum_j \frac{P_j(x)P_j(x_i)w_i}{\gamma_j}\]

with $w_i$ the weights from an integration on the Gauss-Legendre-Lobatto (or Radau) points $x_i$, i.e.,

\[ \int^1_{-1} f(x) d x = \sum_{i} f(x_i)w_i,\]

and $\gamma_j = \sum_k w_k P_j(x_k)P_j(x_k)$ the numerical inner-product. Then, a matrix element can be expressed in integrals over Legendre polynomials rather than Lagrange polynomials, i.e.,

\[ M_{ij} = \int^1_{-1} l_i(x)l_j(x) d x = \sum_{mn} \frac{w_m P_m(x_i) w_n P_n(x_j)}{\gamma_m\gamma_n} \int^{1}_{-1} P_m(x)P_n(x) d x. \]

Defining

\[ A_{mn} = \int^{1}_{-1} P_m(x)P_n(x) d x, \]

we can thus write

\[ M_{ij} = \sum_{mn} \frac{w_m P_m(x_i) w_n P_n(x_j)}{\gamma_m\gamma_n} A_{mn}. \]

We can use a quadrature which yields exact results (to machine precision) to evaluate $A_{mn}$ using fast library functions for the Legendre polynomials, and then carry out the sum $\sum_{mn}$ to obtain exact results (to machine-precision). Here we use a Gauss-Legendre integration quadrature with exact results up to polynomials with order $k_{max} = 4N +1$, with $N=$ngrid and the highest order polynomial product that we integrate is $P_{N-1}(x)P_{N-1}(x)x^2$, which has order $k=2N < k_{max}$.

Result of the inner product of Legendre polynomials of order k.

Formula for Gauss-Legendre-Lobatto differentiation matrix taken from p196 of Chpt The Spectral Elemtent Method' ofComputational Seismology'. Heiner Igel First Edition. Published in 2017 by Oxford University Press. Or https://doc.nektar.info/tutorials/latest/fundamentals/differentiation/fundamentals-differentiationch2.html

D -- differentiation matrix 
x -- Gauss-Legendre-Lobatto points in [-1,1]
ngrid -- number of points per element (incl. boundary points)

Note that D has does not include a scaling factor

Formula for Gauss-Legendre-Radau differentiation matrix taken from https://doc.nektar.info/tutorials/latest/fundamentals/differentiation/fundamentals-differentiationch2.html

D -- differentiation matrix 
x -- Gauss-Legendre-Radau points in [-1,1)
ngrid -- number of points per element (incl. boundary points)

Note that D has does not include a scaling factor

If called for coord.name = vperp elemental matrix KJ on the $i^{th}$ element is

\[ (KJ)_{jk} = -\int^{v_\perp^U}_{v_\perp^L} \frac{\partial\varphi_j(v_\perp)}{\partial v_\perp}\frac{\partial\varphi_k(v_\perp)}{\partial v_\perp} v_\perp^2 d v_\perp = -\int^1_{-1} (c_i + x s_i)^2l_j^\prime(x)l_k^\prime(x) d x /s_i\]

with $c_i$ and $s_i$ the appropriate shift and scale factors, respectively. Otherwise, if called for any other coordinate elemental matrix KJ is the same as LL (see get_LL_local()!).

If called for coord.name = vperp elemental matrix KK on the $i^{th}$ element is

\[ K_{jk} = -\int^{v_\perp^U}_{v_\perp^L} \left(v_\perp \frac{\partial\varphi_j(v_\perp)}{\partial v_\perp} + \varphi_j(v_\perp) \right) \frac{\partial\varphi_k(v_\perp)}{\partial v_\perp} d v_\perp = -\int^1_{-1} ((c_i + x s_i)l_j^\prime(x) + l_j(x))l_k^\prime(x) d x /s_i\]

with $c_i$ and $s_i$ the appropriate shift and scale factors, respectively. Otherwise, if called for any other coordinate elemental matrix KK is the same as LL (see get_LL_local!()). IfexplicitBCterms = true`, boundary terms arising from integration by parts are included at the extreme boundary points.

If called for coord.name = vperp elemental matrix LL on the $i^{th}$ element is

\[ L_{jk} = -\int^{v_\perp^U}_{v_\perp^L} \frac{\partial\varphi_j(v_\perp)}{\partial v_\perp}\frac{\partial\varphi_k(v_\perp)}{\partial v_\perp} v_\perp d v_\perp = -\int^1_{-1} (c_i + x s_i)l_j^\prime(x)l_k^\prime(x) d x /s_i\]

with $c_i$ and $s_i$ the appropriate shift and scale factors, respectively. Otherwise, if called for any other coordinate elemental matrix LL is

\[ L_{jk} = -\int^{v_\|^U}_{v_\|^L} \frac{\partial\varphi_j(v_\|)}{\partial v_\|}\frac{\partial\varphi_k(v_\|)}{\partial v_\|} d v_\| = -\int^1_{-1} l_j^\prime(x)l_k^\prime(x) d x /s_i.\]

If explicit_BC_terms = true, boundary terms arising from integration by parts are included at the extreme boundary points.

If called for coord.name = vperp elemental matrix MM on the $i^{th}$ element is

\[ M_{jk} = \int^{v_\perp^U}_{v_\perp^L} \varphi_j(v_\perp)\varphi_k(v_\perp) v_\perp d v_\perp = \int^1_{-1} (c_i + x s_i)l_j(x)l_k(x) s_i d x \]

with $c_i$ and $s_i$ the appropriate shift and scale factors, respectively. Otherwise, if called for any other coordinate elemental matrix MM is

\[ M_{jk} = \int^{v_\|^U}_{v_\|^L} \varphi_j(v_\|)\varphi_k(v_\|) d v_\| = \int^1_{-1} l_j(x)l_k(x) s_i d x.\]

If called for coord.name = vperp elemental matrix MN on the $i^{th}$ element is

\[ (MN)_{jk} = \int^{v_\perp^U}_{v_\perp^L} \varphi_j(v_\perp)\varphi_k(v_\perp) d v_\perp = \int^1_{-1} l_j(x)l_k(x) s_i d x \]

with $c_i$ and $s_i$ the appropriate shift and scale factors, respectively. Otherwise, if called for any other coordinate elemental matrix MN is the same as MM (see get_MM_local!()).

If called for coord.name = vperp elemental matrix MR on the $i^{th}$ element is

\[ (MR)_{jk} = \int^{v_\perp^U}_{v_\perp^L} \varphi_j(v_\perp)\varphi_k(v_\perp) v_\perp^2 d v_\perp = \int^1_{-1} (c_i + s_i x)^2 l_j(x)l_k(x) s_i d x \]

with $c_i$ and $s_i$ the appropriate shift and scale factors, respectively. Otherwise, if called for any other coordinate elemental matrix MR is the same as MM (see get_MM_local!()).

If called for coord.name = vperp elemental matrix PP on the $i^{th}$ element is

\[ P_{jk} = \int^{v_\perp^U}_{v_\perp^L} \varphi_j(v_\perp)\frac{\partial\varphi_k(v_\perp)}{\partial v_\perp} v_\perp d v_\perp = \int^1_{-1} (c_i + x s_i)l_j(x)l_k^\prime(x) d x \]

with $c_i$ and $s_i$ the appropriate shift and scale factors, respectively. Otherwise, if called for any other coordinate elemental matrix PP is

\[ P_{jk} = \int^{v_\|^U}_{v_\|^L} \varphi_j(v_\|)\frac{\partial\varphi_k(v_\|)}{\partial v_\|} d v_\| = \int^1_{-1} l_j(x)l_k^\prime(x) d x.\]

If called for coord.name = vperp elemental matrix PP on the $i^{th}$ element is

\[ (PU)_{jk} = \int^{v_\perp^U}_{v_\perp^L} \varphi_j(v_\perp)\frac{\partial\varphi_k(v_\perp)}{\partial v_\perp} v_\perp^2 d v_\perp = \int^1_{-1} (c_i + x s_i)^2l_j(x)l_k^\prime(x) d x \]

with $c_i$ and $s_i$ the appropriate shift and scale factors, respectively. Otherwise, if called for any other coordinate elemental matrix PU is the same as PP see get_PP_local!().

Construction function for nonlinear diffusion matrices, only used in the assembly of the collision operator

Construction function to provide the appropriate elemental matrix Q to the global matrix assembly functions.

If called for coord.name = vperp elemental matrix YY0 on the $i^{th}$ element is

\[ (YY0)_{jkm} = \int^{v_\perp^U}_{v_\perp^L} \varphi_j(v_\perp)\varphi_k(v_\perp)\varphi_m(v_\perp) v_\perp d v_\perp = \int^1_{-1} (c_i + x s_i)l_j(x)l_k(x)l_m(x) s_i d x \]

with $c_i$ and $s_i$ the appropriate shift and scale factors, respectively. Otherwise, if called for any other coordinate elemental matrix YY0 is

\[ (YY0)_{jkm} = \int^{v_\|^U}_{v_\|^L} \varphi_j(v_\|)\varphi_k(v_\|)\varphi_m(v_\|) d v_\| = \int^1_{-1} l_j(x)l_k(x)l_m(x) s_i d x.\]

If called for coord.name = vperp elemental matrix YY1 on the $i^{th}$ element is

\[ (YY1)_{jkm} = \int^{v_\perp^U}_{v_\perp^L} \varphi_j(v_\perp)\varphi_k(v_\perp)\frac{\partial\varphi_m(v_\perp)}{\partial v_\perp} v_\perp d v_\perp = \int^1_{-1} (c_i + x s_i)l_j(x)l_k(x)l_m^\prime(x) d x \]

with $c_i$ and $s_i$ the appropriate shift and scale factors, respectively. Otherwise, if called for any other coordinate elemental matrix YY1 is

\[ (YY1)_{jkm} = \int^{v_\|^U}_{v_\|^L} \varphi_j(v_\|)\varphi_k(v_\|)\frac{\partial\varphi_m(v_\|)}{\partial v_\|} d v_\| = \int^1_{-1} l_j(x)l_k(x)l_m^\prime(x) d x.\]

If called for coord.name = vperp elemental matrix YY2 on the $i^{th}$ element is

\[ (YY2)_{jkm} = \int^{v_\perp^U}_{v_\perp^L} \varphi_j(v_\perp)\frac{\partial\varphi_k(v_\|)}{\partial v_\|}\frac{\partial\varphi_m(v_\perp)}{\partial v_\perp} v_\perp d v_\perp = \int^1_{-1} (c_i + x s_i)l_j(x)l_k^\prime(x)l_m^\prime(x) d x/s_i \]

with $c_i$ and $s_i$ the appropriate shift and scale factors, respectively. Otherwise, if called for any other coordinate elemental matrix YY2 is

\[ (YY2)_{jkm} = \int^{v_\|^U}_{v_\|^L} \varphi_j(v_\|)\frac{\partial\varphi_k(v_\|)}{\partial v_\|}\frac{\partial\varphi_m(v_\|)}{\partial v_\|} d v_\| = \int^1_{-1} l_j(x)l_k^\prime(x)l_m^\prime(x) d x /s_i.\]

If called for coord.name = vperp elemental matrix YY3 on the $i^{th}$ element is

\[ (YY3)_{jkm} = \int^{v_\perp^U}_{v_\perp^L} \varphi_j(v_\perp)\frac{\partial\varphi_k(v_\|)}{\partial v_\|}\varphi_m(v_\perp) v_\perp d v_\perp = \int^1_{-1} (c_i + x s_i)l_j(x)l_k^\prime(x)l_m(x) d x \]

with $c_i$ and $s_i$ the appropriate shift and scale factors, respectively. Otherwise, if called for any other coordinate elemental matrix YY3 is

\[ (YY3)_{jkm} = \int^{v_\|^U}_{v_\|^L} \varphi_j(v_\|)\frac{\partial\varphi_k(v_\|)}{\partial v_\|}\varphi_m(v_\|) d v_\| = \int^1_{-1} l_j(x)l_k^\prime(x)l_m(x) d x.\]

Function that fills and n x n array with the values of the identity matrix I.

Function for finding the elemental index in the global distributed-memory grid. Distributed-memory for global finite-element operators is not yet supported.

Function for computing the scale factor on a grid with uniformed spaced element boundaries. Unused.

Function for setting up the full Gauss-Legendre-Lobatto grid and collocation point weights.

Function for setting up the full Gauss-Legendre-Radau grid and collocation point weights.

Function to create gausslegendre_info struct.

Function that creates the gausslegendre_base_info struct for Lobatto points. If collision_operator_dim = true, assign the elemental matrices used to implement the Fokker-Planck collision operator.

Function that creates the gausslegendre_base_info struct for Lobatto points. If collision_operator_dim = true, assign the elemental matrices used to implement the Fokker-Planck collision operator.

A function that assigns the local matrices to a global array QQ_global for later evaluating strong form of required 1D equation.

The 'option' variable is a flag for choosing the type of matrix to be constructed. Currently the function is set up to assemble the elemental matrices without imposing boundary conditions on the first and final rows of the matrix. This means that the operators constructed from this function can only be used for differentiation, and not solving 1D ODEs. The shared points in the element assembly are averaged (instead of simply added) to be consistent with the derivative_elements_to_full_grid!() function in calculus.jl.

A function that assigns the local weak-form matrices to a global array QQ_global for later solving weak form of required 1D equation.

The 'option' variable is a flag for choosing the type of matrix to be constructed. Currently the function is set up to assemble the elemental matrices without imposing boundary conditions on the first and final rows of the matrix by default. This means that the operators constructed from this function can only be used for differentiation, and not solving 1D ODEs. This assembly function assumes that the coordinate is not distributed. To extend this function to support distributed-memory MPI, addition of off-memory matrix elements to the exterior points would be required.

The typical use of this function is to assemble matrixes M and K in

M * d2f = K * f

where M is the mass matrix and K is the stiffness matrix, and we wish to solve for d2f given f. To solve 1D ODEs

K * f = b = M * d2f

for f given boundary data on f with periodic or dirichlet boundary conditions, set

periodic_bc = true, b[end] = 0

or

dirichlet_bc = true, b[1] = f[1] (except for cylindrical coordinates), b[end] = f[end]

in the function call, and create new matrices for this purpose in the gausslegendre_info struct. Currently the Laplacian matrix is supported with boundary conditions.

Function for computing the shift factor on a grid with uniformed spaced element boundaries. Unused.

Function to perform interpolation on a single element.