text:
\subsection{High-order absorbing boundary conditions for planar boundary}
\label{sec:form:HABC}

Let us consider the half-space problem defined on the domain $\Omega = \{\vec{x}\in\mathbb{R}^3: x<0\}$ with the planar boundary $\Gamma = \{\vec{x}\in\mathbb{R}^3: x=0\}$, where $x$ is the coordinate in the Cartesian direction $\vec{e}_x$.
We seek a non-reflective boundary condition to prescribe on $\Gamma$.
For convenience, the transverse component of the position is denoted $\vec{y}$, such that $\vec{x}=(x,\vec{y})$.

The exact non-reflective boundary condition of the half-space problem is well-known.
Using notations borrowed from the pseudo-differential theory, it reads (see \textit{e.g.}~ \cite{engquist1977absorbing,hagstrom1999radiation})
\begin{equation}
  \mathcal{B}^x p = 0, \label{eqn:exactBC}
\end{equation}
with the pseudo-differential operator
\begin{align}
  \mathcal{B}^x &\equiv c\:\partial_x + \partial_t \sqrt{1 - c^2\Delta_\perp^x/\partial_{tt}}, \label{eqn:exactBCoperator}
\end{align}
where $\Delta_\perp^x \equiv \Delta - \partial_{xx}$ is the Laplace-Beltrami operator defined in the plan tangent to the direction $\vec{e}_x$.
Unfortunately, this condition is non-local in both time and space because of the square root, which makes it an impractical boundary treatment.

Local absorbing boundary conditions are classically obtained by approximating the square root $\sqrt{1+X}$ to localize the operator \eqref{eqn:exactBCoperator}.
The features of the obtained conditions depend on the approximation that is used for the square root.
In their seminal work, Engquist and Majda \cite{engquist1977absorbing} showed that Pad\'e approximations lead to stable conditions, while polynomial approximations based on Taylor expansions can lead to unstable conditions.
Other rational approximations have been used to derive one-way wave equations and absorbing boundary conditions with better accuracy for grazing waves, evanescent modes or long-duration simulations (see \textit{e.g.}~ \cite{antoine2006improved,asvadurov2003optimal,hagstrom2009complete,hagstrom2010radiation,halpern1988wide,higdon1986absorbing,ingerman2000optimal,kechroud2005numerical,lu1998complex,milinazzo1997rational}).
In this work, we restrict ourselves to the Pad\'e approximation, which corresponds to an asymptotic case for these rational approximations.
The boundary treatment and the computational procedure described hereafter will be extended to other approximations in the future.

The $(2N+1)^\text{th}$-order Pad\'e approximation of the square root $\sqrt{1+X}$ is classically written as the rational function \cite{bamberger1988higher,guan1985high,lu1998complex}
\begin{align}
  f_N(X) &= 1 + \frac{2}{M} \sum_{n=1}^N \frac{a_n X}{1 + b_n X}, \label{eqn:approxRoot1}
\end{align}
which we rewrite as
\begin{align}
  f_N(X) &= 1 + \frac{2}{M} \sum_{n=1}^N c_n \left(1 - \frac{1+c_n}{1+c_n + X}\right), \label{eqn:approxRoot2}
\end{align}
where $a_n = \sin^2(n\pi/M)$, $b_n = \cos^2(n\pi/M)$, $c_n = \tan^2(n\pi/M)$ and $M=2N+1$.
Using the formula \eqref{eqn:approxRoot2} to approximate the square root in the exact boundary operator \eqref{eqn:exactBCoperator} gives to the approximate boundary condition
\begin{equation}
  \mathcal{L}^x p = 0, \label{eqn:approxBC}
\end{equation}
with the pseudo-differential operator
\begin{align*}
  \mathcal{L}^x &\equiv \partial_t + c \: \partial_x + \frac{2}{M} \sum_{n=1}^{N} c_n \: \partial_t \left(1 - \frac{(1 + c_n) \partial_{tt}}{(1 + c_n) \partial_{tt} - c^2 \Delta_\perp^x}\right).
\end{align*}
We introduce $N$ auxiliary fields $p_n$ defined on the boundary $\Gamma$ as
\begin{align*}
  p_n \equiv \mathcal{M}_n^x\:p, \quad\quad \text{for } n=1,\dots,N,
\end{align*}
with
\begin{align}
  \mathcal{M}_n^x &\equiv \frac{(1 + c_n) \partial_{tt}}{(1 + c_n) \partial_{tt} - c^2 \Delta_\perp^x}, \quad\quad \text{for } n=1,\dots,N.
  \label{eqn:operatorM}
\end{align}
We can then explicitly rewrite the boundary condition \eqref{eqn:approxBC} as
\begin{align}
  \partial_t p + c \: \partial_x p &= \frac{2}{M} \sum_{n=1}^{N} c_n \partial_t (p_n-p), \label{eqn:HABCplan1}
\end{align}
where the auxiliary fields are governed by
\begin{align}
  \left(1+c_n\right) \partial_{tt} \left(p_n - p\right) - c^2 \Delta_\perp^x p_n &= 0, \quad\quad \text{for } n=1,\dots,N. \label{eqn:HABCplan2}
\end{align}
The boundary condition is local and requires the computation of the $N$ auxiliary equations on the boundary $\Gamma$.
Increasing the order of the Pad\'e approximation increases the number of auxiliary equations and the computational cost, but it also improves the accuracy for outgoing traveling waves \cite{collino1993conditions,collino1993high,engquist1977absorbing}.

We note that rational approximations of the square root can be written in several ways, leading to different representations of the boundary conditions.
Collino \cite{collino1993conditions,collino1993high} used the rational representation \eqref{eqn:approxRoot1} and obtained HABCs with auxiliary equations very close to equations \eqref{eqn:HABCplan1}-\eqref{eqn:HABCplan2}.
This representation can also be used with different parameters as soon as they verify some relations \cite{ha1994stability}.
Alternatively, rational functions written as continued fractions lead to HABCs with auxiliary fields governed by coupled equations (see \textit{e.g.}~ \cite{guddati2000continued,hagstrom2004new,asvadurov2003optimal,givoli2003high}).
In this work, we choose the representation \eqref{eqn:approxRoot2} because it leads to compatibility conditions that are easier to incorporate in a discontinuous Galerkin framework than those obtained in previous works.
This aspect is discussed later in the text.

\subsection{Compatibility conditions at edges and corners}
\label{sec:form:corners}

We now extend the boundary treatment to the borders of a cuboidal domain $\Omega=[-L^x,L^x]\times[-L^y,L^y]\times[-L^z,L^z]$ to simulate the infinite space $\mathbb{R}^3$.
The initial conditions and any source are assumed to be compactly supported inside the domain $\Omega$.
Under this assumption, the exact boundary condition \eqref{eqn:exactBC} and its approximation \eqref{eqn:approxBC} can be prescribed on the planes containing each of the six faces of the domain (figure \ref{fig:illustPlanes:1}).
In practice, we would like to prescribe the boundary condition \eqref{eqn:HABCplan1} only on the faces (figure \ref{fig:illustPlanes:2}).
However, computing the auxiliary equations \eqref{eqn:HABCplan2} only on the faces requires boundary conditions for the auxiliary fields on the borders of each face, which are on the edges of the domain.

\begin{figure}[!t]
\centering
\begin{tabular}{c@{\:}c}
\begin{subfigure}[b]{6.5cm}
  \centering
  \caption{HABC defined on planes}
  \includegraphics[width=65mm]{img/cubePlanes1.pdf}
  \label{fig:illustPlanes:1}
\end{subfigure}
&
\begin{subfigure}[b]{6.5cm}
  \centering
  \caption{HABC defined on faces}
  \includegraphics[width=65mm]{img/cubePlanes2.pdf}
  \label{fig:illustPlanes:2}
\end{subfigure}
\end{tabular}
\vspace{-0.5cm}
\caption{Illustration of HABCs defined on the planes containing the faces of the cuboidal domain (a) or defined only on the faces (b). When defining HABCs only on the faces, supplementary boundary conditions must be prescribed for the auxiliary variables on the borders of each face.}
\label{fig:illustPlanes}
\end{figure}

To derive such conditions, we have adapted a strategy proposed by Hagstrom and Warburton \cite{hagstrom2004new} and based on purely algebraic manipulations of the equations.
It gives edge compatibility conditions which preserve the accuracy of the solution, and require the computation of supplementary auxiliary fields governed by auxiliary equations on the edges.
Proceeding in a hierarchical fashion, these new equations require boundary conditions on the borders of each edge, which are at the corners of the domain.
Algebraic manipulations of these equations similarly provide compatibility conditions with supplementary auxiliary fields at the corners.

In this section, we derive compatibility conditions for the edges belonging to the lines $(x,y)=(L^x,L^y)$, $(x,z)=(L^x,L^z)$ and $(y,z)=(L^y,L^z)$ and for the corner $(x,y,z)=(L^x,L^y,L^z)$.
The conditions for the other edges and corners can be obtained straightforwardly by symmetry.

\subsubsection*{Auxiliary fields}

For the sake of clarity, we first define the fields as
\begin{align}
  p_{ijk} &\equiv \mathcal{M}_i^x \mathcal{M}_j^y  \mathcal{M}_k^z \: p,
  \quad\quad \text{for } i,j,k=0,\dots,N,
  \label{eqn:defAddFields}
\end{align}
where the operators $\mathcal{M}_i^x$, $\mathcal{M}_j^y$, $\mathcal{M}_k^z$ are defined using equation \eqref{eqn:operatorM} for $i,j,k>0$ and $\mathcal{M}_0^x=\mathcal{M}_0^y=\mathcal{M}_0^z=1$.
The fields with zero, one, two and three non-zero subscripts are computed on the volume, the faces, the edges and the corners, respectively.
The field $p_{000}$ corresponds to the pressure field $p$.
The fields $p_{i00}$, $p_{0j0}$, $p_{00k}$, with $i,j,k=1,\dots,N$, are defined on the faces $x=L^x$, $y=L^y$, $z=L^z$, respectively.
The fields $p_{ij0}$, $p_{i0k}$, $p_{0jk}$, with $i,j,k=1,\dots,N$, are defined on the edges $(x,y)=(L^x,L^y)$, $(x,z)=(L^x,L^z)$, $(y,z)=(L^y,L^z)$, respectively.
The fields $p_{ijk}$, with $i,j,k=1,\dots,N$, are defined on the corner $(x,y,z)=(L^x,L^y,L^z)$.
Therefore, there are $N$, $N^2$ and $N^3$ auxiliary fields per face, edge and corner, respectively.
Note that, since the initial conditions on $p$ are assumed to be compactly supported inside the domain, the initial conditions on the auxiliary fields are equal to $0$.

The auxiliary fields have two properties that are used when deriving the compatibility conditions.
First, all the auxiliary fields satisfy the wave equation,
\begin{equation}
  \partial_{tt} p_{ijk} - c^2 \Delta p_{ijk} = 0,
  \quad\quad \text{for } i,j,k=0,\dots,N.
  \label{eqn:waveForPijk}
\end{equation}
Indeed, since the pressure field is in the kernel of the wave operator $\left[\partial_{tt}-c^2\Delta\right]$, the auxiliary fields defined by equation \eqref{eqn:defAddFields} are also in this kernel.
Then, the auxiliary fields are related by the following relations
\begin{subequations}
\begin{align}
  p_{ijk} &= \mathcal{M}_i^x \: p_{0jk}, \quad\quad \text{for } i,j,k=0,\dots,N, \label{eqn:relPijk_x} \\
  p_{ijk} &= \mathcal{M}_j^y \: p_{i0k}, \quad\quad \text{for } i,j,k=0,\dots,N, \label{eqn:relPijk_y} \\
  p_{ijk} &= \mathcal{M}_k^z \: p_{ij0}, \quad\quad \text{for } i,j,k=0,\dots,N. \label{eqn:relPijk_z}
\end{align}
\end{subequations}
These relations are trivially obtained by using definition \eqref{eqn:defAddFields} and observing that the operators $\mathcal{M}_i^x$, $\mathcal{M}_j^y$, $\mathcal{M}_k^z$ commute.

\subsubsection*{Boundary conditions on faces and 2D relations}

We write the boundary conditions for the 3D field on the faces $x=L^x$, $y=L^y$ and $z=L^z$.
On these faces, the field $p_{000}$ satisfies HABCs corresponding to half spaces,
\begin{align*}
  \mathcal{L}^x p_{000} &= 0, \quad\quad \text{for } x=L^x, \\
  \mathcal{L}^y p_{000} &= 0, \quad\quad \text{for } y=L^y, \\
  \mathcal{L}^z p_{000} &= 0, \quad\quad \text{for } z=L^z,
\end{align*}
which can be rewritten
\begin{subequations}
\begin{align}
  \partial_t p_{000} + c \: \partial_x p_{000}
    &= \frac{2}{M} \sum_{i=1}^{N} c_i \partial_t (p_{i00}-p_{000}), \quad\quad \text{for } x=L^x, \label{eqn:bc1:x} \\
  \partial_t p_{000} + c \: \partial_y p_{000}
    &= \frac{2}{M} \sum_{j=1}^{N} c_j \partial_t (p_{0j0}-p_{000}), \quad\quad \text{for } y=L^y, \label{eqn:bc1:y} \\
  \partial_t p_{000} + c \: \partial_z p_{000}
    &= \frac{2}{M} \sum_{k=1}^{N} c_k \partial_t (p_{00k}-p_{000}), \quad\quad \text{for } z=L^z. \label{eqn:bc1:z}
\end{align}
\end{subequations}
Using the definition \eqref{eqn:defAddFields}, we have the 2D relations
\begin{subequations}
\begin{align}
  (1+c_i) \partial_{tt} (p_{i00} - p_{000}) - c^2 (\partial_{yy} + \partial_{zz}) p_{i00}
    &= 0, \quad\quad \text{for } x=L^x, \label{eqn:lev1:iX} \\
  (1+c_j) \partial_{tt} (p_{0j0} - p_{000}) - c^2 (\partial_{xx} + \partial_{zz}) p_{0j0}
    &= 0, \quad\quad \text{for } y=L^y, \label{eqn:lev1:jY} \\
  (1+c_k) \partial_{tt} (p_{00k} - p_{000}) - c^2 (\partial_{xx} + \partial_{yy}) p_{00k}
    &= 0, \quad\quad \text{for } z=L^z, \label{eqn:lev1:kZ}
\end{align}
\end{subequations}
for $i,j,k = 1,\dots,N$.

\subsubsection*{Boundary conditions on edges and 1D relations}

On the borders of each face (\textit{i.e.}~ on the edges of the domain), the 2D fields need boundary conditions.
We derive here the compatibility conditions for the 2D fields on the edges $(x,y)=(L^x,L^y)$, $(x,z)=(L^x,L^z)$ and $(y,z)=(L^y,L^z)$.

Because $p_{000}$ belongs to the kernel of $\mathcal{L}^x$ and $\mathcal{L}^y$ on the edge $(x,y)=(L^x,L^y)$ and because these operators commute with $\mathcal{M}^x$ and $\mathcal{M}^y$, the 2D fields $p_{0j0}$ and $p_{i00}$ also belong to the kernel of $\mathcal{L}^x$ and $\mathcal{L}^y$ seen their definition.
This result gives boundary conditions for the 2D fields on the edge $(x,y)=(L^x,L^y)$.
Using similar reasoning, we obtain boundary conditions on the other edges.
We then have
\begin{align*}
  \mathcal{L}^x p_{0j0} &= 0, \quad\quad \text{for } (x,y)=(L^x,L^y), \\
  \mathcal{L}^y p_{i00} &= 0, \quad\quad \text{for } (x,y)=(L^x,L^y), \\
  \mathcal{L}^x p_{00k} &= 0, \quad\quad \text{for } (x,z)=(L^x,L^z), \\
  \mathcal{L}^z p_{i00} &= 0, \quad\quad \text{for } (x,z)=(L^x,L^z), \\
  \mathcal{L}^y p_{00k} &= 0, \quad\quad \text{for } (y,x)=(L^y,L^z), \\
  \mathcal{L}^z p_{0j0} &= 0, \quad\quad \text{for } (y,x)=(L^y,L^z),
\end{align*}
which can be rewritten
\begin{subequations}
\begin{align}
  \partial_t p_{0j0} + c\: \partial_x p_{0j0} &= \frac{2}{M} \sum_{i=1}^{N} c_i \partial_t (p_{ij0}-p_{0j0}),
    \quad\quad \text{for } (x,y)=(L^x,L^y), \label{eqn:bc2:1} \\
  \partial_t p_{i00} + c\: \partial_y p_{i00} &= \frac{2}{M} \sum_{j=1}^{N} c_j \partial_t (p_{ij0}-p_{i00}),
    \quad\quad \text{for } (x,y)=(L^x,L^y), \label{eqn:bc2:2} \\
  \partial_t p_{00k} + c\: \partial_x p_{00k} &= \frac{2}{M} \sum_{i=1}^{N} c_i \partial_t (p_{i0k}-p_{00k}),
    \quad\quad \text{for } (x,z)=(L^x,L^z), \label{eqn:bc2:3} \\
  \partial_t p_{i00} + c\: \partial_z p_{i00} &= \frac{2}{M} \sum_{k=1}^{N} c_k \partial_t (p_{i0k}-p_{i00}),
    \quad\quad \text{for } (x,z)=(L^x,L^z), \label{eqn:bc2:4} \\
  \partial_t p_{00k} + c\: \partial_y p_{00k} &= \frac{2}{M} \sum_{j=1}^{N} c_j \partial_t (p_{0jk}-p_{00k}),
    \quad\quad \text{for } (y,z)=(L^y,L^z), \label{eqn:bc2:5} \\
  \partial_t p_{0j0} + c\: \partial_z p_{0j0} &= \frac{2}{M} \sum_{k=1}^{N} c_k \partial_t (p_{0jk}-p_{0j0}),
    \quad\quad \text{for } (y,z)=(L^y,L^z), \label{eqn:bc2:6}
\end{align}
\end{subequations}
for $i,j,k = 1,\dots,N$.

We next derive governing equations for the 1D fields on the edges.
Using equations \eqref{eqn:relPijk_x}-\eqref{eqn:relPijk_z} to connect the 1D and the 2D fields, one has
\begin{subequations}
\begin{align}
  \partial_{tt} (1+c_i) (p_{ij0} - p_{0j0}) - c^2 (\partial_{yy} + \partial_{zz}) p_{ij0} &= 0,
    \quad\quad \text{for } (x,y)=(L^x,L^y), \label{eqn:lev2:ijX} \\
  \partial_{tt} (1+c_j) (p_{ij0} - p_{i00}) - c^2 (\partial_{xx} + \partial_{zz}) p_{ij0}&= 0,
    \quad\quad \text{for } (x,y)=(L^x,L^y), \label{eqn:lev2:ijY} \\
  \partial_{tt} (1+c_i) (p_{i0k} - p_{00k}) - c^2 (\partial_{yy} + \partial_{zz}) p_{i0k} &= 0,
    \quad\quad \text{for } (x,z)=(L^x,L^z), \label{eqn:lev2:ikX} \\
  \partial_{tt} (1+c_k) (p_{i0k} - p_{i00}) - c^2 (\partial_{xx} + \partial_{yy}) p_{i0k} &= 0,
    \quad\quad \text{for } (x,z)=(L^x,L^z), \label{eqn:lev2:ikZ} \\
  \partial_{tt} (1+c_j) (p_{0jk} - p_{00k}) - c^2 (\partial_{xx} + \partial_{zz}) p_{0jk} &= 0,
    \quad\quad \text{for } (y,z)=(L^y,L^z), \label{eqn:lev2:jkY} \\
  \partial_{tt} (1+c_k) (p_{0jk} - p_{0j0}) - c^2 (\partial_{xx} + \partial_{yy}) p_{0jk} &= 0,
    \quad\quad \text{for } (y,z)=(L^y,L^z), \label{eqn:lev2:jkZ}
\end{align}
\end{subequations}
for $i,j,k = 1,\dots,N$.
Unfortunately, these relations involve spatial derivatives that cannot be computed for 1D fields defined only on the edges.
For instance, only derivatives with respect to $z$ can be computed on the edge $(x,y)=(L^x,L^y)$, while equations \eqref{eqn:lev2:ijX} and \eqref{eqn:lev2:ijY} involve derivatives with respect to $y$ and $x$, respectively.
We then manipulate equations \eqref{eqn:lev2:ijX}-\eqref{eqn:lev2:jkZ} to eliminate such inadmissible derivatives.
Adding the equations corresponding to each edge and using equation \eqref{eqn:waveForPijk} give
\begin{subequations}
\begin{align}
  \partial_{tt} \left[(1+c_i+c_j) p_{ij0} - (1+c_i) p_{0j0} - (1+c_j) p_{i00}\right] - c^2 \partial_{zz} p_{ij0} &= 0,
    \quad\ \text{for } (x,y)=(L^x,L^y), \label{eqn:lev2:ijEdge} \\
  \partial_{tt} \left[(1+c_i+c_k) p_{i0k} - (1+c_i) p_{00k} - (1+c_k) p_{i00}\right] - c^2 \partial_{yy} p_{i0k} &= 0,
    \quad\ \text{for } (x,z)=(L^x,L^z), \label{eqn:lev2:ikEdge} \\
  \partial_{tt} \left[(1+c_j+c_k) p_{0jk} - (1+c_j) p_{00k} - (1+c_k) p_{0j0}\right] - c^2 \partial_{xx} p_{0jk} &= 0,
    \quad\ \text{for } (y,z)=(L^y,L^z), \label{eqn:lev2:jkEdge}
\end{align}
\end{subequations}
for $i,j,k = 1,\dots,N$.
These 1D relations involve only spatial derivatives which are well-defined on edges.

\subsubsection*{Boundary conditions at corners and 0D relations}

On the borders of each edge (\textit{i.e.}~ at the corners of the domain), the 1D field equations require boundary conditions.
We derive the compatibility conditions for the 1D fields at the corner $(x,y,z)=(L^x,L^y,L^z)$ by using a similar strategy than for the edges.

Because $p_{000}$ belongs to the kernel of $\mathcal{L}^x$, $\mathcal{L}^y$ and $\mathcal{L}^z$ and because these operators commute with $\mathcal{M}^x$, $\mathcal{M}^y$ and $\mathcal{M}^z$, the 1D fields $p_{0jk}$, $p_{i0k}$ and $p_{ij0}$ also belong to the kernel of $\mathcal{L}^x$, $\mathcal{L}^y$ and $\mathcal{L}^z$ by definition \eqref{eqn:defAddFields}.
At the corner $(x,y,z)=(L^x,L^y,L^z)$, we then have
\begin{align*}
  \mathcal{L}^x p_{0jk} &= 0, \\
  \mathcal{L}^y p_{i0k} &= 0, \\
  \mathcal{L}^z p_{ij0} &= 0,
\end{align*}
which can be rewritten
\begin{subequations}
\begin{align}
  \partial_t p_{0jk} + c \: \partial_x p_{0jk} &= \frac{2}{M} \sum_{i=1}^N c_i \partial_t (p_{ijk}-p_{0jk}),
  \label{eqn:bc3:1} \\
  \partial_t p_{i0k} + c \: \partial_y p_{i0k} &= \frac{2}{M} \sum_{j=1}^N c_j \partial_t (p_{ijk}-p_{i0k}),
  \label{eqn:bc3:2}  \\
  \partial_t p_{ij0} + c \: \partial_z p_{ij0} &= \frac{2}{M} \sum_{k=1}^N c_k \partial_t (p_{ijk}-p_{ij0}),
  \label{eqn:bc3:3}
\end{align}
\end{subequations}
for $i,j,k = 1,\dots,N$.

We now derive relations for the 0D auxiliary fields at the corner.
Using equations \eqref{eqn:relPijk_x}-\eqref{eqn:relPijk_z} to connect the 0D and the 1D fields, one has
\begin{subequations}
\begin{align}
  (1+c_i) \partial_{tt} (p_{ijk} - p_{0jk}) - c^2 (\partial_{yy} + \partial_{zz}) p_{ijk} &= 0, \label{eqn:lev3:ijkX} \\
  (1+c_j) \partial_{tt} (p_{ijk} - p_{i0k}) - c^2 (\partial_{xx} + \partial_{zz}) p_{ijk} &= 0, \label{eqn:lev3:ijkY} \\
  (1+c_k) \partial_{tt} (p_{ijk} - p_{ij0}) - c^2 (\partial_{xx} + \partial_{yy}) p_{ijk} &= 0, \label{eqn:lev3:ijkZ}
\end{align}
\end{subequations}
for $i,j,k = 1,\dots,N$.
Unfortunately, these relations again involve spatial derivatives which are inadmissible if the 0D fields are defined only on the corner.
Again, we manipulate the equations to eliminate these derivatives.
Adding equations \eqref{eqn:lev3:ijkX}--\eqref{eqn:lev3:ijkZ}, removing the spatial derivatives by using equation \eqref{eqn:waveForPijk} and integrating in time give
\begin{align}
  (1+c_i+c_j+c_k) p_{ijk} - (1+c_i)p_{0jk} - (1+c_j)p_{i0k} - (1+c_k)p_{ij0} &= 0,
  \label{eqn:lev3:ijkCorner}
\end{align}
for $i,j,k = 1,\dots,N$.
These 0D relations do not involve any derivatives, and they can be used to remove the 0D fields $p_{ijk}$ from the boundary conditions \eqref{eqn:bc3:1}-\eqref{eqn:bc3:3}.

\subsection{Formulation with first-order partial differential equations}
\label{sec:form:final}

Since all the auxiliary fields $p_{ijk}(t,\vec{x})$ satisfy the wave equation \eqref{eqn:waveForPijk}, we can define auxiliary velocities $\vec{u}_{ijk}(t,\vec{x})$ such that each pair $(p_{ijk},\vec{u}_{ijk})$ satisfies the pressure-velocity system
\begin{subequations}
\begin{align}
  \partial_t p_{ijk} + \rho c^2 \: \nabla\cdot\vec{u}_{ijk} &= 0, \label{eqn:pressureEqnAux} \\
  \rho \: \partial_t \vec{u}_{ijk} + \nabla p_{ijk} &= 0,         \label{eqn:velocityEqnAux}
\end{align}
\end{subequations}
with $i,j,k=0,\dots,N$.
Introducing these auxiliary velocities and using integration in time on the equations derived in the previous section give a HABC formulation with only algebraic relations as boundary conditions and first-order differential equations as governing equations for the auxiliary fields:
\begin{itemize}

\item on the face $x=L^x$, the boundary condition \eqref{eqn:bc1:x} becomes
\begin{align}
  p_{000} - \rho c u_{000} &= \frac{2}{M} \sum_{i=1}^{N} c_i \: (p_{i00} - p_{000}),
  \label{eqn:HABCsyst:BClev1}
\end{align}
and the 2D fields are governed by
\begin{subequations}
\begin{align}
  (1+c_i) \partial_t p_{i00} + \rho c^2 \big( \partial_y v_{i00} + \partial_z w_{i00} \big) &= (1+c_i) \partial_t p_{000},
    \label{eqn:HABCsyst:GovLev1} \\
  \rho \partial_t v_{i00} + \partial_y p_{i00} &= 0, \label{eqn:HABCsyst:GovLev1b} \\
  \rho \partial_t w_{i00} + \partial_z p_{i00} &= 0, \label{eqn:HABCsyst:GovLev1c}
\end{align}
\end{subequations}
for $i=1,\dots,N$;

\item on the edge $(x,y)=(L^x,L^y)$, the boundary conditions \eqref{eqn:bc2:1}-\eqref{eqn:bc2:2} become
\begin{subequations}
\begin{align}
  p_{0j0} - \rho cu_{0j0} &= \frac{2}{M} \sum_{i=1}^{N} c_i (p_{ij0} - p_{0j0}),
  \label{eqn:HABCsyst:BClev2a} \\
  p_{i00} - \rho cv_{i00} &= \frac{2}{M} \sum_{j=1}^{N} c_j (p_{ij0} - p_{i00}),
  \label{eqn:HABCsyst:BClev2b}
\end{align}
\end{subequations}
and the 1D fields are governed by
\begin{subequations}
\begin{align}
  (1+c_i+c_j) \partial_t p_{ij0} + \rho c^2 \partial_z w_{ij0} &= (1+c_j) \partial_t p_{i00} + (1+c_i) \partial_t p_{0j0}
    \label{eqn:HABCsyst:GovLev2} \\
  \rho \partial_t w_{ij0} + \partial_z p_{ij0} &= 0,
    \label{eqn:HABCsyst:GovLev2b}
\end{align}
\end{subequations}
for $i,j=1,\dots,N$;

\item at the corner $(x,y,z)=(L^x,L^y,L^z)$, the boundary conditions \eqref{eqn:bc3:1}-\eqref{eqn:bc3:3} become
\begin{subequations}
\begin{align}
  p_{0jk} - \rho cu_{0jk} &= \frac{2}{M} \sum_{i=1}^{N} c_i \big( (1-C_{ijk}) p_{0jk} + C_{jki} p_{i0k} + C_{kij} p_{ij0} \big) ,
  \label{eqn:HABCsyst:BClev3a} \\
  p_{i0k} - \rho cv_{i0k} &= \frac{2}{M} \sum_{j=1}^{N} c_j \big( C_{ijk} p_{0jk} + (1-C_{jki}) p_{i0k} + C_{kij} p_{ij0} \big) ,
  \label{eqn:HABCsyst:BClev3b} \\
  p_{ij0} - \rho cv_{ij0} &= \frac{2}{M} \sum_{k=1}^{N} c_k \big( C_{ijk} p_{0jk} + C_{jki} p_{i0k} + (1-C_{kij}) p_{ij0} \big) ,
  \label{eqn:HABCsyst:BClev3c}
\end{align}
\end{subequations}
with $C_{ijk} = (1+c_i)/(1+c_i+c_j+c_k)$, for $i,j,k=1,\dots,N$.

\end{itemize}
Similar relations can be obtained for the other faces, edges and corners.

The boundary treatment presented in this section then consists of a multidimensional solver with equations to solve in the volume, on the faces, on the edges and at the corners of the domain.
The formulation is summarized in this way:
\begin{itemize}
\setlength\itemsep{0em}
\item on the volume, the 3D fields (pressure and velocity) are governed by the classical pressure-velocity system with initial conditions and/or sources which are compactly supported inside the domain;
\item on the faces, one boundary condition is prescribed using the algebraic relation \eqref{eqn:HABCsyst:BClev1} and $N$ sets of 2D auxiliary fields are governed by the first-order differential equations \eqref{eqn:HABCsyst:GovLev1}-\eqref{eqn:HABCsyst:GovLev1c};
\item on the edges, $2N$ boundary conditions ($N$ for each adjacent face) are prescribed using the algebraic relations \eqref{eqn:HABCsyst:BClev2a}-\eqref{eqn:HABCsyst:BClev2b} and $N^2$ sets of 1D auxiliary fields are governed by the first-order differential equations \eqref{eqn:HABCsyst:GovLev2}-\eqref{eqn:HABCsyst:GovLev2b};
\item at the corners, $3N^2$ boundary conditions ($N^2$ for each adjacent edge) are prescribed using the algebraic relations \eqref{eqn:HABCsyst:BClev3a}-\eqref{eqn:HABCsyst:BClev3c}.
\end{itemize}
The equation numbers correspond to the face $x=L^x$, the edge $(x,y)=(L^x,L^y)$ and the corner $(x,y,z)=(L^x,L^y,L^z)$.
By the assumption of the compact support of the initial condition inside the domain, the initial conditions for the auxiliary variable fields are all zero.

\subsection{Extension to mixed boundary conditions}
\label{sec:form:mixed}

We briefly address the case where a homogeneous boundary condition is prescribed on one or more faces of the computational domain.
In exploration geophysics, for instance, the computational cuboidal domain must represent the underground structure.
A HABC can be used on the lateral and bottom faces of the domain, while the so-called \textit{free-surface} boundary condition (which is $p=0$ in the acoustic model) must be prescribed on the upper face to simulate the Earth's surface.

A homogeneous boundary condition on the 3D fields is straightforwardly incorporated into the boundary procedure by using the same condition on the auxiliary 2D and 1D fields.
Indeed, let us consider the homogeneous Dirichlet condition $p=0$ on one face of the domain and a HABC on the adjacent faces.
Seen definition \eqref{eqn:defAddFields}, the 2D auxiliary fields of the adjacent faces must be set to zero on the edges if $p=0$.
Therefore, no auxiliary field must be computed on edges and corners belonging to faces where a homogeneous boundary condition is prescribed.

\subsection{Explicit time-stepping procedure}
\label{sec:scheme:Time}

Because of the coupling between the differential equations and the algebraic relations of the different levels, the formulation described in section \ref{sec:form:final} cannot be straightforwardly solved with an explicit time-stepping procedure.
For instance, the time derivatives of both 3D pressure and 2D pressures appear in equation \eqref{eqn:HABCsyst:GovLev1}, and the time derivatives of both 2D pressures and 1D pressures appear in equation \eqref{eqn:HABCsyst:GovLev2}.
In order to allow an explicit time-stepping procedure, we first reformulate the boundary conditions introducing temporary variables and using characteristic variables (as in \cite{hagstrom2004new}).
In the directions $\vec{e}_x$, $\vec{e}_y$, $\vec{e}_z$, the characteristic variables read
\begin{subequations}
\begin{align}
  r^{-x}_{ijk} &= p_{ijk} - \rho c \: u_{ijk}, &
  r^{+x}_{ijk} &= p_{ijk} + \rho c \: u_{ijk}, \label{eqn:proc:charact1} \\
  r^{-y}_{ijk} &= p_{ijk} - \rho c \: v_{ijk}, &
  r^{+y}_{ijk} &= p_{ijk} + \rho c \: v_{ijk}, \label{eqn:proc:charact2} \\
  r^{-z}_{ijk} &= p_{ijk} - \rho c \: w_{ijk}, &
  r^{+z}_{ijk} &= p_{ijk} + \rho c \: w_{ijk}, \label{eqn:proc:charact3}
\end{align}
\end{subequations}
where $r^{-x}_{ijk}$, $r^{-y}_{ijk}$, $r^{-z}_{ijk}$ and $r^{+x}_{ijk}$, $r^{+y}_{ijk}$, $r^{+z}_{ijk}$ contain information traveling downwardly and upwardly, respectively, along the direction $\vec{e}_x$, $\vec{e}_y$, $\vec{e}_z$.
The comprehensive boundary formulation reads

\begin{itemize}

\item on the face $x=L^x$, temporary fields $p_{i00}^\star$ are defined as
\begin{align}
  p_{i00}^\star &= p_{i00} - (r^{-x}_{000}+r^{+x}_{000}),
  \label{eqn:proc:pi00Star}
\end{align}
the boundary condition gives
\begin{align}
  r^{-x}_{000} = \frac{1}{M} \sum_{i=1}^N c_i \: p_{i00}^\star,
  \label{eqn:proc:faceBC}
\end{align}
and the 2D relations becomes
\begin{subequations}
\begin{align}
  (1+c_i) \  \partial_t p_{i00}^\star + \rho c^2 \big( \partial_y v_{i00} + \partial_z w_{i00} \big) &= 0,
  \label{eqn:proc:2Dpvsyst1} \\
  \rho \partial_t v_{i00} + \partial_y p_{i00} &= 0,
  \label{eqn:proc:2Dpvsyst2} \\
  \rho \partial_t w_{i00} + \partial_z p_{i00} &= 0,
  \label{eqn:proc:2Dpvsyst3}
\end{align}
\end{subequations}
for $i=1,\dots,N$;

\item on the edge $(x,y)=(L^x,L^y)$, temporary fields $p_{ij0}^\star$ are defined as
\begin{align}
  p_{ij0}^\star &= p_{ij0} - \big( C_{ji}  (r^{-y}_{i00}+r^{+y}_{i00}) + C_{ij} (r^{-x}_{0j0}+r^{+x}_{0j0}) \big),
  \label{eqn:proc:pij0Star}
\end{align}
with $C_{ij} = (1+c_i)/(1+c_i+c_j)$, the boundary conditions give
\begin{subequations}
\begin{align}
  r^{-y}_{i00}
    &= \frac{1}{M} \sum_{j=1}^N c_j \: \big(p_{ij0}^\star
      + (C_{ji}-1) \: (r^{-y}_{i00}+r^{+y}_{i00})
      + C_{ij} \: (r^{-x}_{0j0}+r^{+x}_{0j0})\big),
  \label{eqn:proc:edgeBC1} \\
  r^{-x}_{0j0}
    &= \frac{1}{M} \sum_{i=1}^N c_i \: \big(p_{ij0}^\star
      + C_{ji} \: (r^{-y}_{i00}+r^{+y}_{i00})
      + (C_{ij}-1) \: (r^{-x}_{0j0}+r^{+x}_{0j0})\big),
  \label{eqn:proc:edgeBC2}
\end{align}
\end{subequations}
and the 1D relations becomes
\begin{subequations}
\begin{align}
  (1+c_i+c_j) \ \partial_t p_{ij0}^\star + \rho c^2 \partial_z w_{ij0} &= 0,
  \label{eqn:proc:1Dpvsyst1} \\
  \rho \partial_t w_{ij0} + \partial_z p_{ij0} &= 0,
  \label{eqn:proc:1Dpvsyst2}
\end{align}
\end{subequations}
for $i,j=1,\dots,N$;

\item at the corner $(x,y,z)=(L^x,L^y,L^z)$, the boundary conditions give
\begin{subequations}
\begin{align}
  r^{-x}_{0jk}
  = \frac{1}{M} \sum_{i=1}^N c_i \: \big(
        (C_{ijk}-1) (r^{-x}_{0jk}+r^{+x}_{0jk})
      + C_{jki}     (r^{-y}_{i0k}+r^{+y}_{i0k})
      + C_{kij}     (r^{-z}_{ij0}+r^{+z}_{ij0}) \big),
  \label{eqn:proc:cornerBC1}
\\
  r^{-y}_{i0k}
  = \frac{1}{M} \sum_{j=1}^N c_j \: \big(
        C_{ijk}     (r^{-x}_{0jk}+r^{+x}_{0jk})
      + (C_{jki}-1) (r^{-y}_{i0k}+r^{+y}_{i0k})
      + C_{kij}     (r^{-z}_{ij0}+r^{+z}_{ij0}) \big),
  \label{eqn:proc:cornerBC2}
\\
  r^{-z}_{ij0}
  = \frac{1}{M} \sum_{k=1}^N c_k  \: \big(
        C_{ijk}     (r^{-x}_{0jk}+r^{+x}_{0jk})
      + C_{jki}     (r^{-y}_{i0k}+r^{+y}_{i0k})
      + (C_{kij}-1) (r^{-z}_{ij0}+r^{+z}_{ij0}) \big),
  \label{eqn:proc:cornerBC3}
\end{align}
\end{subequations}
for $i,j,k=1,\dots,N$.

\end{itemize}
Similar relations can be obtained for the other faces, edges and corners.

\input{algoHABC2}

The numerical solution of this formulation can be computed with an explicit time-stepping scheme by solving the different levels successively at each time step, starting with the 3D and ending with the 0D.
The complete procedure to update the solution from time $t_1$ to time $t_2=t_1+\Delta t$ is sketched in algorithm \ref{algo:multiSolver}.
The computational load is mainly due to two kinds of operations:
\begin{enumerate}

\item
The 3D, 2D, 1D solvers solve first-order differential systems to update the corresponding pressure and velocity fields.
Since the systems of the faces (equations \eqref{eqn:proc:2Dpvsyst1}-\eqref{eqn:proc:2Dpvsyst3}) and the edges (equations \eqref{eqn:proc:1Dpvsyst1}-\eqref{eqn:proc:1Dpvsyst2}) resemble to pressure-velocity systems, we use the same numerical scheme to solve the systems over each dimension (\textit{i.e.}~ volume, faces and edges).
To update the 3D, 2D, 1D fields at $t=t_2$, the 3D, 2D, 1D solvers use boundary conditions based on incoming characteristics computed at $t=t_1$.
Such characteristic-based boundary conditions are naturally incorporated in discontinuous Galerkin formulations (see section \ref{sec:scheme:DG}).

\item
The 2D, 1D, 0D solvers update respectively the 3D, 2D, 1D incoming characteristics at $t=t_2$, using the fields which are already computed at $t=t_2$.
The computation is straightforward on the faces (equation \eqref{eqn:proc:faceBC}), but it requires the solution of linear systems with $2N$ unknowns on the edges (equations \eqref{eqn:proc:edgeBC1}-\eqref{eqn:proc:edgeBC2}) and $3N^2$ unknowns at the corners (equations \eqref{eqn:proc:cornerBC1}-\eqref{eqn:proc:cornerBC3}).
For instance, for the edge $(x,y)=(L^x,L^y)$, the $2N \times 2N$  system with equations \eqref{eqn:proc:edgeBC1}-\eqref{eqn:proc:edgeBC2} can be written
\begin{align}
  \left[\begin{array}{cc} M\vec{I}-\vec{C}^{(1)} & \vec{C}^{(2)} \\ \vec{C}^{(2)} & M\vec{I}-\vec{C}^{(1)} \end{array}\right]
  \left[\begin{array}{c} \vec{r}^{-x} \\ \vec{r}^{-y} \end{array}\right]
  =
  \left[\begin{array}{cc} \vec{C}^{(1)} & \vec{C}^{(2)} \\ \vec{C}^{(2)} & \vec{C}^{(1)} \end{array}\right]
  \left[\begin{array}{c} \vec{r}^{+x} \\ \vec{r}^{+y} \end{array}\right]
  +
  \left[\begin{array}{c} \vec{b}^{x} \\ \vec{b}^{y} \end{array}\right]
\label{eqn:systEdgeSymmetry}
\end{align}
with
\begin{align*}
  \vec{r}^{-x} &=
  \left[
    r^{-x}_{010} \;\;\;
    r^{-x}_{020} \;\;\;
    \cdots \;\;\;
    r^{-x}_{0N0}
  \right]^T,
&
  \vec{r}^{+x} &=
  \left[
    r^{+x}_{010} \;\;\;
    r^{+x}_{020} \;\;\;
    \cdots \;\;\;
    r^{+x}_{0N0}
  \right]^T,
\\
  \vec{r}^{-y} &=
  \left[
    r^{-y}_{100} \;\;\;
    r^{-y}_{200} \;\;\;
    \cdots \;\;\;
    r^{-y}_{N00}
  \right]^T,
&
  \vec{r}^{-y} &=
  \left[
    r^{+y}_{100} \;\;\;
    r^{+y}_{200} \;\;\;
    \cdots \;\;\;
    r^{+y}_{N00}
  \right]^T,
\end{align*}
\begin{align*}
  \vec{b}^{x} &=
  \Big[ \textstyle
    \sum_i c_i p^\star_{i10} \;\;\;
    \sum_i c_i p^\star_{i20} \;\;\;
    \cdots \;\;\;
    \sum_i c_i p^\star_{iN0}
  \Big]^T,
\\
  \vec{b}^{y} &=
  \Big[ \textstyle
    \sum_j c_j p^\star_{1j0} \;\;\;
    \sum_j c_j p^\star_{2j0} \;\;\;
    \cdots \;\;\;
    \sum_j c_j p^\star_{Nj0}
  \Big]^T,
\end{align*}
where $\mat{I}$ is the $N \times N$ identity matrix and $\mat{C}^{(1)}$ and $\mat{C}^{(2)}$ are a $N \times N$ sparse matrices that only depend on the parameters $c_i$'s.
The vectors of incoming characteristics, $\vec{r}^{-x}$ and $\vec{r}^{-y}$, contain the unknowns of system \eqref{eqn:systEdgeSymmetry}, while $\vec{r}^{+x}$, $\vec{r}^{+y}$, $\vec{b}^x$ and $\vec{b}^y$ are computed using the fields at $t=t_2$.
Because the matrix of the system (\textit{i.e.}~ the matrix in the left-hand side) is the same for all the edges, we precompute and store its inverse, after which we only need to do matrix-vector multiplications at each time step.
The matrix is composed of four $N \times N$ matrices, only two of which are independent.
Since its inverse has the same structure, only two $N \times N$ matrices must then be stored in memory.
Similarly, for the corner conditions, the inverse matrix of the $3N^2 \times 3N^2$ system \eqref{eqn:proc:cornerBC1}-\eqref{eqn:proc:cornerBC3} can be precomputed and stored.
It is composed of nine $N^2 \times N^2$ matrices, only three of which are independent and must be stored.

\end{enumerate}

This procedure is similar to the one described by Hagstrom and Warburton \cite{hagstrom2004new} for 2D cases, but with a different choice for the planar HABC.
Their formulation has been extended to a family of more general HABC, the complete radiation boundary conditions (CRBC) \cite{hagstrom2009complete,hagstrom2010radiation}, which can be accurate for both traveling and evanescent waves.
With a specific choice of parameters, corresponding to the Pad\'e case, these HABCs are equivalent to the one used here, but the formulations are written differently.

With the HABC of Hagstrom and Warburton in the Pad\'e case, fields from all the levels appear in the compatibility condition at corners.
This leads to an inconsistent formulation when using discontinuous Galerkin schemes based on an unstructured mesh.
Indeed, both 3D and 2D fields can have more than one value at the corners of the domain if several tetrahedral or triangular mesh cells touch this corner.
The same inconsistency appears when deriving the compatibility conditions with the HABC proposed by Collino \cite{collino1993conditions,collino1993high}.
By contrast, with our formulation, the 0D fields $p_{ijk}$ are defined only with 1D fields (equation \eqref{eqn:lev3:ijkCorner}), which have only one value at corners.
This observation motivated our choice for the HABC specifically based on the approximate square root represented by equation \eqref{eqn:approxRoot2}.

\subsection{Discontinuous Galerkin time domain scheme}
\label{sec:scheme:DG}

The HABC and compatibility conditions are discretized using a nodal discontinuous Galerkin finite element method with upwind fluxes in space and a low-storage fourth-order Runge-Kutta method in time \cite{hesthaven2007nodal}.
In the complete boundary procedure for cuboidal domains, the pressure-velocity system must be solved on the edges, on the faces and in the volume of the domain, which leads to a multidimensional solver.
We have used the 1D, 2D and 3D versions of the same spatial scheme.

The cuboidal domain is partitioned into a volume mesh of $K^\text{tet}$ non-overlapping tetrahedral cells, $\Omega = \bigcup_k \mathsf{D}_k^\text{tet}$, where $\mathsf{D}_k^\text{tet}$ is the $k^\text{th}$ cell.
Surface and line meshes are built on this volume mesh: the surface mesh is composed of the cell faces belonging to faces of the domain where the HABC is prescribed, while the line mesh is composed of the cell edges belonging to edges of the domain where two HABCs cross.
We denote by $K^\text{tri}$ the number of triangular cells $\mathsf{D}_k^\text{tri}$, by $K^\text{lin}$ the number of line cells $\mathsf{D}_k^\text{lin}$, and by $K^\text{pnt}$ the number of corners where three HABCs cross.
For instance, we have $K^\text{pnt}=8$ if the original problem is defined on the infinite space $\mathbb{R}^3$, and $K^\text{pnt}=4$ if a homogeneous Dirichlet boundary condition is prescribed on one face of the domain.
The volume mesh has $F_\text{bnd}^\text{tet}=K^\text{tri}$ boundary cell faces where a HABC are prescribed.
Since, an a cube, each edge has two neighboring faces and each corner has three neighboring edges, we assume that the surface and line meshes have respectively $F_\text{bnd}^\text{tri}=2K^\text{lin}$ and $F_\text{bnd}^\text{lin}=3K^\text{pnt}$ boundary cell faces where a HABC are prescribed.

For each of the 1D, 2D and 3D solvers, the pressure fields and the Cartesian components of the velocity fields are approximated by piecewise polynomial functions, which are discontinuous at the interface between two cells.
The discrete unknowns correspond to the values of fields at nodes distributed over the boundary and the interior of an element \cite{hesthaven2007nodal}.
In this work, the spatial distribution of nodes in the reference tetrahedron is defined using the \textit{Warp \& Blend} technique \cite{warburton2006explicit}.
The nodes in the reference triangle and on the reference line are chosen to match the face nodes and the edge nodes, respectively, of the tetrahedron.

The spatial scheme is built on a variational form of the equations.
Hereafter, the material properties $\rho$ and $c$ are assumed to be constant over each cell, but potentially discontinuous at the interfaces.
For each line, triangular and tetrahedral cell $\mathsf{D}_k$, we consider the variational form
\begin{align}
     \int_{\mathsf{D}_k} \dvp{p}{t}\:\psi \:d\vec{x}
   + \int_{\mathsf{D}_k} \rho c^2\:(\nabla\cdot\vec{u}) \: \psi \:d\vec{x}
   + \int_{\partial \mathsf{D}_k} (\rho c^2)^{\text{int}} \: \left((\vec{n}\cdot\vec{u})^\star-(\vec{n}\cdot\vec{u})^{\text{int}}\right) \: \psi \: d\vec{x}
  &= 0, \label{eqn:weakDG1}
  \\
     \int_{\mathsf{D}_k} \dvp{\vec{u}}{t}\cdot\boldsymbol{\psi}\:d\vec{x}
   + \int_{\mathsf{D}_k} \frac{1}{\rho} \: (\nabla p)\cdot\boldsymbol{\psi} \:d\vec{x}
   + \int_{\partial \mathsf{D}_k} \frac{1}{\rho^{\text{int}}} \: (p^\star-p^{\text{int}}) \: (\vec{n}\cdot\boldsymbol{\psi}) \:d\vec{x}
  &= 0, \label{eqn:weakDG2}
\end{align}
where $\psi(\vec{x})$ and $\boldsymbol{\psi}(\vec{x})$ are test functions, $\partial\mathsf{D}_k$ is the cell boundary and $\vec{n}$ is the outward unit normal to $\partial\mathsf{D}_k$.
The boundary conditions are prescribed and the solutions at the interface between two cells are coupled by selecting specific values for the numerical fluxes $p^\star$ and $(\vec{n}\cdot\vec{u})^\star$ in the boundary integrals of both equations.
At the interface between two elements, we consider the classical upwind fluxes provided by the exact Riemann solver \cite{wilcox2010high,hesthaven2007nodal,leveque2002finite},
\begin{align*}
  p^\star   &= \frac{\left\{p/(\rho c)\right\} - \vec{n}\cdot\jump{\vec{u}}}{\left\{1/(\rho c)\right\}}, \\
  \left(\vec{n}\cdot\vec{u}\right)^\star &= \frac{\vec{n}\cdot\left\{\rho c \vec{u}\right\} - \jump{p}}{\left\{\rho c\right\}},
\end{align*}
where $\{X\}=(X^{\text{ext}}+X^{\text{int}})/2$ and $\jump{X}=(X^{\text{ext}}-X^{\text{int}})/2$ are the average and the semi-jump, respectively, of any scalar or vector $X$.
The superscripts $\:^{\text{ext}}$ and $\:^{\text{int}}$ denote the exterior and interior values at the interface.
If the medium is homogeneous at the interface, the numerical fluxes can be conveniently rewritten as
\begin{align}
  p^\star &= (r^+)^{\text{int}} + (r^-)^{\text{ext}}, \label{eqn:fluxDG1} \\
  \left(\vec{n}\cdot\vec{u}\right)^\star &= \frac{(r^+)^{\text{int}} - (r^-)^{\text{ext}}}{\rho c}, \label{eqn:fluxDG2}
\end{align}
where $r^+$ and $r^-$ are the outgoing and incoming characteristics, respectively, defined as
\begin{align*}
  r^+ &= p + \rho c \: (\vec{n}\cdot\vec{u}), \\
  r^- &= p - \rho c \: (\vec{n}\cdot\vec{u}).
\end{align*}
At the domain boundary, the basic ABC is straightforwardly incorporated in the formulation by using the numerical fluxes \eqref{eqn:fluxDG1}-\eqref{eqn:fluxDG2} with the incoming characteristic equal to zero.
In the HABC procedure, boundary conditions for the 3D, 2D and 1D solvers are enforced by defining incoming characteristics using the 2D, 1D and 0D solvers, respectively.
Finally, the homogeneous boundary condition $p=0$ is enforced by taking
\begin{align*}
  p^\star &= 0, \\
  \left(\vec{n}\cdot\vec{u}\right)^\star &= \vec{n}\cdot\vec{u}^{\text{int}} + \frac{p^{\text{int}}}{(\rho c)^{\text{int}}}.
\end{align*}

For each element $\mathsf{D}_k$, the semi-discrete equations are obtained by substituting the semi-discrete fields into the variational form \eqref{eqn:weakDG1}-\eqref{eqn:weakDG2}, and using the Lagrange polynomials as test functions \cite{hesthaven2007nodal}.
For each field, this leads to a system that reads
\begin{align*}
  \frac{d\vec{q}_k}{dt} = \vec{r}_k,
\end{align*}
where the vectors $\vec{q}_k$ and $\vec{r}_k$ contain the discrete unknowns and the values of the right-hand side  terms for $\mathsf{D}_k$.
The right-hand side vector can be written as
\begin{equation}
  \vec{r}_{k}
    = \sum_{i=1}^{N_\text{dim}} \sum_{j=1}^{N_\text{dim}}
      g_{k,i,j}^{\text{vol}} \: \vec{D}_{j} \: \vec{f}_{k,j}
    + \sum_{f=1}^{N_\text{faces}}
      g_{k,f}^{\text{sur}} \: \vec{L}_{f} \: \vec{p}_{k,f},
  \label{eqn:dgLocalRhs}
\end{equation}
where $N_\text{dim}$ is the spatial dimension of the element, $N_\text{faces}$ is the number of faces, $\vec{f}_{k,j}$ corresponds to the physical flux in the $x_j$-direction for all the nodes of $\mathsf{D}_k$, and the vector $\vec{p}_{k,f}$ contains the boundary term for all the nodes belonging to the face $f$.
In the right-hand side vector \eqref{eqn:dgLocalRhs}, the first term (called the \textit{volume term}) corresponds to the integrals over the cell $\mathsf{D}_k$, and the second term (called the \textit{surface term}) corresponds to those over its boundary $\partial\mathsf{D}_k$.
The matrices $\mat{D}_{j}$ and $\mat{L}_{f}$ are respectively differentiation and lifting matrices for the reference element, while the geometric factors $g_{k,i,j}^{\text{vol}}$ and $g_{k,f}^{\text{sur}}$ depend on the shape of each element.
The matrices and factors are defined in \cite{hesthaven2007nodal}.
The semi-discrete equations are explicitly derived in \cite{modave2015nodal} for the three-dimensional case.

The low-storage fourth-order Runge-Kutta scheme is used for time discretization.
This scheme has five stages and require the storage of an auxiliary residual vector $\vec{s}_k$.
At each stage $n$ of each time iteration $m$, the residual vector $\vec{s}_k^{m+n/5}$ and the unknown vector $\vec{q}_k^{m+n/5}$ are updated according to
\begin{align}
  \vec{s}_k^{m+n/5} &= a_n \: \vec{s}_k^{m+(n-1)/5} + \vec{r}_k^{m+(n-1)/5}, \label{eqn:timeStepProc1} \\
  \vec{q}_k^{m+n/5} &= \vec{q}_k^{m+(n-1)/5} + b_n \: \vec{s}_k^{m+n/5},     \label{eqn:timeStepProc2}
\end{align}
where $\vec{r}_k^{m+n/5}$ and $\vec{s}_k^{m+n/5}$ correspond to the vector computed at time $t=(m+c_i n)\Delta t$.
The values of the coefficients $a_n$, $b_n$ and $c_n$ can be found in \citep{carpenter1994fourth}.

\subsection{GPU-accelerated computational implementation}
\label{sec:implementation}

We have implemented the boundary procedure in a discontinuous Galerkin code programmed using the \CC~language with the OCCA library \cite{medina2014occa} for GPU computing.
Discontinuous finite element schemes have attractive features for parallel computing on multi-threading devices such as GPU, but a careful implementation is required to optimize the efficiency of the solver (see \textit{e.g.}~ \cite{klockner2009nodal,modave2016gpu,fuhry2014discontinuous}).
In order to improve the computational efficiency, implementation strategies have been studied for advanced discontinuous Galerkin schemes with hybrid meshes \cite{chan2016gpu,chan2016reduced}, Bernstein-Bezier basis functions \cite{chan2015gpu}, multi-rate time-stepping schemes \cite{godel2010gpu,gandham2015gpu} and distributed parallel computing on GPU clusters \cite{modave2015nodal} in several application contexts.
In this work, we propose a single-GPU implementation based on the nodal discontinuous Galerkin method and the  time-stepping scheme presented in the previous section.
We highlight that the implementation strategies used here are compatible with those presented in the above references.

The 3D solver is implemented following strategies described in \cite{modave2015nodal,modave2016gpu}.
We have implemented the 1D and 2D solvers in a similar way, with the supplementary tasks required for the boundary procedure.
A specific implementation has been conceived for the 0D solver, which only solves the compatibility system at corners using the inverse matrix of this system.

\subsubsection*{Memory management}

All the data required for computation are stored in the global memory of the GPU.
For each of the 1D, 2D and 3D solvers, a floating-point array \texttt{q} stores all the discrete unknowns of the solver, while the array \texttt{qf} contains a copy of traces associated to face nodes (\textit{i.e.}~ $p$ and $\vec{n}\cdot\vec{u}$, where $\vec{n}$ is the outward unit normal to the face).
The arrays \texttt{rhs} and \texttt{res} store the right-hand side terms and the residual, respectively, used for the time-stepping procedure \eqref{eqn:timeStepProc1}-\eqref{eqn:timeStepProc2}.
The array \texttt{qb} contains the incoming characteristic variables used as boundary condition at boundary nodes of each mesh.
Arrays are used to store the elemental matrices (\texttt{Drst} and \texttt{Lift}) and the geometric and physical parameters (\texttt{volPar} and \texttt{surPar}) required to compute the right-hand side terms \eqref{eqn:dgLocalRhs}.
Additional arrays used by the 0D, 1D and 2D solvers for the boundary procedure are allocated: the HABC coefficients (\texttt{coefHabc}), the inverse matrices of the compatibility systems (\texttt{matHabcEdge} and \texttt{matHabcCorner}) and a temporary storage \texttt{q0} (its purpose is explained later).

An array of integers, the connectivity array \texttt{map}, is used when computing the numerical fluxes at face nodes in the surface kernel.
This array has one entry for each face node: a positive value gives the address of the corresponding face node on the neighboring cell  in \texttt{qf}, a negative value corresponds to a characteristic-based boundary condition and gives the address of the incoming characteristic in \texttt{qb}, and a zero value corresponds to an homogeneous Dirichlet boundary condition.
Connectivity arrays between the nodes and faces nodes of meshes with different spatial dimensions are also used.

The granularity of storage of all the arrays has been chosen in order to maximize coalescing transfers and data reuse (see \textit{e.g.}~ \cite{klockner2009nodal,modave2016gpu}).
The main parameters of the solvers are defined in table \ref{tab:impHabc:params}.
The sizes and granularity of the main arrays are given in table \ref{tab:impHabc:arrays}.

\begin{table}[!t]
\centering
\small
\begin{tabular}{l|l|c|c|c} \hline
  Definition & Symbol & 1D & 2D & 3D \\ \hline
  Number of elements in the mesh       & $K$              & $K^\text{lin}$ & $K^\text{tri}$ & $K^\text{tet}$ \\
  Number of HABC boundary faces        & $F_\text{bnd}$   & $3K^\text{pnt}$ & $2K^\text{lin}$ & $K^\text{tri}$ \\
  Spatial dimension           & $N_\text{dim}$   & 1 & 2 & 3 \\
  Number of faces per element & $N_\text{faces}$ & 2 & 3 & 4 \\
  Number of nodes per element & $N_p$    &\small $P+1$ &\small $(P+1)(P+2)/2$ &\small $(P+1)(P+2)(P+3)/6$ \\
  Number of nodes per face    & $N_{fp}$ &\small $1$   &\small $P+1$          &\small $(P+1)(P+2)/2$ \\
  Number of (scalar) fields   & $N_\text{fields}$ & 2 & 3 & 4 \\
  Number of traces            & $N_\text{traces}$ & 2 & 2 & 2 \\
  Number of set of fields     & $N_\text{sets}$ & $N^2$ & $N$ & $1$ \\ \hline
\end{tabular}
\caption{Definition and values of the differents parameters used in the multidimensional implementation.
$P$ is the polynomial degree of the basis functions.
$N$ is the number of auxiliary fields in the planar HABC.}
\label{tab:impHabc:params}
\end{table}

\begin{table}[!t]
\centering
\small
\begin{tabular}{l|l|l} \hline
  Definition & Symbol & Size \\ \hline
    Unknown fields at nodes & \texttt{q} & $K \cdot N_\text{sets} \cdot N_\text{fields} \cdot N_p$ \\
    Unknown traces at face nodes & \texttt{qf} & $K \cdot N_\text{sets} \cdot N_\text{traces} \cdot N_\text{faces} \cdot N_{fp} $ \\
    Incoming characteristics at boundary face nodes & \texttt{qb} & $F_\text{bnd} \cdot N_\text{sets} \cdot N_{fp} $ \\
    Right-hand side array & \texttt{rhs} & $K \cdot N_\text{sets} \cdot N_\text{fields} \cdot N_p$ \\
    Residual array & \texttt{res} & $K \cdot N_\text{sets} \cdot N_\text{fields} \cdot N_p$ \\
    Differentiation matrices ($\mat{D}_1$, $\mat{D}_2$, \dots) & \texttt{Drst} & $N_p^2 \cdot N_\text{dim}$ \\
    Lifting matrices ($\mat{L}_1$, $\mat{L}_2$, \dots) & \texttt{Lift} & $N_\text{faces} \cdot N_{fp} \cdot N_p$ \\
    HABC coefficients $c_i$'s &  \texttt{coefHabc} & $N$ \\
    Inverse matrix of the compatibility system at edges & \texttt{matHabcEdge}    & $2\cdot N^2$ \\
    Inverse matrix of the compatibility system at corners & \texttt{matHabcCorner} & $3\cdot N^4$ \\ \hline
\end{tabular}
\caption{Main arrays stored in the global memory of the GPU for the 1D, 2D and 3D solvers.
The last three arrays are used by the 0D, 1D and 2D solvers.
Sizes of arrays are written from the coarsest to the finest granularity of storage.
The symbols are defined in table \ref{tab:impHabc:params}.
}
\label{tab:impHabc:arrays}
\end{table}

\subsubsection*{Kernels}

The computational procedure is decomposed into several subtasks implemented in separate OCCA kernels.
This allows us to optimize each task considering the properties of both the task and the GPU.
Our implementation has three main kernels for each of the 1D, 2D, 3D solvers:
\begin{enumerate}
\setlength\itemsep{0em}
  \item the \textit{volume kernel} computes the first term of the right-hand side vector \eqref{eqn:dgLocalRhs};
  \item the \textit{surface kernel} computes the second term of the right-hand side vector \eqref{eqn:dgLocalRhs};
  \item the \textit{update kernel} performs the time stepping (equations \eqref{eqn:timeStepProc1}-\eqref{eqn:timeStepProc2}) and updates the incoming and outgoing characteristics by performing the operations listed in algorithm \ref{algo:multiSolver}.
\end{enumerate}
For the 0D solver, there is a single \textit{update kernel}, which updates the incoming characteristics used by 1D solver.
All the kernels are called at each stage of each time step in a specific order: first the three volume kernels (in any order), then the three surface kernels (in any order), and finally the four update kernels (starting with the 3D and ending with the 0D, following the procedure in algorithm \ref{algo:multiSolver}).

The volume and surface kernels consist of streaming operations and element-wise matrix-vector multiplications.
In a nutshell, the volume kernels load the values of fields from \texttt{q} for each element, compute the physical fluxes at each node, perform the matrix-vector products using \texttt{Drst}, and store the result in \texttt{rhs}.
The surface kernels load the values of traces and incoming characteristics from \texttt{qf} and \texttt{qb} for each element, compute the numerical fluxes at each face node, perform the matrix-vector products using \texttt{Lift}, and update \texttt{rhs} with the result.

All the volume and surface kernels are written and optimized in a similar way.
In the GPU programming model, a \textit{thread} is the smallest sequence of instructions that are managed independently with their own private memory.
Threads belonging to the same \textit{thread block} run concurrently and can collaborate using shared memory.
Following \cite{klockner2009nodal}, the tasks of the volume and surface kernels are parallelized by associating one thread to the computational work required for one node, and by associating one thread block to several elements.
In the volume and surface kernels, $N_p$ and $\max(N_p, N_{\text{faces}} N_{fp})$ threads are dedicated to one element, and one thread block is dedicated to $K_\text{blkV}$ and $K_\text{blkS}$ elements, respectively.
The parameters $K_\text{blkV}$ and $K_\text{blkS}$ provide a way to tune the occupation of the GPU for each kernel of each solver.
The 3D kernels and further details about the optimization strategies can be found in \cite{modave2015nodal}.
In the 1D and 2D kernels, several sets of fields ($N^2$ and $N$, respectively) are associated to each node, and the operations are performed several times with the different fields, still associating one node per thread.
Since the elemental matrices and parameters are identical for each set of fields, the kernels are written to enable reuse of these data.
Aside from this difference, the 1D and 2D kernels are similar to the 3D kernels.

\input{algoKernelUpdate3D}
\input{algoKernelUpdate2D}
\input{algoKernelUpdate1D}
\input{algoKernelUpdate0D}

The update kernels perform the time-stepping for the fields, and computes the incoming and outgoing characteristics following the boundary procedure in algorithm \ref{algo:multiSolver}.
The operations performed by the update kernels are rather different: the 3D kernel only updates fields (algorithm \ref{algo:updateKernel3}), the 2D kernel updates fields and performs streaming operations (algorithm \ref{algo:updateKernel2}), the 1D kernel updates fields and solves $2N\times2N$ systems (algorithm \ref{algo:updateKernel1}), the 0D kernel only solves $3N^2\times3N^2$ systems (algorithm \ref{algo:updateKernel0}).
The different kernels are written in algorithms \ref{algo:updateKernel3}-\ref{algo:updateKernel0} using pseudo-code to give an overview of the implementation.
In these algorithms, \textbf{parfor} denotes a parallel loop, while \textbf{for} denotes a sequential loop.
The most external parallel loop iterates over thread blocks, while the others iterate over threads.
Since the physical/geometrical factor arrays and the three connectivity arrays are used in a straightforward way, they are not mentioned in the algorithms for the sake of clarity.
We describe hereafter the key aspects of the kernels.

\begin{itemize}

\item
The 3D update kernel (algorithm \ref{algo:updateKernel3}) performs the time-stepping for the 3D fields in three steps: first the residual is updated at nodes with equation \eqref{eqn:timeStepProc1}, then the fields are updated at nodes with equation \eqref{eqn:timeStepProc2}, and finally the traces are computed at face nodes.
As for the volume and surface kernels, each thread deals with the tasks associated to a given node, and each thread block deals with $K_\text{blkU}^\text{tet}$ elements.
There are therefore $K_\text{blkU}^\text{tet}\cdot\max(N_p^\text{tet},N_\text{faces}^\text{tet}\cdot N_{fp}^\text{tet})$ threads per thread block.

\item
The time-stepping of the 2D fields is performed by the 2D update kernel using the same parallelization strategy, with one node per thread and $K_\text{blkU}^\text{tri}$ elements per thread blocks.
Since there are $N$ sets of 2D pressure and velocity fields per node, a sequential loop is used to iterate over the sets for each operation over these fields (algorithm \ref{algo:updateKernel2}).
In addition, the 3D incoming characteristic is computed sequentially by performing the operations described in algorithm \ref{algo:multiSolver}.

In this procedure, both the previous value and the updated value of the 3D outgoing characteristic are needed.
The first is used to compute the temporary 2D fields $p^\star$ at the beginning of the time step, and the second is used when updating the 3D incoming characteristic.
While the updated 3D outgoing characteristic is computed using the updated traces from array $\texttt{qf}^\text{tet}$, we have introduced an additional array $\texttt{q0}^\text{tet}$ that stores the sum of the 3D characteristics to compute $p^\star$ at the beginning of the time step.
This array is updated when the updated 3D incoming characteristic is available.

\item
The role of the 1D update kernel is similar to that of the 2D kernel: updating the 1D fields on the edges of the domain and computing the 2D incoming characteristics for the neighbor faces.
However, the size of data and the type of operations are different: there are $N^2$ 1D fields to update and $2N$ 2D incoming characteristics ($N$ for each neighbor face) to compute by performing matrix-vector products.

For this update kernel, we have modified the parallelization strategy in order to reduce the use of private memory storage and to perform matrix-vector products in parallel.
In the 1D kernel, $N$ threads deals with the tasks associated to one node.
Each of these threads is dedicated to $N$ sets of 1D fields and two 2D incoming characteristics.
In algorithm \ref{algo:updateKernel1}, there are therefore three inner parallel loops (over the elements, the sets and the nodes).
To preserve coalescing memory transfers, the inner most loop processes the nodes since they correspond to the finest granularity of storage for all the arrays.

As explained at the end of section \ref{sec:form:final}, the $2N \times 2N$ system to compute the 2D incoming characteristics can simply be done by performing four matrix-vector products with $N \times N$ matrices, where only two vectors and two matrices are different.
The 1D kernel computes the entries of both vectors (stored in a shared array \texttt{vecX}) and performs the matrix-vector products in a parallel way.
Each thread computes one entry of each vector, and computes one output value of each matrix-vector product.

\item
The only task of the 0D update kernel is computing the 1D outgoing characteristics by solving a $3N^2 \times 3N^2$ linear system for each corner, which the inverse matrix is available.
Using symmetry in the inverse matrix, the task consists in performing nine matrix-vector products with $N^2 \times N^2$ matrices, where only three vectors and three matrices are different.
The 0D kernel computes the entries of the three vectors (stored in a shared array \texttt{vecX}) and performs the matrix-vector products in a parallel fashion.
In our implementation, each thread block deals with one corner, and each thread computes one entry of each vector and one output value of each matrix-vector products.

\end{itemize}

\subsection{Validation benchmark}
\label{sec:resu:valid}

In this section, the accuracy of the HABC formulation is studied with an academic benchmark which the solution is known.
We consider the propagation of a spherical wave in the infinite space $\mathbb{R}^3$.
The wave is generated with a source point at position $\vec{x}_\text{s}\in\Omega$ using the Ricker wavelet $s(t)$ defined as
\begin{align}
  s(t) \equiv \big(1-2\pi^2 f_\text{peak}^2(t-t_\text{s})^2\big) e^{-\pi^2 f_\text{peak}^2(t-t_\text{s})^2},
  \label{eqn:validBench:ricker}
\end{align}
where $f_\text{peak}$ is the peak frequency and $t_\text{s}$ is a time offset.
The source point is incorporated in the pressure-velocity system using a Dirac delta in the pressure equation:
\begin{align*}
  \dvp{p}{t} + \rho c^2 \nabla\cdot\vec{u} &= \delta(\vec{x}-\vec{x}_\text{s}) \:S(t),
\end{align*}
where $S(t)$ is the integral of the wavelet,
\begin{align*}
  S(t) = \int_{-\infty}^t s(t')\:dt' = (t-t_\text{s}) e^{-\pi^2 f_\text{peak}^2(t-t_\text{s})^2}.
\end{align*}
The analytic solution is obtained by taking the convolution of the Green function of the wave equation with the Ricker wavelet.
This solution then reads
\begin{subequations}
\begin{align}
  p^\text{ref}(t,\vec{x}) &= \frac{1}{4\pi r} s(t-r/c),
  \label{eqn:validBench:refSol1} \\
  \vec{u}^\text{ref}(t,\vec{x}) &= \frac{\vec{x}-\vec{x}_\text{s}}{4\pi \rho r^2}
    \left(\frac{1}{r} S(t-r/c) + \frac{1}{c} s(t-r/c) \right),
  \label{eqn:validBench:refSol2}
\end{align}
\end{subequations}
with $r=\|\vec{x}-\vec{x}_\text{s}\|^{1/2}$.

\subsubsection*{Setting}

The numerical simulation is performed on the cuboidal domain $\Omega = [-0.5,0.5]^3$ with a mesh composed of 70895 tetrahedra.
Third-degree polynomial basis functions are used (\textit{i.e.}~ $P=3$).
The time step $\Delta t$ is chosen according to
\begin{align}
  \Delta t = \mathrm{max}_k \frac{1}{c_k(P+1)^2 F_{\text{scale},k}},
  \label{eqn:validBench:timeStep}
\end{align}
where $c_k$ is the wave velocity and $F_{\text{scale},k}$ is the maximum ratio of surface to volume Jacobian of the $k^\text{th}$ mesh cell (see \textit{e.g.}~ \cite{warburton2003constants,chan2016gpu,hesthaven2007nodal}).
For this benchmark, we use dimensionless physical parameters $\rho_k$ and $c_k$ set to $1$.
The duration of the simulation is $t_\text{final}=5$ and the peak of the Ricker wavelet is generated at $t_\text{s}=0.5$ with the peak frequency $f_\text{peak}=2.5$.

We compare the numerical solution obtained with approximate boundary treatments (basic ABC and HABC) to the infinite-space reference solution \eqref{eqn:validBench:refSol1}-\eqref{eqn:validBench:refSol2}.
The difference is quantified with the relative $L^2$-error on the domain $\Omega$ defined as
\begin{align}
  \mathrm{Error}(t)
  = \sqrt{\frac{ \displaystyle
             \int_{\Omega} \left(
                \frac{1}{2\rho c^2} \left(p^{\text{ref}}(t,\vec{x})
                                         -p^{\text{num}}(t,\vec{x})\right)^2
              + \frac{\rho}{2} \left\|\vec{u}^{\text{ref}}(t,\vec{x})
                                     -\vec{u}^{\text{num}}(t,\vec{x})\right\|^2\right) d\vec{x}}{
           \text{Total energy generated by the source}}}.
  \label{eqn:validBench:error}
\end{align}
The total energy generated by the source (used in the denominator) is computed by performing the simulation with the homogeneous Dirichlet condition $p=0$ on $\partial\Omega$, which does not allow outgoing energy flux.
The energy then is evaluated in $\Omega$ when it reaches a constant value, that is when the wavelet is totally generated.
Both this energy and the error are evaluated using a numerical integration with a quadrature rule that is exact for seventh-degree polynomials.
Note that error \eqref{eqn:validBench:error} measures both modeling errors due to the approximate boundary treatment and numerical errors due to the discretization of the problem.

\subsubsection*{Results}

Figure \ref{fig:convBench:fullHabc} shows the time-evolution of the error when the source is placed at the center of the domain ($\vec{x}_\text{s}=(0,0,0)$) and when it is slightly shifted ($\vec{x}_\text{s}=(0.2,0.1,0)$).
The generated spherical wavefront propagates in the domain and the peak reaches the boundary at $t=1$ and $t=0.8$, respectively.
At this instant, the error reaches $\sim 10^{-2}$ with all the boundary treatments and both sources.
In this first period, the total error is dominated by the numerical error.
After, the behavior depends on the boundary treatment.

\begin{figure}[!t]
\centering
\begin{tabular}{cc}
\begin{subfigure}[b]{70mm}
  \centering
  \caption{Source centered $\vec{x}_\text{s}=(0,0,0)$}
  \includegraphics[width=68mm]{img/conv_freeSpace_centPulse.pdf}
  \label{fig:convBench:fullHabc:1}
\end{subfigure}
&
\begin{subfigure}[b]{70mm}
  \centering
  \caption{Source with offset $\vec{x}_\text{s}=(0.2,0.1,0)$}
  \includegraphics[width=68mm]{img/conv_freeSpace_decPulse.pdf}
  \label{fig:convBench:fullHabc:2}
\end{subfigure}
\end{tabular}
\vspace{-0.3cm}
\caption{Time-evolution of the error for the validation benchmark with approximate boundary treatments on all the faces of the domain.
The source is placed at the center of the domain (a) or slightly shifted away from the center (b).
The numbers on the curves indicate the number of additional fields in the boundary treatment ($N=0$ for the basic ABC and $N=2,4,\dots,12$ for HABC).}
\label{fig:convBench:fullHabc}
\end{figure}

With the basic ABC, the error increases until it reaches  $\sim 0.17$ for both sources, which means that approximately $17\%$ of the total energy generated in the domain has been reflected.
The error is clearly dominated by the modeling error: the reflected wavefront propagated in the domain is partially reflected at the boundary.
The error is continuously decreasing as the multiple reflections are absorbed.
With the basic ABC, it can be proved that, for both the continuous model and the numerical scheme, the energy cannot increases once the source is totally generated.
The observed error decay is therefore the expected result.

With the HABC, the error decreases until a minimum is reached between $t=1$ and $t=2.3$.
During this decrease, the error is the same with all the HABCs and is dominated by the numerical error.
The attained minimum depends on the order of the HABC: a larger order provides a smaller error.
We note that this minimum occurs earlier and is larger when the source is shifted (figure \ref{fig:convBench:fullHabc:2}).
Indeed, because the source is closer to a boundary, the reflection occurs earlier, and the amplitudes of both incident and reflected wavefronts are larger.
After the minimum, the error oscillates (for small $N$) or increases (for large $N$) to converge towards the same value for all orders of HABC.
Such phenomenon is well known with Pad\'e-like and Higdon-like boundary conditions (see the numerical results in \textit{e.g.}~ \cite{hagstrom2004new,hagstrom2007local,hagstrom2010radiation}).
It is due to the poor long time error behavior of these conditions, which can be overcome, for instance, with the CRBC \cite{hagstrom2009complete,hagstrom2010radiation}.

In order to validate HABC coupled with a homogeneous boundary condition, we consider a variant of the benchmark where $p=0$ is prescribed on the upper face (\textit{i.e.}~ $z=0.5$), while HABC are used on the other faces.
In this benchmark, a primary wavefront is generated by the source, and a secondary wavefront appears after the reflection of the primary front on the upper boundary.
The reference solution is straightforwardly obtained by using the infinite-space solution \eqref{eqn:validBench:refSol1}-\eqref{eqn:validBench:refSol2} with the method of images.

\begin{figure}[!t]
\centering
\begin{tabular}{cc}
\begin{subfigure}[b]{70mm}
  \centering
  \caption{Source centered $\vec{x}_\text{s}=(0,0,0)$}
  \includegraphics[width=68mm]{img/conv_freeSurf_centPulse.pdf}
  \label{fig:convBench:habcWithDirBC:1}
\end{subfigure}
&
\begin{subfigure}[b]{70mm}
  \centering
  \caption{Source with offset $\vec{x}_\text{s}=(0.2,0.1,0)$}
  \includegraphics[width=68mm]{img/conv_freeSurf_decPulse.pdf}
  \label{fig:convBench:habcWithDirBC:2}
\end{subfigure}
\end{tabular}
\vspace{-0.3cm}
\caption{Time-evolution of the error for the validation benchmark with a free-surface boundary on the upper face of the domain (\textit{i.e.}~ $z=0.5$) and approximate boundary treatments on the other faces.
The source is placed at the center of the domain (a) or slightly shifted away from the center (b).
The numbers on the curves indicate the number of additional fields in the boundary treatment ($N=0$ for the basic ABC and $N=2,4,\dots,12$ for HABC).}
\label{fig:convBench:habcWithDirBC}
\end{figure}

Figure \ref{fig:convBench:habcWithDirBC} shows the time-evolution of the error for the modified benchmark.
The behavior of the error is the same as with the previous benchmark for $t<1$ because the second wavefront has not yet appeared.
In the range $t\in[1,2]$, all the HABC give the same error, which corresponds to the numerical error when both primary and secondary wavefronts are traveling in the domain.
With the basic ABC, the error is clearly dominated by modeling error due the spurious reflection of waves.
The long time error behavior is similar to the previous benchmark.

\subsection{Realistic benchmark}
\label{sec:resu:seam}

In order to test our approach with a more realistic situation, we have built a benchmark based on the \textit{SEAM Phase I} model produced by the \textit{SEG Advanced Modeling Program} \cite{fehler2011seam}.

\subsubsection*{Setting}

The computational domain $\Omega = [0,35~\textrm{km}]\times[0,40~\textrm{km}]\times[0,15~\textrm{km}]$ of the \textit{SEAM Phase I} model is a 3D representation of a deepwater Gulf of Mexico salt domain with a stratigraphy.
The last dimension of the domain corresponds to the vertical direction.
The coordinate $z$ is the depth from the sea level.
The domain is partitioned into an unstructured mesh made of $1,179,989$ tetrahedra.
The mesh has been generated with smaller cells in regions with smaller P-wave velocity in order to accurately represent the spatial oscillations.
Since the physical parameters must be constant over each mesh cell, the mean density $\rho_k$ and the mean P-wave velocity $c_k$ of the SEAM model are taken.
The physical parameters and the mesh are shown on figures \ref{fig:SEAM:model} and \ref{fig:SEAM:mesh}.
A salt body is visible in the middle of the domain.
The upper part of the domain represents the ocean.

\begin{figure}[!t]
\centering
\begin{tabular}{cc}
\begin{subfigure}[b]{8cm}
  \centering
  \caption{Velocity model}
  \includegraphics[width=78mm]{img/seam_modelVelocityB.png}
  \label{fig:SEAM:model:V}
\end{subfigure}
&
\begin{subfigure}[b]{8cm}
  \centering
  \caption{Density model}
  \includegraphics[width=78mm]{img/seam_modelDensityB.png}
  \label{fig:SEAM:model:D}
\end{subfigure}
\end{tabular}
\caption{Density $(a)$ and velocity $(b)$ based on the \textit{SEAM Phase I} model for the realistic benchmark.}
\label{fig:SEAM:model}
\end{figure}

\begin{figure}[!t]
  \centering
  \medskip
  \includegraphics[width=78mm]{img/seam_meshVelocity.png}
\caption{Mesh based on the \textit{SEAM Phase I} model for the realistic benchmark.}
\label{fig:SEAM:mesh}
\end{figure}

A wavefront is generated in the upper part of the domain by using a point source with the Ricker wavelet.
The position of the source is $\vec{x}_\text{s} = (10~\textrm{km},10~\textrm{km},1.5~\textrm{km})$, the time offset is $t_\text{s}=2~\textrm{s}$ and the peak frequency is $f_\text{peak} = 2~\textrm{Hz}$.
The free-surface boundary condition is prescribed on the upper border of the domain (at $z=0$), while an HABC is used on the lateral and bottom borders.
For the boundary procedure, a surface mesh with $23,892$ triangles covers the lateral and bottom faces.
The line mesh composed of $332$ lines is used for the lateral and bottom edges.
Only the $4$ corners at the bottom are considered in the procedure.
The simulations have been performed with third-degree polynomial basis functions, for the duration $15~\textrm{s}$ with the global time-step $\Delta t = 0.124723~\textrm{ms}$, which has been computed using equation \eqref{eqn:validBench:timeStep}.

We solver this benchmark with one single Nvidia K40 GPU, which constrains the size of problems that can be solved.
For realistic applications, the solver can be improved by using strategies for parallel computing on GPU clusters and multi-rate time stepping \cite{modave2015nodal}.
Using these strategies dramatically accelerates the computation, and allows for larger and more refined meshes which support higher frequencies.

\subsubsection*{Qualitative comparison}

We have performed simulations using the basic ABC and HABCs with $N=3$ and $N=6$.
Figure \ref{fig:SEAG:solution} shows snapshots of the solution at different instants for the HABC with $N=3$.
At $t=3.75~\textrm{s}$, we can see the primary wavefront, as well as the secondary wavefront generated after reflection on the free surface of the domain.
The source, represented with a yellow bullet, is in the upper part of the domain, which corresponds to an ocean.
In the remainder of the simulation, both wavefronts are propagated in the ocean and the subsurface, and multiple reflections appear due to geological structures.
Waves travel significantly faster in the subsurface than in the ocean.

\begin{figure}[!p]
\begin{center}
  \includegraphics[width=160mm]{img/seam_sol3_5to5_fin.png}
\end{center}
\vspace{-0.4cm}
\caption{Snapshots of the pressure wavefield at different instants for the realistic benchmark.
  The free-surface boundary condition $(p=0)$ is used for the upper face of the domain,
  while the HABC with $N=3$ is used for the other faces.
  The location of the source point is represented with a yellow bullet on all the figures.}
\label{fig:SEAG:solution}
\end{figure}

\begin{figure}[!tbp]
\centering
\bigskip
\begin{tabular}{c@{\:}c}
\begin{subfigure}[b]{70mm}
  \centering
  \caption{P-wave velocity}
  \includegraphics[width=70mm]{img/seam_cutA_velocity.png}
  \label{fig:SEAM:rezuA:1}
\end{subfigure}
&
\begin{subfigure}[b]{70mm}
  \centering
  \caption{Density}
  \includegraphics[width=70mm]{img/seam_cutA_rho.png}
  \label{fig:SEAM:rezuA:2}
\end{subfigure}
\vspace{0.25cm} \\
\begin{subfigure}[b]{70mm}
  \centering
  \caption{Pressure field for basic ABC}
  \includegraphics[width=70mm]{img/seam_cutA_HABC0_120000.png}
  \label{fig:SEAM:rezuA:3}
\end{subfigure}
&
\begin{subfigure}[b]{70mm}
  \centering
  \caption{Pressure field for HABC with $N=3$}
  \includegraphics[width=70mm]{img/seam_cutA_HABC3_120000.png}
  \label{fig:SEAM:rezuA:4}
\end{subfigure}
\end{tabular}
\vspace{0cm}
\caption{Snapshots of the medium properties and the pressure wavefield in an horizontal planar cut of the sea at depth $z=1.5~\textrm{km}$ and time $t=15~\textrm{s}$.
  The velocity and density models are represented in figures (a) and (b), respectively.
  The wavefield is computed using the basic ABC (c) and the HABC with $N=3$ (d).}
\label{fig:SEAM:rezuA}
\end{figure}

\begin{figure}[!tbp]
\centering
\bigskip
\begin{tabular}{c@{\:}c}
\begin{subfigure}[b]{70mm}
  \centering
  \caption{P-wave velocity}
  \includegraphics[width=70mm]{img/seam_cutB_velocity.png}
  \label{fig:SEAM:rezuB:1}
\end{subfigure}
&
\begin{subfigure}[b]{70mm}
  \centering
  \caption{Density}
  \includegraphics[width=70mm]{img/seam_cutB_rho.png}
  \label{fig:SEAM:rezuB:2}
\end{subfigure}
\vspace{0.25cm} \\
\begin{subfigure}[b]{70mm}
  \centering
  \caption{Pressure field for basic ABC}
  \includegraphics[width=70mm]{img/seam_cutB_HABC0_90000.png}
  \label{fig:SEAM:rezuB:3}
\end{subfigure}
&
\begin{subfigure}[b]{70mm}
  \centering
  \caption{Pressure field for HABC with $N=3$}
  \includegraphics[width=70mm]{img/seam_cutB_HABC3_90000.png}
  \label{fig:SEAM:rezuB:4}
\end{subfigure}
\end{tabular}
\vspace{0cm}
\caption{Snapshots of the medium properties and the pressure wavefield in an horizontal planar cut of the ground at depth $z=5~\textrm{km}$ and time $t=11.25~\textrm{s}$.
  The velocity and density models are represented in figures (a) and (b), respectively.
  The wavefield is computed using the basic ABC (c) and the HABC with $N=3$ (d).}
\label{fig:SEAM:rezuB}
\end{figure}

In order to compare the basic ABC with the HABCs, horizontal planar cuts of the pressure wavefield are shown in figures \ref{fig:SEAM:rezuA} and \ref{fig:SEAM:rezuB} at depths $z=1.5~\textrm{km}$ and $z=5~\textrm{km}$.
HABC results are shown only for $N=3$, because the images obtained with $N=6$ are visually nearly identical.

The horizontal cuts shown on figure \ref{fig:SEAM:rezuA} are taken at the ocean level where the source is placed, at the end of the simulation.
On the left border of the pictures, the medium is heterogeneous because the salt body touches the border (figures \ref{fig:SEAM:rezuA:1} and \ref{fig:SEAM:rezuA:2}).
There is therefore a large change of the medium properties that must be represented by the boundary treatment.
Comparing the results obtained with both boundary treatments, we can clearly observe reflections close to the left and lower sides of figure \ref{fig:SEAM:rezuA:3} (basic ABC) that are not on figure \ref{fig:SEAM:rezuA:4} (HABC).
The larger reflections correspond to oblique waves, while the wavefront corresponding to a normal incidence is not reflected by the basic ABC.
Note that the solution obtained with the HABC does not exhibit any incoherent behavior on the left side of the picture, where the medium is heterogeneous.

Figure \ref{fig:SEAM:rezuB} shows horizontal cuts.
These cuts are taken earlier in the simulation in order to observe eventual reflections of the primary wavefront at the boundary.
Again, we observe reflections when the ABC is used, and not with the HABC.

\subsubsection*{Computational performance}

The total runtime of the realistic simulation using the basic ABC is approximately 8h45 on a Nvidia K40 GPU in single precision.
Using the HABCs increases the runtime by 10 and 22 minutes for $N=3$ and $N=6$, respectively.
As shown on table \ref{table:SEAM:perf}, the runtime rises nearly proportionally to the total number of unknowns.

We have optimized the number of elements per thread block $K_\text{blk}$ for each kernel and each boundary treatment (see section \ref{sec:implementation}).
The optimum values are given in table \ref{table:SEAM:optiKblk}.
Optimizing these parameters $K_\text{blk}$ has a significant impact on the efficiency.
As shown on table \ref{table:SEAM:rezuRuntime}, the speedup achieved after optimization of $K_\text{blk}$ kernel by kernel is between 1.25 and 3.95.
The largest speedups are obtained with the 2D kernels and the 3D volume kernel.
The speedup of the complete implementation after optimization is approximately 1.85.

\begin{table}[!tb]
\begin{center}
\small
\begin{tabular}{|l|r|r|r|} \hline
                            & Basic ABC & HABC \footnotesize$(N=3)$ & HABC \footnotesize$(N=6)$ \\ \hline
  Total number of unknowns  & 94,399,120 & 96,573,304 \small(+2.3\%) & 98,795,296 \small(+4.6\%) \\
  Total runtime             & 8h44 50' & 8h54 45' \small(+1.9\%) & 9h07 43' \small(+4.4\%) \\ \hline
\end{tabular}
\vspace{-0.25cm}
\end{center}
\caption{Performance statistics for the realistic benchmark when using a basic ABC or a HABC as boundary treatment.}
\label{table:SEAM:perf}
\end{table}

\begin{table}[!tb]
\begin{center}
\small
\begin{tabular}{|c|c|c|c|} \hline
  Kernel & Basic ABC & HABC \footnotesize$(N=3)$ & HABC \footnotesize$(N=6)$ \\ \hline
  1D volume  & & 3 & 3 \\
  2D volume  & & 16 & 13 \\
  3D volume  & 8 & 8 & 8 \\
  1D surface & & 3 & 3 \\
  2D surface & & 8 & 5 \\
  3D surface & 2 & 2 & 2 \\
  1D update  & & 8 & 6 \\
  2D update  & & 10 & 8 \\
  3D update  & 3 & 3 & 3 \\ \hline
\end{tabular}
\vspace{-0.25cm}
\end{center}
\caption{Optimum value of the thread block size $K_\text{blk}$ for each kernel and each variant of the realistic benchmark.}
\label{table:SEAM:optiKblk}
\end{table}

\begin{table}[!tb]
\begin{center}
{\small
\begin{tabular}{|c|r|r|r|r|r|r|r|r|r|} \hline
  Kernel &
  \multicolumn{3}{c|}{Basic ABC} &
  \multicolumn{3}{c|}{HABC \footnotesize$(N=3)$} &
  \multicolumn{3}{c|}{HABC \footnotesize$(N=6)$} \\ \hline
             & \footnotesize $K_\text{blk}=1$
             & \footnotesize Opti $K_\text{blk}$
             & \footnotesize SpUp
             & \footnotesize $K_\text{blk}=1$
             & \footnotesize Opti $K_\text{blk}$
             & \footnotesize SpUp
             & \footnotesize $K_\text{blk}=1$
             & \footnotesize Opti $K_\text{blk}$
             & \footnotesize SpUp \\ \hline
  1D volume  &            &          &      & $12.9~\:\mu\textrm{s}$ & $9.52~\:\mu\textrm{s}$ & 1.36 & $27.4~\:\mu\textrm{s}$ & $19.9~\:\mu\textrm{s}$ & 1.38 \\
  2D volume  &            &          &      &  $533~\:\mu\textrm{s}$ &  $135~\:\mu\textrm{s}$ & 3.95 &  $783~\:\mu\textrm{s}$ &  $305~\:\mu\textrm{s}$ & 2.57 \\
  3D volume  & $33.2~\textrm{ms}$ & $13.1~\textrm{ms}$ & 2.53 &  $33.5~\textrm{ms}$ &  $13.1~\textrm{ms}$ & 2.56 &  $33.4~\textrm{ms}$ &  $13.1~\textrm{ms}$ & 2.55 \\
  1D surface &            &          &      & $24.6~\:\mu\textrm{s}$ & $19.3~\:\mu\textrm{s}$ & 1.27 & $68.5~\:\mu\textrm{s}$ & $55.0~\:\mu\textrm{s}$ & 1.25 \\
  2D surface &            &          &      &  $870~\:\mu\textrm{s}$ &  $309~\:\mu\textrm{s}$ & 2.82 &  $1.42~\textrm{ms}$ &  $721~\:\mu\textrm{s}$ & 1.97 \\
  3D surface & $43.5~\textrm{ms}$ & $24.8~\textrm{ms}$ & 1.75 &  $44.0~\textrm{ms}$ &  $25.1~\textrm{ms}$ & 1.75 &  $44.0~\textrm{ms}$ &  $25.1~\textrm{ms}$ & 1.75 \\
  0D update  &            &          &      & $12.7~\:\mu\textrm{s}$ & $12.7~\:\mu\textrm{s}$ &      & $22.7~\:\mu\textrm{s}$ & $22.7~\:\mu\textrm{s}$ &      \\
  1D update  &            &          &      & $29.0~\:\mu\textrm{s}$ & $19.9~\:\mu\textrm{s}$ & 1.46 & $54.4~\:\mu\textrm{s}$ & $37.3~\:\mu\textrm{s}$ & 1.46 \\
  2D update  &            &          &      &  $958~\:\mu\textrm{s}$ &  $526~\:\mu\textrm{s}$ & 1.82 &  $1.53~\textrm{ms}$ &  $875~\:\mu\textrm{s}$ & 1.75 \\
  3D update  & $19.1~\textrm{ms}$ & $14.2~\textrm{ms}$ & 1.35 &  $19.2~\textrm{ms}$ &  $14.2~\textrm{ms}$ & 1.35 &  $19.2~\textrm{ms}$ &  $14.2~\textrm{ms}$ & 1.35 \\ \hline
\end{tabular}
}
\end{center}
\vspace{-0.25cm}
\caption{Average runtime of each kernel per call with and without optimized $K_\text{blk}$, and speedup (SpUp), for each variant of the realistic benchmark.
The 0D update kernel has no tuning parameter.}
\label{table:SEAM:rezuRuntime}
\end{table}

\section*{Acknowledgements}

This work was funded by a grant from TOTAL E\&P Research and Technology USA.
The authors thank TOTAL for permission to publish.
Axel Modave was partially supported by an excellence grant from Wallonie-Bruxelles International (WBI), and was a Postdoctoral Researcher on leave with the F.R.S-FNRS.
The first author thanks Thomas Hagstrom for helpful and informative discussions.

\section{Introduction}
\label{sec:intro}

\input{0_intro}

\section{Non-reflective boundary treatment for cuboidal domain}
\label{sec:HABCform}

\input{1a_conditionPlane}
\input{1b_compatibility}
\input{1c_formulationPV}

\section{Numerical scheme and computational implementation}
\label{sec:scheme}

\input{2a_schemeProcedure}
\input{2b_schemeDG}
\input{2c_schemeImplementation}

\section{Numerical results}
\label{sec:results}

\input{3a_resultsValid}
\input{3b_resultsSEAM}

\newpage
\section{Conclusion}
\label{sec:conclusion}

\input{4_conclusion}

\small
\setlength{\bibsep}{0pt plus 0ex}
\bibliographystyle{abbrvnat}

redpajama_set_name: RedPajamaArXiv


text:
Bob Segar
Shakin Stevens
TODAY IN HISTORY – AUGUST 1st
1939 - Band leader Glenn Miller recorded In the Mood which would later become his theme tune.
1942 - Born this day, Jerry Garcia [Jerome John], in San Francisco, guitarist, banjo player, lyricist, tie designer, lead singer of the Grateful Dead. 1970 UK No.69 and US No.127 album, Workingman's Dead. Died 9 August 1995.
1951 - Born this day, Tommy Bolin, guitarist, Zephyr, The James Gang, joined Deep Purple in 1975. Died 4 December 1976 of a heroin overdose
1953 - Born this day, Robert Cray, in Columbus, Georgia, blues singer, songwriter, guitarist (1987 Grammy), 1987 UK No.50 single Right Next Door, Because Of Me.
1959 - Born this day, Joe Elliott, vocals, Def Leppard, 1987 UK No.6 single Animal, 1987 world wide No.1 album Hysteria, 1988 US No.1 single Love Bites.
1964 - Rod Stewart's group Steam Packet opened two shows for the Rolling Stones at the London Palladium.
1964 - The Beatles scored their fifth US No.1 single in seven months when A Hard Day's Night went to the top of the charts. The group had now spent 17 weeks at the No.1 position in this year.
1964 - Born this day, Adam Duritz, vocals, piano, Counting Crows, 1994 UK No.28 single Mr Jones, 1996 US No.1 album Recovering The Satellites.
1970 - The film 'Performance' featuring Mick Jagger in his acting debut, had its UK premiere in London. The British crime drama film directed by Donald Cammell and Nicolas Roeg also featured James Fox.
1971 - The Concert for Bangladesh was held at Madison Square Garden in New York City. George Harrison, Ringo Starr, Eric Clapton, Bob Dylan, Leon Russell, Ravi Shankar and Billy Preston and members from Badfinger performed. A multi-record set commemorating the event was a super sales success. Together, the concert and the album raised over $11 million to help the starving millions of Bangladesh.
1971 - As a summer replacement, CBS launched The Sonny and Cher Comedy Hour on this date. The unlikely hosting couple, Sonny and Cher Bono, were fresh and funny, and they were embraced by the American public. The show was a huge hit, and was picked up as a permanent CBS fixture. The musical-variety show, with the glamorous, exotic Cher and the farcical, diminutive Sonny poking fun affectionately at each other, lasted many years until the Bonos announced their marriage was over and they were cancelling the show.
1978 - Born this day, Dhani Harrison, George Harrison's first child.
1981 – At 12.01am "MTV" broadcast for the first time. First songs played were "Video Killed the Radio Star" – The Buggles followed by "You Better Run" – Pat Benatar. Most featured artist of the first day was Rod Stewart, first video to play twice was "You Better You Bet" – The Who, Split Enz featured 6 times during the day with 4 different songs including "History Never Repeats", "One Step Ahead", "I Hope I Never" & "I Got You".
1981 - Australian singer Rick Springfield started a two week run at No.1 on the US singles chart with Jessie's Girl, it made No.43 hit in the UK. Rick had played Noah Drake in the TV show General Hospital.
1981 - Shakin' Stevens had his second UK No.1 single with his version of Green Door. It was originally a hit for Frankie Vaughan, a 1956 UK No.2 single.
1987 - Dave Stewart married ex Bananarama member SiobharnFahey.
1987 - Bob Seger scored his first US No.1 single with Shakedown, taken from the film Beverly Hills Cop II.
1987 - Los Lobos were at No.1 on the UK singles chart with their version of the Ritchie Valens song La Bamba. It also made No.1 in the US. The song was the title track from the film based on Ritchie Valens who died in the same plane crash the killed Buddy Holly.
1994 - Lisa Marie Presley confirmed rumours that she had married pop star Michael Jackson on 26 May, in the Dominican Republic. The couple divorced less than two years later.
1999 - Ronan Keating scored his first UK No.1 solo single with When You Say Nothing At All.
2000 - Madonna's forthcoming single Music had it's release date brought forward by two weeks after the track was made available as an illegal MP3 file on the Internet.
2002 - A new book 'Show the Girl the Door' written by a former tour manager disclosed some strange demands by female acts. It revealed that Shania Twain would travel with a sniffer dog in case of bombs. Jennifer Lopez liked her dressing room to be all white, including carpets flowers and furniture. Cher would have high security rooms for her wigs. Janet Jackson would have a full medical team on standby including a doctor nurse and throat specialist and Britney Spears would demand her favourite Gummie Bear soft sweets.
2007 - John Lennon's "granny" sunglasses were snapped up by a British collector at auction. The sunglasses, from one of the last Beatles concerts, were expected to fetch around £1m, but auction bosses refused to say what the actual figure was. Lennon gave the gold-rimmed glasses to his Japanese interpreter in Tokyo in 1966, and the translator removed the lenses when Lennon died.
2013 - The British government was trying to stop American Idol singer Kelly Clarkson from taking a rare turquoise and gold ring once owned by Jane Austen out of the country. The 2002 winner of the Idol TV show bought the jewellery at auction last year for more than £150,000, but Culture minister Ed Vaizey had put a temporary export bar on it saying he wanted the "national treasure" to be "saved for the nation".

redpajama_set_name: RedPajamaCommonCrawl


Technical research is aimed at solutions and design: translating knowledge into things. It comes to making tangible materials, products or processes. This way, it distinguishes itself from scientific research (science). Within the Centre of Applied Explore of the Pogled u Plavo, lecturers, researchers and students cooperate with practitioners from the industry in research projects and research labs such as the SEFLab, the Innovation Lab, OpenBIM-lab and the Logistics Lab.
Demonstrators and prototypes of future products, like the H2A, the hydrogen-car of the Pogled u Plavo which won the first price in the hydrogenclass of the Shell EcoMarathon Europe with a fuel consumption of one liter per 3042 km.
Articles in professional journals, professional-websites and scientific magazines.
Publications in the publication series of the Pogled u Plavo Centre of Applied Explore Technology.
Contributions on seminars and conferences.


text:
import express from 'express';
import * as request from 'request';
import * as util from 'util';
import * as bparser from 'body-parser';
import { createHmac } from 'crypto';
import * as http from 'http';
import * as https from 'https';
import * as oauth from './oauth';
import * as ssl from './ssl';
import * as company from './company';
import debug from 'debug';

// Debug log
const log = debug('watsonwork-companyinfo-app');

// Look for company entities mentioned in messages and automatically retrieve
// and post information about these companies
export const companyInfo = (appId, token, fruserId, frkey) => (req, res) => {
  // Respond to the Webhook right away, as the response message will
  // be sent asynchronously
  res.status(201).end();

  // Only handle message-created Webhook events, and ignore the app's
  // own messages
  if(req.body.type !== 'message-created' || req.body.userId === appId)
    return;

  log('Got a message %o', req.body);
  const msg = req.body;

  // React to mentions of company names in the message and send company
  // information back to the conversation in the originating space
  if(msg.content)

    // Attempt to recognize mentioned company entities
    company.entities(msg.content,
      msg.spaceId, msg.messageId, msg.time,
      msg.userId, msg.userName,
      fruserId, frkey, (err, entities) => {

        (entities || []).map((entity) => {
          if(entity.score < 60)
            return;

          // Retrieve each company metadata
          company.metadata(entity.id, fruserId, frkey, (err, info) => {
            if(info)
              // Send a message with the company metadata / info
              send(req.body.spaceId,
                util.format(
                  '*Company*\n%s\n' +
                  '*Industries*\n%s\n' +
                  '*Sectors*\n%s\n' +
                  '*Segments*\n%s\n',
                  info.name,
                  info.industries.map((i) => i.name).join(','),
                  info.sectors.map((s) => s.name).join(', '),
                  info.segments.map((s) => s.name).join(', ')),
                token(),
                (err, res) => {
                  if(!err)
                    log('Sent message to space %s', req.body.spaceId);
                });
          });
        });
      });
};

// Send an app message to the conversation in a space
const send = (spaceId, text, tok, cb) => {
  request.post(
    'https://api.watsonwork.ibm.com/v1/spaces/' + spaceId + '/messages', {
      headers: {
        Authorization: 'Bearer ' + tok
      },
      json: true,
      // An App message can specify a color, a title, markdown text and
      // an 'actor' useful to show where the message is coming from
      body: {
        type: 'appMessage',
        version: 1.0,
        annotations: [{
          type: 'generic',
          version: 1.0,

          color: '#6CB7FB',
          text: text
        }]
      }
    }, (err, res) => {
      if(err || res.statusCode !== 201) {
        log('Error sending message %o', err || res.statusCode);
        cb(err || new Error(res.statusCode));
        return;
      }
      log('Send result %d, %o', res.statusCode, res.body);
      cb(null, res.body);
    });
};

// Verify Watson Work request signature
export const verify = (wsecret) => (req, res, buf, encoding) => {
  if(req.get('X-OUTBOUND-TOKEN') !==
    createHmac('sha256', wsecret).update(buf).digest('hex')) {
    log('Invalid request signature');
    const err = new Error('Invalid request signature');
    err.status = 401;
    throw err;
  }
};

// Handle Watson Work Webhook challenge requests
export const challenge = (wsecret) => (req, res, next) => {
  if(req.body.type === 'verification') {
    log('Got Webhook verification challenge %o', req.body);
    const body = JSON.stringify({
      response: req.body.challenge
    });
    res.set('X-OUTBOUND-TOKEN',
      createHmac('sha256', wsecret).update(body).digest('hex'));
    res.type('json').send(body);
    return;
  }
  next();
};

// Create Express Web app
export const webapp = (appId, secret, wsecret, fruserId, frkey, cb) => {
  // Authenticate the app and get an OAuth token
  oauth.run(appId, secret, (err, token) => {
    if(err) {
      cb(err);
      return;
    }

    // Return the Express Web app
    cb(null, express()

      // Configure Express route for the app Webhook
      .post('/companyinfo',

        // Verify Watson Work request signature and parse request body
        bparser.json({
          type: '*/*',
          verify: verify(wsecret)
        }),

        // Handle Watson Work Webhook challenge requests
        challenge(wsecret),

        // Handle Watson Work messages
        companyInfo(appId, token, fruserId, frkey)));
  });
};

// App main entry point
const main = (argv, env, cb) => {
  // Create Express Web app
  webapp(
    env.COMPANYINFO_APP_ID, env.COMPANYINFO_APP_SECRET,
    env.COMPANYINFO_WEBHOOK_SECRET,
    env.COMPANYINFO_FR_USER_ID, env.COMPANYINFO_FR_KEY, (err, app) => {
      if(err) {
        cb(err);
        return;
      }

      if(env.PORT) {
        // In a hosting environment like Bluemix for example, HTTPS is
        // handled by a reverse proxy in front of the app, just listen
        // on the configured HTTP port
        log('HTTP server listening on port %d', env.PORT);
        http.createServer(app).listen(env.PORT, cb);
      }

      else
        // Listen on the configured HTTPS port, default to 443
        ssl.conf(env, (err, conf) => {
          if(err) {
            cb(err);
            return;
          }
          const port = env.SSLPORT || 443;
          log('HTTPS server listening on port %d', port);
          https.createServer(conf, app).listen(port, cb);
        });
    });
};

if (require.main === module)
  main(process.argv, process.env, (err) => {
    if(err) {
      console.log('Error starting app:', err);
      return;
    }
    log('App started');
  });

redpajama_set_name: RedPajamaGithub