# Analysis of tensor-product discontinous Galerkin operators for Vlasov-
Poisson simulations and GPU implementation on Python

D.W. Crews

###### Abstract

The discontinuous Galerkin (DG) finite element method is conservative, lends
itself well to parallelization, and is high-order accurate due to its close
affinity with the theory of quadrature and orthogonal polynomials. When
applied with an orthogonal discretization (i.e. a rectilinear grid) the DG
method may be efficiently implemented on a GPU in just a few lines of high-
level language such as Python. This work demonstrates such an implementation
by writing the DG semi-discrete equation in a tensor-product form and then
computing the products using open source GPU libraries. The results are
illustrated by simulating a problem in plasma physics, namely an instability
in the magnetized Vlasov-Poisson system. Further, as DG is closely related to
spectral methods through its orthogonal basis it is possible to calculate a
transformation to an alternative set of global eigenfunctions for purposes of
analysis or to perform additional operations. This transformation is also
posed as a tensor product and may be GPU-accelerated. In this work a Fourier
series is computed for example (although this does not beat discrete Fourier
transform), and is used to solve the Poisson part of the Vlasov-Poisson system
to $\mathcal{O}(\Delta x^{n+1/2})$-accuracy.

##  1 Introduction

In science and engineering a typical problem is to model the change of some
quantities in time, such as pressure and density in a fluid flow or the
distribution of temperature in a heated object. When the underlying physics
result in linear equations then solutions can be found. Yet in applications
these models usually take the form of nonlinear partial differential equations
(PDEs). To name a few there are the reaction-diffusion systems of chemistry;
the famous Navier-Stokes equations; the Hamilton-Jacobi equation of mechanics
and optimal control; reduced wave equations of optics, seismology, and
oceanography; field equations like Maxwell’s and Einstein’s; Boltzmann’s
kinetic equation; quantum theory’s Schrödinger equation; etc. So it is useful
to solve nonlinear PDEs.

Lacking a well-understood technique to solve nonlinear equations, PDEs are
often simulated numerically by representing the solution on some set of
points. These systems evolve in many dimensions $d$, typically of $d\geq 3$,
so in solving them one experiences the “curse of dimensionality.” This means
that if $N$ points are used to represent the solution in any one direction
then the total number of points is $N^{d}$. To remain computationally
tractable one is limited to a certain total $M$, so the approximately
$M^{1/d}$ points of each dimension should be chosen judiciously to obtain the
most accurate solution, provided by e.g. quadrature nodes.

The discontinuous Galerkin (DG) finite element method combines a prudent
choice of points with a highly parallelizable approach [1]. An excellent
review of DG in a historical context is given in [2]. Many expositions are
precisely formulated and, rightly so, are heavy with mathematical formalisms
such as function spaces and various error estimates [3] [4]. Additionally, the
DG literature pushes the boundaries of the method with hybridizable [5], semi-
Lagrangian [6], superconvergent [7], and space-time [8] innovations. Studies
show impressive performance and scaling of DG-type methods on GPUs [9]. Yet
there seems to be a gap in the recent literature, namely an easy-to-understand
description of a vanilla DG method on a GPU. This article intends to fill this
gap with a cogent explanation of an implementation easily performed on a
desktop with a graphics card. Ideally this will complement the existing
literature for beginners with the method.

The DG method is based on a discontinuous piecewise interpolation through
points in a collection of elements. The PDE is written in an integral (weak)
form and then projected onto the polynomials interpolated by these points.
Being integrated, the points are chosen as the nodes of Gauss-Lobatto
quadrature for a high-order accurate approximation. This means the approximate
solution is on an element-wise orthogonal basis which furthers an intepretable
analytical approach. The projection leads to a first-order ODE in time, the
semi-discrete equation, whose right-hand side is made up of tensor products of
the discretized gradient matrices with the PDE’s flux. When the finite
elements are arranged rectilinearly then these tensor products decouple into
simpler products.

By casting an algorithm as tensor operations, an implementation may be
accelerated by between 10-100 times with parallelization techniques [10].
Parallelization is increasingly indispensable for scientific computing as its
technology evolves towards increased memory and processing power.
Supercomputing clusters use message-passing interface (MPI) or multiple
graphics processing units (GPUs), and most of today’s desktop computers have a
GPU. On clusters, MPI-based simulations using parallelized tensor products
have been run with a million processors [11] and research is examining MPI
scaling versus parallelism on multiple GPUs [12]. Yet even consumer-grade GPUs
are attaining sufficient on-board memory for solution of relatively large
PDEs, with up to 24 GB of dedicated memory on the recently-released NVIDIA RTX
3000’s. This is an exciting development as it allows anyone to solve PDEs in
parallel on a simple desktop machine.

As this article concerns rectilinear grids, it is also worthwhile to
thoroughly review the most commonly used basis functions. The DG method is
often applied on unstructured meshes with simplex elements and the integrals
needed to obtain the gradient discretization matrices are challenging, so they
are generally computed numerically [1][13]. It was recently discovered in the
context of astrophysics calculations on a curvilinear grid that their
integrals may be found in explicit forms for the Legendre-Gauss-Lobatto (LGL)
quadrature rule [14][15]. The LGL interpolation polynomials are found along
the way to be naturally expanded in Legendre series.

This Legendre spectral expansion connects the DG method to the theory of
spherical harmonics, and the connection is useful in this work for Fourier
transformation. To understand this, first note that Gauss-Legendre quadrature
is of unit weighting because the measure of a spherical surface projects
uniformly onto an interval (its axis) due to a theorem of Archimedes [16]. As
Legendre polynomials are those harmonics invariant under rotations of the
polar axis, the quadrature shares their completeness properties, i.e. the
relation [17]

$\sum_{s=0}^{\infty}\Big{(}s+\frac{1}{2}\Big{)}P_{s}(x)P_{s}(y)=\delta(y-x),\quad\quad-1\leq
x,y\leq 1.$ (1)

Specifically, the interpolation polynomials $\ell_{j}(\xi)$ through $n$ Gauss-
Legendre nodes satisfy

$\ell_{j}(\xi)=w_{j}\sum_{s=0}^{n-1}\Big{(}s+\frac{1}{2}\Big{)}P_{s}(\xi_{j})P_{s}(\xi),\quad\text{
and }\quad\ell_{j}(\xi_{i})=\delta_{ij}$ (2)

with $\\{w_{j}\\}_{j=0}^{n-1}$, $\\{\xi_{j}\\}_{j=0}^{n-1}$ the quadrature
weights and nodes respectively111Considering the summation
$\lim_{j\to\infty}\sum_{j=0}^{n-1}\ell_{j}(\xi)f(\xi)$ suggests that the the
interpolation polynomial plays a role in discrete integration as if it were a
$\delta$-sequence on $[-1,1]$ in the sense of distributions [18].. The LGL
interpolant is identical to this form with a reweighted final eigenvalue
$(s+1/2)\to(s/2)$ [14]. With this relation the DG gradient matrices in the LGL
basis may be found explicitly [14][15].

This article is structured as follows. Section 2 reviews the DG projection for
PDEs and demonstrates that for a rectilinear grid the elements may be
represented as tensor products of one-dimensional elements. Consequently only
the one-dimensional discretization matrices are needed. Using these results,
Section 3 demonstrates a simple Eulerian implementation of computing the semi-
discrete equation in $n$-dimensions with GPU acceleration using upwind
numerical fluxes for example. Following this, Section 4 introduces an example
problem by posing a problem in plasma physics, a magnetized plasma instability
present in the Vlasov-Poisson system. Similar problems in plasma theory are
often approached by solving the Vlasov equation with DG method [19] [20] [21]
[22]. It should be noted that the approach advocated here will work for
systems describing flows in any (reasonable) number of dimensions. For
instance by changing variables, flux, and boundary conditions in the example
one can solve the Euler equations in three dimensions.

Having posed the example problem, Section 5 then reviews the DG basis matrices
in the Lobatto basis. The basis analysis is then extended in Section 5.3.2 to
a Fourier spectral method used to solve the Poisson part of Vlasov-Poisson.
Simulation results are discussed in Section 6 with an emphasis on the
computational performance achieved using an RTX 3090 desktop GPU. Finally,
Appendix 4.1 contains a short discussion of the initial condition used for the
example problem. All code used to produce the work shown in this article can
be found on the author’s GitHub at https://github.com/crewsdw/Vlasov1D2V.

##  2 The DG method for rectilinear (tensor-product) grids

Consider the conservation law for a scalar $u=u(t,x_{i})$ and its flux
$F^{j}=F^{j}(x_{i},u)$,

$\partial_{t}u+\partial_{j}F^{j}=0.$ (3)

Throughout this section the up/down summation notation is used. Equation 3
also describes vector-valued $u$ and higher-order derivatives reduced to
first-order systems. The DG method splits the domain $\Omega$ into finite
elements $\Omega_{\alpha}$ with boundary $\partial\Omega_{\alpha}$, and
projects the weak or variational form of Eqn. 3 onto a polynomial basis of
each element, resulting in an approximate solution $\widetilde{u}^{\alpha}$. A
set of nodes is chosen $\\{\xi_{j}\\}_{j=0}^{n-1}$ within each element which
are interpolated by the Lagrange polynomials

$\ell_{j}(\xi)=\prod_{\begin{subarray}{c}k=0\\\ k\neq
j\end{subarray}}^{n-1}\frac{\xi-\xi_{k}}{\xi_{j}-\xi_{k}}.$ (4)

The element basis is defined as these polynomials $\ell_{j}(\xi)$. Rectilinear
$d$-dimensional elements $\Omega$ are built up by tensor products of the one-
dimensional line element $\mathcal{L}$, so that
$\Omega=\mathcal{L}_{d}\otimes\mathcal{L}_{d-1}\otimes\cdots\otimes\mathcal{L}_{1}$.
The node set of $\Omega$ is correspondingly a tensor product of each line’s
node set, so the $d$-dimensional Lagrange functions can be factorized into
products of the Lagrange polynomials in each dimension. To avoid a
proliferation of indices, define an index $\beta_{i}$ to range through the
nodes $\\{\xi_{j}\\}_{j=0}^{n-1}$ in the line element $\mathcal{L}_{i}$, and
then let $\beta=\\{\beta_{0},\beta_{1},\cdots,\beta_{d-1}\\}$ be the multi-
index collecting the nodal index of each dimension. Additional discussion on
ordering of a multi-index for tensor-product geometry constructions can be
found in [23]. The Lagrange polynomial of $\Omega$ is then

$L_{\beta}(x_{0},x_{1},\cdots,x_{d})=\ell_{\beta_{0}}(x_{0})\ell_{\beta_{1}}(x_{1})\cdots\ell_{\beta_{d}}(x_{d})$
(5)

with $x_{i}$ the i’th coordinate. In the projection of the approximation onto
the basis functions $L_{\beta}$ the expansion coefficients $u^{\alpha,\beta}$
may be labelled by the same multi-index,

$\widetilde{u}^{\alpha}(t,x)=\sum_{\beta}u^{\alpha,\beta}(t)L_{\beta}(x).$ (6)

Galerkin projection of Eq. 3 leads to inner products
$\langle\cdot|\cdot\rangle_{\Omega_{\alpha}}$ over the basis functions.
Specifically, Eq. 3 is integrated on the domain against the basis functions
$L_{\gamma}$,

$\langle\partial_{t}u|L_{\gamma}\rangle_{\Omega_{\alpha}}+\langle\partial_{j}F^{j}|L_{\gamma}\rangle_{\Omega_{\alpha}}=0,$
(7)

and the basis expansion Eq. 6 is substituted. The resulting projected integral
form is

$\frac{du^{\alpha,\beta}}{dt}\langle
L_{\beta}|L_{\gamma}\rangle_{\Omega_{\alpha}}+\langle\partial_{j}F^{\alpha,j}|L_{\gamma}\rangle_{\Omega_{\alpha}}=0.$
(8)

Information must be passed between elements $\Omega_{\alpha}$ in the form of a
flux, or else the elements would decouple and the discretization would not be
consistent. This is done by integrating the flux term
$\langle\partial_{j}F^{\alpha,j}|L_{\gamma}\rangle_{\Omega_{\alpha}}$ by parts
and taking the boundary term to be a function of both the local state and that
of the neighbor, i.e.
$F^{\alpha}|_{\partial\Omega_{\alpha}}\equiv\mathcal{F}(\partial\Omega_{\alpha}^{-},\partial\Omega_{\alpha}^{+})$
where $\partial\Omega_{\alpha}^{-}$ is the state interior to the element and
$\partial\Omega_{\alpha}^{+}$ the state exterior. The function $\mathcal{F}$,
termed numerical flux, is chosen so that Eq. 3 is discretized consistently.
The requirement for consistency means $\mathcal{F}$ is related to the nature
of information propagation in the system. The single integration by parts
results in the DG weak form, and one further integration by parts in the
strong form, given respectively by

$\displaystyle\frac{du^{\alpha,\beta}}{dt}\langle
L_{\beta}|L_{\gamma}\rangle_{\Omega_{\alpha}}$
$\displaystyle=F^{\alpha,\beta,j}\langle
L_{\beta}|\partial_{j}L_{\gamma}\rangle_{\Omega_{\alpha}}-\mathcal{F}^{\alpha,\beta,j}\langle
L_{\beta}|L_{\gamma}\rangle_{\partial\Omega^{j}_{\alpha}},$ (9)
$\displaystyle\frac{du^{\alpha,\beta}}{dt}\langle
L_{\beta}|L_{\gamma}\rangle_{\Omega_{\alpha}}+F^{\alpha,\beta,j}\langle\partial_{j}L_{\beta}|L_{\gamma}\rangle_{\Omega_{\alpha}}$
$\displaystyle=(F^{\alpha,\beta,j}-\mathcal{F}^{\alpha,\beta,j})\langle
L_{\beta}|L_{\gamma}\rangle_{\partial\Omega^{j}_{\alpha}}$ (10)

These inner product matrices are termed the mass $M$, face mass $\Gamma$,
advection $\mathcal{A}$, and stiffness $S$ matrices by analogy with continuum
mechanics, and their short-hand definitions are

$\displaystyle M^{\alpha}_{\beta\gamma}$ $\displaystyle\equiv\langle
L_{\beta}|L_{\gamma}\rangle_{\Omega_{\alpha}},\quad\quad\Gamma^{\alpha,j}_{\beta\gamma}\equiv\langle
L_{\beta}|L_{\gamma}\rangle_{\partial\Omega^{j}_{\alpha}},$ (11)
$\displaystyle\mathcal{A}_{\beta\gamma}^{\alpha,j}$
$\displaystyle\equiv\langle
L_{\beta}|\partial_{j}L_{\gamma}\rangle_{\Omega_{\alpha}},\quad
S_{\beta\gamma}^{\alpha,j}\equiv\langle\partial_{j}L_{\beta}|L_{\gamma}\rangle_{\Omega_{\alpha}}.$
(12)

Now solving for the expansion coefficients results in element-wise operators
composed of the products

$\displaystyle\Upsilon^{\alpha,\beta}_{\gamma,j}\equiv(M^{-1})^{\alpha,\beta\epsilon}\mathcal{A}_{\gamma\epsilon,j}^{\alpha},\quad\quad\Xi^{\alpha,\beta}_{\gamma,j}\equiv(M^{-1})^{\alpha,\beta\epsilon}\Gamma^{\alpha}_{\gamma\epsilon,j},\quad\quad
D^{\alpha,\beta}_{\gamma,j}\equiv(M^{-1})^{\alpha,\beta\epsilon}S_{\gamma\epsilon,j},$
(13)

resulting in the weak and strong form semi-discrete equations per element,

$\displaystyle\frac{d}{dt}u^{\alpha,\beta}$
$\displaystyle=\Upsilon^{\alpha,\beta}_{\gamma,j}F^{\alpha,\gamma,j}-\Xi^{\alpha,\beta}_{\gamma,j}\mathcal{F}^{\alpha,\gamma,j}$
(14)
$\displaystyle\frac{d}{dt}u^{\alpha,\beta}+D^{\alpha,\beta}_{\gamma,j}F^{\alpha,\gamma,j}$
$\displaystyle=\Xi^{\alpha,\beta}_{\gamma,j}(F^{\alpha,\gamma,j}-\mathcal{F}^{\alpha,\gamma,j}).$
(15)

Considering the weak form, Eqn. 14, the first term represents fluxes due to
internal degrees of freedom within an element, and the second term boundary
fluxes. The strong form operator $D^{\alpha,\beta}_{\gamma,\mu}$ is a direct
gradient discretization called the derivative matrix, while the operator
$\Upsilon^{j,\alpha}_{k,\mu}$ approximates the gradient in integral form. In
both cases $\Xi^{\alpha,\beta}_{\gamma,\mu}$ discretizes the surface integrals
between elements. Explicitly, the multi-index $\alpha$ denotes the elements
$\Omega_{\alpha}$, $\beta$ the element’s nodes, and $j$ the coordinates.

The basis of element $\Omega_{\alpha}$ is related to a reference element by an
isoparametric transform and its Jacobian $J^{\alpha}_{ij}$, so that only one
set of matrices
$\\{\Upsilon^{\beta}_{\gamma,j},\Xi^{\beta}_{\gamma,j},D^{\beta}_{\gamma,j}\\}$
need be determined. These reference operators are then related to element-wise
ones by
$\Upsilon^{\alpha,\beta}_{\gamma,\nu}=J^{\mu,\alpha}_{\nu}\Upsilon^{\beta}_{\gamma,\mu}$.
The following will consider operations with only the reference operators by
taking elements to have an identical Jacobian. Equations 14 & 15 represent
ODEs whose right-hand side is given by a tensor product contracting the multi-
index $\gamma$ and the coordinates $j$. The following section shows that for
rectilinear elements these contractions simplify into a sum over products with
the one-dimensional matrices $\Upsilon,D,\Xi$.

###  2.1 Constituent reduction for rectilinear (tensor-product) elements

Now consider rectilinear $n$-cube elements specifically, with the multi-index
$(\alpha,\beta)$ as shown in Fig. 1. Substitution of the factorized basis
polynomials defined by Eq. 5 into Eqs. 11-12 shows the resulting object to be
composed of products of the lower-dimensional constituents. Specifically, the
mass matrix $M_{\beta,\gamma}$ with sub-element nodal multi-indices
$\beta,\gamma$ decomposes into tensor products as

$M_{\beta\gamma}=M_{\beta_{d-1},\gamma_{d-1}}\otimes
M_{\beta_{d-2},\gamma_{d-2}}\otimes\cdots\otimes M_{\beta_{0},\gamma_{0}}.$
(16)

so that it may be termed the mass tensor for $d>1$. It follows that the
inverse mass tensor is
$M^{-1}_{\beta\gamma}=M^{-1}_{\beta_{d-1},\gamma_{d-1}}\otimes
M^{-1}_{\beta_{d-2},\gamma_{d-2}}\otimes\cdots\otimes
M^{-1}_{\beta_{0},\gamma_{0}}$. Similarly, the advection tensor
$\mathcal{A}_{\beta\gamma}^{j}=\langle\partial_{j}L_{\beta}|L_{\gamma}\rangle_{\Omega_{\alpha}}$
has $d$ components (ranging over index $j$) each given by the tensor products

$\mathcal{A}_{\beta\gamma}^{j}=\begin{bmatrix}M_{\beta_{d-1},\gamma_{d-1}}\otimes
M_{\beta_{d-2},\gamma_{d-2}}\otimes\cdots\otimes\mathcal{A}_{\beta_{0},\gamma_{0}}\\\
\vdots\quad\quad\quad\quad\vdots\quad\quad\quad\vdots\\\
M_{\beta_{d-1},\gamma_{d-1}}\otimes\mathcal{A}_{\beta_{d-2},\gamma_{d-2}}\otimes\cdots\otimes
M_{\beta_{0},\gamma_{0}}\\\ \mathcal{A}_{\beta_{d-1},\gamma_{d-1}}\otimes
M_{\beta_{d-2},\gamma_{d-2}}\otimes\cdots\otimes
M_{\beta_{0},\gamma_{0}}\end{bmatrix}$ (17)

with analogous results for the face mass $\Gamma$ and stiffness matrices $S$.
The components corresponding to the sub-indices of the multi-indices
$\beta,\gamma$ are seen to be the lower-dimensional matrices themselves. Due
to the tensor product mixed-product property $(A\otimes B)(C\otimes
D)=(AC)\otimes(BD)$, the internal flux tensor
$\Upsilon_{\mathcal{D}}^{\mu}=M^{-1}_{\mathcal{D}}A^{\mu}_{\mathcal{D}}$ has
the $d$ components,

$\Upsilon^{j}_{\beta\gamma}=\begin{bmatrix}I_{\beta_{d-1},\gamma_{d-1}}\otimes
I_{\beta_{d-2},\gamma_{d-2}}\otimes\cdots\otimes\Upsilon_{\beta_{0},\gamma_{0}}\\\
\vdots\quad\quad\quad\quad\vdots\quad\quad\quad\vdots\\\
I_{\beta_{d-1},\gamma_{d-1}}\otimes\Upsilon_{\beta_{d-2},\gamma_{d-2}}\otimes\cdots\otimes
I_{\beta_{0},\gamma_{0}}\\\ \Upsilon_{\beta_{d-1},\gamma_{d-1}}\otimes
I_{\beta_{d-2},\gamma_{d-2}}\otimes\cdots\otimes
I_{\beta_{0},\gamma_{0}}\end{bmatrix},$ (18)

each consisting of $(d-1)$ tensor products with the identity. The numerical
flux tensor $\Xi^{j}_{\beta\gamma}$ follows an identical pattern, with
component $k$ as $I_{\beta_{d-1},\gamma_{d-1}}\otimes
I_{\beta_{d-2},\gamma_{d-2}}\otimes\cdots\otimes\Xi_{\beta_{k},\gamma_{k}}\otimes\cdots\otimes
I_{\beta_{0},\gamma_{0}}$.

Equation 18 means that the flux of direction $j$ only needs to be contracted
with the corresponding one-dimensional matrix of direction $j$. To illustrate,
consider the case of a $d=3$ flux $F^{j}$ with a rectilinear discretization.
Separate the element multi-index $\alpha$ into directional indices $a,b,c$ and
$\beta$ into sub-element nodal indices $p,q,r$ so that
$F^{\alpha}_{j,\gamma}\equiv F^{abc}_{j,pqr}$. The internal flux product is
then

$\displaystyle\Upsilon^{\gamma,j}_{\beta}F^{\alpha}_{j,\gamma}$
$\displaystyle\equiv(I_{2})^{p}_{n}(I_{1})^{q}_{m}(\Upsilon_{0})^{r}_{\ell}F^{abc}_{0,pqr}+(I_{2})^{p}_{n}(\Upsilon_{1})^{q}_{m}(I_{0})^{r}_{\ell}F^{abc}_{1,pqr}+(\Upsilon_{2})^{p}_{n}(I_{1})^{q}_{m}(I_{0})^{r}_{\ell}F^{abc}_{0,pqr}$
(19)
$\displaystyle=(\Upsilon_{0})^{r}_{\ell}F^{abc}_{0,nmr}+(\Upsilon_{1})^{q}_{m}F^{abc}_{1,nq\ell}+(\Upsilon_{2})^{p}_{n}F^{abc}_{2,pm\ell}$
(20)

by carrying through each product with the identity. By similar reasoning,
having formed a numerical flux array $\mathcal{F}$ for the two faces of each
sub-element axis, the boundary flux product is simply

$\Xi^{j,\gamma}_{\beta}\mathcal{F}^{\alpha}_{j,\gamma}=(\Xi_{0})^{\gamma_{r}}_{\ell}\mathcal{F}^{abc}_{0,nm\gamma_{r}}+(\Xi_{1})^{\gamma_{q}}_{m}\mathcal{F}^{abc}_{1,n\gamma_{q}r}+(\Xi_{2})^{\gamma_{p}}_{n}\mathcal{F}^{abc}_{2,\gamma_{p}m\ell}$
(21)

where $\gamma_{i}=\delta_{0i}+\delta_{i(n-1)}$ picks out the boundary nodes of
the sub-element index.

Figure 1: Rectilinear elements have a natural tensor product structure which
is illustrated by the use of multi-indexes. Shown is a two-dimensional example
with an eight-node Legendre-Gauss-Lobatto (LGL) node set, where the elements
in directions $(0,1)$ are indexed by $\alpha_{0},\alpha_{1}$ and the nodal
sub-element positions by $\beta_{0},\beta_{1}$. Any given node’s index
consists of its tuple $(\alpha_{0},\beta_{0},\alpha_{1},\beta_{1})$.

Equation 20 demonstrates that orthogonal elements reduce the number of
calculations needed to evaluate the right-hand side of the semi-discrete
equation. The gain increases with dimension versus using full basis tensors.
With full tensors, $\Upsilon$ is a $\mathbb{R}^{d\times N\times N}$ array for
$N\approx n^{d}$ the total number of nodes in the element and $n$ the nodes
per sub-element. In contrast, with a rectilinear discretization only the $d$
one-dimensional matrices $\Upsilon\in\mathbb{R}^{n\times n}$ are needed.
Unless a non-orthogonal mesh is required for particular domain geometry, an
orthogonal grid provides significant simplification.

In some cases non-orthogonal elements are desired in certain dimensions of a
problem while other dimensions are free to use orthogonal discretization [24].
For example the collisionless Boltzmann equation is an advection equation in
$d=(D+V)$-dimensional phase space consisting of $D$ configurational (spatial)
dimensions and $V$ velocity dimensions. As irregular boundaries will occur
only in spatial dimensions, one is free to discretize velocity space
orthogonally. In this case a blended element basis $E$ may be built up from
tensor products of the simplex element basis $\Delta$ and $V$ copies of the
line element $\mathcal{L}_{i}$, i.e. for $D=V=2$ one has
$E\equiv\Delta\otimes\mathcal{L}_{1}\otimes\mathcal{L}_{2}$. In this case

$\Upsilon_{E}^{j}=\begin{bmatrix}I_{2}\otimes I_{1}\otimes\Upsilon_{\Delta}\\\
I_{2}\otimes\Upsilon_{1}\otimes I_{\Delta}\\\ \Upsilon_{2}\otimes I_{1}\otimes
I_{\Delta}\end{bmatrix}$ (22)

and similarly for the numerical flux array. The multi-index is naturally
composed of unstructured spatial element information, velocity-space grid
indices $a,b$, and sub-element nodal indices $p,q$. Equation 20 then
generalizes to the blended case.

##  3 GPU implementation with a high-level language

This section details an efficient implementation of the DG method using GPU-
accelerated Python libraries for CUDA devices. Any language which can
efficiently compute tensor products of large arrays can implement DG well;
this example utilizes Python and the CuPy library [25]. This library mimics
NumPy data array structures and may be programmed in only a few lines of code
to launch optimized kernels at execution for the tensor products appearing in
the DG semi-discrete equations. Whichever library is used, its operations
should occur entirely on device memory so that CPU-GPU transfers are only
needed to saving data. Such a GPU-centered DG implementation is shown
schematically in Fig 2. Finally, the DG basis matrices are calculated via the
Eqs. 54 and 57 in Section 5 using i) a scientific library with Legendre
polynomials such as SciPy, and ii) tables of Legendre-Gauss-Lobatto quadrature
nodes and weights available online, for example at [26].

Figure 2: For performance of a parallelized DG implementation it’s critical
that all operations in the main time-stepping loop occur on the device (GPU).
Device-host memory transfer should be minimized to the greatest extent
possible. Transferring the array is only necessary at whichever time
increments data is to be saved to the disk.

###  3.1 Evaluating the semi-discrete equation with tensor products

As discussed in Section 2, when solving the semi-discrete equation on an
orthogonal grid one need only contract the flux $F^{\alpha,\beta,j}$ of each
direction $j$ with its corresponding sub-element index. With a NumPy-like data
array structure this is accomplished by structuring the array data
$u^{\alpha,\beta}$ by its natural tensor-product structure. As in Fig. 1, each
pair of indices for elements $\alpha_{j}$ and nodes $\beta_{j}$ corresponds to
the direction $j$, so a natural array ordering is u[$\alpha_{0}$, $\beta_{0}$,
$\alpha_{1}$, $\beta_{1}$, $\cdots$, $\alpha_{d-1}$, $\beta_{d-1}$] and fluxes
$F^{\alpha,\beta,j}=F^{j}(u)$ index in the same manner per dimension $j$.

Structuring the data array by both element and sub-element indices allows the
semi-discrete equation to be elegantly programmed using GPU-accelerated tensor
product operations such as cupy.tensordot(a, b, axes=), which compiles a just-
in-time GPU kernel to compute a tensor product with contraction on specified
axes like NumPy’s tensordot function. For example, consider evaluating the
internal flux product $\Upsilon^{\beta}_{\gamma,j}F^{\alpha,\gamma,j}$ in the
DG weak-form of Eqn. 14. As the product decouples in an orthogonal
discretization (Eqn. 20) it may be computed as a sum of tensor products for
each direction using the basis_product function,

        def basis_product(flux, basis_matrix, axis, permutation):
            return cupy.transpose(cupy.tensordot(flux, basis_matrix,
                                  axes=([axis], [0])), axes=permutation)

where axis gives the sub-element index of the direction of flux, and the
function transpose returns a view to the original index order according to the
tuple permutation as tensordot concatenates the tensor indices remaining after
contraction of the specified axes.

####  3.1.1 Computing numerical fluxes

The function basis_product is also used to globally compute the numerical flux
product $\Xi^{\beta}_{\gamma,j}\mathcal{F}^{\alpha,\gamma,j}$, by first
arranging a numerical flux array $\mathcal{F}^{\alpha,\gamma,j}$ according to
the tensor-product index ordering. The form of the function $\mathcal{F}$
depends on the chosen numerical flux scheme. The upwind method, where
$\mathcal{F}$ is a function of only one side of the interface, is well-suited
for hyperbolic problems as it solves the underlying interface, or Riemann,
problem [27]. As a unidirectional numerical flux the upwind method is also
simple to adapt for the alternating fluxes used for parabolic problems [28].

The idea behind upwind flux is for information to flow across an element’s
boundaries in accordance with the PDE’s dispersion relation. For example, in
the problem $u_{t}+Vu_{x}=0$ the dispersion relation is $\omega-Vk=0$ and
information flows to the right. Therefore information should be sent from left
to right, and numerical flux would always be chosen as depending on values
from the left.

Let $R$ stand for right, and $L$ stand for left. Defining auxiliary variables
$a^{R}=\text{max}(a,0)$ and $a^{L}=\text{min}(a,0)$ only one of which is
nonzero, the upwind flux for $F=au$ with advection speed $a$ is

$\mathcal{F}(u^{R},u^{L})=a^{R}u^{L}+a^{L}u^{R}.$ (23)

An element-wise implementation of the upwind method on GPU, a bit more lengthy
than the internal flux product, is discussed in Appendix A.

##  4 An example problem: the Vlasov-Poisson system

A representative problem to illustrate the DG implementation is the Vlasov-
Poisson system of plasma physics. These equations describe a probability
distribution of charged particles evolving under a phase space transport
equation, self-consistently coupled through its zeroth moment to Gauss’s law
for the electric field as represented in potential form by Poisson’s equation.
The system is widely used in the study of low-frequency waves in hot plasmas
[29] and is often the starting point for derivation of reduced models
describing turbulent plasma dynamics such as the gyrokinetic and quasilinear
approximations [30]. For a distribution of electrons $f(\bm{x},\bm{v},t)$
interacting through the electric potential $\Phi(x)$ in the presence of a
fixed, neutralizing background, the system consists of

$\displaystyle\partial_{t}f$
$\displaystyle+\bm{v}\cdot\nabla_{x}f+\frac{e}{m}(\nabla_{x}\Phi-\bm{v}\times\bm{B}_{\text{ext}})\cdot\nabla_{v}f=0$
(24) $\displaystyle\nabla^{2}\Phi$
$\displaystyle=\frac{en_{0}}{\epsilon_{0}}\Big{(}\int_{-\infty}^{\infty}fdv-1\Big{)}$
(25)

where $e,m$ represent electron charge and mass respectively, $\bm{B}_{ext}$ is
an external magnetic field, and $n_{0}$ is a reference particle density
normalizing the distribution $f$.

A useful model problem to benchmark an implementation of the system in one
spatial and two velocity dimensions, i.e. $f(\bm{x},\bm{v})=f(x,u,v)$, is the
instability of a loss-cone distribution to perpendicular-propagating
cyclotron-harmonic waves [20]. By normalizing length to the Debye length
$\lambda_{D}$, time to the plasma frequency $\omega_{p}^{-1}$, velocities to
the thermal velocity $v_{th}=\lambda_{D}\omega_{p}$, and the fields by setting
$E_{0}=v_{th}B_{0}$ such that the reference plasma-cyclotron frequency ratio
is unity $(\omega_{p}/\omega_{c})_{0}=1$, the resulting normalized system is
given by the equations

$\displaystyle\partial_{t}f$ $\displaystyle+F^{j}\partial_{j}f=0,$ (26)
$\displaystyle\partial_{xx}\Phi$ $\displaystyle=\int_{-\infty}^{\infty}fdv-1,$
(27) $\displaystyle F$ $\displaystyle=\begin{bmatrix}u,&\partial_{x}\Phi-
vB_{\text{ext}},&uB_{\text{ext}}\end{bmatrix}^{T}.$ (28)

Here the normalized external magnetic field $B_{\text{ext}}$ is set such that
the desired plasma-cyclotron frequency ratio is
$\omega_{p}/\omega_{c}=B_{\text{ext}}^{-1}$ in normalized units.

###  4.1 Initial condition for the Vlasov-Poisson system

The following subsection is particular to plasma physics, and the uninitiated
reader should skip to the end. To determine an appropriate initial condition,
the fundamental modes of Eqs. 24 & 25 are determined by their linearization
about zero-order cyclotron orbits followed by Fourier transformation in the
polar velocity coordinates $(u,v)\to(v_{\perp},\varphi)$ to yield the
linearized solution [29] [31],

$f_{1}(\omega,k,v_{\perp},\varphi)=\Phi(\omega,k)e^{ikv_{\perp}\sin(\varphi)}\Big{(}\sum_{n=-\infty}^{\infty}\frac{n}{n-\omega}J_{n}(kv)e^{-in\varphi}\Big{)}\frac{1}{v_{\perp}}\frac{\partial
f_{0}}{\partial v_{\perp}}$ (29)

with frequency $\omega$ normalized to the cyclotron frequency $\omega_{c}$ and
wavenumber to the Larmor radius $r_{L}=v_{t}/\omega_{c}$. Combination of the
zeroth moment $\int_{-\infty}^{\infty}f_{1}dv$ together with the Fourier-
transformed Poisson equation yields the Harris dispersion relation, here
written in the integral form [20] [32]

$\epsilon(\omega,k)\equiv
1+\Big{(}\frac{\omega_{p}}{\omega_{c}}\Big{)}^{2}\frac{1}{\sin(\pi\omega)}\int_{0}^{\pi}\sin(\theta)\sin(\omega\theta)\mathcal{H}_{0}[f_{0}](\lambda(\theta))d\theta=0$
(30)

where $\mathcal{H}_{0}[f](q)$ is the zero-order Hankel transform of $f(v)$ and
the parameter $\lambda\equiv 2k\cos(\frac{1}{2}\theta)$. A zero-order function
commonly used to model loss-cone-like probability distributions is the
function

$f_{0}(x,v_{\perp},\varphi)=\frac{1}{2\pi\alpha^{2}j!}\Big{(}\frac{v_{\perp}^{2}}{2\alpha^{2}}\Big{)}^{j}\exp\Big{(}-\frac{v_{\perp}^{2}}{2\alpha^{2}}\Big{)}$
(31)

with ring parameter $j$ and radially-normalized gradient,

$v_{\perp}^{-1}\partial_{v_{\perp}}f_{0}=\frac{1}{v_{\perp}^{2}}\Big{(}j-\frac{v_{\perp}^{2}}{2\alpha^{2}}\Big{)}f_{0}.$
(32)

Equation 31 describes a ring distribution, and its thermal properties are
summarized in Appendix LABEL:sec:therm_appendix. From the Fourier-multiplier
property
$\mathcal{H}_{0}[(v^{2})^{n}f(v)](q)=(-\nabla^{2}_{q})^{n}\mathcal{H}_{0}[f(v)](q)$,
and the fact that the coefficients of successive radial Laplacians of the
polar Gaussian are the Laguerre polynomials $L_{j}(x)$, it follows that the
transform of Eqn. 31 is given by

$\mathcal{H}_{0}[f_{0}](\lambda)=L_{j}\Big{(}\frac{1}{2}\alpha^{2}\lambda^{2}\Big{)}\exp\Big{(}-\frac{1}{2}\alpha^{2}\lambda^{2}\Big{)}.$
(33)

The spatio-temporal modes of the linearized Vlasov-Poisson system are then
solutions of

$\epsilon(\omega,k)=1+\Big{(}\frac{\omega_{p}}{\omega_{c}}\Big{)}^{2}\frac{1}{\sin(\pi\omega)}\int_{0}^{\pi}\sin(\omega\theta)\sin(\theta)L_{j}(\beta)e^{-\beta}d\theta=0$
(34)

where $\beta=2k^{2}\cos^{2}(\frac{1}{2}\theta)$. The dispersion function
$\epsilon(\omega,k)$ may be written in closed form using hypergeometric
functions of the form ${}_{2}F_{2}$. Yet computing them requires a power
series and the trigonometric integral may be accurately computed using a
quadrature method. A fifty-point Gauss-Legendre method was used for this work.

Equation 34 has many solutions representing the cyclotron harmonic modes of
the plasma [32]. With ring parameter $j=6$ and frequency ratio
$\omega_{p}/\omega_{c}=10$, an unstable solution of maximum growth-rate occurs
at $(k_{0},\omega_{0})\sim(0.886,0.349i)$, for which the corresponding phase
space mode is found by inverse-transforming Eq.29 with
$\Phi(k,\omega)=\Phi_{0}\delta(k-k_{0})\delta(\omega-\omega_{0})$. The
resulting mode can be written as a Fourier series in the cylindrical angle
$\varphi$,

$\displaystyle f_{1}(x,v_{\perp},\varphi)$
$\displaystyle\sim\frac{1}{v_{\perp}}\frac{\partial f_{0}}{\partial
v_{\perp}}\text{Re}(\psi_{0})$ (35)
$\displaystyle\psi_{0}(x,v_{\perp},\varphi)$ $\displaystyle\equiv
e^{i(k_{0}x+k_{0}v_{\perp}\sin(\varphi))}\sum_{n=-\infty}^{\infty}\frac{n}{n-\omega_{0}}J_{n}(k_{0}v_{\perp})e^{-in\varphi}$
(36)

where a partial sum truncated near the approximate harmonic mode number
results in a good approximation. For a well-converged sum this work uses an
$n_{\text{max}}=20$-term approximation.

The zeroth moment of Eq. 35 is a single harmonic $(\delta
n)\sin(k_{0}x+\theta)$ of amplitude $\delta n$ with some phase shift $\theta$
(though if desired a centered mode may be found by also considering the
conjugate wavenumber’s solution at $k_{0}^{\prime}=-k_{0}$). The initial
condition $f_{1}$ should be multiplied by a scaling factor so that the
perturbation amplitude $A$ is small enough for a valid linearization. Based on
Eq. 26, an estimate of the validity condition is $E_{1}\ll
v_{th}B_{\text{ext}}$. Combining with Poisson’s equation, it follows that the
condition for a valid linearization is

$\frac{\delta
n}{n_{0}}\ll\Big{(}\frac{\omega_{c}}{\omega_{p}}\Big{)}(k\lambda_{D})\sim\Big{(}\frac{\omega_{c}}{\omega_{p}}\Big{)}^{2}(kr_{L}).$
(37)

This condition can be quite restrictive. For instance, the initial condition
described above must satisfy $\delta n/n_{0}\ll 0.0089$, or else a faster
growing two stream-like instability will be observed.

##  5 Element-wise operators in the Lobatto basis

This section reviews the theory of DG matrices in the semi-discrete equation
according to the method of [14] [15]. A review of this theory is useful as the
solution is both explicit and interpretable to arbitrarily high order. A
novice reader may wish to skip to the results given in Eqs. 54, 57, and 58.

###  5.1 Basic patterns in the Gauss-Legendre basis

It’s useful to first review the result for full Gauss-Legendre quadrature as
the steps are simple and the patterns reoccur in the Lobatto basis. The
situation is simple in full Gaussian quadrature, but easy evaluation of the
numerical fluxes is greatly simplified by Lobatto quadrature by including the
endpoints [1]. Denoting the node locations by the $n$ quadrature points
$\\{\xi_{j}\\}_{j=0}^{n-1}$, in addition to continuous orthogonality the
Legendre polynomials are discretely orthogonal [33] [34],

$\sum_{s=0}^{n-1}w_{s}P_{k}(\xi_{s})P_{\ell}(\xi_{s})=\frac{1}{k+\frac{1}{2}}\delta_{k\ell}$
(38)

with $w_{s}$ the quadrature weights. This follows from the order $2n-1$
quadrature property as the sum corresponds to the continuous integral
$\int_{-1}^{1}P_{k}(x)P_{\ell}(x)dx$. Expansions in the Lagrange basis of Eqn.
4 and in the Legendre basis are linearly related due to the interpolation
property $\ell_{j}(\xi_{i})=\delta_{ij}$,

$f(\xi)=\sum_{j=0}^{n-1}f_{j}\ell_{j}(\xi)=\sum_{k=0}^{n-1}c_{k}P_{k}(\xi),\quad\implies\quad
f_{j}=\mathcal{V}_{j}^{k}c_{k}$ (39)

where $\mathcal{V}^{k}_{j}=P_{k}(\xi_{j})$ is termed the (generalized)
Vandermonde matrix. The inverse transform follows directly from Eqn. 38 as
$(\mathcal{V}^{-1})_{j}^{k}=w_{j}P_{k}(\xi_{j})(k+\frac{1}{2})$. In
particular, the Lagrange function itself is

$\ell_{j}(\xi)=w_{j}\sum_{s=0}^{n-1}\Big{(}s+\frac{1}{2}\Big{)}P_{s}(\xi_{j})P_{s}(\xi)\approx
w_{j}\delta(\xi-\xi_{j})$ (40)

making the interpolation polynomial a weighted partial sum of the completeness
relation. As the $P_{j}(\xi)$ are orthogonal, the mass matrix
$M_{ij}=\int_{-1}^{1}\ell_{i}(\xi)\ell_{j}(\xi)d\xi=\omega_{i}\delta_{ij}$ is
diagonal, with inverse $M^{-1}_{ij}=\omega_{i}^{-1}\delta_{ij}$. The advection
and stiffness $S=\mathcal{A}^{T}$ matrices follow from evaluating the
integrals by quadrature, e.g.

$\mathcal{A}_{ij}=\omega_{i}\ell_{j}^{\prime}(\xi_{i})=\omega_{i}\omega_{j}\sum_{s=0}^{n-1}\Big{(}s+\frac{1}{2}\Big{)}P_{s}(\xi_{i})P_{s}^{\prime}(\xi_{j}).$
(41)

Then the weak and strong form internal flux operators are partial sums of the
delta derivative,

$\displaystyle\Upsilon^{i}_{j}$
$\displaystyle=w_{j}\sum_{s=0}^{n-1}\Big{(}s+\frac{1}{2}\Big{)}P_{s}(\xi_{i})P_{s}^{\prime}(\xi_{j})\approx
w_{j}\delta^{\prime}(\xi-\xi_{i})\big{|}_{\xi=\xi_{j}},$ (42) $\displaystyle
D^{i}_{j}$
$\displaystyle=w_{j}\sum_{s=0}^{n-1}\Big{(}s+\frac{1}{2}\Big{)}P_{s}^{\prime}(\xi_{i})P_{s}(\xi_{j})\approx
w_{j}\delta^{\prime}(\xi-\xi_{j})\big{|}_{\xi=\xi_{i}},$ (43)

in the distributional sense. For example, the product $\Upsilon^{i}_{j}F^{j}$
integrates the weak form’s flux term

$\Upsilon^{i}_{j}F^{j}=\sum_{j=0}^{n-1}w_{j}F_{j}\sum_{s=0}^{n-1}\Big{(}s+\frac{1}{2}\Big{)}P_{s}(\xi_{i})P_{s}^{\prime}(\xi_{j})\approx\sum_{j=0}^{n-1}w_{j}\nabla(F(\xi-\xi_{i}))\Big{|}_{\xi=\xi_{j}}\approx\int_{-1}^{1}(\nabla
F^{i})d\xi$ (44)

by quadrature within the element.

###  5.2 DG matrices in the Lobatto basis

Lobatto type quadrature includes interval end-points, so that boundary
information is localized to a single interpolation basis function. The
Legendre-Gauss-Lobatto (LGL) quadrature scheme is expressed by, for
$\int_{-1}^{1}f(\xi)d\xi\approx\sum_{j=0}^{n-1}w_{j}f(\xi_{j})$, the nodes
$\\{\xi_{j}\\}_{j=0}^{n-1}$ and weights $\\{w_{j}\\}_{j=0}^{n-1}$ according to

$\\{\xi_{j}|(1-\xi_{j}^{2})P_{n-1}^{\prime}(\xi_{j})=0\\},\quad\quad
w_{j}=\frac{2}{n(n-1)}\frac{1}{(P_{n-1}(\xi_{j}))^{2}}.$ (45)

The above rule integrates polynomials of degree $\leq 2n-3$ [34].

Discrete orthogonality of classical orthogonal polynomials under Gaussian
quadrature follows from that of the eigenvectors of the defining recurrence
relation represented as a tri-diagonal Jacobi matrix [33]. According to the
article of Gautschi [35], the quadrature nodes and weights were often
constructed via the eigenvalues of the matrix for the Legendre recurrence
relation with a modified final row and column, particularly when using Jacobi
polynomials, rather than the form of Eq. 45. This was originally proposed in a
classic work by Golub [36]. However, a consequence of this modified Jacobi
matrix eigenvalue problem is that the Legendre polynomials remain discretely
orthogonal under LGL quadrature,

$\sum_{s=0}^{n-1}w_{s}P_{k}(\xi_{s})P_{\ell}(\xi_{s})=\gamma_{k}\delta_{k\ell},\quad\quad\gamma_{k}\equiv\begin{cases}\frac{2}{2k+1}&k<n-1,\\\
\frac{2}{n-1}&k=n-1\end{cases}$ (46)

with a modified final eigenvalue from $(k+\frac{1}{2})^{-1}$ to $(k/2)^{-1}$.
This also follows by direction calculation as in [14], where $k<n-1$ follows
by quadrature and $\gamma_{n-1}$ by the boundary weights $w_{0},w_{n-1}$.

Just as for Gauss-Legendre quadrature, the discrete orthogonality relation
expresses the coefficients of the Vandermonde matrix and its inverse. That is,
with the same expansion of Eqn. 39 the matrix elements are
$\mathcal{V}_{j}^{k}\equiv P_{k}(\xi_{j})$ with inverse components
$(\mathcal{V}^{-1})_{j}^{k}=w_{j}\gamma_{k}^{-1}P_{k}(\xi_{j})$ following
directly from Eqn. 46. The spectral transform from the Lagrange to Legendre
spectral basis is identical in form to full Gauss-Legendre quadrature. In
particular the modal form of the Lagrange basis functions themselves is of the
same form,

$\ell_{j}(\xi)=(\mathcal{V}^{-1})_{j}^{k}P_{k}(\xi)=w_{j}\sum_{s=0}^{n-1}\gamma_{k}^{-1}P_{s}(\xi_{j})P_{s}(\xi)\approx\omega_{j}\delta(\xi-\xi_{j})$
(47)

and expresses a re-weighted partial summation of the completeness theorem due
to the Lobatto weighting of the last term in the series, $\gamma_{n-1}$.

As discovered in [14], the Lagrange function spectral form of Eqn. 47 reveals
the Lobatto basis mass matrix to be diagonal with a rank-one update. This
makes it explicitly invertible. To review the quoted result, let
$A_{m}\equiv\int_{-1}^{1}P_{m}(\xi)P_{m}(\xi)d\xi$. Noting that
$A_{m}=\frac{2}{2m+1}=\gamma_{m}$ of Eqn. 46 for $m<n-1$, the mass matrix
integrates to

$\displaystyle M_{ij}=\int_{-1}^{1}\ell_{i}(\xi)\ell_{j}(\xi)d\xi$
$\displaystyle=\sum_{k,m=0}^{n-1}\frac{w_{i}w_{j}}{\gamma_{k}\gamma_{m}}P_{k}(\xi_{i})P_{m}(\xi_{j})\int_{-1}^{1}P_{k}(\xi)P_{m}(\xi)d\xi$
(48)
$\displaystyle=w_{i}\sum_{k=0}^{n-1}\frac{w_{j}}{\gamma_{k}}P_{k}(\xi_{i})P_{k}(\xi_{j})+\frac{A_{n-1}-\gamma_{n-1}}{\gamma_{n-1}^{2}}w_{i}w_{j}P_{n-1}(\xi_{i})P_{n-1}(\xi_{j})$
(49) $\displaystyle=w_{i}(\delta_{ij}-\alpha
w_{j}P_{n-1}(\xi_{i})P_{n-1}(\xi_{j}))$ (50)

as the first sum equals $\ell_{i}(\xi_{j})$ and where
$\alpha=(\gamma_{n-1}-A_{n-1})\gamma_{n-1}^{-2}=\frac{n(n-1)}{2(2n-1)}$. The
mass matrix is full, but differs from the diagonal Gauss-Legendre result by a
rank-one matrix. Applying the identity
$(I-uv^{T})^{-1}=I+\frac{uv^{T}}{1-v^{T}u}$ inverts the mass matrix as

$M^{-1}_{ij}=\frac{1}{w_{i}}\delta_{ij}+\frac{n}{2}P_{n-1}(\xi_{i})P_{n-1}(\xi_{j}).$
(51)

The DG matrices in the Lobatto basis follow from Eqn. 51. Fortunately, the
face mass matrix $\Gamma$, advection matrix $A$ and stiffness matrix $S$
follow more easily than the mass matrix. For the interval $[-1,1]$ the face
mass matrix $\Gamma$ is a Kronecker delta picking out the boundary nodes,
while the advection matrix $\mathcal{A}$ is maximal order $2n-3$, so by
quadrature

$\mathcal{A}_{ij}=\langle\ell_{i}(\xi)|\partial_{\xi}\ell_{j}(\xi)\rangle=w_{i}\ell_{j}^{\prime}(\xi_{i}),$
(52)

with the stiffness matrix its transpose
$S_{ij}=\mathcal{A}_{ij}^{T}=w_{j}\ell_{i}^{\prime}(\xi_{j})$. Using the
projection matrices $M^{-1}$, $\mathcal{A}$, and $S$, and the Lagrange
function of Eqn. 47, the DG operators follow as

$\displaystyle\Upsilon^{j}_{k}$
$\displaystyle=\frac{w_{k}}{w_{j}}\ell_{j}^{\prime}(\xi_{k})+\frac{n}{2}w_{k}P_{n-1}(\xi_{j})P_{n-1}^{\prime}(\xi_{k}),$
(53) $\displaystyle\Xi^{j}_{\gamma^{k}}$
$\displaystyle=\frac{1}{w_{j}}\delta_{j\gamma^{k}}+\frac{n}{2}P_{n-1}(\xi_{j})P_{n-1}(\xi_{\gamma^{k}}),$
(54) $\displaystyle D^{j}_{k}$ $\displaystyle=\ell_{k}^{\prime}(\xi_{j})$ (55)

Here, the numerical flux array $\Xi$ is simply the first and last columns of
the inverse mass matrix where the index
$\gamma^{k}=\delta_{k0}+\delta_{k(n-1)}$ denotes the interval end-points. The
term $P_{n-1}^{\prime}(\xi_{k})$ in $\Upsilon^{j}_{k}$ expressing the rank-one
difference from the Gauss-Legendre result arises from the identity

$\sum_{s=0}^{n-1}P_{n-1}(\xi_{s})\ell_{s}^{\prime}(\xi_{j})=P_{n-1}^{\prime}(\xi_{j})$
(56)

which follows from the orthogonality relation Eqn. 46. The result for
$D^{j}_{k}$ is found by accounting for boundary terms after a discrete
integration by parts and canceling the remainder. By then applying Eqn. 47 the
weak form operator $\Upsilon$ becomes simply

$\Upsilon_{k}^{j}=w_{k}\sum_{s=0}^{n-1}\Big{(}s+\frac{1}{2}\Big{)}P_{s}(\xi_{j})P_{s}^{\prime}(\xi_{k}).$
(57)

This form of the internal flux operator clearly demonstrates it as a flux
discretization, that is, as a partial derivative of the Legendre completeness
relation. In the limit of many quadrature nodes the internal flux operator
approaches $\sim\omega_{k}\delta^{\prime}(\xi-\xi_{j})\big{|}_{\xi=\xi_{k}}$,
picking out the flux gradient as
$\Upsilon^{j}_{k}F^{k}\sim\sum_{k}w_{k}\delta^{\prime}(\xi-\xi_{j})|_{\xi=\xi_{k}}F^{k}\approx\nabla
F^{j}$, just as Gauss-Legendre but of lower accuracy due to the reweighted
value $\gamma_{n-1}$. Similarly, the derivative matrix is seen to be

$D^{j}_{k}=w_{k}\sum_{s=0}^{n-1}\frac{1}{\gamma_{s}}P_{s}^{\prime}(\xi_{j})P_{s}(\xi_{k})$
(58)

so that, through integration by parts to return to the strong form, the
components differ from the completeness relation by the last summation term
[14]. This difference alters only the boundary values of $D^{j}_{k}$ as
$P_{n-1}^{\prime}(\xi_{j})=0$ except for $\xi_{0}$ and $\xi_{n-1}$.

###  5.3 Spectral methods on a discontinuous Lobatto grid: Fourier analysis

This section discusses a complementary approach to the DG method using its
discontinuous interpolation polynomial basis. Problems can be solved by
transforming to a global spectral basis, i.e. a complete set of orthogonal
functions for the entire domain. Since the underlying node set remains the LGL
nodes this approach naturally complements the DG method. The transformation is
summarized by a tensor whose components are called connection coefficients,
because they connect different families of basis functions [23].

Global spectral methods work well for elliptic problems, as illustrated here
by solving Poisson’s equation. This makes the method of interest for coupled
hyperbolic-elliptic systems, with the hyperbolic part solved using DG method
and the elliptic part by the method of this section. Other applications of
this high-order Fourier transform are as a global filter for anti-aliasing,
accurate calculation of convolutions, and post-processing spectral analysis.

Note that in DG method functions are approximated by an $N$-element piecewise
interpolation,

$f(x)=\bigoplus_{m=0}^{N-1}\sum_{j=0}^{n-1}f_{mj}\ell_{j}(x)$ (59)

with the direct sum $\bigoplus$ over elements and $\ell_{j}(x)$ the Lagrange
polynomial of the nodes $\\{\xi_{j}\\}_{j=0}^{n-1}$. Now consider the
determination of Eqn. 59’s Fourier coefficients. As the approximation
satisfies the Dirichlet conditions, it is expandable in Fourier series as [37]

$f(x)=\sum_{p=-\infty}^{\infty}c_{p}e^{ik_{p}x},\quad\quad
c_{p}=\frac{1}{L}\int_{0}^{L}f(x)e^{-ik_{p}x}dx$ (60)

where $L$ is the domain length and $k_{p}=\frac{2\pi}{L}p$ the wavenumbers.
Because the interpolation polynomial is an approximate identity
$\ell_{j}(\xi)\approx w_{j}\delta(\xi-\xi_{j})$, its Fourier transform
approximates the mode of that node, $w_{j}e^{-ik_{p}\xi_{j}}$. This is seen by
carrying out the Fourier integral using the interpolation polynomial expansion
(Eqn. 47) having applied the affine transformation

$\xi=J_{m}(x-\bar{x}_{m})$ (61)

per element with $\bar{x}_{m}$ the mid-point of the $m$’th element,
$J_{m}=2/(\Delta x)_{m}$ the element Jacobian, and $(\Delta x)_{m}$ its width.
The function $\ell_{j}(\xi)$ is a combination of spherical harmonics
$P_{j}(\xi)$, so its Fourier transform is a combination of spherical waves as
the transformation reduces to the integral [23]

$\int_{-1}^{1}P_{s}(\xi)e^{-ik_{p}J_{m}^{-1}\xi}d\xi=2(-i)^{s}j_{s}(J_{m}^{-1}k_{p})$
(62)

with $j_{s}(\zeta)$ the spherical Bessel function. The solution $c_{p}$ is
cleanly expressed by defining

$\displaystyle\mathbb{T}_{p}^{mj}\equiv$ $\displaystyle~{}\frac{(\Delta
x)_{m}}{L}e^{-ik_{p}\bar{x}_{m}}w_{j}\sum_{s=0}^{n-1}\gamma_{s}^{-1}(-i)^{s}P_{s}(\xi_{j})j_{s}(J_{m}^{-1}k_{p})$
(63) $\displaystyle\approx$
$\displaystyle~{}w_{j}e^{-ik_{p}x_{mj}}\frac{2(\Delta x)_{m}}{L}$ (64)

such that $c_{p}=\mathbb{T}_{p}^{mj}f_{mj}$. Equation 63 is an $n$-term sum of
the expansion [17],

$e^{-ikx}=2\sum_{s=0}^{\infty}\Big{(}s+\frac{1}{2}\Big{)}(-i)^{s}P_{s}(x)j_{s}(k).$
(65)

with a Lobatto-weighted coefficient $\gamma_{n-1}$. The transform encoded by
$\mathbb{T}_{p}^{mj}$ depends on the grid only.

####  5.3.1 Approximation capability of the transformation

As the polynomial order of the basis $\ell_{j}(\xi)$ is increased, the
spectrum is approximated by successively higher order quadratures. In practice
the inverse-transformation of Eqn. 60 should be summed up to the highest mode
$p_{\text{max}}$ prior to aliasing. To diagnose the aliasing phenomenon,
consider approximating the elliptic cosine $\text{cn}(x|m)$ normalized to
$x\in[-1/2,1/2]$. Its Fourier expansion is

$\text{cn}(4Kx|m)=\frac{\pi}{K\sqrt{m}}\sum_{n=0}^{\infty}\text{sech}(\pi(1+2n)K^{\prime}/K)\cos(2\pi(1+2n)x)$
(66)

where $K(m)$ is the quarter-period (complete elliptic integral) and
$K^{\prime}(m)=K(1-m)$. Examining the spectral error associated with a
discontinuous interpolation as in Fig. 3 (with $\sqrt{m}=1-10^{-6}$) shows the
number of accurate modes prior to the aliasing limit to go beyond the usual
Nyquist frequency (that is, half the sampling rate) associated with the usual
discrete Fourier transform. At high enough order ($n\gtrsim 6$), the aliasing
point occurs at greater than twice the Nyquist limit of a Fourier series on
equispaced points (or piecewise constant interpolating polynomials).

Figure 3: Normalized coefficient error
$\epsilon=|\widetilde{c}_{p}-c_{p}|/|c_{p}|$ of the elliptic cosine on
$[-1/2,1/2]$ using the transform $\widetilde{c}_{p}=\mathbb{T}_{p}^{mj}f_{mj}$
with $N=20$-elements. Even-numbered modes, not shown, float at $10^{-16}$.
This example explores the location of the aliasing point based on the first
local error minimum. The usual discrete Fourier transform’s aliasing point
(the Nyquist frequency) is at half the sampling frequency, here the tenth mode
$p=10$. The transform by $\mathbb{T}_{p}^{mj}$ for $n>1$ shifts the aliasing
point to greater $p$. Based on these observations, for $n\gtrsim 6$ the
Nyquist mode can be taken as roughly the number of elements, here
$p_{\text{max}}\sim 20$.

Note that if all elements are of equal width, then the indices corresponding
to nodes $j$ and elements $m$ factorize, and the transform can be written much
as a classic Fourier transform matrix on equidistant data. The piecewise-
constant polynomial $n=1$ corresponds with the discrete Fourier transform when
the inverse-transform sum is truncated at the Nyquist mode, as
$j_{0}(k_{p})=\text{sinc}(k_{p})$ factors from the transform components and
the matrix may then be inverted by the discrete orthogonality of Fourier modes
on equidistant points. For $n>1$ the inverse transformation is not exact yet
increases in accuracy with $n$, illustrated by solving Poisson’s equation. In
this sense the scheme constitutes a kind of high-order discrete Fourier
transform.

####  5.3.2 An application: solving Poisson’s equation

For an application of the spectral method, consider the problem

$\frac{d^{2}\Phi}{dx^{2}}=f(x)$ (67)

with periodic boundary conditions $\Phi(x+L)=\Phi(x)$ and $x\in[0,L]$. A
spectral solution is obtained by sampling the source $f(x)$ on the piecewise
quadrature nodes and expanding it in the form of Eqn. 59. Having obtained the
Fourier coefficients $c_{p}=\mathbb{T}_{p}^{mj}f_{mj}$ by application of the
spectral transform tensor (Eqn. 63), the solution in the spectral domain
follows from the transform of Eq. 67,

$-k_{p}^{2}\widetilde{\Phi}(k_{p})=c_{p}$ (68)

for a spectral solution
$\widetilde{\Phi}(k)=\sum_{p=-\infty}^{\infty}-k_{p}^{-2}c_{p}\delta(k-k_{p})$.
Summing modes up to $p_{\text{max}}$, the solution obtained is

$\displaystyle\Phi(x)$
$\displaystyle=\bigoplus_{m=0}^{N-1}\sum_{j=0}^{n-1}\Phi_{mj}\ell_{j}(x)$ (69)
$\displaystyle\Phi_{mj}$
$\displaystyle=\sum_{p=-p_{\text{max}}}^{p_{\text{max}}}-k_{p}^{-2}c_{p}e^{ik_{p}x_{mj}}$
(70)

where $x_{mj}$ is the location of the j’th node in the m’th element. The
electric field $E(x)=-\frac{d\Phi}{dx}$ is obtained in the same manner,

$E_{mj}=\sum_{p=-p_{\text{max}}}^{p_{\text{max}}}-((ik_{p})^{-1}c_{p})e^{ik_{p}x_{mj}}.$
(71)

Often only the field $E(x)$ is needed, like in the Vlasov-Poisson system. The
solution obtained by this high-order spectral method is observed in Fig. 4 to
be $\mathcal{O}((\Delta x)^{n+1/2})$-accurate for a simple benchmark problem
where the error is calculated using the broken $\ell_{2}$-norm,

$\varepsilon=||u-u_{\text{exact}}||_{2}=\frac{1}{N}\sum_{m=1}^{N}||u_{m}-u_{\text{exact}}||_{2,m},\quad||u||_{2,m}\equiv\sqrt{\int_{x_{m}}^{x_{m}+(\Delta
x)_{m}}u(x)^{2}dx}.$ (72)

(a) Error in the potential $\Phi(x)$.

(b) Error in the field $E(x)=-\frac{d\Phi}{dx}$.

Figure 4: Convergence with element width $h$ and polynomial order $p$, or
$hp$-refinement, for the spectral solution of Poisson’s equation
$\Phi^{\prime\prime}=\sin(2\pi x)$ using LGL nodes of $n=3$-$8$ and $N=10$,
$15$, $20$, and $25$ elements, where the legends indicate the line of best
fit’s slope sequentially for $n=3$-$8$. An $h$-convergence of
$\mathcal{O}((\Delta x)^{n+1/2})$ is observed, the same as obtained by solving
the DG strong form with stabilized central numerical fluxes as in [1].

Figure 5 studies the case of a source density $f(x)$ with discontinuities on
element boundaries, a situation admitted by the DG scheme and encountered in
practice when solving coupled hyperbolic-elliptic equations, by introducing
large-amplitude random noise. As the Fourier series is single-valued on
element boundaries, the field and potential solutions via the spectral method
are $\mathcal{C}^{0}$.

(a) Noisy source function $f(x)$.

(b) Potential solution $\Phi(x)$.

(c) Field solution $E(x)=-\frac{d\Phi}{dx}$.

Figure 5: Comparison of solution to Poisson’s equation $\nabla^{2}\Phi=f(x)$
with $x\in[-\pi,\pi]$ between the LGL-Fourier spectral method and a matrix-
based method using stabilized central numerical fluxes as in [1], where the
source density $f(x)=\sin(x)+\mathcal{N}(0,1)$ and $\mathcal{N}(0,1)$ is a
normally distributed random variable of mean $0$ and variance $1$. The domain
consists of $N=10$ elements of $n=8$ nodes each, and Figs. 5(b), 5(c) plot the
solutions using the basis functions $\bigoplus_{m}\sum_{j}y_{mj}\ell_{j}(x)$.

##  6 Example results: simulating the Vlasov-Poisson system

This section explores numerical solutions to the Vlasov-Poisson problem posed
in Section 4 in order to illustrate the DG method on GPU. As in Section 4.1,
the chosen initial condition consists of the velocity distribution
$f(x,u,v)=f_{0}(x,u,v)+\epsilon f_{1}(x,u,v)$ with $\epsilon$ a small
parameter. For the equilibrium distribution $f_{0}$, the loss cone
distribution Eq. 31 with ring parameter $j=6$ and thermal velocity $\alpha=1$
is chosen, while the perturbation mode $f_{1}$ is given by Eq. 35. The
parameter $\epsilon$ is chosen small enough to satisfy the linearization
condition of Eq. 37, though $\epsilon$ is not otherwise not important. For
these examples $\epsilon$ is chosen such that the perturbed density
$\int_{-\infty}^{\infty}fdv=0.002$. Further, the normalized magnetic field is
taken as $B_{\text{ext}}=0.1$.

To summarize the numerical method, Eq. 26 is a conservation law and so
discretized as in Section 2. As the Vlasov equation is a first-order
hyperbolic problem the upwind numerical fluxes of Section 3.1.1 are used,
which are also described in [1] [27]. For time integration, the Shu-Osher
third-order explicit SSP-RK method [38] with spatial-order dependent CFL
numbers given in [39] is used to evolve the semi-discrete equation in the weak
form of Eq. 14. Poisson’s equation, Eq. 27, is solved at each RK stage via the
method of Section 5.3.2, i.e. the DG Fourier-spectral method.

###  6.1 Two cyclotron-harmonic instability case studies

Two example simulations $A$ and $B$ are performed using as perturbation
eigenvalues the pairs $(k_{A},\omega_{A})\sim(0.886,0.349i)$ and
$(k_{B},\omega_{B})\sim(1.4,1.182+0.131i)$, each a solution of Eq. 30, the
dispersion function. Directly exciting a mode is a general method; any
solution $\\{f_{1},(k,\omega)\\}$ of the linearized equations can be simulated
for comparison to linear theory. The modes are chosen with $\omega_{i}>0$ in
order to observe the transition from linear to nonlinear amplitudes with a
verifiable growth rate. Lastly, the domain length is taken as $L=2\pi/k$ and
the domain velocity limits are set to $u_{\text{max}},v_{\text{max}}\sim\pm
8.5$. Figure 6 shows the perturbations $f_{1}(x,u,v)$ at $t=0$.

(a) Mode of simulation $A$ with $\omega_{r}=0$. The mode changes sign in
velocity space where $\partial_{v_{\perp}}f_{0}\sim 0$.

(b) Simulation $B$’s perturbation corresponding to a mode with $\omega_{r}\neq
0$.

Figure 6: The perturbations $f_{1}(x,u,v)$ consist of twisting islands in the
phase space, capturing the combined physics of translation due to particle
momentum, acceleration by the electric field, and rotation by the magnetic
field. Shown here are the perturbations $f_{1}\equiv f-f_{0}$ at time $t=0$
for simulations $A$ and $B$, each with iso-surfaces at $30\%$ of the minimum
(green) and maximum (yellow). The plot axes $(X,Y,Z)$ correspond to the phase
space coordinates $(x,u,v)$.

The problem is run with an element resolution $(N_{x},N_{u},N_{v})=(25,50,50)$
and $n=8$ nodal basis per dimension for a total of $32$ million nodes. These
instabilities grow on a slow time-scale relative to the plasma frequency; that
is, they grow at a fraction of the cyclotron time-scale $\omega_{c}^{-1}$,
while time $t$ is normalized to the plasma frequency $\omega_{p}^{-1}$.
Further, the run conditions have $\omega_{p}=10\omega_{c}$. This means that
the instabilities reach their nonlinear saturation phase after many hundreds
of plasma periods. Simulation $A$ reaches saturation around $t=100$ but was
run to $t=175$, while simulation $B$ saturates at around $t=175$ and was
stopped at $t=200$. Each case requires approximately $50000$ three-stage time
steps to reach the stop time, with a machine run-time of several hours on an
RTX 3090 GPU. For a sense of magnitude, an equivalent single-threaded
implementation on CPU, at least thirty times slower, would require at least
one week of calculation time.

Three-dimensional isosurface plots were produced using PyVista, a Python
package for VTK. To prepare the data, an average was first taken for nodes
lying on element boundaries for smoothness, and the $8$-nodes per element were
resampled to $25$ linearly spaced points per axis and per element on the basis
functions of Section 5. These iso-contours are shown for simulations $A$ and
$B$ in Figs. 7 and 8 respectively. Both cases result in phase space structures
with fine features, a phenomenon in self-consistent kinetic dynamics called
filamentation [40]. These filaments develop shortly into the saturated state,
showing the importance of high-resolution and high-order techniques in
Eulerian simulation of Vlasov-Poisson systems.

Figure 9 shows the electric potentials $\Phi(x)$ in simulations $A$ and $B$.
In simulation $A$ the wave potential $\Phi(x)$ is stationary with a weakly
fluctuating boundary, so that part of the density $f(x,v)$ within the
potential well executes trapped orbits. This results in a trapping structure
with orbits tracing a nonlinear potential similar to that seen in
electrostatic two-stream instability simulations where the potential is like
that of a pendulum, i.e. $\Phi(x)=\sin(x)$, with its characteristic separatrix
structure. In this case particles also execute cyclotron motion so that the
separatrix of the potential in simulation $A$ is similar in form to that shown
by the isosurfaces.

(a) Initial condition, $t=0$.

(b) Developed trapping, $t=80$.

(c) Further filamentation, $t=120$.

Figure 7: Phase space view in $(X,Y,Z)\equiv(x,u,v)$ of simulation $A$ focused
on $(x,u)$-plane as iso-contours at $15\%$ of $\text{max}(f)$ (yellow). Only
the inner domain is shown with $u,v\in(-4,4)$, so the outer edges are visible.
A trapped particle structure develops consisting of a ribbon of density
winding around a a separatrix, while the outer ring consists of untrapped
(passing) particles.

On the other hand, the saturated wave potential $\Phi(x)$ of simulation $B$
translates with positive phase velocity $v_{\varphi}\approx\omega_{r}/k$. The
region of particle interaction translates along with the wave potential and
forms a vortex structure with the appearance of a kink in the phase space
density $f(x,v)$. The center of this kink continues to tighten as the
simulation progresses, leading to progressively finer structures just as in
simulation $A$. This effect is in agreement with the filamentation phenomenon
and introduces simulation error as the structures lead to large gradients on
the grid scale where discreteness produces dispersion error. For this reason
the simulation is stopped at $t=200$.

(a) Developing mode, $t=100$.

(b) Developed vortex, $t=160$.

(c) Translating vortex, $t=180$.

Figure 8: Phase space view looking on $(-x,v)$ \- plane of simulation $B$ at
$15\%$ iso-contours of $\text{max}(f)$ (yellow). The mode is seen to be a
growing, translating potential $\Phi(x)$ of positive phase velocity
$v_{\varphi}=\omega_{r}/k$ with an underlying phase space vortex structure
centered at $(u,v)=0$. The vortex shape is explained by considering the
trajectory of a test particle in the wave. That is, particles with a velocity
close to that of the wave see a stationary potential and are accelerated to a
high $u$-velocity. They then translate towards positive $x$ while their
velocity vector is rotated by the Lorentz force to $-u$ at a rate close to the
wave frequency (as $\omega_{r}\approx 1.2\omega_{c}$). The particle then
advects opposite the direction of wave propagation, before repeating the
cycle.

(a) Simulation A at $t=120$.

(b) Simulation B at $t=180$

Figure 9: Electric potentials $\Phi(x)$ at saturation of the two studied
cases. The potential of $A$ is stationary while that of $B$ is translating to
the right. The negative of the potential $-\Phi(x)$ is shown in order to
account for the electron’s negative charge. In both cases particle trapping
structures develop in the potential wells, or minimum regions, in $-\Phi(x)$.

(a) Simulation A, with $\omega_{A}=0.349i$.

(b) Simulation B, with $\omega_{B}=1.18+0.131i$

Figure 10: Domain-integrated field energy growth profiles for simulations $A$
and $B$. The thermal energy of the zero-order distribution is $0.25$ per unit
length, according to Appendix B with $\alpha=1$. In simulation $A$ this
corresponds to a domain-integrated thermal energy of $E_{A}\approx 17.8$ and
in simulation $B$ to $E_{B}\approx 11.2$. In both cases the instability
saturates with an electric energy a few percent of the thermal energy,
approximately $6\%$ in $A$ and $2.5\%$ in $B$.

##  7 Conclusion

This study is a self-contained description of how to implement the
discontinuous Galerkin method on the GPU, i.e. a technique to solve PDEs using
many finite elements, where the approximation in each element is supported by
interpolation polynomials. When these basis polynomials interpolate the Gauss-
Lobatto quadrature points of the one-dimensional interval
$\mathcal{I}\equiv[-1,1]$ then explicit and interpretable forms of the
discontinuous Galerkin operators can be found. This is the case due to the
affinity of quadrature on $\mathcal{I}$ to the spherical harmonics through the
theory of Legendre polynomials. Recognizing this affinity allows the
calculation of connection coefficients to alternative basis sets; for example
in this work the Fourier spectrum of the interpolation polynomials is utilized
to solve the Poisson equation to $\mathcal{O}((\Delta x)^{n+1/2})$ accuracy.

An important finding is that orthogonal discretization of a computational
domain results in significant computational savings over simplex elements, and
increasingly so per dimension. That is, only one-dimensional basis matrices
need be constructed for general orthogonal discretizations using $n$-cube
elements built up from tensor products of $\mathcal{I}$, while full tensors
must be constructed for simplex elements. Finally, these results are combined
to solve an example problem, namely the Vlasov-Poisson system of coupled
hyperbolic-elliptic PDEs, using GPU-accelerated libraries. The semi-discrete
equation is evaluated in only a few lines of code by utilizing CUDA wrappers
with NumPy-like data arrays, allowing tensor-product index ordering in a
simple routine. This approach does not outperform a custom implementation with
CUDA code [9], yet it has the advantage for beginners of simplicity.

## References

  * [1] J.S. Hesthaven and T. Warburton.  Nodal discontinuous Galerkin methods.  Springer, 2008.
  * [2] B. Cockburn, G.E. Karniadakis, and C.W. Shu.  Discontinuous Galerkin methods: theory, computation, and applications.  Springer Science and Business Media, 2012.
  * [3] V. Dolejší and M. Feistauer.  Discontinuous Galerkin Method: Analysis and Applications to Compressible Flow.  Springer, 2015.
  * [4] C. Johnson and J. Pitkäranta.  An analysis of the discontinuous Galerkin method for a scalar hyperbolic equation.  Math. Comp., 46:1–26, 1986.
  * [5] S. Rhebergen and G.N. Wells.  A Hybridizable Discontinuous Galerkin Method for the Navier–Stokes Equations with Pointwise Divergence-Free Velocity Field.  J. Sci. Comp., 76:1484–1501, 2018.
  * [6] L. Dongmi, W. Huang, and J. Qui.  A quasi-Lagrangian moving mesh discontinuous Galerkin method for hyperbolic conservation laws.  J. Comp. Phys., 396:544–578, 2019.
  * [7] P. Castillo and S. Gòmez.  Conservative super-convergent and hybrid discontinuous Galerkin methods applied to nonlinear Schrödinger equations.  Appl. Math. Comp., 371, 2020.
  * [8] M. Tavelli and M. Dumbser.  A pressure-based semi-implicit space–time discontinuous Galerkin method on staggered unstructured meshes for the solution of the compressible Navier–Stokes equations at all Mach numbers.  J. Comp. Phys., 341:341–376, 2017.
  * [9] L. Einkemmer.  Semi-Lagrangian Vlasov Simulation on GPUs.  Comput. Phys. Comm., 254, 2020.
  * [10] K. Świrydowicz, N. Chalmers, A. Karakus, and T. Warburton.  Acceleration of tensor-product operations for high-order finite element methods.  Int. J. High Perform. Comput. Appl., 33(4), 2019.
  * [11] P. Fischer, K. Heisey, and M. Min.  Scaling limits for PDE-based simulation.  Proceedings of 22nd AIAA computational fluid dynamics conference, Dallas, United States, 22 June 2015.
  * [12] P. Fischer, M. Min, T. Rathnayake, et al.  Scalability of high-performance PDE solvers.  Int. J. High Perform. Comput. Appl., 34(5), 2020.
  * [13] G.M. Karniadakis and S.J. Sherwin.  Spectral/hp Element Methods for CFD.  Oxford University Press, 2005.
  * [14] S.A. Teukolsky.  Short note on the mass matrix for Gauss-Lobatto grid points.  J. Comp. Phys., 312(333), 2016.
  * [15] S.A. Teukolsky.  Formulation of discontinuous Galerkin methods for relativistic astrophysics.  J. Comp. Phys., 312, 2016.
  * [16] G. Kuperberg.  Numerical Cubature from Archimedes’ Hat-box Theorem.  SIAM J. Numer. Anal., 44, 2006.
  * [17] S. Hassani.  Mathematical Physics: A Modern Introduction to its Foundations.  Springer, Second edition, 2013.
  * [18] I. Stakgold.  Boundary Value Problems of Mathematical Physics, volume I.  The Macmillan Company, 1967.
  * [19] J.A. Rossmanith and D.C. Seal.  A positivity-preserving high-order semi-lagrangian discontinuous Galerkin scheme for the Vlasov–Poisson equations.  J. Comp. Phys., 230:6203–6232, 2011.
  * [20] G.V. Vogman, P. Colella, and U. Shumlak.  Dory-Guest-Harris instability as a benchmark for continuum kinetic Vlasov-Poisson simulations of magnetized plasmas.  J. Comp. Phys., 277, 2014.
  * [21] J. Juno, A. Hakim, J. TenBarge, E. Shi, and W. Dorland.  Discontinuous Galerkin algorithms for fully kinetic plasmas.  J. Comp. Phys., 353, 2018.
  * [22] A. Hakim, M. Francisquez, J. Juno, and G. Hammett.  Conservative discontinuous Galerkin schemes for nonlinear Dougherty-Fokker-Planck collision operators.  J. Plasma Phys., 86(4), 2020.
  * [23] S. Olver, R.M. Slevinsky, and A. Townsend.  Fast algorithms using orthogonal polynomials.  Acta Numerica, 29:573–699, 2020.
  * [24] A. Ho, I.A.M. Datta, and U. Shumlak.  Physics-based-adaptive plasma model for high-fidelity numerical simulations.  Front. Phys., 2018.
  * [25] R. Okuta, Y. Unno, D. Nishino, S. Hido, and C. Loomis.  CuPy: A NumPy-Compatible Library for NVIDIA GPU calculations.  In Proceedings of Workshop on Machine Learning Systems (LearningSys) in The Thirty-first Annual Conference on Neural Information Processing Systems (NIPS), 2017.
  * [26] Nodes and Weights of Gauss-Lobatto Calculator.  https://kesian.casio.com/, May 2021.
  * [27] R.J. Leveque.  Finite volume methods for hyperbolic problems.  Cambridge Texts in Applied Mathematics, 2002.
  * [28] B. Cockburn and C.W. Shu.  The Local Discontinuous Galerkin Method for Time-Dependent Convection-Diffusion Systems.  SIAM J. Numer. Anal., 35:2440–2463.
  * [29] T.H. Stix.  Waves in Plasmas.  Springer-Verlag, New York, 1992.
  * [30] R. Balescu.  Aspects of anomalous transport in plasmas.  Institute of Physics, Bristol, 2005.
  * [31] D.A. Gurnett and A. Bhattacharjee.  Introduction to plasma physics.  Cambridge University Press, Second edition, 2017.
  * [32] J.A. Tataronis and F.W. Crawford.  Cyclotron harmonic wave propagation and instabilities: I. Perpendicular propagation.  J. Plasma Phys., 4, 1970.
  * [33] C. Dunkl and Y. Xu.  Orthogonal polynomials of several variables.  Cambridge University Press, Second edition, 2014.
  * [34] R. Bulirsch and J. Stoer.  Introduction to numerical analysis.  Springer, New York, 1991.
  * [35] W. Gautschi.  High-order Gauss-Lobatto formulae.  Numerical Algorithms, 25:213–222, 2000.
  * [36] G. Golub.  Some modified matrix eigenvalue problems.  SIAM Rev., 15:318–334, 1973.
  * [37] C. Lanczos.  Discourse on Fourier Series.  SIAM Classics in Applied Mathematics, 2016.
  * [38] S. Gottlieb.  On High Order Strong Stability Preserving Runge-Kutta and Multi Step Time Discretizations.  J. Sci. Comp., 25(112), 2005.
  * [39] B. Cockburn and C.W. Shu.  Runge-Kutta Discontinuous Galerkin Methods for Convection-Dominated Problems.  J. Sci. Comp., 16(3), 2001.
  * [40] C.Z. Cheng and G. Knorr.  The integration of the Vlasov equation in configuration space.  J. Comp. Phys., 22:330–351, 1976.

##  Appendix A Implemention of upwind fluxes in Python with broadcasting

This appendix details how the upwind numerical flux may be prepared for the
basis product function described in Section 3.1. The function cupy.roll(a,
shift) (viz. numpy.roll) returns a view into a shifted axis of the tensor
array, and is ideal to compute the numerical flux using GPU-accelerated
element-wise operations. In each direction $j$ the numerical flux array
$\mathcal{F}^{\alpha,\beta,j}$ is to be constructed of shape
[$N_{0},n_{0},N_{1},n_{1},\cdots,N_{j},2,\cdots,N_{d-1},n_{d-1}$] for $N_{j}$
the elements in $j$ and $n_{d}$ its sub-element nodes. The advection speed
array, of size [$N_{j},n_{j}$], has its sign measured by the two arrays

        one_negatives = cp.where(condition=speed < 0, x=1, y=0)
        one_positives = cp.where(condition=speed > 0, x=1, y=0)

which calculates $a^{R}$ and $a^{L}$. The numerical flux array for each
direction can then be built with broadcasting. To do this in a dimension-
independent fashion, define boundary_slices for element boundaries and
advection_slices for the flux vector components.

In two dimensions, for example, boundary_slices is defined as the list of
tuples

        e0, e1 = slice(elements[0]), slice(elements[1])
        n0, n1 = slice(nodes[0]), slice(nodes[1])
        self.boundary_slices = [
                # x-directed face slices [(left), (right)]
                [(e0, 0, e1, n1), (e0, -1, e1, n1)],
                # y-directed face slices [(left), (right)]
                [(e0, n0, e1, 0), (e0, n0, e1, -1)] ]

On the other hand, advection_slices is an array to prepare broadcasting with
one of two directions of the e.g. shape [elements[0], nodes[0], elements[1],
2] numerical flux array $\mathcal{F}^{\alpha,\beta,j}$ for direction 1. For
example, if this 2D flux were for a rotation
$F=\begin{bmatrix}-y,x\end{bmatrix}^{T}$, the advection slices are

        self.advection_slices = [(None, slice(elements[1]), slice(nodes[1])),
                                 (slice(elements[0]), slice(nodes[0]), None)]

Having set up this infrastructure, the numerical flux in direction dim is
constructed by the function:

        def flux_in_dim(self, flux, one_negatives, one_positives, basis_matrix, dim):
            # allocate
            num_flux = cp.zeros(self.num_flux_sizes[dim])

            # Upwind flux, left face (with surface normal vector -1.0)
            num_flux[self.boundary_slices[dim][0]] =
                -1.0 * (cp.multiply(cp.roll(flux[self.boundary_slices[dim][1]],
                                            shift=1, axis=self.grid_axis[dim]),
                                    one_positives[self.advection_slices[dim]]) +
                        cp.multiply(flux[self.boundary_slices[dim][0]],
                                    one_negatives[self.advection_slices[dim]]))

            # Upwind flux, right face (with surface normal vector +1.0)
            num_flux[self.boundary_slices[dim][1]] =
                (cp.multiply(flux[self.boundary_slices[dim][1]],
                             one_positives[self.advection_slices[dim]]) +
                 cp.multiply(cp.roll(flux[self.boundary_slices[dim][0]], shift=-1,
                                     axis=self.grid_axis[dim]),
                             one_negatives[self.advection_slices[dim]]))

            return basis_product(flux=num_flux, basis_matrix=basis.xi,
                                 axis=self.sub_element_axis[dim],
                                 permutation=self.permutations[dim])

In the case of a hyperbolic problem, if the sign of the advection speed is
constant during a problem then the sign arrays should be computed prior to the
main loop.

##  Appendix B Thermal properties of the ring distribution

This appendix discusses the first and second moments of the ring distribution

$f_{0}(v)=\frac{1}{2\pi\alpha^{2}j!}\Big{(}\frac{v^{2}}{2\alpha^{2}}\Big{)}^{j}\exp\Big{(}-\frac{v^{2}}{2\alpha^{2}}\Big{)}$
(73)

where $v$ is the radial part in polar coordinates. The first radial moment
$\langle v\rangle$ measures the average velocity of the particles in the ring.
The moment is given by

$\displaystyle\langle v\rangle$
$\displaystyle\equiv\int_{0}^{\infty}vf_{0}(v)(2\pi v)dv$ (74)
$\displaystyle=\frac{2}{j!}\int_{0}^{\infty}\Big{(}\frac{v^{2}}{2\alpha^{2}}\Big{)}^{j+1}\exp\Big{(}-\frac{v^{2}}{2\alpha^{2}}\Big{)}dv$
(75) $\displaystyle=\sqrt{2}\alpha\frac{\Gamma(j+\frac{3}{2})}{\Gamma(j+1)}$
(76)

by using Euler’s integral, where $\Gamma(j+1)=j!$ is the Gamma or factorial
function. Expanding the Gamma ratio in Laurent series about $j=\infty$ yields

$\frac{\Gamma(j+\frac{3}{2})}{\Gamma(j+1)}=\sqrt{j}\Big{(}1+\frac{3}{8j}+\mathcal{O}(j^{-2})\Big{)}\approx\sqrt{j}$
(77)

so that $\langle v\rangle\approx\sqrt{2j}\alpha$, which corresponds with the
peak of the distribution. Likewise, the centered second moment
$\frac{1}{2}\langle(v-\langle v\rangle)^{2}\rangle$ measures the average
energy of the particles. It is given by

$\displaystyle\frac{1}{2}\langle(v-\langle v\rangle)^{2}\rangle$
$\displaystyle=\frac{1}{2}\int_{0}^{\infty}(v-\langle
v\rangle)^{2}f_{0}(v)(2\pi v)dv,$ (78)
$\displaystyle=\frac{1}{2}\Big{(}\Big{[}\int_{0}^{\infty}v^{2}f_{0}(v)(2\pi
v)dv\Big{]}-2\langle v\rangle^{2}+\langle v\rangle^{2}\Big{)},$ (79)
$\displaystyle=\alpha^{2}(j+1-\Big{(}\frac{\Gamma(j+\frac{3}{2})}{\Gamma(j+1)}\Big{)}^{2}).$
(80)

Expanding this result about $j=\infty$ in Laurent series, one has

$j+1-\Big{(}\frac{\Gamma(j+\frac{3}{2})}{\Gamma(j+1)}\Big{)}^{2}=\frac{1}{4}\Big{(}1-\frac{1}{8j}+\mathcal{O}(j^{-2})\Big{)}\approx\frac{1}{4}.$
(81)

The average velocity and thermal energy of the ring distribution can be taken
to a fair degree of accuracy as (with leading-order accuracy of $3/(8j)$ and
$1/(8j)$ respectively),

$\displaystyle\langle v\rangle$ $\displaystyle\approx\sqrt{2j}\alpha,$ (82)
$\displaystyle\frac{1}{2}\langle(v-\langle v\rangle)^{2}\rangle$
$\displaystyle\approx\Big{(}\frac{\alpha}{2}\Big{)}^{2}.$ (83)