remarkRemark hypothesisHypothesis claimClaim Scalable Local Timestepping on
Octree GridsMilinda Fernando, Hari Sundar

# Scalable Local Timestepping on Octree Grids††thanks: Submitted to the
editors August 17, 2020. This work was supported by the National Science
Foundation(NSF) grant PHY-1912930, and the National Aeronautics and Space
Administration (NASA) grant 80NSSC20K0528. This work used computing resources
of the Extreme Science and Engineering Discovery Environment (XSEDE)
allocation TG-PHY180054 and Frontera pathway allocation PHY20033.

Milinda Fernando School of Computing, University of Utah. ().
<EMAIL_ADDRESS>Hari Sundar School of Computing, University of Utah. ().
<EMAIL_ADDRESS>
###### Abstract

Numerical solutions of hyperbolic partial differential equations(PDEs) are
ubiquitous in science and engineering. Method of lines is a popular approach
to discretize PDEs defined in spacetime, where space and time are discretized
independently. When using explicit timesteppers on adaptive grids, the use of
a global timestep-size dictated by the finest grid-spacing leads to
inefficiencies in the coarser regions. Even though adaptive space
discretizations are widely used in computational sciences, temporal adaptivity
is less common due to its sophisticated nature. In this paper, we present
highly scalable algorithms to enable local timestepping (LTS) for explicit
timestepping schemes on fully adaptive octrees. We demonstrate the accuracy of
our methods as well as the scalability of our framework across 16K cores in
TACC’s Frontera. We also present a speed up estimation model for LTS, which
predicts the speedup compared to global timestepping(GTS) with an average of
$0.1$ relative error.

###### keywords:

spacetime adaptivity, PDEs, local timestepping, octree grids

68U01, 68W10, 68U05,65M06

##  1 Introduction

The numerical solution of hyperbolic Partial differential equations (PDEs)
plays an important role in science in engineering, with a wide range of
applications from modeling earthquakes [8] to simulating gravitational waves
[24, 22, 33]. These show up commonly as initial value problems that are
typically solved using the method of lines by first discretizing in space and
then solving the resulting set of ordinary differential equations (ODEs). This
is commonly done using timestepping schemes, with explicit timestepping
methods such as the Runge-Kutta methods [30] being more common for evolving
hyperbolic systems. Additionally, these systems are characterized by the need
for high levels of spatial adaptivity [43, 22, 24]. High levels of adaptivity
impose severe stability restrictions for the explicit timestepping schemes
popular for solving such systems. Therefore, it is common—especially for large
scale distributed memory codes—to use a global (everywhere in space) timestep,
which is dictated by the smallest grid resolution in space [21]. This is
highly inefficient for large systems with several orders of magnitude
difference between the finest and coarsest grid resolution, as the coarser
regions are forced to take extremely small timesteps than would be needed for
stability [29, 38]. Local timestepping schemes can greatly speed up such codes
by ensuring that the adaptivity in space is matched by a corresponding
adaptivity in time. In this work, we develop an efficient, scalable local
timestepping scheme (LTS) for explicit single and multi-stage timestepping
algorithms on octree-adaptive spatial grids. We demonstrate the efficacy of
our scheme using linear and non-linear wave equations on up to $16K$ processes
on TACC’s Frontera. We also present a model to estimate the expected speedup
from using our scheme with an average of $10\%$ relative error in the
estimation of speedup. This provides a reliable way to determine in which
cases the use of local timesteping can be beneficial.

Our framework allows space adaptivity via the use of octrees (quadtrees in
$2D$) and uses high-order finite difference (FD) methods for space
discretization. Octree based adaptive space discretizations [18, 31, 27, 39,
46] are popular in large scale simulations because of their quasi-structured
nature, allowing for efficient and scalable data-structures and algorithms.
Our framework targets applications in computational relativity[24] and uses a
3+1 decomposition of spacetime operators, to compute a space slice where time
is constant and uses standard time integration methods (explicit) to evolve
the space slice forward in time. The time integration method of choice has
been the Runge-Kutta methods because of the larger stability region and the
availability of low-storage versions. For the specific problem of simulating
binary black hole mergers to estimate gravitational waves, the spatial grid is
characterized by $12$-$19$ levels of adaptivity and by the need to evolve for
extremely long times. These results are simulations that need to run for
months on thousands of processes. The use of LTS for these simulations can
provide 10-70x speed up, which can greatly reduce the time and cost of
obtaining gravitational waveforms.

While the theoretical aspects of LTS have been an area of active research in
recent years [36, 37, 32, 42, 5], performing LTS in a distributed computing
environment comes with additional challenges. A central bottleneck to
scalability with LTS is the variability in computational loads for different
regions of space based on their spatial adaptivity, as finely refined regions
take exponentially more timesteps than coarser regions. This requires
partitioning approaches that can account for such variable workloads.
Additionally, since in our target applications, the meshes are dynamic,
frequent re-partitioning is required, requiring fast partitioning algorithms
that are able to adapt to the variable computational load. While the use of
graph partitioning approaches [42] is likely to produce superior partitions,
the cost of partitioning our meshes in parallel makes the approach infeasible.
The key contributions presented in the paper can be summarized as follows.

  * •

Scalable LTS on octrees: We present a scalable LTS framework for multi-stage
explicit methods, on 2:1 balanced octree grids. While several community octree
frameworks are available[16, 49, 7, 48], to the best of our knowledge, they
are limited to space adaptivity.

  * •

Load balancing in LTS: The number of local timesteps needed to reach the
coarsest time on the grid depends on the spatial adaptivity. This leads to
load balancing issues in LTS. In order to resolve this, we propose an space
filling curve (SFC)-based weighted partitioning scheme. Compared to
traditional SFC-based partitioning, we compute the weighted length of the
curve to achieve a balanced load for LTSpartitions.

  * •

Low overhead of LTS compared to GTS: We present both strong and weak scaling
results for LTS and GTS approaches on octrees. These results demonstrate that
local block synchronization in LTS followed by time interpolations have a
lower cost compared to global block synchronization present in GTS approach.

  * •

Accuracy of LTS: We conduct numerical experiments to demonstrate the
correctness of the implemented LTS framework. The presented numerical results
demonstrate the accuracy of the LTS framework for both linear and non-linear
problems.

  * •

LTS performance model: We present an analytical performance modelto estimate
the speed up of LTS over the GTS scheme. The analytical model is extended to
compute a theoretical upper bound for the speed up that can be achieved for a
given spatial adaptivity structure.

Organization of the paper: The rest of the paper is organized as follows. In
§2, we give a brief motivation on the importance of LTS and a quick overview
of the existing state-of-the-art approaches in the field. In §3, we present
the algorithms and methods developed in detail to compare its efficiency to
GTS. In §4, we discuss the experimental setup, demonstrate strong and weak
scalability of our approach, and the accuracy of our scheme. In §5, we
conclude with directions for future work.

##  2 Background

In comparison to spatial adaptivity, local time adaptivity is less frequently
used by large-scale applications. Local timestepping requires additional
corrections in mismatching regions in time using interpolations and/or
extrapolations. Depending on the nature of the differential operator, these
operations can lead to problems in stability[28] of the numerical scheme.
Recent LTS methods are influenced by the split Runge-Kutta(RK) methods[41],
where two ODE systems are integrated using different step sizes (one called
active with the smaller timestep and the other called latent with the larger
timestep size) on the grid. Corrections using interpolations are performed for
the interface between the two grid regions. A complete numerical analysis of
spacetime adaptive timestepping methods is complicated, but early work by
Berger provided the first mathematical analysis for adaptive schemes for the
wave equation [11, 12]. Algorithms presented in these papers, discuss two main
approaches, interpolation based and coarse-mesh approximations. In the
interpolation methods, the solution at the coarse mesh is used to interpolate
the values needed for the finer mesh. In the coarse mesh approximation, the
coarse mesh is used to take a pseudo timestep that is used by the finer mesh.
In [20], the authors present methods for energy-conserving corrections in time
for Maxwell’s equations. There is a rich literature of LTS for discontinuous
Galerkin methods with special focus on energy-conserving time correction
operators [45, 34, 38, 35] which are important for complicated non-linear
spacetime differential operators. As mentioned previously, to enable LTS in a
distributed parallel setting requires specialized partitioning methods to
ensure load balance. Dynamic load balancing for adaptive mesh refinement is an
active research area [9, 23, 19, 26]. Sophisticated hypergraph partitioning
techniques have been used for LTS for the wave equation [42], but the cost of
partitioning makes it prohibitively expensive for AMR applications requiring
frequent re-meshing and therefore re-partitioning.

AMR in space and spacetime is an area of active research area. Here we present
a brief overview of AMR packages that focus on both space and spacetime.Block-
structured or patch-based AMR is widely used in the astrophysics and
computational fluid dynamics (CFD) communities. In block-structured AMR, the
adaptivity structure is predetermined and evolved during the simulation
appropriately. Some codes support LTS with block AMR[15, 22, 6], primarily
based on the Berger-Oliger AMR criteria[13]. Berger-Oliger AMR criteria
supports adaptivity in space and time, but requires strong constraints on the
structure of the adaptivity, such as all grids at level $l+1$(child grids)
should be entirely contained within the grids at level $l$ (parent grids)
while grid at the same level may overlap. There exist other block-AMR
codes[2], which only supports space adaptivity and no adaptivity in time.

Another commonly used approach for large-scale AMR is octree-based AMR[16, 49,
7, 48]. In octree-based AMR, the adaptive grid is represented using quadtrees
and octrees. Unlike block-based AMR, octree-based AMR has relaxed constraints
on refinement, providing highly adaptive quasi-structured (point-structured)
grids in space. To the best of our knowledge, currently available octree-AMR
codes are limited only to space adaptivity with no support to enable
adaptivity in time. In this paper, our main focus is to perform large scale
LTS on adaptive octree grids, we choose the simple time interpolation methods
presented in the papers [36, 37] which form the mathematical basis for the
methods presented in this work. Additional details and analysis can be found
in these papers on the LTS scheme used in this paper.

\̃\

##  3 Methodology

###  3.1 $3+1$ decomposition of PDEs

In this paper, we focus on differential operators defined on the traditional
spacetime that is a 4d manifold (3 space + 1 time dimension). Let
$\mathcal{L}$ be a differential operator (linear or non-linear),
$\mathcal{L}:X\rightarrow Y$ where $X,Y$ are appropriate functional spaces
which $\mathcal{L}$ acts upon. For example, when
$\mathcal{L}\equiv\partial_{t}-\partial^{2}_{xx}$, we get the heat operator or
when $\mathcal{L}\equiv\partial^{2}_{tt}-\partial^{2}_{xx}$ one attains a
linear wave operator. Throughout this paper, we focus only on the operators
$\mathcal{L}$, which can be transformed into evolution equation of the form
(1), which we refer to as $3+1$ decomposition of operator $\mathcal{L}$

(1) $\partial_{t}u=F(t,u(t))$

where, for $T\in\mathcal{R}^{+}$, $F:[0,T]\times W\rightarrow X$, and
$W\subset X$. For a given $s\in[0,T)$ and $\phi\in W$, the solution or the
integral curve of $F$ with respect to $s,\text{ and }W$ is a map where the
range is $C^{0}([s,T]\times W)\cap C^{1}([s,T]\times X)$ that satisfies (1) on
$[s,T]$ and $u(s)=\phi$. Analysis of the well-posed nature of these integral
curves, are out of the scope of this paper, hence we assume these integral
curves are well-posed, and can be computed with numerical timestepping, with
the appropriate necessary stability constraints.

###  3.2 Adaptivity and parallelism in space

While this paper is centered on adaptivity in time, it is closely related to
our realization of adaptivity in space using octrees. We give a brief overview
of our spatial adaptivity framework in this section. The framework is freely
available via an MIT license[25] and additional details on our algorithms can
be found in [23, 24, 47].

For a given spatial domain $\Omega=[0,2^{L}]^{3}$ where $L$ is the maximum
depth parameter, we use the octree data structure to represent spatial
discretization. Octrees are widely used [44, 17, 49, 10] in computational
sciences for its simplicity, efficient data-structures, ease of partitioning,
and parallel scalability. For a given function
$f:\Omega\rightarrow\mathcal{R}$, we use axis-aligned octrees to generate our
space discretization. The adaptive mesh refinement (AMR) criteria can be
application specific, in this paper, we focus on the wavelet based AMR schemes
described in [24]. To perform numerical computations on adaptive octrees, we
need to have neighborhood information for elements (octants) as well as at the
nodal level.

####  3.2.1 Octree partitioning

The problem size or the load varies significantly during octree construction,
balancing, and meshing during the simulation. This necessitates the efficient
partitioning of the octree. As is common, we use a Hilbert curve based
partitioning scheme similar to [23], but with minor modifications.

####  3.2.2 Octree Construction and Refinement

The octree construction is based on expanding user-specified functions in
accordance with the specified AMR criteria, but for the work presented in this
paper, we use the wavelet transform and truncate the expansion (i.e. stopping
the refinement at that level) when the coefficients are smaller than a user-
specified tolerance $\epsilon>0$. Intuitively, the wavelet coefficient
measures the failure of the field to be interpolated from the coarser level.
In distributed memory, all processes start from the root and refine until at
least $p^{2}$ octants are produced. These are equally partitioned across all
processes. Subsequent refinements happen in an element-local fashion and are
embarrassingly parallel. A re-partition is performed at the end of
construction to ensure load balance.

####  3.2.3 2:1 Balancing

Following the octree construction, we enforce a 2:1 balance condition, i.e.,
any octant can have a neighbor that is either the same size, half as big, or
twice as big. This makes subsequent operations simpler without affecting the
adaptive properties. Our balancing algorithm is similar to existing approaches
for balancing octrees [14, 48] with minor changes in the choice of data
structures, and process-local balancing algorithms.

####  3.2.4 Mesh generation

By meshing, we refer to the construction of the data structures required to
perform $d^{th}$ order numerical computations octree data. One of the key
steps of the mesh generation phase is to construct neighborhood information
for octants. Primarily, there are two maps that are produced. The first map
o2o to determine the neighboring octants of a given octant and a map o2n to
compute the nodes corresponding to a given octant. o2o map is generated, by
performing parallel searches similar to approaches described in [23, 16] and
optimized as per the methods described in [24]. Assuming we have $n$ octants
per partition, these search operations and building required o2o and o2n data
structures can be performed in $\mathcal{O}(n\log(n))$ and $\mathcal{O}(n)$
complexity, respectively. Computation of the o2n map is a local operation,
where we start with octant local nodes( with duplicate nodes for shared octant
boundaries), and eliminating duplicate nodes with a globally well defined
criteria such as space filling curve (SFC) ordering operator(see figure 1).

Figure 1:  A $2D$ example of octant local nodes (in the center) and shared
octant nodes (the rightmost figure) nodal representation (for octant order of
$2$) of the adaptive quadtree shown in the leftmost figure. Note that in
octant local nodes, representation nodes are local to each octant and contain
duplicate nodes. By removing all the duplicate and hanging nodes by the rule
of nodal ownership, we get the shared octant nodes representation. Note that
the nodes are color-coded based on the octant level.

####  3.2.5 Octree to block decomposition

For a given octree $\mathcal{T}$, we compute a compressed octree of
$\mathcal{T}$, denoted as $\mathcal{B}$ (blocks) where each leaf node in
$\mathcal{B}$ is a node in $\mathcal{T}$ with uniformly refined sub-octree.
Octree to block decomposition allows, to perform finite difference (FD)
computations on adaptive octrees, as well as simplifies the computation of
data structures required to perform LTS (see figure 2).

octree to blocks $(b_{1},0)$$(b_{2},1)$$(b_{3},0)$$(b_{4},1)$

Figure 2:  A simplistic example of octree ($\mathcal{T}$) to block
decomposition. The left figure shows the considering adaptive octree, and its
block decomposition is shown in the right. Each block $(b_{k},l)$ is
associated with uniform grid level parameter $l$, where $l$ denotes the level
of refinement of the sub-tree rooted at $b_{k}$ node in $\mathcal{T}$.

###  3.3 Finite difference computations on adaptive octrees.

We use the finite difference method (FDM) to discretize differential operators
in space. Performing FD computations on adaptive octrees require additional
processing, to ensure that each refinement level has neighboring points
available at the same resolution. This is achieved by performing octree to
block decomposition where for each block, we compute a padding region of a
specific width corresponding to the FD stencil. In order to compute the padded
blocks, we use the computed octant to octant (o2o) and octant to nodal (o2n)
information with the corresponding space interpolations. For example, between
finer and coarser grid blocks, we use coarser to finer interpolation to
compute the padding region for the finer block, while using finer to coarser
injections for the padding region of the coarser block. Note that in the
paper, the octant shared node representation is also referred to as zipped
representation (see figure 1) and block with padding computed referred to as
the unzipped representation. All the FD stencils are applied at the unzipped
representation, while just prior to the communication, we perform zip
operation so the inter-process communication happens in the efficient more
compact form (see figure 3).

Figure 3:  A simplistic example of octree to block decomposition and unzip
operation. The leftmost figure shows the considering adaptive octree with
shared octant nodes, and its block decomposition is shown in the middle. Note
that the given octree is decomposed into four regular blocks of different
sizes. The rightmost figure shows the decomposed blocks padded with values
coming from neighboring octants with interpolation if needed. In order to
perform unzip operation, both o2o and o2n mappings are used.

###  3.4 Explicit timestepping schemes

Explicit timestepping (ETS) is a class of numerical schemes that compute the
solution curve for (1). In explicit methods, the solution at time $t_{n+1}$,
$u^{n+1}\equiv u^{n+1}(t_{n+1},\cdot)$ is computed directly from the solution
at the previous timestep $u^{n}$ and does not require a linear solve. In order
to numerically evolve (1), we discretize $F$, i.e., discretization in space
first. The resulting set of ODEs are discretized using explicit time
integration. Depending on the properties on $\mathcal{L}$, there can be
additional constraints on $\Delta t,\Delta x$. For most hyperbolic operators,
the Courant–Friedrichs–Lewy (CFL) condition [21] is a necessary condition for
stability for numerical time evolution. The CFL condition $\frac{\Delta
t}{\Delta x}<C$, where C is a constant that specifies a necessary condition
for stability. Intuitively, it imposes the constraint that we cannot propagate
spatial information in time, faster than the speed of information propagation
defined by operator $\mathcal{L}$. Runge-Kutta (RK) [40] schemes are widely
used explicit timesteppers (see equation (2)), that will be our main focus for
spatially adaptive local timestepping.

$\displaystyle k_{1}$ $\displaystyle=F(u^{n})$ $\displaystyle k_{2}$
$\displaystyle=F(u^{n}+a_{2,1}k_{1})$ $\displaystyle\;\;\vdots$ $\displaystyle
k_{p}$ $\displaystyle=F(u^{n}+a_{p,1}k_{1}+...+a_{p,p-1}k_{p-1})$ (2)
$\displaystyle u^{n+1}$ $\displaystyle=u^{n}+\Delta
t(\sum_{i=1}^{p}w_{i}k_{i})$

Being one-step methods, they do not require a starting procedure, have a large
stability region, and from a large-scale computational perspective can easily
be converted into low-storage versions.

unzipGTS $t$ loopRK stagesglobal blk_syncLTS $t$ loopblock level $l$ loopRK
stages $B_{l}\subset B$={ $k\geq l_{max}$ to $l$}local blk_synczip

Figure 4: This figure illustrates the overview of time evolution using GTS and
LTS for multi-stage explicit timestepping, on octree grids. In order to
perform FD computations, the adaptive grid is decomposed into uniform block
patches with appropriate padding and spatial derivatives are evaluated on
equispaced block representation(unzipped) computed using previous timestep
solution $u_{n}$. The unzip operation results in a sequence of that which are
used to compute the solution on the internal block ($\blacksquare$), using the
padding values at the block boundary ($\blacksquare$). After time evolution,
the next timestep $u_{n+1}$ is projected back to sparse grid (zip)
representation. In GTS, for each explicit stage, we evolve all the blocks
using $\Delta t_{finest}$, followed by a global block synchronization
operation. This global synchronization operation consists of projection of
block local solution for zipped representation (shared octant nodes ),
followed by inter-process communication and project back to unzipped
representation. Note that shared octant nodes are used for inter-process
communication, since it is compact and does not contain node duplicate values.
In LTS, we have a block level loop which selects subset of blocks $B_{l}$
which are eligible to evolve, followed by the explicit stage loop. Once
$B_{l}$ is evolved, we perform block synchronization for $B_{l}$, hence this
is a local block synchronization. For this local block synchronization, for
inter-process communication we use octant local nodes representation
(unzipped) without padding region). Unlike GTS we cannot use shared octant
nodes representation, since blocks evolved are at different times.

###  3.5 GTS: Global timestepping

In this section, we describe how we perform global timestepping (GTS) on
adaptive octrees with FD computations (see figure 4). In GTS, the finest space
resolution will determine the timestep size for the entire grid, and the
entire domain will march in a synchronized fashion. In most cases, this is not
efficient especially when the overall percentage of finer octants across the
grid is low, which is the case for many computational science applications [1,
43, 24]. Let us assume we have the solution $U^{n}$ defined on the adaptive
grid. Then for each explicit stage, we unzip the intermediate timestep, loop
over all the blocks to compute the block internal using FD stencils, perform
zip operation and compute the next intermediate timestep(see algorithm 1). In
order to synchronize the block padding regions we need inter-process
communication, therefore before each unzip, we perform data exchange with
neighboring processes using the compact zipped representation. Note that all
blocks do not depend on the values from neighboring processes. We label these
as independent blocks, and those requiring values from other processes as
dependent blocks. We exploit this property to overlap the unzip computation
with communication, i.e., we initiate the communication before starting the
computation of padding for independent blocks, and upon receiving the data
from neighboring processes perform the unzip operation for dependent blocks.

Algorithm 1 Global timestepping (GTS)

1:$U^{n}$ previous timestep, $\Delta t$, $B$ sequence of blocks

2:$U^{n+1}=U(t^{n}\Delta t,\cdot)$

3:$u\\_unzip\leftarrow unzip(U^{n})$

4:for $s=1:k$ do

5: for $b\in B$ do

6: $k[b,s]\leftarrow compute\\_stage(u\\_unzip,k[1,..,s-1])$

7: $k[s]\leftarrow blk\\_sync(k[s])$

8:$U^{n+1}\leftarrow compute\\_step(U^{n},k)$

###  3.6 LTS: Local timestepping

We now describe how to enable LTS for explicit single and multi-stage
timesteppers on 2:1 balanced octree grids. Since $\Delta t=c\Delta x$, where
$c$ is the CFL constant, we know that timestep sizes between coarser ($\Delta
t_{c}$) and finer ($\Delta t_{f}$) grids are also 2:1 balanced, therefore,
$\Delta t_{c}=2\Delta t_{f}$.

####  3.6.1 Single stage explicit schemes

For single stage explicit schemes such as forward Euler, we can enable the
correction between coarser ($b_{c}$) and finer ($b_{f}$) blocks, by making
$b_{c}$ takes a pseudo timestep, for block $b_{f}$ (see figure 5). This
approach is simple for single stage timesteppers, since only a single layer of
coarser blocks that needs to take a pseudo timestep, as a correction for the
finer blocks. Extending this to multi-step schemes is complicated and loses
the local adaptive timestepping, since $N_{s}$ stage timestepper require
$N_{s}$ layer of coarser blocks to take pseudo timesteps to make stage
corrections. Therefore, we approach the multi-stage LTS differently (see
figure 4).

$t_{n}$$b_{1}$$b_{2}$$b_{3}$$b_{4}$$t_{n+\frac{1}{4}}$$t_{n+\frac{1}{2}}$$t_{n+\frac{3}{4}}$$t_{n+1}$

Figure 5: Simple illustration of local timestepping for single-stage schemes
on 2:1 balanced adaptive grids. Here we have a refined block $b_{1}$ with
neighbor $b_{2}$ that can be at most twice as large. We also consider block
$b_{3}$ that doesn’t have any neighbors smaller than itself. Note that 2:1
balancing ensures that these are the only cases that can exist for any
adaptive mesh. Block $b_{1}$ and $b_{3}$ take timesteps corresponding to the
size of the blocks. Block $b_{2}$ however takes $2\times$ as many timesteps
compared to $b_{3}$, first a half-step to help its neighbor $b_{1}$ take its
second timestep, and then a full size timestep to reach $t_{n+1}$.

####  3.6.2 Multi-stage explicit schemes

Algorithm 2 Local timestepping (LTS)

1:$U^{n}$ previous timestep, $\Delta t$, $B$ sequence of blocks

2:$U^{n+1}=U(t^{n}\Delta t,\cdot)$

3:$V\leftarrow U^{n}$

4:$u\\_unzip\leftarrow unzip(U^{n})$

5:for $l=l_{max}$ to $l_{min}$ do

6: $B_{l}\leftarrow compute\\_blk\\_subset(B,l)$

7: for $b\in B_{l}$ do

8: for $s=1:k$ do

9: $k[b,s]\leftarrow compute\\_stage(u\\_unzip,k[1,..,s-1])$

10: $k[s]\leftarrow blk\\_sync\\_local(k[s],B_{l})$

11: $V\leftarrow compute\\_step\\_partial(V,k,B_{l})$

12:$U(t^{n}+\Delta t_{coarsest},\cdot)\leftarrow V$

In this section, we present how we perform LTS for multi-stage explicit
schemes (see algorithm 2), on 2:1 balanced distributed octrees. The main
challenge extending the single-stage timestepping approach for multi-stage
schemes, is that the layers of blocks which have to perform pseudo timestep
increases with the number of stages in the timestepping scheme. In order to
avoid this, we need to decouple the stages $k_{i}$ from the timestep size used
to evolve the solution. The linear relation (see equation 3) between stages
$k_{i}$ and the $\partial^{i}_{t}U$ can be derived as described in [37, 36],
where $C$ is a lower triangular coefficient matrix, coefficients are derived
from $a_{ij}$ coefficients of the explicit scheme.

(3) $K=\begin{bmatrix}k_{1}\\\ k_{2}\\\ \vdots\\\ k_{p}\end{bmatrix}=P_{\Delta
t}\times C\times\begin{bmatrix}\partial_{t}U_{t=t^{n}}\\\
\partial_{t}^{2}U_{t=t^{n}}\\\ \vdots\\\
\partial_{t}^{p}U_{t=t^{n}}\end{bmatrix}$

where $P_{\Delta t}$ and $C$ are matrices defined as,

(4) $\small P_{\Delta t}=\begin{bmatrix}1&0&\cdots&0\\\ 0&\Delta t&\cdots&0\\\
\vdots&\vdots&\ddots&\vdots\\\ 0&0&\cdots&\Delta
t^{p-1}\end{bmatrix},C=\begin{bmatrix}1&0&\cdots&0\\\ 1&c_{2,2}&\cdots&0\\\
\vdots&\vdots&\ddots&\vdots\\\ 1&c_{p,2}&\cdots&c_{p,p}\end{bmatrix}$

Due to the 2:1 balancing, any given coarser and finer octant interfaces can
differ at most by 1 refinement level. Hence timestep size ratio between
coarser and finer levels are 2:1 balanced as well. Initially, we assume that
all the blocks are synchronized in time. Let’s consider adjacent finer block
$b_{f}$ and coarser block $b_{c}$ (see Figure 6). Initially, both blocks have
the solution at $t^{n}$. Both blocks compute the stage $k_{i}$ using their
corresponding timestep size $\Delta t_{f}$ and $\Delta t_{c}$, but we need to
apply a correction between the finer $K_{f}$ and coarser $K_{c}$ stages as per
(5), where $B_{dt}$ denotes an upper triangular matrix containing the
coefficients computed from the Taylor expansion of each stage [37].

(5) $\small K_{c|t=t^{n}}=P_{\Delta t_{c}}CB_{0}P^{-1}_{\Delta
t_{f}}C^{-1}K_{f|t=t^{n}}=M^{(1)}_{fc}K_{f|t=t^{n}}$

Note that during the first step of the $b_{f}$ correction is applied to both
coarse and finer blocks, after the first step of $b_{f}$, $b_{c}$ has already
reached the coarser timestep, hence for the second step of $b_{f}$ the
correction comes only from coarser grid (see equation (6) and Figure 6).

(6) $\small K_{f|t=t^{n}+1/2}=P_{\Delta t_{f}}CB_{\Delta_{f}}P^{-1}_{\Delta
t_{c}}C^{-1}K_{c|t=t^{n}}=M^{(2)}_{cf}K_{c|t=t^{n}}$

More details on the error and stability analysis of these correction operators
can be found in [37, 36].

$t$$\Delta t_{c}=2\Delta t_{f}$$\Delta
t_{f}$$b_{f}$$b_{c}$$K_{f|t=t^{n}}$$K_{c|t=t^{n}}$$K_{f|t=t^{n+1/2}}$$M^{(1)}_{fc}$
& $[M^{(1)}_{fc}]^{-1}$$[M^{(2)}_{cf}]$

Figure 6: Example of finer ($b_{f}$) to coarser ($b_{c}$) and coarser($b_{c}$)
to finer ($b_{f}$) grid corrections. $M^{(1)}_{fc}$ is finer to coarser
correction between fine step 1 and coarse step. $M^{(1)}_{fc}$ applied to
finer grid we get the stage correction to coarser grids, and
$[M^{(1)}_{fc}]^{-1}$ applied to coarser grid stages, we get the finer step 1
corrections. For the step 2 of the finer grid, we need the stage correction
from coarser which is computed by applying, $[M^{(2)}_{cf}]$ to coarser step.

####  3.6.3 Synchronization between blocks

The computed correction operators allows us to relax the data dependency
between adjacent blocks. The error caused by correction operators can be
written as $\mathcal{O}(\Delta t)^{p}$, where $p$ is the order of accuracy in
the explicit timestepping method. Therefore it is important to make sure that
these correction operators are computed between blocks where the time
difference is minimal. For a given octree $T$ let $l_{max}$ and $l_{min}$ be
the finest and the coarsest levels of refinement of the blocks. Then we
evolve, the finest block with frequency of 1, next finest with frequency of
$\frac{1}{2}$ and coarsest with the frequency of $\frac{1}{2^{\Delta L}}$
where $\Delta L=l_{max}-l_{min}$, until all the blocks reached the next
coarsest time (see figure 7). The above with 2:1 balancing the maximum time
gap between two adjacent blocks ($b_{f}$, $b_{f}$) is bounded by $\Delta
t_{f}$, the step-size of the finer block. Once all the blocks that can be
processed have been computed and stored in the octant local nodes
representation (see figure 1), we perform data exchange between partitions to
update ghosted values. Note that unlike GTS, we cannot use the shared octant
nodes representation for communication since different grid points can be at a
different time during LTS. Therefore we need duplicate degrees of freedoms at
octant level to perform LTS ghost synchronization. Following ghost
synchronization, we perform the block-wise correction as per (5) and (6).
These are applied and copied to the block padding region (see figure 6) of the
neighboring blocks.

#####  3.6.3.1 Partial block synchronization

At a given time level $l$, we evolve a subset of blocks $B_{l}$, hence at the
time of the inter-process synchronization (ghost exchange), we only need to
communicate data corresponding to $B_{l}$. We pre-compute maps for each time
level to perform partial ghost synchronization, to reduce the amount of data
communicated between processors. This is more efficient than the global block
synchronization performed in GTS, and in addition to the reduction of work by
reducing the number of timesteps taken by coarser regions, it also reduces the
amount of data synchronization needed.

time$1\Delta t$$2\Delta t$$3\Delta t$$4\Delta
t$$b_{h}$$b_{2h}$$b_{4h}$$l=0$$l=0$$l=0$$l=1$$l=2$$l=2$$l=3$

Figure 7: A simple illustration of block set $B=\\{b_{h},b_{2h},b_{4h}\\}$
that is evolving at each level $l$ of the time level loop. Assume in the
beginning, all the blocks are synchronized in time. At level $l=0$,
$B_{0}=\\{b_{h},b_{2h},b_{4h}\\}$, perform timestep of its corresponding
timestep size $\\{\Delta t,2\Delta t,4\Delta t\\}$. Note now the block
$b_{4h}$ already reached the coarsest time over $B$. Similarly, for $l=1$,
$B_{1}=\\{b_{h}\\}$ , for $l=2$, $B_{2}=\\{b_{h},b_{2h}\\}$ and for $l=3$,
$B_{3}=\\{b_{h}\\}$, are evolved with corresponding timestep sizes. After
$B_{3}$, all the blocks have reached the next coarsest timestep size over $B$.
Note that, appropriate time interpolation and corrections are applied across
different refinement regions (see figure 6).

####  3.6.4 Weighted partitioning

We use SFC based partitioning scheme to distribute work among the processors.
For each local partition $\tau_{k}$, the amount of work that each block has to
perform to reach to the global coarsest timestep depends on the block
refinement level. For example a block of level $l$ has to perform 2x timesteps
compared to its adjacent coarser block at level $l-1$. A partitioning scheme
with equal weights assigned to each octant will result in load imbalanced
partitions for LTS timestepping. In order to overcome that, we perform
weighted partitioning, where the relative weight of block increases with the
refinement level (see Figure 8). We modified our SFC-based partitioning scheme
to account for the specified weights of the blocks. A 3d SFC curve can be
considered to be an injective mapping between 1d domain to 3d octree domain.
Using SFC ordering, we can sort the octants which results in a linear ordering
of the octants, while ensuring spatial locality. Once ordered according to the
SFC, partitioning the domain reduces to partitioning a $1D$ curve. In weighted
SFC-partitioning, we use weighted length of the curve i.e. we aim for equal
$\sum_{e\in\tau_{k}}w_{e}$, where $w_{e}$ denotes the weight of the octant,
$\tau_{k}$ a given partition. Note that for GTS, we use $w_{e}=1,\forall
e\in\tau_{k}$.

$0$$5$$10$$15$$20$$25$$30$$0$$100$$200$$300$$400$processor id (rank)
$\rightarrow$number of octants $\rightarrow$octants$0$$5\cdot 10^{9}$$1\cdot
10^{10}$$1.5\cdot 10^{10}$$2\cdot 10^{10}$weight $\rightarrow$partition weight
Figure 8: A simple illustration of our weighted SFC-based partitioning scheme
for an adaptive octree partitioned across 32 processors. The bar plot shows
the number of octants each partition has, and the line plot represents the
total weight of each partition. Note that partitions with low octant counts
have highly refined regions; hence they need to perform a larger number of
timesteps to reach the coarser level timestep and vice versa, but the total
weight of each partition is roughly equal, i.e., all partitions perform
roughly the same amount of work.

###  3.7 GTS Vs. LTS: Approximating the speedup

Here we present theoretical bounds for the work performed by the GTS and LTS
timesteping approaches, on 2:1 balanced octrees. Let $L=(l_{max}-l_{min})$ be
the difference between the maximum and minimum refinement levels for a given
octree. Let $W=\\{\alpha_{0},...,\alpha_{L}\\}$, be the corresponding work for
a block sequence $B=\\{b_{l}\\}_{l=l_{min}}^{l_{max}}$, then for matching the
finest timestep sizes between GTS and LTS finest levels, the work performed by
the GTS and LTS schemes can be written as (8) and (7). Assuming constant time
to perform a work unit, for speedup $S$, $1/S$ can be written as (9). Since,
$\frac{\alpha_{0}}{|W|}<\frac{1}{S}$, maximum theoretical speedup that can
achieved, for a given block distribution can be written as
$S<\frac{|W|}{\alpha_{0}}$.

(7) $\displaystyle W_{lts}$ $\displaystyle=\sum_{l=0}^{L}2^{L-l}\alpha_{l}$
(8) $\displaystyle W_{gts}$ $\displaystyle=2^{L}\sum_{l=0}^{L}\alpha_{l}$ (9)
$\displaystyle\frac{W_{lts}}{W_{gts}}$
$\displaystyle=\frac{1}{|W|}(\alpha_{0}+\frac{\alpha_{1}}{2}+...+\frac{\alpha_{L}}{2^{L}})<1\text{
Where, }$ (10) $\displaystyle|W|$ $\displaystyle=\sum_{l=0}^{L}\alpha_{l}$

##  4 Results

###  4.1 Experimental setup

The large scalability experiments reported in this paper were performed on
TACC’s Frontera supercomputer. Frontera is an Intel supercomputer at Texas
advanced computing center (TACC) with a total of 8,008 nodes, each consisting
of a Xeon Platinum 8280 (”Cascade Lake”) processor, with a total of 448,448
cores. Each node has 192GB of memory. The interconnect is based on Mellanox
HDR technology with full HDR (200 Gb/s) connectivity between the switches and
HDR100 (100 Gb/s) connectivity to the compute nodes.

###  4.2 Non-linear and linear wave propagation

In this section, we introduce a simple model to demonstrate LTS, the classical
wave equation. We write the classical wave equation in a form with first
derivatives in time and second derivatives in space. This allows us to easily
apply the specified, explicit schemes.

For a scalar function $\chi(t,x^{i})$, the classical wave equation in
Cartesian coordinates $(t,x,y,z)$ with a non-linear source term can be written
as,

(11) $\frac{\partial^{2}\chi}{\partial
t^{2}}-\left(\frac{\partial^{2}}{\partial x^{2}}+\frac{\partial^{2}}{\partial
y^{2}}+\frac{\partial^{2}}{\partial
z^{2}}\right)\chi=-c\frac{\sin(2\chi)}{r^{2}},$

where $r=\sqrt{x^{2}+y^{2}+z^{2}}$. We write the equation as first order in
time system by introducing the variable $\phi$ as

(12) $\displaystyle\frac{\partial\chi}{\partial t}$ $\displaystyle=\phi$ (13)
$\displaystyle\frac{\partial\phi}{\partial t}$
$\displaystyle=\left(\frac{\partial^{2}}{\partial
x^{2}}+\frac{\partial^{2}}{\partial y^{2}}+\frac{\partial^{2}}{\partial
z^{2}}\right)\chi-c\frac{\sin(2\chi)}{r^{2}}.$

For non-linear wave propagation results presented in this paper, we used
$c=1$. We choose outgoing radiative boundary conditions for this system [4].
We assume that the variables $\chi$ and $\phi$ approach the form of spherical
waves as $r\to\infty$, which decay as $1/r^{k}$. The radiative boundary
conditions then have the form

(14) $\frac{\partial f}{\partial t}=\frac{1}{r}\left(x\frac{\partial
f}{\partial x}+y\frac{\partial f}{\partial y}+z\frac{\partial f}{\partial
z}\right)-k(f-f_{0}),$

where $f$ represents the functions $\chi$ and $\phi$, and $f_{0}$ is an
asymptotic value. We assume $k=1$ for $\chi$ and $k=2$ for $\phi$.

For the linear wave propagation results presented, we simply zero out the non-
linear source term (i.e. set $c=0$). The analytical solution for the 1d linear
wave operator in (11) with zero source term can be written as,

(15) $\chi(t,x)=\frac{f(x-t)+f(x+t)}{2},\text{ where }\chi(0,x)=f(x)$

Figure 9: Plots for linear wave propagation with LTS with a velocity vector
$(1,1,1)$ using a Gaussian pulse centered at $(0,0,0)$ as the initial
condition. Images shown from left to right and top to bottom in increasing
order of simulation time.

###  4.3 Accuracy

We conduct numerical experiments using linear and non-linear wave propagation
to test the accuracy of our methods and implementation. For linear wave
propagation, we compute the analytical solution for wave propagation in the
$x$ direction and compare the analytical solution with the computed solution
using global and local timestepping. For the above experiments, we used a
third-order RK scheme with increasing spacetime resolution. Figure 10 shows
numerical error for LTS and GTS approaches with increase resolution.

$0$$0.5$$1$$1.5$$2$$0$$2$$4$$6$$\cdot 10^{-7}$timeerror
$l_{2}$GTS(d=8)LTS(d=8)GTS(d=10)LTS(d=10) Figure 10: Discrete $l_{2}$ error
compared to the analytical solution for the 1d wave operator in 3d for GTS and
LTS timestepping. For this experiment we used maxdepth of 8 and 10 with
refinement trigger tolerance of $10^{-5}$.

Since the computation of the analytical solution for the 1d wave operator with
non-linear source term is complicated, we compare the $l_{\infty}$ norm
computed on the numerical difference between global and local evolved
timesteps. In table 1 we present the difference between the solution $\chi$
evolved using GTS and LTS for increasing maximum refinement level(maxdepth) 8
and 10. Again, both GTS and LTS are in agreement to machine precision. This
demonstrates that we can use LTS in lieu of GTS without sacrificing accuracy
or stability for both linear as well as non-linear problems.

Time(s) | $||\chi_{lts}-\chi_{gts}||_{\infty}$
---|---
|  maxdepth=8 |  maxdepth=10
0.000000 |  0 |  0
0.130208 |  9.93E-39 |  2.22E-16
0.260417 |  6.92E-34 |  1.11E-15
0.390625 |  4.86E-30 |  1.50E-15
0.520833 |  6.20E-27 |  2.08E-15
0.651042 |  1.39E-17 |  5.11E-15
0.781250 |  5.55E-17 |  1.09E-14
0.911458 |  1.39E-16 |  9.69E-15
1.041670 |  2.09E-16 |  9.21E-15
Table 1: $l_{\infty}$ difference between GTS and LTS at corresponding
timesteps for non-linear wave propagation. For this experiment we used
maxdepth of 8 and 10 with refinement trigger tolerance of $10^{-5}$. Note that
for maxdepth 8, 8 global timesteps and for maxdepth 10, 32 global timesteps
were equivalent for a single LTS step.
8163264128256512102420484096819216384$0$$2$$4$$6$number of cores
$\rightarrow$time (s)
$\rightarrow$blk_sync(LTS)time_interp.rhs(LTS)comm(LTS)blk_sync(GTS)rhs(GTS)comm(GTS)
Figure 11: Weak scaling results on TACC’s Frontera for RK3 timestepping using
LTS(left) and GTS(right) approaches. For this experiment, the maximum and
minimum refinement levels are $l_{max}=8$ and $l_{min}=3$, hence $\Delta L=5$.
Therefore a single LTS step is equivalent to $2^{\Delta L}=32$ global
timesteps. For LTS, the plot shows the overall cost breakdown between block
synchronization (blk_sync), applying time interpolations between blocks
(LTS_interp.), computing the right hand side (rhs) and communication costs
(comm). For GTS, we show the cost breakdown between blk_sync, rhs and comm.
Note that for GTS time interpolations are not required. Note the significant
difference of blk_sync cost between LTS and GTS. For GTS blk_sync is a global
operation, while in LTS blk_sync is a local operation, where synchronization
performed only on the subset of blocks, which are currently being evolved.
These weak scaling results were performed using a grain size of $\sim 100K$
unknowns per core, with the number of cores ranges from $8$ to $16,384$ cores.
The largest problem recorded had $1.6\times 10^{9}$ unknowns. The above
results are generated for radial wave propagation with a maxdepth 10 and a
refinement tolerance of $10^{-5}$.

###  4.4 LTS efficiency and space adaptivity

As mentioned in §3.7, we can approximate the speed up $S$ between LTS and GTS
for a given octant distribution. Since we can end up with meshes where the use
of LTS will not provide significant advantages over GTS, we can selectively
use LTS based on the expected speed up. To evaluate our speed up model (9),
and to assess the overhead of applying the LTS correction operators, we
computed the actual speed up reported for the linear wave propagation with
increasing maxdepth. The estimated and reported speed up values are presented
in the table 2. As can be seen, the estimated speed up values are sufficiently
close to the predicted ones, allowing applications to determine when it is
beneficial to use LTS.

d | $l_{min}$ | $l_{max}$ |  GTS(s) |  LTS(s) | est. speed up | reported speed up
---|---|---|---|---|---|---
9 | 2 | 7 | 2.91 | 0.88 | 3.31 | 3.30
10 | 3 | 8 | 26.69 | 9.26 | 3.31 | 2.87
11 | 2 | 9 | 114.74 | 36.76 | 3.49 | 3.12
12 | 2 | 10 | 238.05 | 73.94 | 3.56 | 3.21

Table 2: The estimated vs. reported speed up for LTS over GTS for linear wave
propagation with increasing maxdepth for adaptive octrees.

Figure 12: Example octree grids generated for black hole binaries with mass
ratio 1,10,100.

Adaptivity in spacetime can be vital for some computational applications,
especially when spacetime adaptivity is necessary to make these simulations
feasible even on leadership architectures. Here we present estimated speed up
by using LTS, for the simulation of binary black hole mergers and the
computation of the resulting gravitational waves[1, 3, 33, 24]. The
computational cost of these simulations increase significantly when the mass
ratio $q$ of the two black holes increases. Assuming we need $n$ grid points
in 1d to capture the larger black hole, to capture the smaller black hole in
the presence of the larger black hole we need $qn$ in 1d, hence, increase of
$q^{3}$ points in 3d. This makes the large mass ratio gravitational wave
simulations infeasible at the time. We use these large mass ratio binary black
hole grids(see figure 12) to estimate the speed up that can be enabled by the
time adaptivity (see table 3). For mass ratios of 10, one can expect up to 70x
speed up, which is a significant reduction in the cost (runtime and energy) of
estimating the gravitational waves and can reduce the burden on supercomputing
resources.

mass ratio | $l_{min}$ | $l_{max}$ | est. speed up
---|---|---|---
1 | 3 | 15 | 9.82
2 | 3 | 16 | 18.105
5 | 3 | 18 | 15.035
10 | 3 | 22 | 71.8302
Table 3: Estimated speed up for binary black hole grid with increasing mass
ratios from 1 to 10.  641282565121024204840968192$0$$100$$200$$300$number of
cores $\rightarrow$time (s)
$\rightarrow$blk_sync(LTS)time_interp.rhs(LTS)comm(LTS)blk_sync(GTS)rhs(GTS)comm(GTS)
Figure 13: Strong scaling results on TACC’s Frontera for LTS(left) and
GTS(right) timestepping using the RK3 explicit scheme. For this experiment,
the refinement levels are $l_{max}=10$ and $l_{min}=3$, hence $\Delta L=7$.
Therefore a single LTS is equivalent to $2^{\Delta L}=128$ global timesteps.
For LTS, the plot shows the overall cost breakdown between block
synchronization (blk_sync), applying time interpolations between blocks
(LTS_interp.), computing the right hand side (rhs) and the communication
costs(comm). For GTS, we show the cost breakdown between blk_sync, rhs and
comm (time interpolations are not required for GTS). The significant
difference of blk_sync cost between LTS and GTS. For GTS blk_sync is a global
operation, while in LTS blk_sync is a local operation, where synchronization
performed only on the subset of blocks, which are currently being
evolved.Presented strong carried out for a fixed problem size of $262M$
unknowns where the number of cores ranging from $64$ to $8192$ cores. Note
that for strong scaling results re-meshing is disabled in order to keep the
problem size fixed and unchanged during evolution.
816326412825651210242048$0$$0.5$$1$number of cores $\rightarrow$time (s)
$\rightarrow$SFC_Wpartbal.(2:1)mesh_gen Figure 14: Weak scaling results in
TACC’s Frontera to perform SFC-based weighted partitioning (SFC_Wpart) for
Gaussian octant distribution centered at (0,0,0) followed by $2:1$ balancing
(bal.2:1) of the octants, which is used as the input for the mesh
generation(mesh_gen). For this experiment, we used 1.6M grid points per core,
using $7^{3}$ points per octant, with the number of cores varying from $8$ to
$2048$. The largest problem reported had a total of $3.3B$ grid points, where
the mesh generation completed under 2s. 326412825651210242048$0$$10$$20$number
of cores $\rightarrow$time (s) $\rightarrow$SFC_Wpartbal.(2:1)mesh_gen Figure
15: Strong scaling results in TACC’s Frontera to perform SFC-based weighted
partitioning (SFC_Wpart) for Gaussian octant distribution centered at (0,0,0)
followed by $2:1$ balancing (bal.2:1) of the octants, which is used as the
input for the mesh generation(mesh_gen). For the depicted strong scaling, we
keep the problem size fixed at $3.3B$ grid points with the number of cores
increasing from $32$ to $2048$.

###  4.5 Weak and strong scaling

Parallel scalability of timesteppers is essential when dealing with large
scale simulations. In this section, we present weak and strong scalability
results for the linear wave propagation problem, on octree meshes using global
and local timestepping. For both weak and strong scaling, we used $7^{3}$ grid
points per octant. For weak scaling, we set the computational domain to
$[-10,10]^{3}$, and used maxdepth 10\. For the maxdepth 10 grid generated
$l_{min}=3$ and $l_{max}=8$, hence 32 global timesteps is equivalent to a
single LTS timestep. Therefore, in order to make the GTS and LTS results
comparable, in the following scaling results, we present timing for 32 steps
in GTS, 1 step (32 partial steps) in LTS. In weak scaling, we increase the
grid size, such that keeping the degrees of freedom per core roughly constant
(100K unknowns per core). Weak scaling results for LTS and GTS are presented
in figure 11. Each bar presents the corresponding evolution time between LTS
and GTS. For each scheme, we present the overall cost breakdown between,
computation of the right-hand side(rhs) of the PDE, and block padding
synchronization (blk_sync) and inter-process communication (comm). The
blk_sync cost consists of space interpolation, which is common for both LTS
and GTS schemes due to space adaptivity. We present an extra fraction of time
interpolation cost between blocks only present in the LTS. For GTS the
blk_sync operation is a global synchronization, i.e., all the blocks are
evolved and need to synchronize the padding regions for the next rhs
computation. In contrast, LTS scheme the blk_sync operation is a
local(partial) synchronization limited to the blocks evolved at the current
partial step followed by time interpolation to correct the padding regions
between blocks. The weak scaling plot shows that the partial synchronization
with appropriate time interpolation is efficient than the global
synchronization in GTS scheme.

To perform strong scaling (see figure 13 ), we use the maxdepth 12 and
recorded $l_{min}$ and $l_{max}$ were 3 and 10, respectively. Therefore, 128
global timesteps are equivalent to a single LTS timestep. For strong scaling
tests, we keep the total grid size fixed $262M$ unknowns and increase the
number of cores from $64$ to $8192$. The strong scaling plot shows, the same
cost breakdown described above. The recorded average parallel efficiencies
between LTS and GTS schemes were $87\%$ and $74\%$, respectively. The low
overhead of blk_sync operation, allows LTS to demonstrate superior weak and
strong scalability compared to GTS scheme.

###  4.6 Weighted partitioning and mesh-generation

The performance of
data partitioning and mesh generation is crucial for AMR applications,
especially when the computational grid changes frequently. We refer to this
process as re-meshing, which require re-partitioning of the data (since the
refinement change may have caused load-imbalance), enforcing 2:1 balancing,
and mesh data structure generation. The performance of the re-meshing is
crucial, but it is not the main focus of this paper. In the current
implementation, we trigger refinement in LTS when all the blocks are
synchronized in time. Figure 9 shows how the grid changes as the wave
propagate radially outwards. Figures 14 and 15 show the weak and strong
scalability of the operations related to re-meshing. The above experiments
show that mesh generation has a relatively high computational cost, compared
to SFC weighted partitioning and 2:1 balancing of octrees. This is mainly
because mesh generation performs a large number of search operations on the
octree to build the neighborhood data structures, which are essential to
perform numerical computations.

##  5 Conclusions

In this paper, we presented methods to enable time adaptivity for solving PDEs
numerically on spatially adaptive grids. We presented experimental results for
the accuracy and scalability of the presented approaches. We show that for
some highly adaptive octrees with high levels of refinement, time adaptivity
can be essential to reduce overall time to solution. As future work, we will
explore more sophisticated energy-conserving time interpolation operators,
with increased accuracy. Currently, re-partitioning, and re-meshing are
triggered as global operations when all the blocks are synchronized in time.
If only a localized portion of the mesh needs to be refined, then ideally one
should only require an update (remesh) of those regions. This is complicated
by the fact that different regions will no longer be syncronized in time.
While the proposed LTS scheme supports this, performing this efficiently
requires new partial remeshing algorithms, that will be the focus of future
work.

## References

  * [1] Einstein at Home.  http://einstein.phys.uwm.edu.
  * [2] FLASH home page.  http://flash.uchicago.edu/website/home/, 2007.
  * [3] A. Abramovici et. al., Science, 256 (1992), p. 325.
  * [4] M. Alcubierre, Introduction to 3+1 numerical relativity, International series of monographs on physics, Oxford Univ. Press, Oxford, 2008\.
  * [5] M. Almquist and M. Mehlin, Multilevel local time-stepping methods of runge–kutta-type for wave equations, SIAM Journal on Scientific Computing, 39 (2017), pp. A2020–A2048, https://doi.org/10.1137/16M1084407, https://doi.org/10.1137/16M1084407, https://arxiv.org/abs/https://doi.org/10.1137/16M1084407.
  * [6] L. B. N. L. Applied Numerical Algorithms Group, Chombo –infrastructure for adaptive mesh refinement.  http://seesar.lbl.gov/anag/chombo/, 2006.
  * [7] W. Bangerth, R. Hartmann, and G. Kanschat, deal. ii—a general-purpose object-oriented finite element library, ACM Transactions on Mathematical Software (TOMS), 33 (2007), pp. 24–es.
  * [8] H. Bao, J. Bielak, O. Ghattas, L. F. Kallivokas, D. R. O’Hallaron, J. R. Shewchuk, and J. Xu, Large-scale simulation of elastic wave propagation in heterogeneous media on parallel computers, Computer methods in applied mechanics and engineering, 152 (1998), pp. 85–102.
  * [9] P. Bastian, Load balancing for adaptive multigrid methods, SIAM Journal on Scientific Computing, 19 (1998), pp. 1303–1321.
  * [10] J. Bédorf, E. Gaburov, and S. Portegies Zwart, Bonsai: A GPU Tree-Code, in Advances in Computational Astrophysics: Methods, Tools, and Outcome, R. Capuzzo-Dolcetta, M. Limongi, and A. Tornambè, eds., vol. 453 of Astronomical Society of the Pacific Conference Series, July 2012, p. 325, https://arxiv.org/abs/1204.2280.
  * [11] M. J. Berger, Adaptive mesh refinement for hyperbolic partial differential equations, tech. report, Stanford Univ Ca. Dept. of Computer Science, 1982.
  * [12] M. J. Berger, Stability of interfaces with mesh refinement, Mathematics of computation, 45 (1985), pp. 301–318.
  * [13] M. J. Berger and J. Oliger, Adaptive mesh refinement for hyperbolic partial differential equations, Journal of Computational Physics, 53 (1984), pp. 484 – 512, https://doi.org/https://doi.org/10.1016/0021-9991(84)90073-1, http://www.sciencedirect.com/science/article/pii/0021999184900731.
  * [14] M. Bern, D. Eppstein, and S.-H. Teng, Parallel construction of quadtrees and quality triangulations, International Journal of Computational Geometry & Applications, 9 (1999), pp. 517–532.
  * [15] G. L. Bryan, M. L. Norman, B. W. O'Shea, T. Abel, J. H. Wise, M. J. Turk, D. R. Reynolds, D. C. Collins, P. Wang, S. W. Skillman, B. Smith, R. P. Harkness, J. Bordner, J. hoon Kim, M. Kuhlen, H. Xu, N. Goldbaum, C. Hummels, A. G. Kritsuk, E. Tasker, S. Skory, C. M. Simpson, O. Hahn, J. S. Oishi, G. C. So, F. Zhao, R. Cen, and Y. L. and, ENZO: An Adaptive mesh refinement code for astrophysics, The Astrophysical Journal Supplement Series, 211 (2014), p. 19, https://doi.org/10.1088/0067-0049/211/2/19, https://doi.org/10.1088%2F0067-0049%2F211%2F2%2F19.
  * [16] C. Burstedde, L. C. Wilcox, and O. Ghattas, p4est: Scalable algorithms for parallel adaptive mesh refinement on forests of octrees, SIAM Journal on Scientific Computing, 33 (2011), pp. 1103–1133, https://doi.org/10.1137/100791634.
  * [17] C. Burstedde, L. C. Wilcox, and O. Ghattas, p4est: Scalable algorithms for parallel adaptive mesh refinement on forests of octrees, SIAM Journal on Scientific Computing, 33 (2011), pp. 1103–1133, https://doi.org/10.1137/100791634.
  * [18] J. Carrier, L. Greengard, and V. Rokhlin, A fast adaptive multipole algorithm for particle simulations, SIAM journal on scientific and statistical computing, 9 (1988), pp. 669–686.
  * [19] U. Catalyurek, E. Boman, K. Devine, D. Bozdag, R. Heaphy, and L. Riesen, Hypergraph-based dynamic load balancing for adaptive scientific computations, in Proc. of 21st International Parallel and Distributed Processing Symposium (IPDPS’07), IEEE, 2007.  Best Algorithms Paper Award.
  * [20] F. Collino, T. Fouquet, and P. Joly, Conservative space-time mesh refinement methods for the fdtd solution of maxwell’s equations, Journal of Computational Physics, 211 (2006), pp. 9–35.
  * [21] R. Courant, K. Friedrichs, and H. Lewy, On the partial difference equations of mathematical physics, IBM Journal of Research and Development, 11 (1967), pp. 215–234, https://doi.org/10.1147/rd.112.0215.
  * [22] Einstein Toolkit.  http://einsteintoolkit.org.
  * [23] M. Fernando, D. Duplyakin, and H. Sundar, Machine and application aware partitioning for adaptive mesh refinement applications, in Proceedings of the 26th International Symposium on High-Performance Parallel and Distributed Computing, HPDC ’17, New York, NY, USA, 2017, ACM, pp. 231–242, https://doi.org/10.1145/3078597.3078610, http://doi.acm.org/10.1145/3078597.3078610.
  * [24] M. Fernando, D. Neilsen, H. Lim, E. Hirschmann, and H. Sundar, Massively parallel simulations of binary black hole intermediate-mass-ratio inspirals, SIAM Journal on Scientific Computing, 41 (2019), pp. C97–C138, https://doi.org/10.1137/18M1196972, https://doi.org/10.1137/18M1196972, https://arxiv.org/abs/https://doi.org/10.1137/18M1196972.
  * [25] M. S. Fernando and H. Sundar, paralab/Dendro-5.01: Local timestepping on octree grids, June 2020, https://doi.org/10.5281/zenodo.3879315, https://doi.org/10.5281/zenodo.3879315.
  * [26] J. E. Flaherty, R. M. Loy, C. Özturan, M. S. Shephard, B. K. Szymanski, J. D. Teresco, and L. H. Ziantz, Parallel structures and dynamic load balancing for adaptive finite element computation, Applied Numerical Mathematics, 26 (1998), pp. 241–263.
  * [27] D. Fuster, A. Bagué, T. Boeck, L. Le Moyne, A. Leboissetier, S. Popinet, P. Ray, R. Scardovelli, and S. Zaleski, Simulation of primary atomization with an octree adaptive mesh refinement and vof method, International Journal of Multiphase Flow, 35 (2009), pp. 550–565.
  * [28] M. J. Gander and L. Halpern, Techniques for locally adaptive time stepping developed over the last two decades, in Domain decomposition methods in science and engineering XX, Springer, 2013, pp. 377–385.
  * [29] C. Gear and D. Wells, Multirate linear multistep methods, (1984).
  * [30] S. Gottlieb and C.-W. Shu, Total variation diminishing runge-kutta schemes, Mathematics of computation, 67 (1998), pp. 73–85.
  * [31] R. Grauer, C. Marliani, and K. Germaschewski, Adaptive mesh refinement for singular solutions of the incompressible euler equations, Physical review letters, 80 (1998), p. 4177.
  * [32] M. J. Grote, M. Mehlin, and T. Mitkova, Runge–kutta-based explicit local time-stepping methods for wave propagation, SIAM Journal on Scientific Computing, 37 (2015), pp. A747–A775, https://doi.org/10.1137/140958293, https://doi.org/10.1137/140958293, https://arxiv.org/abs/https://doi.org/10.1137/140958293.
  * [33] had home page.  http://had.liu.edu, 2009.
  * [34] P. Joly and J. Rodriguez, An error analysis of conservative space-time mesh refinement methods for the one-dimensional wave equation, SIAM journal on numerical analysis, 43 (2005), pp. 825–859.
  * [35] A. Kanevsky, M. H. Carpenter, D. Gottlieb, and J. S. Hesthaven, Application of implicit–explicit high order runge–kutta methods to discontinuous-galerkin schemes, Journal of Computational Physics, 225 (2007), pp. 1753–1781.
  * [36] L. Liu, X. Li, and F. Q. Hu, Nonuniform time-step runge–kutta discontinuous galerkin method for computational aeroacoustics, Journal of Computational Physics, 229 (2010), pp. 6874–6897.
  * [37] L. Liu, X. Li, and F. Q. Hu, Nonuniform-time-step explicit runge–kutta scheme for high-order finite difference method, Computers & Fluids, 105 (2014), pp. 166–178.
  * [38] G. Marcus and M. Teodora, High-order explicit local time-stepping methods for damped wave equations, Journal of Computational and Applied Mathematics, 239 (2013), pp. 270 – 289, https://doi.org/https://doi.org/10.1016/j.cam.2012.09.046, http://www.sciencedirect.com/science/article/pii/S0377042712004190.
  * [39] K. D. Nikitin, M. A. Olshanskii, K. M. Terekhov, and Y. V. Vassilevski, A numerical method for the simulation of free surface flows of viscoplastic fluid in 3d, Journal of Computational Mathematics, (2011), pp. 605–622.
  * [40] W. H. Press, B. P. Flannery, S. A. Teukolosky, and W. T. Vetterling, Numerical Recipes in C, Cambridge, 1992, pp. 350, 379.
  * [41] J. R. Rice, Split runge-kutta method for simultaneous, Journal of Research of the National Bureau of Standards: Mathematics and mathematical physics. B, 64 (1960), p. 151.
  * [42] M. Rietmann, D. Peter, O. Schenk, B. Uçar, and M. Grote, Load-balanced local time stepping for large-scale wave propagation, in 2015 IEEE International Parallel and Distributed Processing Symposium, 2015, pp. 925–935.
  * [43] J. Rudi, A. C. I. Malossi, T. Isaac, G. Stadler, M. Gurnis, P. W. J. Staar, Y. Ineichen, C. Bekas, A. Curioni, and O. Ghattas, An extreme-scale implicit solver for complex pdes: Highly heterogeneous flow in earth’s mantle, in Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, SC ’15, New York, NY, USA, 2015, ACM, pp. 5:1–5:12, https://doi.org/10.1145/2807591.2807675, http://doi.acm.org/10.1145/2807591.2807675.
  * [44] R. Sampath, H. Sundar, S. S. Adavani, I. Lashuk, and G. Biros, Dendro home page, 2008.  http://www.seas.upenn.edu/csela/dendro.
  * [45] A. Sandu and E. M. Constantinescu, Multirate explicit adams methods for time integration of conservation laws, Journal of Scientific Computing, 38 (2009), pp. 229–249.
  * [46] H. Sundar, G. Biros, C. Burstedde, J. Rudi, O. Ghattas, and G. Stadler, Parallel geometric-algebraic multigrid on unstructured forests of octrees, in SC’12: Proceedings of the International Conference on High Performance Computing, Networking, Storage and Analysis, IEEE, 2012, pp. 1–11.
  * [47] H. Sundar, R. Sampath, and G. Biros, Bottom-up construction and 2:1 balance refinement of linear octrees in parallel, SIAM Journal on Scientific Computing, 30 (2008), pp. 2675–2708, https://doi.org/10.1137/070681727.
  * [48] H. Sundar, R. S. Sampath, S. S. Adavani, C. Davatzikos, and G. Biros, Low-constant parallel algorithms for finite element simulations using linear octrees, in SC’07: Proceedings of the International Conference for High Performance Computing, Networking, Storage, and Analysis, ACM/IEEE, 2007\.
  * [49] T. Weinzierl et al., Peano—a Framework for PDE Solvers on Spacetree Grids, 201x, http://www.peano-framework.org.