∎

11institutetext:  Andrew J. Christlieb 22institutetext: Michigan State
University
Department of Computational Mathematics, Science and Engineering
428 S. Shaw Lane
East Lansing, Michigan 48824, USA
22email<EMAIL_ADDRESS>33institutetext: Pierson T. Guthrey 44institutetext:
Weapons and Complex Integration
Lawrence Livermore National Laboratory
Livermore, CA 94550, USA
Tel.: +337-781-5574
44email<EMAIL_ADDRESS>55institutetext: James A. Rossmanith
66institutetext: Iowa State University
Department of Mathematics
411 Morrill Road
Ames, Iowa 50011, USA
66email<EMAIL_ADDRESS>
# Parallel Scaling of the Regionally-Implicit Discontinuous Galerkin Method
with Quasi-Quadrature-Free Matrix Assembly ††thanks:  This work was performed
under the auspices of the U.S. Department of Energy by Lawrence Livermore
National Laboratory under Contract DE-AC52-07NA27344. Scaling studies were
performed at Institute for Cyber-Enabled Research at Michigan State University
and at Livermore Computing at Lawrence Livermore National Laboratory. This
work was funded in part by ONR High Order Scalable Solvers grant
N0014-19-1-2476, and AFOSR Computational Non-Ideal Plasma Physics grant
FA9550-17-1-0394, and DoE SciDAC TEAMS grant DE-SC0017955. JAR was supported
in part by NSF Grants DMS–1620128 and DMS–2012699.

Andrew J. Christlieb  Pierson T. Guthrey111Corresponding Author James A.
Rossmanith

(Received: date / Accepted: date)

###### Abstract

In this work we investigate the parallel scalability of the numerical method
developed in Guthrey and Rossmanith [The regionally implicit discontinuous
Galerkin method: Improving the stability of DG-FEM, SIAM J. Numer. Anal.
(2019)]. We develop an implementation of the regionally-implicit discontinuous
Galerkin (RIDG) method in DoGPack, which is an open source C++ software
package for discontinuous Galerkin methods. Specifically, we develop and test
a hybrid OpenMP and MPI parallelized implementation of DoGPack with the goal
of exploring the efficiency and scalability of RIDG in comparison to the
popular strong stability-preserving Runge-Kutta discontinuous Galerkin (SSP-
RKDG) method. We demonstrate that RIDG methods are able to hide communication
latency associated with distributed memory parallelism, due to the fact that
almost all of the work involved in the method is highly localized to each
element, producing a localized prediction for each region. We demonstrate the
enhanced efficiency and scalability of the of the RIDG method and compare it
to SSP-RKDG methods and show extensibility to very high order schemes. The
two-dimensional scaling study is performed on machines at the Institute for
Cyber-Enabled Research at Michigan State University, using up to 1440 total
cores on Intel(R) Xeon(R) Gold 6148 CPU @ 2.40GHz CPUs. The three-dimensional
scaling study is performed on Livermore Computing clusters at at Lawrence
Livermore National Laboratory, using up to 28672 total cores on Intel Xeon
CLX-8276L CPUs with Omni-Path interconnects.

###### Keywords:

discontinuous Galerkin hyperbolic conservation laws Courant-Friedrichs-Lewy
condition time-stepping numerical stability strong scaling high performance
computing domain decomposition quadrature-free

###### MSC:

65M12 65M60 65Y20 35L03

††journal: BIT

##  1 Introduction

In high-performance computing, parallel scaling is an important measure of
algorithm efficiency, since it provides a quantifiable measure of how much
computational benefit is gained by moving to more and more computing resources
(i.e., more processors/cores/threads). There are two standard approaches
measure scalability: (1) weak scaling, where the problem size (i.e., the
degrees of freedom) increases with the number of cores such that the problem
size per processor/core/thread is fixed; and (2) strong scaling, where the
number of processors/cores/threads is increased while the problem size is
fixed. For problems that are CPU-bound, i.e., problems where memory is not the
primary concern, but computing the solution on a single processor takes a long
time, the relevant scaling measure is strong scaling.

For the numerical solution of partial differential equations (PDEs), the
standard strategy for strong scaling is to subdivide a problem of size $P$
over $n$ physical compute cores until a minimal problem size per physical
compute core is reached (e.g., see Fischer, Heisey, and Min Fischer ). The
ratio of problem size to number of compute cores is known as the granularity,
$\eta$, and can be used to predict the ratio of MPI-communication and OpenMP
overhead to parallelized work. In the limit of minimal granularity $\eta\to
1$, parallelization overhead is maximized, reducing efficiency.

In this work we are concerned with strong scaling via a hybrid OpenMP openmp
\+ MPI Gabriel2004  approach for a specific class of problems: hyperbolic
conservation laws, with a specific type of spatial discretization:
discontinuous Galerkin (DG) finite element methods (e.g., see Cockburn and Shu
cockshu5 ). The most common technique to time-advance DG spatial
discretizations is the strong-stability-preserving Runge-Kutta (SSP-RK)
scheme. This approach is explicit in time, and thus very efficient in terms of
CPU and memory resources per time-step. However, a drawback is that SSP-RK
schemes applied to high-order DG methods suffer from very small allowable
time-steps, which then results in an overall reduction in efficiency (i.e.,
many time-steps to reach a fixed final time). In order to overcome these small
time-step restrictions, Guthrey and Rossmanith RIDG_paper_2019  developed a
new time-stepping approach for DG that was dubbed the regionally-implicit DG
(RIDG) scheme. The primary goal of the current work is to develop, implement,
and test a parallelized version of the regionally-implicit discontinuous
Galerkin (RIDG) method.

###  1.1 Scalability challenges

The algorithms discussed in this paper are high-order discontinuous Galerkin
(DG) methods with explicit time-stepping for hyperbolic conservation laws.
Explicit time-stepping algorithms for DG methods are usually associated with
small time-steps, typically with Courant-Friedrichs-Lewy (CFL) numbers much
less than $1$ (and indeed inversely proportional to the method order of
accuracy). The benefit of such schemes is that they have a very low
computational cost associated with the (typically) nearest neighbor
stenciling. Thus, there are two challenges associated with attempts to scale
such schemes to multiple processors. First, as we divide the computational
domain into several pieces, each piece must synchronize with its neighbors
across the domain-decomposition pseudo-boundaries, a process also known as
synchronization of the halo regions Consortium2017 . This must occur at least
once per time-step (typically several times for multi-stage methods), and thus
must be performed many times over the course of forming a solution for some
fixed final time. Secondly, for each time-step all processes must synchronize
and compute the global maximum wave speed, so that each process is using the
same time-step size given by the CFL-restriction. This synchronization takes
the form of an all-to-all reduction. This can be potentially avoided by using
local-time stepping, but we do not discuss this in this paper.

In the face of the challenges, discontinuous Galerkin methods enjoy enhanced
scalability properties due to the nearest-neighbor stenciling typical of such
methods and the highly element-localized compute intensity. Adaptive mesh
refinement has been shown to work very well when applied to DG methods at
scale Dumbser2013a . The arbitrary high order schemes using derivatives (ADER-
DG) are predictor-corrector schemes that avoid repeated communication for each
time-step by being fully discrete Dumbser2018 . That is, there is a space-time
basis underlying the time-step updates as opposed to a method-of-lines update
of a purely spatial basis. This approach shows improved MPI scaling efficiency
compared to RKDG schemes, resulting in up to $23\%$ smaller $L_{2}$ errors
with $1-2\%$ longer runtimes. One goal of this paper is to perform a similar
comparison of RIDG and RKDG methods in this paper.

For discontinuous Galerkin methods, the granularity $\eta$ is usually
discussed in terms of the total number of elements $P$, as opposed to the
total number of problem degrees of freedom. This is because domain-
decomposition is not typically implemented below the level of the elements.
That is, individual elements are not typically broken up as any part of a
domain decomposition strategy.

Domain decomposition for strong scaling such a problem involves dividing the
physical domain into a number of sub-domains equal to the number of MPI tasks
or compute nodes, each with a number of physical compute cores. This gives
rise to domain decomposition pseudo-boundaries, in addition to the physical
boundaries. In this paper we consider periodic boundary conditions, where
physical boundaries become pseudo-boundaries, as they separate problem
subdomains.

As discussed above, dynamic explicit time-stepping relies on an all-to-all
reduction to ensure that each time-step taken uses a global estimate for the
maximum wave speed in the computational domain. That is, each subdomain
assigned to an MPI task must compute its “task-local” maximum wave speed, and
a “task-global” maximum wave speed must be computed via an all-to-all
reduction to ensure that the Courant-Friedrichs-Lewy (CFL) condition is
satisfied by all of the sub-domains. Since this procedure merely involves
comparing wave speeds to find the maximum wave speed among all processes, the
communication latency of this procedure is dominated by the total number of
processes as $\log_{2}P$, not the granularity $\eta$.

Every subdomain must perform halo-region communications with all subdomains
with which it shares a pseudo-boundary. This implies that periodic boundary
conditions are more difficult to scale than Dirichlet or outflow conditions,
since they inherently involve extra communications, specifically the wrap-
around communications (e.g., the east side of the mesh communicates with the
west side of the mesh). As the granularity increases, the time spent copying
data to MPI buffers relative to the work in the rest of the problem increases,
potentially degrading efficiency. However, this is a memory operation and thus
can be relatively quick, even at low granularity.

Explicit time-stepping methods usually exhibit restrictive CFL conditions
(even more-so with high-order DG spatial discretizations), which increase the
total number of time steps, which in turn increases the total amount of
communications per simulation time. Thus, the scalability of an explicit time-
stepping scheme depends on the scheme’s ability to overcome communication
overheads and hide latency as $\eta\to 1$.

###  1.2 Scope of this work

The purpose of the current work is to develop, implement, and study an OpenMP
+ MPI parallelization of the regionally implicit discontinuous Galerkin scheme
(RIDG) Guthrey2017 ; RIDG_paper_2019 . In particular, we are concerned with
comparing it to parallel implementations of the most commonly used variant of
the discontinuous Galerkin method: the strong-stability-preserving Runge-Kutta
DG (SSP-RKDG) scheme cockshu5 ; article:GoShu98 ; gottliebShuTadmor01 . In
section 2 we will briefly review the regionally-implicit method, and in
particular, the two key parts of the time-stepping strategy: (1) regionally-
implicit prediction and (2) explicit correction. In section 3 we develop a
quasi-quadrature-free approach for the Jacobian matrix assembly that is
required in the regionally-implicit prediction step. Next we briefly define
the terminology needed to quantify the accuracy and efficiency of our parallel
implementations in section 4. A brief discussion and explanation of the
proposed parallel implementation is provided in section 5. In section 6.1,
section 6.2, and section 6.3 we report our findings for examples in one, two,
and three spatial dimensions, respectively. The two-dimensional scaling study
presented in section 6.2 is performed on the Institute for Cyber-Enabled
Research at Michigan State University’s Intel(R) Xeon(R) Gold 6148 CPU @
2.40GHz CPUs, using up to 1440 total cores. The three-dimensional scaling
study presented in section 6.3 is performed on Livermore Computing clusters at
at Lawrence Livermore National Laboratory, using up to 28672 total cores on
Intel Xeon CLX-8276L CPUs with Omni-Path interconnects. In these sections we
demonstrate the benefits of using RIDG over SSP-RKDG, especially as the number
of processors increase and as the order of the method increases. Finally, we
conclude our findings in section 7.

##  2 Regionally-implicit discontinuous Galerkin (RIDG) method

$y$$x$${\underline{\Psi}}^{T}\,{\underline{\underline{\widehat{W}}}}_{\,i-1\,j-1}$${\underline{\Psi}}^{T}\,{\underline{\underline{\widehat{W}}}}_{\,i\,j-1}$${\underline{\Psi}}^{T}\,{\underline{\underline{\widehat{W}}}}_{\,i+1\,j-1}$${\underline{\Psi}}^{T}\,{\underline{\underline{\widehat{W}}}}_{\,i-1\,j}$${\underline{\Psi}}^{T}\,{\underline{\underline{W}}}_{ij}$${\underline{\Psi}}^{T}\,{\underline{\underline{\widehat{W}}}}_{\,i+1\,j}$${\underline{\Psi}}^{T}\,{\underline{\underline{\widehat{W}}}}_{\,i-1\,j+1}$${\underline{\Psi}}^{T}\,{\underline{\underline{\widehat{W}}}}_{\,i\,j+1}$${\underline{\Psi}}^{T}\,{\underline{\underline{\widehat{W}}}}_{\,i+1\,j+1}$interior
fluxproper upwind flux Figure 1: Shown is the 2D Cartesian region
${\mathcal{R}}_{i}$ over which the RIDG prediction step is carried out. In the
RIDG prediction step, all of the states, excepting only the one belonging to
the middle element, are only temporary variables and will be discarded once
the predicted solution in element $ij$ has been computed – to make note of
this we place hats over the temporary variables.

The regionally implicit discontinuous Galerkin (RIDG) method was developed in
recent work by Guthrey and Rossmanith RIDG_paper_2019 . The goal of that work
was to improve the linear stability of the the ADER-DG scheme
article:Dumbser2006 ; article:GasDumHinMun2011 ; article:Zanotti2015 . In
particular, the origin of the RIDG method is tied to the method of Gassner et
al. article:GasDumHinMun2011 , where it was shown that Lax-Wendroff
discontinuous Galerkin Qiu2005a  can be formulated as a predictor-corrector
method. The predictor is a local version of a spacetime DG method
article:KlaVegVen2006 ; article:Sudirham2006  (i.e., the predictor is
something like a block-Jacobi update for a fully implicit spacetime DG
method), and the corrector is an explicit method that uses the spacetime
reconstructed solution from the predictor step. The drawback of the locally-
implicit approach is that the resulting ADER-DG scheme suffers from small
maximum allowable time-steps that decrease with increasing order of accuracy.
The RIDG method instead uses a regionally-implicit prediction step, which has
the benefit of resulting in a DG method with maximum allowable time-steps that
are much closer to the optimal limit; and, in particular, the maximum allowed
time-step does not degrade as the order of accuracy is increased. The drawback
of RIDG is that per time-step it is more computational expensive than standard
ADER-DG; however, this penalty is more than offset by the significantly
increased stable time-step, especially for very high-order schemes. In this
section we briefly review the key concepts of the RIDG approach. For more
details we refer the reader to RIDG_paper_2019 .

###  2.1 General setup

Consider hyperbolic conservation laws of the form

${\underline{q}}_{,t}+{\underline{\nabla}}\cdot{\underline{\underline{F}}}\left({\underline{q}}\right)={\underline{0}},$
(1)

where
${\underline{q}}\left(t,{\underline{x}}\right):\mathbb{R}^{+}\times\mathbb{R}^{M_{\text{dim}}}\mapsto\mathbb{R}^{M_{\text{eqn}}}$
is the vector of conserved variables,
${\underline{\underline{F}}}\left({\underline{q}}\right):\mathbb{R}^{M_{\text{eqn}}}\mapsto\mathbb{R}^{{M_{\text{eqn}}}\times{M_{\text{dim}}}}$
is the flux function, ${M_{\text{dim}}}$ is the number of spatial dimensions,
and ${M_{\text{eqn}}}$ is the number of conserved variables. We define
$\Omega\subset\mathbb{R}^{M_{\text{dim}}}$ to be a polygonal domain with
boundary $\partial\Omega$, and discretize $\Omega$ using a finite set of non-
overlapping elements, ${\mathcal{T}}_{i}$, such that
$\cup_{i=1}^{M_{\text{elem}}}{\mathcal{T}}_{i}=\Omega$, where
${M_{\text{elem}}}$ is the total number of elements. Let
${\mathbb{Q}}\left({M_{\text{deg}}},{M_{\text{dim}}}\right)$ denote the set of
polynomials from $\mathbb{R}^{M_{\text{dim}}}$ to $\mathbb{R}$ with maximal
polynomial degree ${M_{\text{deg}}}$. On the mesh of ${M_{\text{elem}}}$
elements we define the broken finite element space:

${\mathcal{W}}^{h}:=\left\\{{\underline{w}}^{h}\in\left[L^{\infty}(\Omega)\right]^{{M_{\text{eqn}}}}:\,{\underline{w}}^{h}\bigl{|}_{{\mathcal{T}}_{i}}\in\left[{\mathbb{Q}}\left({M_{\text{deg}}},{M_{\text{dim}}}\right)\right]^{{M_{\text{eqn}}}}\,\,\forall{\mathcal{T}}_{i}\right\\},$

where $h$ is the grid spacing.

Let $\Phi_{k}\left({\underline{x}}\right)$ for $k=1,\ldots,M_{\text{basis}}$
be a basis that spans
${\mathbb{Q}}\left({M_{\text{deg}}},{M_{\text{dim}}}\right)$ over
${\mathcal{T}}_{i}$ and is orthonormal:

$\frac{1}{|{\mathcal{T}}_{i}|}\int_{{\mathcal{T}}_{i}}\Phi_{k}\left({\underline{x}}\right)\Phi_{\ell}\left({\underline{x}}\right)\,d{\underline{x}}=\delta_{k\ell},$

where $\delta_{k\ell}$ is the Kronecker delta and $|{\mathcal{T}}_{i}|$ is the
volume of element ${\mathcal{T}}_{i}$.

In order to setup the RIDG method we assume that on each element the solution
is of the following form:

${\underline{q}}^{h}\left(t^{n},{\underline{x}}\right)\Bigl{|}_{{\mathcal{T}}_{i}}={\underline{\Phi}}\left({\underline{x}}\right)^{T}{\underline{\underline{Q}}}_{i}=\sum_{\ell=1}^{M_{\text{basis}}}{\underline{Q}}_{i}^{\ell}\left(t^{n}\right)\,\Phi_{\ell}\left({\underline{x}}\right),$
(2)

where ${\underline{Q}}_{i}^{\ell}(t^{n})$ are the unknown degrees of freedom
at time level $t^{n}$. Using this ansatz, multiplying eq. 1 by $\varphi_{k}$,
integrating over the element ${\mathcal{T}}_{i}$, using integration-by-parts
in space, and integrating over the time slab $\left[t^{n},t^{n+1}\right]$,
results in the following equation:

$\begin{split}{\underline{Q}}_{i}^{k}\left(t^{n+1}\right)&={\underline{Q}}_{i}^{k}\left(t^{n}\right)+\frac{1}{|{\mathcal{T}}_{i}|}\int_{t^{n}}^{t^{n+1}}\int_{{\mathcal{T}}_{i}}{\underline{\nabla}}\Phi_{k}\cdot{\underline{\underline{F}}}\left(\,{\underline{q}}^{h}\right)\,d{\underline{x}}\,dt\\\
&-\frac{1}{|{\mathcal{T}}_{i}|}\int_{t^{n}}^{t^{n+1}}\oint_{\partial{\mathcal{T}}_{i}}\Phi_{k}\,{\underline{{\mathcal{F}}}}\left({\underline{q}}^{h}_{+},{\underline{q}}^{h}_{-};{\underline{n}}\right)\,d{\underline{s}}\,dt,\end{split}$
(3)

where ${\underline{n}}$ is an outward-pointing normal vector to
$\partial{\mathcal{T}}_{i}$, ${\underline{q}}^{h}_{+}$ and
${\underline{q}}^{h}_{-}$ are the states on either side of the boundary
$\partial{\mathcal{T}}_{i}$, and ${\underline{{\mathcal{F}}}}$ is some
appropriate numerical flux. Equation eq. 3 has the look and feel of a fully-
discrete numerical method; however, the above volume and surface integral
cannot be evaluated directly since they require knowledge of the solution over
the entire time-slab: $\left[t^{n},t^{n+1}\right]$.

The idea of schemes such as ADER-DG article:Dumbser2006 ;
article:GasDumHinMun2011 ; article:Zanotti2015 , its close cousin Lax-Wendroff
DG Qiu2005a , and indeed RIDG RIDG_paper_2019 , is to form a prediction for
solution ${\underline{q}}^{h}$ over the time-slab:
$\left[t^{n},t^{n+1}\right]$, and then to insert this prediction into equation
eq. 3 to obtain a fully discrete update formula. The main difference between
ADER-DG, Lax-Wendroff DG, and RIDG is only in how the prediction step is
formed. We briefly review the details for RIDG below.

###  2.2 RIDG prediction step

In order to describe the prediction step, we define spacetime elements and
regions:

${\mathcal{S}}_{i}:=\left[t^{n},\,t^{n+1}\right]\times{\mathcal{T}}_{i}\qquad\text{and}\qquad{\mathcal{R}}_{i}:=\bigcup_{j:{\mathcal{T}}_{i}\cap{\mathcal{T}}_{j}\neq\emptyset}{\mathcal{S}}_{j},$
(4)

respectively. Note that definition eq. 4 includes all of the elements
${\mathcal{T}}_{j}$ that share either a face or a vertex with element
${\mathcal{T}}_{i}$ – we refer to all of the elements in ${\mathcal{R}}_{i}$
as vertex-neighbors of element ${\mathcal{T}}_{i}$; a 2D version is depicted
in fig. 1 . Also note that on a Cartesian mesh in ${M_{\text{dim}}}$
dimensions, region ${\mathcal{R}}_{i}$ contains exactly $3^{{M_{\text{dim}}}}$
spacetime elements. On each spacetime element we write the predicted solution
as

${\underline{q}}\left(t,{\underline{x}}\right)\Bigl{|}_{{\mathcal{S}}_{i}}\approx
w_{i}:={\underline{\Psi}}^{T}{\underline{\underline{W}}}_{i},$ (5)

where
${\underline{\underline{W}}}_{i}\in\mathbb{R}^{{M_{\text{P}}}\times{M_{\text{eqn}}}}$,
${\underline{\Psi}}\in\mathbb{R}^{{M_{\text{P}}}}$,

${M_{\text{P}}}:=\left({M_{\text{deg}}}+1\right)^{{M_{\text{dim}}}+1},$

and ${\underline{\Psi}}$ is an orthonormal basis on the spacetime element:

$\frac{1}{\bigl{|}{\mathcal{S}}_{i}\bigr{|}}\int_{S_{i}}{\underline{\Psi}}\,{\underline{\Psi}}^{T}\,d\tau\,d\xi={\underline{\underline{\mathbb{I}}}}\in\mathbb{R}^{{M_{\text{P}}}\times{M_{\text{P}}}}.$

In order to arrive at an algebraic system of equations to solve for the
unknown spacetime coefficients in region ${\mathcal{R}}_{i}$, we multiply eq.
1 by test functions ${\underline{\Psi}}$ and then integrate over each
${\mathcal{S}}_{j}\in{\mathcal{R}}_{i}$. We apply integration-by-parts in both
space and time. The integration-by-parts in time connects the current time-
slab $[t^{n},t^{n+1}]$ to the solution in the previous time-slab (i.e., this
enforces causality). In space we treat element boundaries differently
depending on whether or not those boundaries are (1) strictly internal to
region ${\mathcal{R}}_{i}$ or (2) on the boundary of region
${\mathcal{R}}_{i}$. For strictly internal boundaries we evaluate the
resulting surface integral through standard upwind fluxes (i.e., Rusanov
fluxes) that utilize the solution on both sides of the interface. However, on
the boundary of ${\mathcal{R}}_{i}$ we evaluate the resulting surface
integrals only using the trace of the solution that is internal to
${\mathcal{R}}_{i}$. The effect of these choices is to insulate region
${\mathcal{R}}_{i}$ from all elements exterior to ${\mathcal{R}}_{i}$. Solving
the algebraic equations in region ${\mathcal{R}}_{i}$ resulting from the above
described integration-by-parts, gives us the regionally-implicit prediction.
This setup is depicted in the 2D setting in fig. 1 .

Algebraically, the system that must be solved in each region has the following
form for each $j:{\mathcal{T}}_{j}\in{\mathcal{R}}_{i}$ and for each
$k=1,2,\ldots,{M_{\text{P}}}$:

$\displaystyle\begin{split}{\underline{R}}_{jk}:=&\underbrace{\left[\int_{{\mathcal{T}}_{j}}\left(\Psi_{k}{\underline{\Psi}}^{T}\right)\Bigl{|}_{t^{n+1}}\,d{\underline{x}}-\int_{t^{n}}^{t^{n+1}}\int_{{\mathcal{T}}_{j}}\Psi_{k,t}\,{\underline{\Psi}}^{T}\,d{\underline{x}}\,dt\right]\,{\underline{\underline{W}}}_{j}}_{\text{(time
term)}}\\\
-&\underbrace{\int_{t^{n}}^{t^{n+1}}\int_{{\mathcal{T}}_{j}}{\underline{\nabla}}\Psi_{k}\cdot{\underline{\underline{F}}}\left({\underline{\Psi}}^{T}\,{\underline{\underline{W}}}_{j}\right)\,d{\underline{x}}\,dt}_{\text{(internal
flux term)}}\\\
+&\underbrace{\int_{t^{n}}^{t^{n+1}}\oint_{\partial{\mathcal{T}}^{\star}_{j}}\Psi_{k}\,{\underline{{\mathcal{F}}}}\left({\underline{w}}^{h}_{+},{\underline{w}}^{h}_{-};{\underline{n}}\right)\,d{\underline{s}}\,dt}_{\text{(surface
flux
term)}}-\underbrace{\left[\int_{{\mathcal{T}}_{j}}\left(\Psi_{k}{\underline{\Phi}}^{T}\right)\Bigl{|}_{t^{n}}\,d{\underline{x}}\right]\,{\underline{\underline{Q}}}^{n}_{j}}_{\text{(causal
source)}}={\underline{0}},\end{split}$ (6)

where $\partial{\mathcal{T}}^{\star}_{j}$ is the part of the boundary of
${\mathcal{T}}_{j}$ that is interior to ${\mathcal{R}}_{i}$ and
${\underline{n}}$ is an outward pointing normal to on
$\partial{\mathcal{T}}^{\star}_{j}$. The states ${\underline{w}}^{h}_{\pm}$
are the solution on either sides of the boundary
$\partial{\mathcal{T}}^{\star}_{j}$. In this work, the numerical flux,
${\underline{{\mathcal{F}}}}$, is taken to be the Rusanov (sometimes called
the local Lax-Friedrichs) flux article:Ru61 .

On a Cartesian mesh in ${M_{\text{dim}}}$ dimensions, equation eq. 6
represents a nonlinear algebraic system of size
${M_{\text{eqn}}}\cdot{M_{\text{P}}}\cdot 3^{M_{\text{dim}}}$. In fact, the
only portion of the solution to eq. 6 that we actually need to retain is the
solution on the central element: ${\mathcal{T}}_{i}$; the solution on the
remaining elements in ${\mathcal{R}}_{i}$ will be discarded (again, see fig. 1
).

This prediction step portion of the RIDG scheme is by far the most expensive.
Furthermore, the solution to eq. 6 by itself is quite useless: this solution
is not even consistent with the original PDE due to the fact that region
${\mathcal{R}}_{i}$ has been insulated from all other regions in the
computational domain. However, using this solution as a prediction inside of
an appropriate correction step produces a scheme that is both high-order
accurate and stable up to CFL numbers that significantly exceed Lax-Wendroff
DG and SSP-RKDG RIDG_paper_2019 .

###  2.3 Correction step for RIDG

For the correction step we assume a solution of the form eq. 2 and use of a
version of equation eq. 3 where the solution inside each of the spacetime
integrals is replaced by the predicted solution eq. 5 :

$\begin{split}{\underline{Q}}_{i}^{k}\left(t^{n+1}\right)&={\underline{Q}}_{i}^{k}\left(t^{n}\right)+\frac{1}{|{\mathcal{T}}_{i}|}\int_{t^{n}}^{t^{n+1}}\int_{{\mathcal{T}}_{i}}{\underline{\nabla}}\Phi_{k}\cdot{\underline{\underline{F}}}\left({\underline{\Psi}}^{T}\,{\underline{\underline{W}}}_{i}\right)\,d{\underline{x}}\,dt\\\
&-\frac{1}{|{\mathcal{T}}_{i}|}\int_{t^{n}}^{t^{n+1}}\oint_{\partial{\mathcal{T}}_{i}}\Phi_{k}\,{\underline{{\mathcal{F}}}}\left({\underline{w}}^{h}_{+},{\underline{w}}^{h}_{-};{\underline{n}}\right)\,d{\underline{s}}\,dt.\end{split}$

In the above equatoion ${\underline{n}}$ is an outward-pointing normal vector
to $\partial{\mathcal{T}}_{i}$, ${\underline{w}}^{h}_{+}$ and
${\underline{w}}^{h}_{-}$ are the states on either side of the boundary
$\partial{\mathcal{T}}_{i}$, and ${\underline{{\mathcal{F}}}}$ is the
numerical flux. In this work we always make use of the Rusanov (sometimes
called the local Lax-Friedrichs) flux article:Ru61 . This portion of the
regionally-implicit DG update is fully explicit and computationally
inexpensive relative to the prediction step.

##  3 Efficient implementation of RIDG via quasi-quadrature-free Jacobian
matrix assembly

As described in section 2.2 , the prediction step in the regionally-implicit
DG method is the most computationally expensive, since it requires the
solution of a nonlinear system of algebraic equations over each region
${\mathcal{R}}_{i}$. Due to the computational complexity of space-time DG
methods, efficient matrix assembly is required if high order RIDG methods are
to be viable. Therefore, in order to make this step as efficient as possible,
we develop in this section a quasi-quadrature-free Jacobian matrix assembly.

Quadrature-free schemes have provided efficient alternatives to quadrature in
DG schemes for many years Atkins1998 . Depending on the basis choice, they can
benefit from reduced computational complexity for high-dimensional problems
without sacrificing accuracy. That is, the quadrature free schemes can be
designed to inherit the same assumptions about an underlying integrand as do
discrete quadrature schemes. Thus, quadrature-free integrals are efficient
replacements for standard quadrature schemes, especially when orthogonal bases
are used. Matrix assembly for nodal elements can have additional exploits
Engsig-Karup2016 . With some work, quarature-free schemes can also be
accurately extended to unstructured geometries chan2016 . In the current work
we extend this previous work to high-order space-time Jacobian matrix
assembly, as required by the Newton method for solving the nonlinear systems
inherent to the RIDG prediction step Guthrey2017 ; RIDG_paper_2019 .

The prediction step algebraic system given by eq. 6 is typically solved via a
Newton’s method, which can be written as

${\underline{W}}^{(m+1)}={\underline{W}}^{(m)}-\left[{{\underline{\underline{J}}}}\left({\underline{W}}^{(m)}\right)\right]^{-1}{\underline{R}}\left({\underline{W}}^{(m)}\right),$

where ${\underline{R}}$ is shorthand for the residual defined in eq. 6 written
as a column vector, ${\underline{W}}$ is shorthand for all of the unknown
cofficients in region ${\mathcal{R}}_{i}$ written as a column vector, $m$ is
the Newton iteration counter, and ${\underline{\underline{J}}}$ is the
Jacobian of ${\underline{R}}$ with respect to ${\underline{W}}$.

For ease of discussion, let us focus our description of the Jacobian-free
implementation on only the internal flux term in the residual and let us
consider only a scalar case (${M_{\text{eqn}}}=1$) – these restrictions can
easily be removed:

$\begin{split}{R}_{jk}&=(\text{time
term})-\int_{t^{n}}^{t^{n+1}}\int_{{\mathcal{T}}_{j}}{\underline{\nabla}}\Psi_{k}\cdot{\underline{F}}\left({\underline{\Psi}}^{T}\,{\underline{W}}_{j}\right)\,d{\underline{x}}\,dt\\\
&+(\text{surface flux term})-(\text{causal source}),\end{split}$ (7)

for each $j:{\mathcal{T}}_{j}\in{\mathcal{R}}_{i}$ and for each
$k=1,2,\ldots,{M_{\text{P}}}$, where
${\underline{W}}_{j}\in\mathbb{R}^{{M_{\text{P}}}}$ are the unknown
coefficients in spacetime element ${\mathcal{S}}_{j}$. The Jacobian of this
residual with respect to the degrees of freedom ${\underline{W}}$ is

$\begin{split}{J}_{jk\ell}&=\frac{\partial}{\partial{W_{\ell}}}(\text{time
term})-\int_{t^{n}}^{t^{n+1}}\int_{{\mathcal{T}}_{j}}{\underline{\nabla}}\Psi_{k}\cdot{\underline{F}}^{\prime}\left({\underline{\Psi}}^{T}\,{\underline{W}}_{j}\right)\,{\Psi}_{\ell}\,d{\underline{x}}\,dt\\\
&+\frac{\partial}{\partial{W_{\ell}}}(\text{surface flux term}),\end{split}$
(8)

where ${\underline{F}}^{\prime}(q)\cdot{\underline{n}}$ is the flux Jacobian
of the hyperbolic conservation law in direction ${\underline{n}}$.

In order to assemble the above Jacobian matrix for the Newton iteration,
several methods are available.

  * •

Quadrature: We could perform the integrals in eq. 8 via quadrature for every
$k,\ell\leq\theta$, where $\theta$ is the number of terms in our space-time
Legendre basis. This method is able to take advantage of the sparsity pattern
of our block-stencil. The computational complexity of such a method is

$\text{timing}\approx\mathcal{O}(\theta^{2}({M_{\text{deg}}}+1)^{{M_{\text{dim}}}+1})\approx\mathcal{O}(({M_{\text{deg}}}+1)^{3{M_{\text{dim}}}+3}),$
(9)

which will become incredibly expensive for high-order methods in three
dimensions.

  * •

Perturbation: One alternative to direct quadrature is to approximate the
Jacobian ${J}_{jk\ell}$ via either finite differences or Gateaux derivatives
of the residual eq. 7 ; this is also known as a perturbation method. This
method is very simple to implement, as you only need to define the residual
for your Newton iteration. However, since the residual needs to be recomputed
a number of times equal to the number of degrees of freedom, the computational
complexity is roughly the same as eq. 9 .

  * •

Quadrature-free: If the flux function is sufficiently simple (i.e., a
polynomial function of the solution), then an alternative to the quadrature
and perturbation approaches is to perform exact integrations of the terms
needed in the Jacobian; this was the method employed in Guthrey and Rossmanith
RIDG_paper_2019 . Unfortunately, this method does not generalize well, even to
relatively simple rational flux functions such as those seen in the
compressible Euler equations.

In this work we consider an alternative to all of these approaches; it is
similar to the quadrature-free approach, but applies to general fluxes (i.e.,
non-polynomial). The idea is this: we project components of the flux Jacobian,
${\underline{F}}^{\prime}(w)$, onto the polynomial basis ${\underline{\Psi}}$:

${\underline{F}}^{\prime}_{\,p}=\left\\{\int_{t^{n}}^{t^{n+1}}\int_{{\mathcal{T}}_{j}}\Psi_{p}\,{\underline{F}}^{\prime}(w)\,d{\underline{x}}\,dt\right\\}\biggl{/}\left\\{\int_{t^{n}}^{t^{n+1}}\int_{{\mathcal{T}}_{j}}\left(\Psi_{p}\right)^{2}\,d{\underline{x}}\,dt\right\\}.$
(10)

Using the expansion of ${\underline{F}}^{\prime}(w)$ in the expression for the
residual eq. 8 , we obtain

$\begin{split}{J}_{jk\ell}&=\frac{\partial}{\partial{W_{\ell}}}(\text{time
term})-\sum_{p=1}^{\infty}\left(\int_{t^{n}}^{t^{n+1}}\int_{{\mathcal{T}}_{j}}{\underline{\nabla}}\Psi_{k}\,\Psi_{p}\,{\Psi}_{\ell}\,d{\underline{x}}\,dt\right)\cdot{\underline{F}}^{\prime}_{\,p}\\\
&+\frac{\partial}{\partial{W_{\ell}}}(\text{surface flux term}).\end{split}$
(11)

Due to the orthogonality of our basis, there exists some $L$ such that for all
$p>L$ and all needed $j$, $k$, and $\ell$, the integrals in the above
expression all vanish (see Gupta and Narasimhan Gupta2007 ). That is, there
are a limited amount of the integrals in eq. 11 that are actually nonzero.
Furthermore, these integral expressions can be precomputed exactly and stored,
and quadrature is only needed for the projection eq. 10 . Thus, for each
$j,k,\ell$ we may form a list $\alpha_{jk\ell}$ of the indices $p$ such that
the associated integral was nonzero and a list ${\underline{\beta}}_{jk\ell
p}$ of the precomputed results of the nonzero integrals. Finally, we can then
simplify eq. 11 to the following:

$\begin{split}{J}_{jk\ell}&=\frac{\partial}{\partial{W_{\ell}}}(\text{time
term})-\sum_{p=1}^{\alpha_{jk\ell}}\left({\underline{\beta}}_{jk\ell
p}\cdot{\underline{F}}^{\prime}_{\,p}\right)+\frac{\partial}{\partial{W_{\ell}}}(\text{surface
flux term}).\end{split}$

23456$10^{1}$$10^{3}$$10^{5}$$10^{7}$$10^{9}$$3.73$$5.59$$7.04$$8.15$$9.18$$2.1$$4.05$$5.5$$6.6$$7.54$$1.58$$2.94$$4.03$$4.94$$5.73$method
order $M_{deg}+1$log timings ($\mu$s)PerturbationQuadratureQuasi-quadrature-
free Figure 2: Newton iteration Jacobian matrix assembly time for one 3D+1
space-time prediction for a nonlinear problem, using (1, blue) a theoretical
estimate of the cost of a perturbation method obtained by multiplying the
timings of the residual computation by the number of degrees of freedom in the
Newton iteration (2, red) traditional quadrature routines to compute the
integrals in eq. 8, and (3, tan) the quasi-quadrature-free routine discussed
in this appendix. Runtimes are $\log_{10}$ of averages of 10 trial runs in
microseconds. We see that the quasi-quadrature-free method is the most
efficient, and scales better as we increase ${M_{\text{deg}}}$. For
${M_{\text{deg}}}+1=6$, the quasi-quadrature-free method is almost two orders
of magnitude faster than the traditional quadrature scheme.

This quasi-quadrature-free method maintains the same order of accuracy as the
quadrature or perturbation strategies described above, but in practice is
orders of magnitude faster, because the space-time quadrature is performed
once per degree of freedom per element as opposed to once per square of the
degrees of freedom per element. We compare timings of these three methods with
${M_{\text{dim}}}=3$ in Table fig. 2 . At this time we do not have an a priori
estimate for the computational complexity of the quasi-quadrature-free method,
but we are able to compute a posteriori estimates for each method:

$\displaystyle\text{perturbation timings}\approx$
$\displaystyle\;\mathcal{O}(({M_{\text{deg}}}+1)^{11.9})$
$\displaystyle\text{quadrature timings}\approx$
$\displaystyle\;\mathcal{O}(({M_{\text{deg}}}+1)^{11.6})$
$\displaystyle\text{quasi-quadrature-free timings}\approx$
$\displaystyle\;\mathcal{O}(({M_{\text{deg}}}+1)^{9.3})$

We see that the perturbation and traditional quadrature routines exhibit the
expected computational complexity described in eq. 9, but that the quasi-
quadrature-free method exhibits a significantly lower computational
complexity.

For brevity and clarity we have omitted the discussion of the temporal
derivative terms and surface flux terms. However, we note that the time terms
of the matrix assembly of $J_{jk\ell}$ can be precomputed and stored, and no
projection is needed during the matrix assembly. Furthermore, the surface
terms in general involve a numerical flux, for which a numerical flux Jacobian
must be generated or approximated before it is projected. These types of terms
affect several blocks of the Jacobian matrix, exactly in accordance to the
block-stencil derived from our numerical flux.

To review, we have efficiently extended quadrature-free schemes to our space-
time prediction procedure needed for our regionally-implicit DG method. This
implementation offers incredible speedups to a very computationally expensive
method in a way that is extensible to arbitrary problems and relatively simple
to implement. We will use this strategy in the remainder of this paper.

##  4 Basic terminology

The goal of this work is to extend the convergence studies performed in
Guthrey and Rossmanith RIDG_paper_2019  to the case of strong scaling in an
HPC setting. Before reporting results in subsequent sections, we describe here
some of the key terms needed to analyze the efficiency and accuracy of
parallelized DG schemes. We also note that here and throughout the remainder
of the paper we consider only the scalar: ${M_{\text{eqn}}}=1$. All the
findings in this work generalize to the more general case, but for brevity and
clarity we assume for the remainder of this current work that
${M_{\text{eqn}}}=1$.

We define the following key terms.

mesh:

The description of the number of elements. All methods considered in this
paper use a uniform Cartesian mesh.

•

dof : The total number of degrees of freedom. We compute this by multiplying
basis element size ($\theta$) by the number of elements (mesh):

$\text{\bf dof}(\theta,\text{\bf mesh})=\theta\times\text{\bf mesh}.$

efom:

High-order DG methods can be compared to first-order finite volume methods by
(among other ways) comparing the equivalent number of the total degrees of
freedom. Thus, efom, equivalent first order mesh, indicates what size of a
Cartesian mesh with one degree of freedom per cell would have the equivalent
number of degrees of freedom for a given mesh and method order. The efom for a
method with a basis of size $\theta$ for a ${M_{\text{dim}}}$-dimensional
problem is computed using

$\text{\bf efom}(\theta,\text{\bf mesh})=\left[\left(\text{\bf
dof}\right)^{\frac{1}{{M_{\text{dim}}}}}\right]^{{M_{\text{dim}}}}=\left[\left(\theta\times\text{\bf
mesh}\right)^{\frac{1}{{M_{\text{dim}}}}}\right]^{{M_{\text{dim}}}}.$ (12)

For example, a $10\times 10$ mesh in 2D with a basis of size $16$ has an efom
of $40\times 40$ since it has the equivalent number of degrees of freedom as a
$40\times 40$ first order mesh.

error:

The approximate relative error in the $L^{2}$ norm of our numerical solution.
This error is obtain by considering an exact projection $q_{\text{exact}}$
onto an infinite basis, ${\underline{\Phi}}$, on each element:
${\underline{X}}_{i}:=\underset{{\underline{\Phi}}^{\infty}}{\operatorname{proj}}\,q_{\text{exact}}$.
Comparing this exact expansion to the numerical solution gives the following
element-wise error definition:

$e({\underline{x}})\Bigl{|}_{{\mathcal{T}}_{i}}:=\left|\underbrace{\sum\limits_{\ell=1}^{\infty}\Phi_{\ell}({\underline{x}})\,X^{\ell}_{i}}_{\text{exact}}-\underbrace{\sum\limits_{\ell=1}^{M_{\text{basis}}}\Phi_{\ell}({\underline{x}})\,Q^{\ell}_{i}}_{\text{numerical}}\right|.$

Using the orthonormality of the basis functions gives the following $L^{2}$
error:

${\left\|e\right\|}_{L^{2}}^{2}=\sum_{i}\left[\underbrace{\sum\limits_{\ell=1}^{M_{\text{basis}}}\left(X_{i}^{\ell}-Q^{\ell}_{i}\right)^{2}}_{\text{error
in
coeffs}}+\underbrace{\sum\limits_{\ell=M_{\text{basis}}+1}^{M_{\text{basis}}^{+}}\left(X_{i}^{\ell}\right)^{2}}_{\text{dominant
trunc.
error}}+\underbrace{\sum\limits_{\ell=M_{\text{basis}}^{+}+1}^{\infty}\left(X_{i}^{\ell}\right)^{2}}_{\text{high-
order trunc. error}}\right],$

where $M_{\text{basis}}$ and $M_{\text{basis}}^{+}$ are the number of basis
functions in ${\mathbb{Q}}\left({M_{\text{deg}}},{M_{\text{dim}}}\right)$ and
${\mathbb{Q}}\left({M_{\text{deg}}}+1,{M_{\text{dim}}}\right)$, respectively.
To sufficient accuracy the error can be approximated by discarding all the
terms past $\ell=M_{\text{basis}}^{+}$, resulting in the following approximate
relative error:

$\text{ \bf
error}=\frac{{\left\|e\right\|}_{L^{2}}}{{\left\|q_{\text{exact}}\right\|}_{L^{2}}}\approx\sqrt{\frac{\sum\limits_{i=1}\left[\sum\limits_{\ell=1}^{M_{\text{basis}}}\left(X_{i}^{\ell}-Q^{\ell}_{i}\right)^{2}+\sum\limits_{\ell=M_{\text{basis}}+1}^{M_{\text{basis}}^{+}}\left(X_{i}^{\ell}\right)^{2}\right]}{\sum\limits_{i=1}\left[\sum\limits_{\ell=1}^{M_{\text{basis}}^{+}}\left(X_{i}^{\ell}\right)^{2}\right]}}.$
(13)

•

aproximate order of accuracy: For the convergence studies performed in this
paper, we approximate the convergence rate $M$ using

$\text{\bf
error}(h)=ch^{M}+{\mathcal{O}}\left(h^{M+1}\right)\quad\Longrightarrow\quad
M\approx\frac{\log\left({\text{\bf error}(h_{1})}/{\text{\bf
error}(h_{2})}\right)}{\log\left({h_{1}}/{h_{2}}\right)},$ (14)

where $h$ is the mesh spacing.

runtime:

For each experiment we compute the wall clock runtime of each method in
seconds. For RIDG methods, the computational effort is dominated by small
dense matrix inverses in the prediction step, and thus for a method that
requires $N_{t}$ time-steps with a mesh of size $N^{M_{\text{dim}}}$ and a
space-time basis of size $\theta_{T}$, we expect the runtime to scale as

$\text{\bf runtime}\approx\mathcal{O}\left(N_{t}\cdot
N^{M_{\text{dim}}}\cdot(3^{M_{\text{dim}}}\theta_{T})^{3}\right).$ (15)

If we use the $\mathcal{Q}$ spacetime basis for the prediction step as
discussed in RIDG_paper_2019  and we let $M={M_{\text{deg}}}+1$, then
$\theta_{T}=M^{{M_{\text{dim}}}+1}$ and so for explicit time-stepping where
$\Delta t\approx\nu\Delta x\approx\nu/N$, where $\nu$ is the CFL number, we
get that

$\text{\bf
runtime}\approx\mathcal{O}\left(27^{{M_{\text{dim}}}}\cdot\nu^{-1}\cdot
N^{{M_{\text{dim}}}+1}\cdot M^{3{M_{\text{dim}}}+3}\right).$ (16)

RIDG methods have the fortunate property that $\nu^{-1}\approx\mathcal{O}(1)$,
whereas methods such as SSP-RKDG experience
$\nu^{-1}\approx\mathcal{O}(M^{2})$.

quality:

When comparing various methods in terms of efficiency, the two most prevalent
metrics are runtime and error. One could consider how fast various methods
reach a fixed error (time to solution), or vice versa. In this paper we
combine these two metrics into a third metric we call quality. We define the
quality of the solution as

$\text{\bf quality}=-\log\left(\text{\bf error}\times\text{\bf
runtime}\right).$ (17)

We note that this metric is on a logarithmic scale, so changes of $\pm 1$ for
this metric are quite significant. If we assume that RIDG methods of order
$M={M_{\text{deg}}}+1$, have an error convergence scaling of
$\mathcal{O}(N^{-M})$, then using eq. 16 we expect the quality to scale as

$\text{\bf
quality}\approx\mathcal{O}\left(\nu+(M-{M_{\text{dim}}}-1)\log(N)-(3{M_{\text{dim}}}+3)\log(M)\right).$
(18)

This provides the intuition that high order RIDG methods will have a high
quality, offset by increased runtime costs associated with the dense matrix
inverses. However, as we consider higher mesh resolutions, the quality of a
method of fixed order will be driven higher by the shrinking error. For
smaller mesh resolutions, the quality will be adversely dominated by the high
runtime cost of the methods.

dof/c:

This is the number of degrees of freedom per compute core, and is computed as:

$\text{\bf dof/c}=\frac{\theta\times{\bf mesh}}{\bf
cores}=\frac{\theta\times{\bf mesh}}{\bf tasks\times cores/task},$ (19)

where $\theta$ is the size of the DG basis on each element. This quantity is
one way to measure the relative compute intensity per core for a given
experiment.

speedup over single task runtimes:

We measure the ratios of runtimes of the single task run versus the multi-task
runs to compare how much faster the latter case runs. Ideally we expect that
the speedup is equal to the total number of tasks being used. It can be
described with the formula

$\text{\bf strong speedup}(\text{\bf tasks}=i)=\frac{\text{\bf
runtime}(\text{\bf tasks}=1)}{\text{\bf runtime}(\text{\bf tasks}=i)}.$ (20)

Note that this definition of speedup compares runtimes for a given set of
tasks as opposed to cores. That is, we do not compare runtimes against serial
execution.

strong efficiency:

This is a measure of how our speedup is to the ideal case. It is simply the
ratio of actual speedup versus ideal speedup. As we add tasks, the parallel
efficiency for a strong scaling test is:

$\text{\bf strong efficiency}(\text{\bf tasks}=i)=\frac{\text{\bf
speedup}(\text{\bf tasks}=i)-1}{i-1}\times 100\%.$ (21)

An efficiency of $100\%$ corresponds to perfect linear speedup scaling with
the number of tasks (or cores). An efficiency of $0\%$ means that as we used
more tasks, the runtime remained unchanged or increased (we map negative
efficiencies to $0\%$).

comms:

This is an estimate for the total number of MPI communications. It is
estimated using the formula

$\text{\bf comms}={\bf tasks}\times{\bf timesteps}\times{\bf
stages/timestep}\times{\bf comms/stage}.$ (22)

We approximate the number of timesteps via the CFL relation

$\Delta t=\nu\Delta x\quad\implies\quad\text{\bf timesteps}=\nu^{-1}N,$ (23)

where $N$ is the number of mesh elements in each direction. We note that RKDG
methods of $M_{deg}=3$ have 9 stages per timestep, and that all RIDG methods
have 2 stages (predictor and corrector). There are 8 communications per stage
in 2D and 26 in 3D.

##  5 Parallel implementation of RIDG

We now briefly discuss the domain decomposition strategy implemented for
achieving high performance computing with RKDG and RIDG methods. We consider a
Cartesian mesh evenly subdivided into a $N\times N$ (2D) or $N\times N\times
N$ (3D) grid of Cartesian submeshes, where $N^{2}$ (2D) or $N^{3}$ (3D) is the
number of compute nodes used. Each node (i.e., submesh) is associated with a
single MPI task. OpenMP is used for shared-memory parallelization of
operations on each Cartesian submesh. MPI is used to perform data
communications across the Cartesian submesh interfaces, often called the ghost
zone or halo region communications. Submeshes that share interfaces will need
to communicate multiple times per timestep, as detailed in the next two
sections.

###  5.1 Domain decomposition strategy for RKDG

In order to compute the needed fluxes for each stage of RKDG, each Cartesian
submesh requires information from face-neighbor submeshes. This means that for
RKDG the submeshes must communicate with neighboring submeshes in each stage.
The communication latency for this operation can be hidden by using non-
blocking MPI routines for the intra-face communications before we perform the
volume integrals over the submesh. Each volume integral has a relatively high
arithmetic intensity and can be completed using information completely local
to each element. Once these volume integrals are computed, we simply wait
until intra-face communications are completed, which ideally has already
occurred. Then, the fluxes and thus the update for the RK stage can be
computed. Lastly, each timestep requires an all-to-all communication to
enforce the CFL time-step restriction by communicating the CFL number used by
each submesh.

###  5.2 Domain decomposition strategy for RIDG

For RIDG the prediction step requires that each given Cartesian submesh has
information from vertex-neighbor submeshes, since forming a region for a cell
involve the cell’s vertex-neighbors. Latency associated with this
communication can be hidden by using non-blocking MPI routines for the
interface communications. After these communications are started, we perform
the prediction step for elements not on the submesh boundary, as data for
their vertex neighbors is located in shared memory. Once the boundary
communications are complete, we can and then compute the predictions for
boundary elements. We must communicate these boundary predictions as they are
needed for the correction step. Again we may use non-blocking routines and
finish computing the non-boundary predictions during these communications.
Once all predictions are formed and the communications are complete, we
continue to the correction step. Just as with RKDG, each timestep requires an
all-to-all communication to enforce the CFL time-step restriction by
communicating the CFL number used by each submesh. This all-to-all
communication can also be performed using a non-blocking MPI routine while the
correction update is computed. For nonlinear problems, the prediction step and
its associated communications must be repeated until the region residual is
driven down to some tolerance.

##  6 Numerical results

In order validate the proposed implementation of the regionally-implicit
scheme we consider in this section examples 1D, 2D, and 3D. The 1D code runs
sufficiently fast that no parallelization is required; we provide the 1D
results mostly to demonstrate the efficiency gains of high-order RIDG.
Subsequently we consider 2D and 3D examples, all of which are implemented in
parallel.

###  6.1 1D results

We consider the following 1D linear advection equation:

$\begin{cases}q_{,t}+q_{,x}=0&\text{for}\quad(t,x)\in\left[0,T\right]\times\left[0,1\right],\\\
q(t=0,x)=e^{((x-c)^{2}-\omega^{2})^{-1}}&\text{if}\quad(x-c)^{2}<\omega^{2},\\\
q(t=0,x)=0&\text{if}\quad(x-c)^{2}\geq\omega^{2},\end{cases}$ (24)

with $c=\frac{1}{2}$ and $\omega=\frac{1}{3}$. Note that the initial condition
is clearly $C^{\infty}$ for $\left|x-c\right|<\omega$; furthermore, it can be
shown that all derivatives for this function vanish as
$(x-c)^{2}\to\omega^{2}$ from the left, and hence the initial condition is
actually $C^{\infty}$ over all of $\mathbb{R}$. This level of smoothness is
required to test the convergence properties of arbitrarily high order methods.

1D RKDG: (${M_{\text{deg}}}=3$, $\nu=0.1$)
---
mesh | dof | $L^{2}$ error (13) |  order (14) | runtime (s) |  quality (17)
$50$ | $200$ | ${3.07}\mathrm{e}{-}4$ | – | ${5.16}\mathrm{e}{-}1$ | $3.8$
$70$ | $280$ | ${7.96}\mathrm{e}{-}5$ | $4.0$ | ${8.95}\mathrm{e}{-}1$ | $4.1$
$120$ | $480$ | ${9.23}\mathrm{e}{-}6$ | $4.0$ | ${2.34}\mathrm{e}{+}0$ | $4.7$
$240$ | $960$ | ${5.77}\mathrm{e}{-}7$ | $4.0$ | ${7.94}\mathrm{e}{+}0$ | $5.3$
1D RIDG: (${M_{\text{deg}}}=3$, $\nu=0.9$)
mesh | dof | $L^{2}$ error (13) |  order (14) | runtime (s) |  quality (17)
$50$ | $200$ | ${6.48}\mathrm{e}{-}4$ | – | ${1.09}\mathrm{e}{+}0$ | $3.2$
$70$ | $280$ | ${1.47}\mathrm{e}{-}4$ | $4.4$ | ${1.07}\mathrm{e}{+}0$ | $3.8$
$120$ | $480$ | ${1.85}\mathrm{e}{-}5$ | $3.8$ | ${2.96}\mathrm{e}{+}0$ | $4.3$
$240$ | $960$ | ${8.31}\mathrm{e}{-}7$ | $4.5$ | ${1.24}\mathrm{e}{+}1$ | $5.0$
1D RIDG: (${M_{\text{deg}}}=5$, $\nu=0.9$)
mesh | dof | $L^{2}$ error (13) |  order (14) | runtime (s) |  quality (17)
$30$ | $180$ | ${1.57}\mathrm{e}{-}5$ | – | ${8.29}\mathrm{e}{-}1$ | $4.9$
$50$ | $300$ | ${8.62}\mathrm{e}{-}7$ | $5.7$ | ${4.21}\mathrm{e}{+}0$ | $5.4$
$70$ | $420$ | ${1.27}\mathrm{e}{-}7$ | $5.7$ | ${4.17}\mathrm{e}{+}0$ | $6.3$
$120$ | $720$ | ${3.18}\mathrm{e}{-}9$ | $6.8$ | ${1.24}\mathrm{e}{+}1$ | $7.4$
1D RIDG: (${M_{\text{deg}}}=7$, $\nu=0.9$)
mesh | dof | $L^{2}$ error (13) |  order (14) | runtime (s) |  quality (17)
$20$ | $160$ | ${1.32}\mathrm{e}{-}6$ | – | ${1.25}\mathrm{e}{+}0$ | $5.8$
$30$ | $240$ | ${4.77}\mathrm{e}{-}8$ | $8.2$ | ${2.86}\mathrm{e}{+}0$ | $6.9$
$50$ | $400$ | ${6.65}\mathrm{e}{-}10$ | $8.4$ | ${1.54}\mathrm{e}{+}1$ | $8.0$
$70$ | $560$ | ${6.63}\mathrm{e}{-}11$ | $6.9$ | ${1.48}\mathrm{e}{+}1$ | $9.0$
1D RIDG: (${M_{\text{deg}}}=9$, $\nu=0.9$)
mesh | dof | $L^{2}$ error (13) |  order (14) | runtime (s) |  quality (17)
$10$ | $100$ | ${4.78}\mathrm{e}{-}6$ | – | ${1.05}\mathrm{e}{+}0$ | $5.3$
$20$ | $200$ | ${4.62}\mathrm{e}{-}9$ | $10.0$ | ${4.22}\mathrm{e}{+}0$ | $7.7$
$30$ | $300$ | ${8.25}\mathrm{e}{-}11$ | $9.9$ | ${9.57}\mathrm{e}{+}0$ | $9.1$
$50$ | $500$ | ${6.11}\mathrm{e}{-}13$ | $9.6$ | ${4.83}\mathrm{e}{+}1$ | $10.5$
1D RIDG: (${M_{\text{deg}}}=11$, $\nu=0.9$)
mesh | dof | $L^{2}$ error (13) |  order (14) | runtime (s) |  quality (17)
$5$ | $60$ | ${1.82}\mathrm{e}{-}4$ | – | ${8.72}\mathrm{e}{-}1$ | $3.8$
$10$ | $120$ | ${5.88}\mathrm{e}{-}8$ | $11.6$ | ${3.52}\mathrm{e}{+}0$ | $6.7$
$20$ | $240$ | ${1.40}\mathrm{e}{-}11$ | $12.0$ | ${1.28}\mathrm{e}{+}1$ | $9.7$
$30$ | $360$ | ${1.21}\mathrm{e}{-}13$ | $11.7$ | ${3.08}\mathrm{e}{+}1$ | $11.4$
Table 1: Convergence and runtime study for 1D SSP-RKDG ${M_{\text{deg}}}=3$
and RIDG methods of various ${M_{\text{deg}}}$ for the problem defined by eq.
24. We list the mesh size, total number of degrees of freedom, approximate
$L^{2}$ error compared to the exact solution, the approximate observed order
of convergence, the wall clock runtime in seconds, and the quality for each
experiment. We see that the RIDG method can use the same CFL $\nu=0.9$ for all
values of ${M_{\text{deg}}}$. The highest quality solutions are those produced
by the highest order RIDG methods.

The results of these convergence tests are shown in table 1 . We notice that
the maximum allowable CFL restriction for the RIDG methods is bounded from
below by some $\nu_{\text{min}}$; the key point is that $\nu_{\text{min}}$ for
RIDG is independent of ${M_{\text{deg}}}$, meaning that the maximum allowable
time-step does not degrade as ${M_{\text{deg}}}$ increases, which is in
contrast to other explicit DG time-stepping approaches (i.e., RKDG and Lax-
Wendroff DG). We observe the expected order of convergence for all methods. We
also notice that as we consider higher order RIDG methods, the maximum quality
increases. Furthermore, we notice that the ${M_{\text{deg}}}=3$ RKDG method
has a higher quality than the ${M_{\text{deg}}}=3$ RIDG method, however higher
order RIDG methods each show a continual increase in the solution quality of
1-2 per 2 increase of ${M_{\text{deg}}}$. Note that each increase of 1 in the
solution quality indicates a tenfold decrease in the product of error and
runtime (see eq. 17 ).

We can conclude from these results that while the ${M_{\text{deg}}}=3$ RKDG
method has a superior quality metric than the RIDG method of the same order,
the RIDG methods of ${M_{\text{deg}}}\geq 5$ have a superior quality metric
compared to the RKDG method. Furthermore, higher order RIDG methods have a
superior quality per degree of freedom, since as we increase
${M_{\text{deg}}}$ we are able to obtain higher solution qualities with
smaller meshes, thus leading to fewer total degrees of freedom.

###  6.2 2D results

Next we consider the RIDG scheme in two spatial dimensions. In particular, we
test the full parallel implementation of RIDG by running strong scaling
studies to test the manycore capabilities of RIDG as compared to RKDG. Strong
scaling is achieved by fixing a problem size and efficiently increasing the
amount of compute cores used to solve the same problem. For explicit
timestepping methods, this translates to considering a fixed mesh size and
subdivide the mesh into smaller pieces, where the number of subdivisions is
equal to the number of tasks, as depicted in table 2 . A full strong scaling
study subdivides the problem among compute cores until the number of elements
(also known as cells or zones) per compute core is minimized (ideally to
unity). This causes the the ratio of MPI communications to work/computation to
increase to some maximum value. The 2D strong scaling studies provided in this
section are performed at the Institute for Cyber Enabled Research at Michigan
State University, using Intel(R) Xeon(R) Gold 6148 2.40GHz CPUs.

2D Strong Scaling Study
---
tasks | mesh/task | subdomain | total cores | elements/core
$1$ | $60^{2}$ | $\left[0.000,1.000\right]^{2}$ | $40$ | $90.0$
$4$ | $30^{2}$ | $\left[0.000,0.500\right]^{2}$ | $160$ | $22.5$
$9$ | $20^{2}$ | $\left[0.000,0.33\bar{3}\right]^{2}$ | $360$ | $10.0$
$16$ | $15^{2}$ | $\left[0.000,0.250\right]^{2}$ | $640$ | $5.6$
$25$ | $12^{2}$ | $\left[0.000,0.200\right]^{2}$ | $1000$ | $3.6$
$36$ | $10^{2}$ | $\left[0.000,0.16\bar{6}\right]^{2}$ | $1440$ | $2.5$
Table 2: Strong scaling study for 2D RKDG and RIDG. We scale the methods from
a fixed physical domain of $\left[0,1\right]^{2}$ to 36 equally sized
subdomains. Each task, which corresponds to one node and thus several
processors, solves one subdomain. Above we present the total number of tasks,
the mesh size for each subdivision, the physical size of each subdomain, the
total number of cores, and the number of mesh elements per core. We see that
this experiment drives the number of mesh elements per core down to nearly
unity. This scaling study was performed at the Institute for Cyber-Enabled
Research at Michigan State University, using Intel(R) Xeon(R) Gold 6148 CPU @
2.40GHz CPUs.

####  6.2.1 2D advection

For our first test of the scalability of RKDG and RIDG methods in two
dimensions, we consider the the 2D scalar advection equation:

$\begin{cases}q_{,t}+q_{,x}+q_{,y}=0&(t,{\underline{x}})\in\left[0,T\right]\times\left[0,1\right]^{2},\\\
q(t=0,{\underline{x}})=e^{\left(\|{\underline{x}}-{\underline{c}}\|^{2}-\omega^{2}\right)^{-1}}&\text{if}\quad\|{\underline{x}}-{\underline{c}}\|^{2}<\omega^{2},\\\
q(t=0,{\underline{x}})=0&\text{if}\quad\|{\underline{x}}-{\underline{c}}\|^{2}\geq\omega^{2},\end{cases}$
(25)

where ${\underline{c}}=\frac{1}{2}(1,1)^{T}$ and $\omega=\frac{1}{3}$ and we
enforce periodic boundary conditions. We perform the strong scaling study
defined by table 2 and report the results in table 3.

2D RKDG: (${M_{\text{deg}}}=3$, $\nu=0.05$)
---
cores | dof/c | efom | comms | time | eq. 20 | eq. 21 | eq. 17
$40$ | $900.0$ | $190^{2}$ | $9.60{e4}$ | $3.50{e1}$ | $--$ | $--$ | $5.4$
$160$ | $225.0$ | $95^{2}$ | $3.84{e5}$ | $2.63{e1}$ | $1.33\times$ | $11.1\%$ | $5.6$
$360$ | $100.0$ | $63^{2}$ | $8.64{e5}$ | $1.69{e1}$ | $2.07\times$ | $13.4\%$ | $5.8$
$1000$ | $36.0$ | $38^{2}$ | $2.40{e6}$ | $5.01{e1}$ | $0.70\times$ | $0.0\%$ | $5.3$
$1440$ | $25.0$ | $32^{2}$ | $3.46{e6}$ | $6.99{e1}$ | $0.50\times$ | $0.0\%$ | $5.1$
2D RIDG: (${M_{\text{deg}}}=3$, $\nu=0.7$)
cores | dof/c | efom | comms | time | eq. 20 | eq. 21 | eq. 17
$40$ | $1440.0$ | $240^{2}$ | $1.20{e3}$ | $5.40{e1}$ | $1.00\times$ |  $--$, | $5.5$
$160$ | $360.0$ | $120^{2}$ | $4.80{e3}$ | $1.56{e1}$ | $3.45\times$ | $81.7\%$ | $6.1$
$360$ | $160.0$ | $80^{2}$ | $1.08{e4}$ | $1.37{e1}$ | $3.94\times$ | $36.8\%$ | $6.1$
$1000$ | $57.6$ | $48^{2}$ | $3.00{e4}$ | $6.89{e0}$ | $7.83\times$ | $28.5\%$ | $6.4$
$1440$ | $40.0$ | $40^{2}$ | $4.32{e4}$ | $6.07{e0}$ | $8.89\times$ | $22.5\%$ | $6.5$
2D RIDG: (${M_{\text{deg}}}=5$, $\nu=0.7$)
cores | dof/c | efom | comms | time | eq. 20 | eq. 21 | eq. 17
$40$ | $3240.0$ | $360^{2}$ | $1.23{e3}$ | $8.71{e2}$ | $1.00\times$ | $--$ | $7.6$
$160$ | $810.0$ | $180^{2}$ | $4.92{e3}$ | $2.82{e2}$ | $3.09\times$ | $69.5\%$ | $8.1$
$360$ | $360.0$ | $120^{2}$ | $1.11{e4}$ | $2.37{e2}$ | $3.68\times$ | $33.5\%$ | $8.2$
$1000$ | $129.6$ | $72^{2}$ | $3.08{e4}$ | $7.55{e1}$ | $11.54\times$ | $43.9\%$ | $8.7$
$1440$ | $90.0$ | $60^{2}$ | $4.43{e4}$ | $4.44{e1}$ | $19.60\times$ | $53.2\%$ | $8.9$
Table 3: Strong scaling study for the 2D RIDG methods with
${M_{\text{deg}}}=3$ and ${M_{\text{deg}}}=5$ on the 2D advection eq. 25. We
list the total number of cores, the number of degrees of freedom per core, the
equivalent first order mesh as described in eq. 12, the estimated total number
of MPI communications as computed by eq. 22, the runtimes in seconds, the
speedup over single node performance as computed by eq. 20, the scaling
efficiency from a single node eq. 21, and the solution eq. 17.

We notice that the total runtime of the RKDG method is initially less than
that of the RIDG method, but does not scale well with granularity. After 1000
cores, the runtime of RKDG increases as we use more compute cores, indicating
runtime scaling breakdown. That is, the many-node RKDG method exhibits low
speedups/efficiencies compared to the single node runtime. The fastest
solution was obtained using 360 cores. The RIDG method, while initially more
computationally expensive, scales well with the number of added compute
resources. As we add more cores, the RIDG methods of ${M_{\text{deg}}}=3,5$
both exhibit monotonic decreases in the overall method runtime, indicating
successful runtime scaling to high granularity. As we consider higher
granularity, the runtime scaling efficiency for the RIDG method
${M_{\text{deg}}}=3$ begins to drop off. However, the RIDG method
${M_{\text{deg}}}=5$ maintains decent runtime scaling efficiency even at very
high granularity. We see that the maximum solution quality obtained by the
RKDG method is 5.8, while for RIDG ${M_{\text{deg}}}=3$ the maximum quality is
6.5. For RIDG ${M_{\text{deg}}}=5$ the maximum quality is 8.9. Recalling that
solution quality is a logarithmic scale, this means that the RIDG
${M_{\text{deg}}}=5$ solution exhibits a far superior efficiency. Although the
RIDG method has a very hefty computational cost compared to RKDG, it is able
to scale to much higher levels of granularity. Thus, the RIDG method is able
to provide a superior solution quality at scale.

####  6.2.2 2D Burgers

To demonstrate the capability of this method to efficiently solve nonlinear
problems, we consider here is the 2D scalar Burgers equation:

$\begin{cases}q_{,t}+\left(\frac{1}{2}q^{2}\right)_{,x}+\left(\frac{1}{2}q^{2}\right)_{,y}=0&(t,{\underline{x}})\in\left[0,T\right]\times\left[0,1\right]^{2},\\\
q(t=0,{\underline{x}})=\frac{1}{4}(1-\cos(x))(1-\cos(y))&(x,y)\in\left[0,1\right]^{2},\end{cases}$
(26)

where we enforce periodic boundary conditions. We perform the same strong
scaling study as described in table 2 . Our results are listed in table 4 .
These studies are also performed at the Institute for Cyber Enabled Research
at Michigan State University, using Intel(R) Xeon(R) Gold 6148 2.40GHz CPUs.
The entries of the Jacobian eq. 8 are calculated using the quasi-quadrature-
free approach sketched in section 3 .

2D RKDG: (${M_{\text{deg}}}=3$, $\nu=0.05$)
---
cores | dof/c | efom | comms | time | eq. 20 | eq. 21 | eq. 17
$40$ | $900.0$ | $190^{2}$ | ${9.60}\mathrm{e}{+}4$ | ${1.90}\mathrm{e}{+}1$ | $--$ | $--$ | $5.5$
$160$ | $225.0$ | $95^{2}$ | ${3.84}\mathrm{e}{+}5$ | ${7.36}\mathrm{e}{+}0$ | $2.58\times$ | $52.7\%$ | $5.9$
$360$ | $100.0$ | $63^{2}$ | ${8.64}\mathrm{e}{+}5$ | ${5.51}\mathrm{e}{+}0$ | $3.45\times$ | $30.6\%$ | $6.0$
$1000$ | $36.0$ | $38^{2}$ | ${2.40}\mathrm{e}{+}6$ | ${1.95}\mathrm{e}{+}1$ | $0.98\times$ | $0.0\%$ | $5.5$
$1440$ | $25.0$ | $32^{2}$ | ${3.46}\mathrm{e}{+}6$ | ${2.91}\mathrm{e}{+}1$ | $0.65\times$ | $0.0\%$ | $5.3$
2D RIDG: (${M_{\text{deg}}}=3$, $\nu=0.7$)
cores | dof/c | efom | comms | time | eq. 20 | eq. 21 | eq. 17
$40$ | $1440.0$ | $240^{2}$ | ${5.48}\mathrm{e}{+}3$ | ${8.89}\mathrm{e}{+}1$ | $1.00\times$ | $--$ | $5.9$
$160$ | $360.0$ | $120^{2}$ | ${2.20}\mathrm{e}{+}4$ | ${2.35}\mathrm{e}{+}1$ | $3.78\times$ | $92.6\%$ | $6.4$
$360$ | $160.0$ | $80^{2}$ | ${4.92}\mathrm{e}{+}4$ | ${1.26}\mathrm{e}{+}1$ | $7.07\times$ | $75.9\%$ | $6.7$
$1000$ | $57.6$ | $48^{2}$ | ${1.37}\mathrm{e}{+}5$ | ${7.94}\mathrm{e}{+}0$ | $11.20\times$ | $42.5\%$ | $6.9$
$1440$ | $40.0$ | $40^{2}$ | ${1.98}\mathrm{e}{+}5$ | ${5.81}\mathrm{e}{+}0$ | $15.31\times$ | $40.9\%$ | $7.1$
Table 4: Strong scaling study for Burgers eq. 26 using the RIDG methods with
${M_{\text{deg}}}=3$ and RKDG methods ${M_{\text{deg}}}=3$, as defined by
table 2. We list the total number of cores, the number of degrees of freedom
per core, the equivalent first order mesh as described in eq. 12, the
estimated total number of MPI communications as computed by eq. 22, the
runtimes in seconds, the speedup over single node performance as computed by
eq. 20, the scaling efficiency from a single node eq. 21, and the solution
quality eq. 17.

We notice that the RKDG method requires many more total halo-region MPI
communications than the RIDG method. We also notice that the total runtime of
the RKDG method is initially less than that of the RIDG method, but does not
scale well with granularity. After 1000 cores, the runtime of RKDG increases
as we use more compute cores, indicating runtime scaling breakdown. The
fastest solution was obtained using 360 cores. The RIDG method, while
initially more computationally expensive, scales well with the number of added
compute resources. As we add more cores, the RIDG method exhibits monotonic
decreases in the overall method runtime, indicating successful runtime scaling
to high granularity. The fastest runtime was obtained at 1440 cores. The many-
node RKDG method exhibits low speedups/efficiencies compared to the single
node runtime. We notice that at 1000 cores, the speedups drop below $1$ and
the efficiencies drop to $0\%$, indicating runtime scaling breakdown. The
many-node RIDG method exhibits much larger speedups, and thus higher scaling
efficiencies. As we consider higher granularity, the runtime scaling
efficiency for the RIDG method maintains decent runtime scaling efficiency
even at very high granularity. As we consider various levels of granularity,
the maximum solution quality obtained by the RKDG method is 6.0, For RIDG the
maximum quality is 7.1.

We conclude that the RKDG method does not experience runtime scaling to high
granularity for this nonlinear problem in two dimensions, but the RIDG method
${M_{\text{deg}}}=3$ experiences good runtime scaling to high granularity,
providing an improvement in the maximum possible solution quality. This is
despite the additional iterations that the RIDG method must perform for
nonlinear problems greatly increase the computational cost of the method and
the total number of MPI communications.

###  6.3 3D results

In this section we test the suitability of the RIDG method at scale for
solving 3D problems. The method efficiency is again compared to the Runge-
Kutta DG methods. In particular we consider the following 3D linear advection
equation:

$\begin{cases}q_{,t}+q_{,x}+q_{,y}+q_{,z}=0&\text{for}\quad\left(t,{\underline{x}}\right)\in\left[0,T\right]\times\left[0,1\right]^{3},\\\
q(t=0,{\underline{x}})=e^{(\|{\underline{x}}-{\underline{c}}\|^{2}-\omega^{2})^{-1}}&\text{if}\quad\|{\underline{x}}-{\underline{c}}\|^{2}<\omega^{2},\\\
q(t=0,{\underline{x}})=0&\text{if}\quad\|{\underline{x}}-{\underline{c}}\|^{2}\geq\omega^{2},\end{cases}$
(27)

with ${\underline{c}}=\frac{1}{2}(1,1,1)^{T}$ and $\omega=\frac{1}{3}$. Note
that the initial condition is clearly $C^{\infty}$ for
$\left|x-c\right|<\omega$; this is the multidimensional analog of eq. 24 .

3D Strong Scaling Study
---
tasks | mesh/task | subdomain | total cores | elements/core
$1$ | $48\times 48\times 48$ | $\left[0.000,1.000\right]^{2}$ | $56$ | $1974.9$
$27$ | $16\times 16\times 16$ | $\left[0.000,0.33\bar{3}\right]^{2}$ | $1512$ | $73.1$
$216$ | $8\times 8\times 8$ | $\left[0.000,0.16\bar{6}\right]^{2}$ | $12096$ | $9.1$
$512$ | $6\times 6\times 6$ | $\left[0.000,0.125\right]^{2}$ | $28672$ | $3.9$
Table 5: Strong scaling study for 3D RKDG and RIDG. We scale the methods from
a fixed physical domain of $\left[0,1\right]^{3}$ to 512 equally sized
subdomains. Each task, which corresponds to one node, solves one subdomain.
Above we present the total number of tasks, the mesh size for each
subdivision, the physical size of each subdomain, the total number of cores,
and the number of mesh elements per core. We see that this experiment scales
the problem granularity down to about 4 elements per core. We also display the
errors for each method scaled. These runs were performed on Livermore
Computing clusters at at Lawrence Livermore National Laboratory using Intel
Xeon CLX-8276L CPUs with Omni-Path interconnects.  3D RKDG:
(${M_{\text{deg}}}=3$, $\nu=0.05$,  $L^{2}\,\text{\bf
error}\,\eqref{eqn:l2relerror}={1.36}\mathrm{e}{-}07$)
---
cores | dof/c |  efom eq. 12 | comms | time | eq. 20 | eq. 21 | eq. 17
$56$ | $126390.9$ | $192^{3}$ | ${2.59}\mathrm{e}{+}5$ | ${4.92}\mathrm{e}{+}3$ | $--$ | $--$ | $4.90$
$1512$ | $4681.1$ | $64^{3}$ | ${7.00}\mathrm{e}{+}6$ | ${2.15}\mathrm{e}{+}2$ | $22.86\times$ | $84.1\%$ | $6.26$
$12096$ | $585.1$ | $32^{3}$ | ${5.60}\mathrm{e}{+}7$ | ${2.89}\mathrm{e}{+}1$ | $170.53\times$ | $78.8\%$ | $7.13$
$28672$ | $246.9$ | $24^{3}$ | ${1.33}\mathrm{e}{+}8$ | ${1.54}\mathrm{e}{+}1$ | $319.98\times$ | $62.4\%$ | $7.41$
3D RIDG: (${M_{\text{deg}}}=3$, $\nu=0.6$,  $L^{2}\,\text{\bf
error}\,\eqref{eqn:l2relerror}={1.53}\mathrm{e}{-}07$)
cores | dof/c |  efom eq. 12 | comms | time | eq. 20 | eq. 21 | eq. 17
$56$ | $126390.9$ | $192^{3}$ | ${4.32}\mathrm{e}{+}3$ | ${8.24}\mathrm{e}{+}3$ | $--$ | $--$ | $4.72$
$1512$ | $4681.1$ | $64^{3}$ | ${1.17}\mathrm{e}{+}5$ | ${3.28}\mathrm{e}{+}2$ | $25.09\times$ | $92.7\%$ | $6.13$
$12096$ | $585.1$ | $32^{3}$ | ${9.33}\mathrm{e}{+}5$ | ${4.43}\mathrm{e}{+}1$ | $186.17\times$ | $86.1\%$ | $7.00$
$28672$ | $246.9$ | $24^{3}$ | ${2.21}\mathrm{e}{+}6$ | ${2.19}\mathrm{e}{+}1$ | $377.10\times$ | $73.6\%$ | $7.31$
3D RIDG: (${M_{\text{deg}}}=5$, $\nu=0.6$,  $L^{2}\,\text{\bf
error}\,\eqref{eqn:l2relerror}={4.97}\mathrm{e}{-}11$)
cores | dof/c |  efom eq. 12 | comms | time | eq. 20 | eq. 21 | eq. 17
$1512$ | $15798.9$ | $96^{3}$ | $1.17{e5}$ | $2.76{e4}$ | $--$ | $--$ | $7.64$
$12096$ | $1974.9$ | $48^{3}$ | $9.33{e5}$ | $3.76{e3}$ | $7.34\times$ |  $90.6\%$* | $8.51$
$28672$ | $833.1$ | $36^{3}$ | $2.21{e6}$ | $1.57{e3}$ | $17.54\times$ |  $92.1\%$* | $8.88$
Table 6: Strong scaling study for the 3D advection eq. 27 as defined by table
5. We list the total number of cores, the number of degrees of freedom per
core, the equivalent first order mesh as described in eq. 12, the estimated
total number of MPI communications as computed by eq. 22, the runtimes in
seconds, the speedup over single node performance as computed by eq. 20, the
scaling efficiency from a single node eq. 21, and the solution eq. 17. In this
work we are most interested in the quality metric, as this combines error and
runtime into a single metric. In this table, we see that the high order RIDG
methods exhibit the highest quality. Note* the scaling efficiencies eq. 21 for
the ${M_{\text{deg}}}=5$ RIDG scheme are computed relative to the 27 node run
as opposed to the single node case.

To test the scalability of the RIDG method for such a problem, we consider the
strong scaling study defined in table 5 . As per our scaling studies in the
earlier sections, our goal is to shrink the number of degrees of freedom per
core. This maximizes the ratio of communication overhead to task-localized
work. In three spatial dimensions, the RIDG prediction stage is a four-
dimensional problem and thus there is plenty of task-local work to hide
latency, as higher dimensionality exponentially increases the number of
elements in each DG basis. For RIDG, this greatly increases the cost of the
region matrix inversions. Despite this our predicted quality metric eq. 18
increases as we consider higher order methods. This is demonstrated in the
results shown in table 6 . These scaling studies were performed at the
Livermore Computing Center at Lawrence Livermore National Laboratory using
Intel Xeon E5-2695 v4 2.1GHz CPUs with Omni-Path interconnects. The entries of
the Jacobian eq. 8 are calculated using the quasi-quadrature-free approach
sketched in section 3.

First we notice that both the RKDG and RIDG methods exhibit better scaling
efficiency for this 3D problem than they do for the 2D problem. The
computational complexity of DG increases the ratio of MPI-communication and
OpenMP overhead to parallelized work, which also increases the the
computational work available to hide the communication latencies, despite the
increased number of communications per time-step. This is indicated by the
“degrees of freedom per core" metric eq. 19 shown in table 6 compared to the
same metric in table 4 . Even so, the RKDG method requires many more total
halo-region MPI communications than the RIDG method of any order. This is due
to the fact that the RIDG methods of ${M_{\text{deg}}}=3,5$ exhibit the
enhanced CFL restriction $\nu=0.6$, and thus all can take a relatively large
timestep compared to the RKDG method. Furthermore, we notice that the total
runtime of the RKDG method is half of that of the RIDG method of the same
${M_{\text{deg}}}$. Since their errors as per table 6 are similar in
magnitude, this implies the RKDG method ${M_{\text{deg}}}=3$ has a higher
quality than the RIDG method ${M_{\text{deg}}}=3$. For the higher order RIDG
method, ${M_{\text{deg}}}=5$, we observe that the runtimes increase by a
factor of 2 orders of magnitudes when compared to the methods of
${M_{\text{deg}}}=3$. However, the errors as seen in table 5 are 3 orders of
magnitude smaller for ${M_{\text{deg}}}=5$. Thus, we find that the higher
order methods have a higher quality metric than the methods of
${M_{\text{deg}}}=3$, and are thus more efficient. This is achieved in part by
the ability of the RIDG method of ${M_{\text{deg}}}=5$ to maintain a strong
scaling efficiency of around $90\%$ above 28,000 cores, as the lower order
schemes begin to drop off in efficiency. Note that in 3D the quality metric
for ADER-DG methods are approximately $+0.12$ higher than the quality metric
for RKDG Dumbser2018 . We conclude that due to their high quality metric, the
high order RIDG method with ${M_{\text{deg}}}=5$ scales very efficiently for
3D problems.

##  7 Conclusion

In this paper we have explored the scalability of the regionally implicit
discontinuous Galerkin (RIDG) method, which was previously shown to have
favorable convergence properties for solving hyperbolic conservation laws
RIDG_paper_2019 . We compared the results to that of the extremely popular
strong-stability-preserving Runge-Kutta DG (SSP-RKDG) method. We performed
these comparisons for a two-dimensional linear problem, a two-dimensional
nonlinear problem, and a three-dimensional linear problem which served as a
toy problem related to the Relativistic Vlasov Maxwell system. We demonstrated
efficient strong-scaling properties of the RIDG DG methods, we are able to
maintain vertex-neighbor stenciling and communications while taking highly
enhanced time-step sizes. We also demonstrated the incredible boost in
efficiency of using the quasi-quadrature free strategy for space-time matrix
assembly. This strategy reduces the computational complexity of space-time
matrix assembly down to what we would expect from purely spatial matrix-
assembly. This is critical to the viability of RIDG methods.

The work used to hide the aforementioned communication and collective
latencies itself benefits from intra-node (shared memory) parellelism via
OpenMP openmp ; Architecture2018 . The RIDG method achieves decent intra-node
efficiency simply by distributing the “regions” over the number of shared
memory compute cores, similar to the ADER-DG parellelism strategy Fambri2018 .
The result is that in the limit of minimum granularity (minimum cells per
core), the work done by each core is dominated by formation of the space-time
region Jacobian and small linear system solves. In contrast, the RKDG method
distributes the computation of the flux quadrature and the volume integral
quadrature to form a residual. In the the limit of minimum granularity, this
means that each CPU core computes the spatial component of the residual of a
cell. Thus, the RIDG method distributes a far greater load of work compared to
the RKDG method, which leads to superior intra-node scaling, at the cost of
incredibly increased computational cost per time step.

Further, the superior stability properties of the RIDG method allowed larger
time steps than what the RKDG method was able to take. The predictor-corrector
strategy of the RIDG method means that the method communicates twice per time-
step compared to once per stage (many times per time-step). Thus, the RIDG
method communicates across domain decomposition pseudo-boundaries (halo
regions) nearly two orders of magnitudes fewer times than the RKDG method. In
addition, the amount of work used to hide the latency of the communications is
much greater for the RIDG method than for the RKDG method. In the nonlinear
case, this advantage is diminished as the RIDG method must communicate once
per Newton iteration during the prediction step, but the RIDG method still
exhibits a superior ability to hide communication latency behind the Newton
iteration matrix solves.

Furthermore, the reduced number of time-steps means that there are fewer all-
to-all collective communications of the maximum wavespeed, an operation
required of explicit time integrators. The RIDG method is able to compute the
maximum wavespeed after the prediction step, and is thus able to hide the
latency associated with this collective using the corrector step.

Lastly, reduction in parallel scaling efficiency is inevitable, as a growing
percentage of each subdomain must be copied into contiguous MPI buffers and a
growing number of compute resources must be synchronized for maximum
efficiency. However, if the ratio of time spent updating the solution over the
time spent performing these copies is maximized, as is with RIDG, then higher
parallel scaling can be maintainted.

Due to excellent inter-node and intra-node parallelism, the RIDG method offers
excellent scalability for explicit time-stepping of hyperbolic conservation
laws, and is able to demonstrate incredible performance with very high-order
methods at high dimensionality.

###### Acknowledgements.

This work was performed under the auspices of the U.S. Department of Energy by
Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344. The
authors would like to thank Philipp Grete and Forrest Glines for their helpful
conversations.

## Disclaimer

This document was prepared as an account of work sponsored by an agency of the
United States government. Neither the United States government nor Lawrence
Livermore National Security, LLC, nor any of their employees makes any
warranty, expressed or implied, or assumes any legal liability or
responsibility for the accuracy, completeness, or usefulness of any
information, apparatus, product, or process disclosed, or represents that its
use would not infringe privately owned rights. Reference herein to any
specific commercial product, process, or service by trade name, trademark,
manufacturer, or otherwise does not necessarily constitute or imply its
endorsement, recommendation, or favoring by the United States government or
Lawrence Livermore National Security, LLC. The views and opinions of authors
expressed herein do not necessarily state or reflect those of the United
States government or Lawrence Livermore National Security, LLC, and shall not
be used for advertising or product endorsement purposes.

## References

  * (1) Atkins, H.L., Shu, C.W.: Quadrature-free implementation of discontinuous Galerkin method for hyperbolic equations.  AIAA Journal 36(5), 775–782 (1998).  DOI 10.2514/3.13891.  URL http://arc.aiaa.org/doi/10.2514/2.436
  * (2) Chan, J., Wang, Z., Modave, A., Remacle, J.F., Warburton, T.: GPU-accelerated discontinuous Galerkin methods on hybrid meshes.  J. Comput. Physics 318(div), 142–168 (2016).  DOI 10.1016/j.jcp.2016.04.003
  * (3) Cockburn, B., Shu, C.W.: The Runge–Kutta discontinuous Galerkin method for conservation laws V: Multidimensional systems.  J. Comput. Physics 141(2), 199–224 (1998)
  * (4) Dagum, Leonardo and Menon, R.: OpenMP: an industry standard API for shared-memory programming.  Computational Science & Engineering, IEEE 5, 46–55 (1998)
  * (5) Dumbser, M., Fambri, F., Tavelli, M., Bader, M., Weinzierl, T.: Efficient implementation of ADER discontinuous Galerkin schemes for a scalable hyperbolic PDE engine.  Axioms pp. 1–26 (2018).  DOI 10.3390/axioms7030063.  URL http://arxiv.org/abs/1808.03788
  * (6) Dumbser, M., Munz, C.D.: Building blocks for arbitrary high order discontinuous Galerkin schemes.  J. Sci. Comput. 27, 215–230 (2006)
  * (7) Dumbser, M., Zanotti, O., Hidalgo, A., Balsara, D.S.: ADER-WENO finite volume schemes with space-time adaptive mesh refinement.  J. Comput. Physics 248, 257–286 (2013).  DOI 10.1016/j.jcp.2013.04.017.  URL http://dx.doi.org/10.1016/j.jcp.2013.04.017
  * (8) Engsig-Karup, A.P., Eskilsson, C., Bigoni, D.: A stabilised nodal spectral element method for fully nonlinear water waves.  J. Comput. Physics 318, 1–21 (2016).  DOI 10.1016/j.jcp.2016.04.060
  * (9) Fambri, F., Dumbser, M., Köppel, S., Rezzolla, L., Zanotti, O.: ADER discontinuous Galerkin schemes for general-relativistic ideal magnetohydrodynamics.  Mon. Not. R. Astron. Soc. 477(4), 4543–4564 (2018).  DOI 10.1093/mnras/sty734
  * (10) Fischer, P.F., Heisey, K., Min, M.: Scaling Limits for PDE-Based Simulation.  In: 22nd AIAA Computational Fluid Dynamics Conference, 22-26 June 2015, Dallas, TX, pp. 1–10 (2015)
  * (11) Gabriel, E., Fagg, G.E., Bosilca, G., Angskun, T., Dongarra, J.J., Squyres, J.M., Sahay, V., Kambadur, P., Barrett, B., Lumsdaine, A., Castain, R.H., Daniel, D.J., Graham, R.L., Woodall, T.S.: Open MPI: Goals, concept, and design of a next generation MPI implementation.  Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics) 3241, 97–104 (2004).  DOI 10.1007/978-3-540-30218-6-19
  * (12) Gassner, G., Dumbser, M., Hindenlang, F., Munz, C.D.: Explicit one-step time discretizations for discontinuous Galerkin and finite volume schemes based on local predictors.  J. Comput. Physics 230, 4232–4247 (2011)
  * (13) Gottlieb, S., Shu, C.W.: Total variation diminshing Runge-Kutta schemes.  Math. of Comput. 67, 73–85 (1998)
  * (14) Gottlieb, S., Shu, C.W., Tadmor, E.: Strong stability-preserving high-order time discretization methods.  SIAM Rev. 43(1), 89–112 (2001)
  * (15) Gupta, M., Narasimhan, S.G.: Legendre polynomials Triple Product Integral and lower-degree approximation of polynomials using Chebyshev polynomials.  Tech. Rep. CMU-RI-TR-07-22, Carnegie Mellon University (2007)
  * (16) Guthrey, P.: Regionally implicit discontinuous Galerkin methods for solving the relativistic Vlasov-Maxwell system.  Ph.D. thesis, Iowa State University (2017)
  * (17) Guthrey, P., Rossmanith, J.: The regionally implicit discontinuous Galerkin method: Improving the stability of DG-FEM.  SIAM J. Numer. Analysis 57(3), 1263–1288 (2019).  DOI 10.1137/17M1156174
  * (18) INTERTWinE Consortium: Best Practice Guide to Hybrid MPI + OpenMP Programming (2017).  URL http://www.intertwine-project.eu/sites/default/files/images/INTERTWinE_Best_Practice_Guide_MPI%2BOpenMP_1.1.pdf
  * (19) Klaij, C., van Der Vegt, J., van Der Ven, H.: Space-time discontinuous Galerkin method for the compressible Navier-Stokes equations.  J. Comput. Physics 217(2), 589–611 (2006)
  * (20) OpenMP Architecture Review Board: OpenMP Application Programming Interface (2018).  URL https://www.openmp.org/wp-content/uploads/OpenMP-API-Specification-5.0.pdf
  * (21) Qiu, J., Dumbser, M., Shu, C.W.: The discontinuous Galerkin method with Lax-Wendroff type time discretizations.  Comput. Methods Appl. Mech. Eng. 194, 4528–4543 (2005).  DOI 10.1016/j.cma.2004.11.007
  * (22) Rusanov, V.: Calculation of interaction of non-steady shock waves with obstacles.  J. Comp. Math. Phys. USSR 1, 267–279 (1961)
  * (23) Sudirham, J., van Der Vegt, J., van Damme, R.: Space-time discontinuous Galerkin method for advection-diffusion problems on time-dependent domains.  Appl. Numer. Math. 56(12), 1491–1518 (2006)
  * (24) Zanotti, O., Fambri, F., Dumbser, M., Hidalgo, A.: Space–time adaptive ADER discontinuous Galerkin finite element schemes with a posteriori sub-cell finite volume limiting.  Computers & Fluids 118, 204 – 224 (2015).  DOI https://doi.org/10.1016/j.compfluid.2015.06.020