High-order asynchrony-tolerant finite difference schemes for partial differential equations

TL;DR

This paper develops high-order finite difference schemes that remain accurate despite asynchronous processing in parallel PDE simulations, improving scalability without sacrificing precision.

Contribution

It introduces a general methodology for deriving high-order asynchrony-tolerant finite difference schemes that preserve accuracy under relaxed synchronization.

Findings

Schemes maintain stability and accuracy with asynchrony

Numerical error analyzed within a statistical framework

Validated performance through numerical experiments

Abstract

Synchronizations of processing elements (PEs) in massively parallel simulations, which arise due to communication or load imbalances between PEs, significantly affect the scalability of scientific applications. We have recently proposed a method based on finite-difference schemes to solve partial differential equations in an asynchronous fashion -- synchronization between PEs is relaxed at a mathematical level. While standard schemes can maintain their stability in the presence of asynchrony, their accuracy is drastically affected. In this work, we present a general methodology to derive asynchrony-tolerant (AT) finite difference schemes of arbitrary order of accuracy, which can maintain their accuracy when synchronizations are relaxed. We show that there are several choices available in selecting a stencil to derive these schemes and discuss their effect on numerical and computational performance. We provide a simple classification of schemes based on the stencil and derive schemes that are representative of different classes. Their numerical error is rigorously analyzed within a statistical framework to obtain the overall accuracy of the solution. Results from numerical experiments are used to validate the performance of the schemes.

Tables

Table 1. Table 1: Summary of classification of asynchrony-tolerant (AT) schemes.

Feature	Classification
Layout of	symmetric	asymmetric
grid points	$J_1=J_2$	$J_1\ne J_2$
Nature of delay	unconstrained	uniform delay
at buffer points		${\tilde{k}}_{i+j}=\tilde{K}\forall i+j\in B$
Artificial delay	zero delay	non-zero delay
at interior point	${\tilde{k}}_{i+j}=0\forall i+j\in I$	${\tilde{k}}_{i+j}\le 0\forall i+j\in I$
Coefficients	symmetric	asymmetric
	$\left|{\tilde{c}}_{-j}^0\right|=\left|{\tilde{c}}_j^0\right|$	$\left|{\tilde{c}}_{-j}^0\right|\ne \left|{\tilde{c}}_j^0\right|$

Table 2. Table 2: Parameters of numerical experiments. In the table: AD - linear advection-diffusion equation, D - diffusion equation, VB - non-linear viscous Burgers’ equation; Eul - first order Euler scheme, AB2 and AB3 - second and third order Adams-Bashforth schemes; CD2, CD4 and CD6 - second, fourth and sixth order central difference schemes; AT schemes referred according to ( d , a , r ) 𝑑 𝑎 𝑟 (d,a,r) notation in Tables  LABEL:tab:atschemes-left and LABEL:tab:atschemes-right .

Case	Equation	Time derivative	Space derivatives
			Synchronous	Asynchrony-tolerant
1	AD	Eul	CD2	$(1,2,2)b$, $(2,2,2)b$
2	D	Eul	CD2	$(2,1,2)$
3	AD	Eul	CD2	$(1,2,2)a$, $(2,2,2)a$
4	AD	AB2	CD4	$(1,4,2)$, $(2,4,2)$
5	D	AB3	CD6	$(2,6,2)$
6	VB	AB2	CD4	$(1,4,2)$, $(2,4,2)$

Equations

$${\frac{{\partial }^{d}u}{\partial x^d}|}_i^n=\sum _{j=-J_1}^{J_2}{c}_j{u}_{i+j}^n+\mathcal{O}(\Delta x^a),$$

$${\frac{{\partial }^{d}u}{\partial x^d}|}_i^n\approx \sum _{j=-J_1}^{J_2}{c}_j{u}_{i+j}^{n-l},$$

$${\frac{{\partial }^{d}u}{\partial x^d}|}_i^n\approx \sum _{j=-J_1}^{J_2}\sum _{l=0}^{L-1}{\tilde{c}}_j^l{u}_{i+j}^{n-l},$$

$$u_{i+j}^{n-l}=\sum _{\eta =0}^{\infty }\sum _{\zeta =0}^{\infty }{u}^{(\eta ,\zeta )}\frac{(j\Delta x{)}^{\eta }(-l\Delta t{)}^{\zeta }}{\eta !\zeta !},$$

$$u_{i+j}^n=\sum _{\eta =0}^{\infty }{u}^{(\eta ,0)}\frac{(j\Delta x{)}^{\eta }}{\eta !}$$

$$\sum _{j=-J_1}^{J_2}\sum _{l=0}^{L-1}{\tilde{c}}_j^l{u}_{i+j}^{n-l}$$

$$\sum _{j=-J_1}^{J_2}\sum _{l=0}^{L-1}{\tilde{c}}_j^l\frac{(j\Delta x{)}^{\eta }(-l\Delta t{)}^{\zeta }}{\eta !\zeta !}=\{\begin{array}{ll}1& \text{for }(\eta ,\zeta )=(d,0)\\ 0& \text{for }\eta +r\zeta <d+a;(\eta ,\zeta )\ne (d,0).\end{array}$$

$${\frac{{\partial }^{d}u}{\partial x^d}|}_i^n\approx \sum _{j=-J_1}^{J_2}\sum _{l={\tilde{l}}_1(j)}^{{\tilde{l}}_2(j)}{\tilde{c}}_j^l{u}_{i+j}^{n-l},$$

$$\sum _{j=-J_1}^{J_2}\sum _{l={\tilde{l}}_1(j)}^{{\tilde{l}}_2(j)}{\tilde{c}}_j^l{j}^{\eta }{l}^{\zeta }=0.$$

$$\sum _{i+j\in B}{\tilde{c}}_j^{\tilde{k}}{j}^{\eta }=0,$$

$${\frac{{\partial }^{d}u}{\partial x^d}|}_i^n\approx \sum _{j=0}^J{\varphi }_j\left(u_{i+j}^n+(-1{)}^d{u}_{i-j}^n\right),$$

$$u_{i+j}^n-u_{i-j}^n=2\left[u^{(1,0)}\frac{(j\Delta x)}{1!}+u^{(3,0)}\frac{(j\Delta x{)}^3}{3!}+u^{(5,0)}\frac{(j\Delta x{)}^5}{5!}+\dots \right]$$

$$u_{i+j}^{n-l}-u_{i-j}^{n-l}$$

$${\frac{{\partial }^{d}u}{\partial x^d}|}_i^n\approx \sum _{j=0}^J\sum _{l={\tilde{l}}_1(-j)}^{{\tilde{l}}_2(-j)}{\tilde{\varphi }}_j^l\left(u_{i+j}^n+(-1{)}^d{u}_{i-j}^{n-l}\right)$$

$$\left[\begin{array}[]{cccc}1&1&1&1\\ -\Delta x&0&\Delta x&2\Delta x\\ \frac{\Delta x^{2}}{2}&0&\frac{\Delta x^{2}}{2}&2\Delta x^{2}\\ 0&0&-{\tilde{k}}_{i+1}\Delta t&-{\tilde{k}}_{i+2}\Delta t\end{array}\right]\left[\begin{array}[]{c}{\tilde{c}}_{i-1}^{0}\\ {\tilde{c}}_{i}^{0}\\ {\tilde{c}}_{i+1}^{{\tilde{k}}_{i+1}}\\ {\tilde{c}}_{i+2}^{{\tilde{k}}_{i+2}}\\ \end{array}\right]=\left[\begin{array}[]{c}0\\ 1\\ 0\\ 0\\ \end{array}\right].$$

$${\frac{\partial u}{\partial x}|}_i^n$$

$${\frac{\partial u}{\partial x}|}_i^n=\frac{-3u_{i-1}^n+3u_i^n-u_{i+1}^{n-\tilde{k}}+u_{i+2}^{n-\tilde{k}}}{4\Delta x}+\mathcal{O}\left(\Delta x^2,\tilde{k}\Delta t\right),$$

$${\frac{{\partial }^{2}u}{\partial x^2}|}_i^n$$

$$(\eta ,\zeta )$$

$$\left[\begin{array}[]{ccc}2&2&2\\ 0&\Delta x^{2}&4\Delta x^{2}\\ 0&\frac{\Delta x^{4}}{12}&\frac{4\Delta x^{4}}{3}\\ \end{array}\right]\left[\begin{array}[]{c}{\tilde{\phi}}_{0}^{l}\\ {\tilde{\phi}}_{1}^{l}\\ {\tilde{\phi}}_{2}^{l}\\ \end{array}\right]=\left[\begin{array}[]{c}0\\ 1\\ 0\\ \end{array}\right]$$

$${\frac{{\partial }^{2}u}{\partial x^2}|}_i^n=\frac{-u_{i+2}^n+16u_{i+1}^n-30u_i^n+16u_{i-1}^{n-l}-u_{i-2}^{n-l}}{12\Delta x^2}+\mathcal{O}(\Delta x^4,l\Delta t,l\Delta t/\Delta x^2).$$

$$\left[\begin{array}[]{ccc}-{\tilde{k}}\frac{5\Delta t}{4\Delta x^{2}}&-({\tilde{k}}+1)\frac{5\Delta t}{4\Delta x^{2}}&-({\tilde{k}}+2)\frac{5\Delta t}{4\Delta x^{2}}\\ {\tilde{k}}^{2}\frac{5\Delta t^{2}}{8\Delta x^{2}}&({\tilde{k}}+1)^{2}\frac{5\Delta t^{2}}{8\Delta x^{2}}&({\tilde{k}}+2)^{2}\frac{5\Delta t^{2}}{8\Delta x^{2}}\\ 1&1&1\\ \end{array}\right]\left[\begin{array}[]{c}{\tilde{\phi}}_{0}^{{\tilde{k}}}\\ {\tilde{\phi}}_{1}^{{\tilde{k}}+1}\\ {\tilde{\phi}}_{2}^{{\tilde{k}}+2}\\ \end{array}\right]=\left[\begin{array}[]{c}0\\ 0\\ 1\\ \end{array}\right]$$

$${\frac{{\partial }^{2}u}{\partial x^2}|}_i^n$$

$$\sum _{k[i]=0}^{L-1}{p}_{k[i]}=1.$$

$$⟨\overline{E}⟩=\frac{1}{N}\sum _{i=1,N}\overline{E_i^n},$$

$$⟨\overline{E}⟩=\frac{1}{N}\left[\sum _{i\in I_I}{E}_i^n+\sum _{i\in I_B}\overline{{\tilde{E}}_i^n}\right]$$

$$⟨\overline{E}⟩=\frac{1}{N}\left[\sum _{i\in I_I}{E}_i^n+\sum _{i\in I_B}{E}_i^n{|}_s+\sum _{i\in I_B}\overline{{\tilde{E}}_i^n}{|}_a\right],$$

$$u_i^{n+1}=u_i^n+\Delta t\sum _{m=0}^{T-1}{\beta }_m{f}_i^{n-m},$$

$$u_i^{n+1}=u_i^n+\Delta t\left(\frac{3}{2}{f}_i^n-\frac{1}{2}{f}_i^{n-1}\right),$$

$$\frac{\partial u}{\partial t}=\alpha \frac{{\partial }^{2}u}{\partial x^2},$$

Full text

1 Introduction

Numerical simulations are an important tool in understanding complex problems in physics and engineering systems. Many of these phenomena are multi-scale in nature, and are governed by nonlinear partial differential equations (PDEs). With a wide range of scales at realistic conditions, like turbulence phenomena in fluid flows, the numerical solution of these equations becomes computationally very expensive. Advances in computing technology have made it possible to carry out intensive simulations on massively parallel computers. Currently, state-of-the-art simulations are routinely being done on tens or hundreds of thousands of processing elements (PEs) [1, 2, 3, 4].

It is known, at extreme scale, that data communication as well as synchronization between PEs pose a major challenge in the scalability of scientific applications [5]. In the case of PDE solvers, where the parallelism is typically realized by decomposing the computational domain among PEs, communications that affect the scalability arise due to the computation of spatial derivatives in order to propagate the physical information across the domain. The problem becomes more acute in simulations of transient phenomena, where spatial derivatives are evaluated at each time step over an integration of large number of steps. Another issue concerning the scalability is related to the performance variations across the PEs in a parallel system. In this case, sub-optimal performance of even a few PEs may lead to idling of others, as dictated by the data dependencies involved in the computations. It is likely that in future Exascale computing systems, which will have an extremely large PE count, communication and synchronization will be a major bottleneck. It is thus not surprising that there is a substantial increased interest in developing numerical methods that minimize communications and relax data synchronizations at the mathematical level [6, 7].

An early effort in solving PDEs in an asynchronous fashion has been presented in [8, 9]. Their method is based on finite-difference schemes and is restricted to the solution of parabolic PDEs with at most second order accuracy. More recent work [10, 11], again based on finite-difference method, has suggested that due to the randomness in the arrival of messages at different PEs, the resulting algebraic difference equations are stochastic in nature. In that work, a statistical framework to analyze such systems was developed to study the numerical properties of commonly used schemes in the presence of asynchrony. Furthermore, they show that though the stability and consistency of the schemes can be maintained, their accuracy is significantly degraded. They also proposed the possibility of deriving schemes that are tolerant to communication data asynchrony. A follow up of this work to a simple specific equation and numerical scheme has be presented in [12]. Although the authors were able to maintain second order accuracy for their chosen scheme when asynchrony is present, one can show using Taylor series that they are severely limited to low order of accuracy. However, as mentioned earlier, a number of natural and engineering systems are multi-scale in nature and will require higher order accurate schemes. In this work, we present a general methodology to generate different classes of high-order asynchrony-tolerant (AT) schemes. This is the main objective of this paper.

The rest of the paper is organized as follows. We first briefly review the concept of asynchronous computing for PDEs in section 2. A general method to derive AT schemes, the choices in stencil available in arriving at these schemes and their classification are presented in section 3. In section 4, we show a statistical framework to analyze the overall accuracy of a numerical method when AT schemes are used. Numerical experiments to validate the performance of AT schemes are shown in section 5. Conclusions and further discussions are presented in section 6.

2 Concept

Let $u(x,t)$ be a function of spatial coordinate $x$ and time $t$, which is governed by a time-dependent PDE in a one-dimensional domain. Fig. 1 illustrates the discretized domain which is decomposed into $P$ number of PEs. Let $i$ and $n$ represent an arbitrary grid point in the domain and time level such that $u(x_{i},t_{n})={u_{i}^{n}}$. For clarity in the exposition, we assume that the grid points are uniformly distributed in the domain with a spacing $\Delta x$. A finite-difference to approximate a spatial derivative at point $i$ and time level $n$ can be expressed, in the most general case, as

[TABLE]

where $d$ is the order of the derivative, $J_{1}$ and $J_{2}$ are the number of points to the left and right of point $i$ in the stencil, and $c_{j}$ is the appropriate coefficient or weight of ${u_{i+j}^{n}}$ such that the scheme is accurate to an order $a$ in space. The term $\mathcal{O}(\Delta x^{a})$ represents the truncation error of the scheme.

Usually, the numerical solution of a time-dependent PDE is obtained by advancing an initial condition according to an algebraic finite-difference equation in small steps of time $\Delta t$. During each time advancement, say, marching from a time level $n$ to $n+1$, spatial derivatives are computed at each grid point using Eq. (1). In general, these computations are trivial to implement in a serial code, as the value of the function at all the grid points will be locally available in the memory of the PE. However, if the domain is decomposed into multiple PEs, computations at points near PE boundaries may need values of the function at stencil points that are computed in the neighboring PEs. Such values are commonly communicated into buffer or ghost points, as shown in Fig. 1. Note that the number of values communicated across the left and right PE boundaries is equal to $J_{1}$ and $J_{2}$, respectively.

Let $I$ represent the set of physical grid points in the domain and $B$ represent the set of buffer points. For convenience we divide the set $I$ further such that $I={I_{I}}\cup{I_{B}}$. The set of grid points near PE boundaries whose computations need data from neighboring PEs will be denoted by ${I_{B}}$. The complementary set of interior points, whose computations are independent of communication between PEs is denoted by ${I_{I}}$. In commonly used parallel algorithms, computations at a point $i\in{I_{B}}$ cannot be advanced until the communication between PEs is complete. This is typically ensured by enforcing communication synchronization after messages are issued from one PE to another. As mentioned earlier, with a large number of PEs such synchronizations become expensive and result in poor scalability of codes at extreme scales. We refer to this as synchronous computing.

In the case of asynchronous computing, communication between PEs is initiated at each time step, however, the data synchronization is not enforced. This means, we cannot ensure that the time level of the function at buffer points is $n$. It can be $n$, $n-1$, $n-2$, … depending on the status of messages from successive time advancements. Due to the random nature of the arrival of messages at different PEs [13], the availability of a particular time level at a buffer point is also random. Let ${\tilde{n}}=n-{\tilde{k}}_{j}$ be the latest available time level at a buffer point $j$, where ${\tilde{k}}_{j}$ is the corresponding random delay at that point222In order to distinguish random from deterministic variables, we will use a tilde ( $\tilde{}$ ) over the variable for the former.. Note that ${\tilde{n}}$ can be different at different locations and time levels. If we restrict the maximum allowable delay levels to $L$, then ${\tilde{n}}\in\left\{n,n-1,...,n-L+1\right\}$ and ${\tilde{k}}_{j}\in\left\{0,1,...,L-1\right\}$. The scheme in Eq. (1), when asynchrony is allowed, can be rewritten as

[TABLE]

where $l=0$ for $i+j\in I$ and $l={\tilde{k}}_{i+j}$ for $i+j\in B$. Unlike the scheme in Eq. (1) which contains a single time level, this scheme uses multiple time levels when some of the points in the stencil belong to the set $B$. It has been shown in [11] that the accuracy of common finite-differences used in such an asynchronous fashion is significantly affected. In particular, accuracy drops to first order regardless of the original finite difference used. Thus, the need to derive AT schemes that maintain accuracy even when there is a communication delay. We do this next.

3 Asynchrony-tolerant (AT) schemes

3.1 General methodology

Taylor series and the method of undetermined coefficients provide a systematic procedure to derive finite-difference schemes. As we show momentarily, this approach can also be used to construct AT schemes to approximate spatial derivatives. Let ${u_{i+j}^{n-l}}$ represent the function at a generic point $i+j$ in the stencil with an arbitrary delay of $l$ levels to compute a spatial derivative at a point $i$ and time level $n$. Using the $L$ possible time delays, we can express an AT scheme as

[TABLE]

where ${\tilde{c}}_{j}^{l}$, for the range of $j$ and $l$, are the appropriate coefficients that have to be determined. Note that this scheme represents the most general case with the function at all possible time levels at each point in the stencil. However, depending on the delay at each grid point and time step, which is given by ${\tilde{k}}_{i+j}$, only one or few time levels may be used in approximating the derivative. The random nature of ${\tilde{k}}_{i+j}$ is, now, embedded into ${\tilde{c}}_{j}^{l}$. The merits of using older time levels not just at buffer points, but also at interior points will be discussed later.

The coefficients in the scheme expressed in Eq. (3) can be obtained by imposing constraints on different terms of the Taylor series, upon expansion of the function at each combination of point and time level in the stencil. Let us consider the Taylor series of ${u_{i+j}^{n-l}}$ about the point $i$ and time level $n$. The series is an expansion in two variables, namely $\Delta x$ and $\Delta t$, which is given by

[TABLE]

where $u^{(\eta,\zeta)}$ denotes the $\eta$-th and $\zeta$-th partial derivative in space and time, respectively, of $u$ evaluated at $i$ and $n$. When $l=0$, the function corresponds to a synchronous value of $u$. This makes the terms in the series a function of $\Delta x$ only:

[TABLE]

To obtain the constraints that will assure a given order of accuracy, we substitute the Taylor series of $u$ in the right hand side of Eq. (3).

[TABLE]

The linear combination of the function values, in the above equation, represents a scheme when ($i$) the coefficient of the $d$-th derivative of $u$ in space is unity, and ($ii$) low-order terms are eliminated according to the desired accuracy of the scheme.

Let $a$ be the desired order of accuracy in space. This means that the leading order term in the truncation error should vary with the grid spacing as $\Delta x^{a}$. In the Taylor series of the function for synchronous schemes, as in Eq. (5), the lower order terms can be readily identified as the ones with the power of $\Delta x$ less than $d+a$. However, when asynchrony is present this is not obvious. The terms in the series can now be a function of either or both $\Delta x$ and $\Delta t$, which are usually not independent. In order to identify the lower order terms, let us assume the relation $\Delta t\sim\Delta x^{r}$. Such a relation is often obtained from analysis of the scheme’s numerical stability or other constraints posed by the physics of the problem. Using this relation, we can arrive at the condition to identify lower order terms that need to be eliminated to obtain a scheme of order $a$. This expression is: $\eta+r\zeta<d+a$. Using Eq. (6), we can then summarize the constraints as

[TABLE]

Clearly, the first condition in the above equation makes the coefficient of the term corresponding to $d$-th spatial derivative on the right hand side of Eq. (6) unity. The second condition will set to zero all the necessary lower order terms to obtain an overall accuracy $a$. For a given stencil, these conditions give rise to a system of linear equations. The number of equations in the system is one more than the number of lower order terms that have to be eliminated from Eq. (6). Let $\boldsymbol{A}\boldsymbol{\tilde{c}}=\boldsymbol{b}$ represent this system, where $\boldsymbol{A}$ is the coefficient matrix whose elements are a function of $j$ and $l$, $\boldsymbol{\tilde{c}}$ is the vector of variables that contains coefficients in the scheme and $\boldsymbol{b}$ is the vector with zero elements except for the row corresponding to the order of the derivative to be approximated. The solution to this system determines the coefficients of the scheme.

Before getting into the discussion on the choice of stencil, we make a few observations regarding the linear system when asynchrony is present. To aid the discussion we express the terms in the Taylor series of the function at the generic stencil point, ${u_{i+j}^{n-l}}$, in a matrix format as shown in Fig. 2. This provides a simple format to visualize different terms in the series and help us easily identify the terms on which the conditions in Eq. (7) have to be imposed.

In this graphical representation, we omit constants in each term for the sake of clarity.

In words, Eq. (7) implies constraints on the term containing the derivative of order $(d,0)$ and on all the terms that satisfy the inequality $\eta+r\zeta<d+a$, that is, all terms above the $\eta+r\zeta=d+a$ line in Fig. 2. With this representation, we can easily separate terms that need to be eliminated from those that do not. To illustrate this, let us choose $d=1$, $a=2$ and $r=1$, which corresponds to a second-order approximation of the first derivative, using a convective-type CFL condition such that $\Delta t\sim\Delta x$. For these parameters, conditions are imposed on the terms with $(\eta,\zeta)=\{(0,0),(1,0),(2,0),(0,1),(1,1),(0,2))\}$, which are the terms above the line $A$ in the figure. If asynchrony is absent, that is, $l=0$, the only terms that are non-zero in the table belong to the first column. This shows that, for a given accuracy, the number of terms on which conditions are imposed is larger when asynchrony is present, which thus results in a larger linear system.

The increase in the number of equations also depends on $r$, which relates $\Delta t$ and $\Delta x$. For example, the situation for $r=2$ is also shown in Fig. 2 with line $B$. The number of terms above the line $B$ (4 terms) is less than $A$ (6 terms), which implies that a higher $r$ will reduce the number of lower order terms due to asynchrony for a given accuracy.

The other aspect is the increase in stencil size with increase in accuracy. In commonly used synchronous schemes ($l=0$), a successive increase in the order of accuracy will impose a new condition on one more term in the Taylor series, which adds an additional equation to the linear system. The linear system is then solved by adding one more grid point to the stencil. However, in deriving AT schemes more than one additional equation may be added to the system (compare the number of terms above lines A and C). Thus, we expect the stencil of AT schemes to grow larger than commonly used synchronous schemes when the desired accuracy is increased.

3.2 Choice of stencil

In principle one can choose a stencil that consists of different grid points and time levels to approximate spatial derivatives. However, the stencil of commonly used synchronous schemes are constructed exclusively with spatial grid points. This has some advantages. First, the function at the synchronous time level is available for spatial derivative evaluation at all points in the domain. Second, as argued in [11], and elaborated in section 3.1 above, this choice avoids the additional terms that will appear in Taylor series when the stencil consists of delayed time levels.

When asynchrony is present, on the other hand, as is clear from Eq. (3), the function can belong to multiple time levels. Thus, we can take advantage of using the function at delayed time levels in deriving AT schemes. The choice of stencil should be made according to the nature of terms in Taylor series on which conditions in Eq. (7) are imposed. To understand this let us recall the tabular representation of the Taylor series of ${u_{i+j}^{n-l}}$, as shown in Fig. 3.

We can classify terms into four groups, as represented by the four different colors in the figure. Terms in blue are a function of $\Delta x$ alone. These terms will appear in the Taylor series of the function when $j\neq 0$. Similarly, terms that are a function of $\Delta t$ only are shown in green and they appear when $l\neq 0$. Terms in red are a function of both $\Delta x$ and $\Delta t$, and these appear when $j\neq 0$ and $l\neq 0$. The term $u$ in black is a function of neither $\Delta x$ nor $\Delta t$ and is present in the Taylor series of the function at any point and time level. In order to eliminate specific terms in the truncation error, it is apparent that we cannot arbitrarily choose the points and time levels in a stencil. They have to be selected according to the number of terms in each of these groups. For example, if a linear system consists of three equations that correspond to condition on terms that belong to the red group, then the scheme would need the function evaluated at a minimum of three combinations of $j$ and $l$ such that $j\neq 0$ and $l\neq 0$. If not, the linear system may not have a solution or may have a solution which correspond to stencils completely biased towards the synchronous side of the stencil, like forward and backward differences.

The choice of stencil has consequences also in terms of the performance of simulation codes on parallel machines. Expanding the stencil in space will lead to larger message sizes to be sent over the network, which may be too expensive at extreme scales. Using multiple levels in time will keep the messages relatively smaller, but will increase the memory requirements in each PE. This choice, thus, would require information on the specific computing system to be used for the simulation.

The rectangular box in Fig. 4 illustrates the layout of the stencil used in expressing the general scheme in Eq. (3). However, as mentioned earlier, not all the time levels at all points are required to approximate the derivative. Instead, one can limit the number of time levels at each grid point in such a way to introduce the exact number of coefficients that would make the linear system solvable. Eq. (3), then, becomes:

[TABLE]

where $\tilde{l}_{1}(j)$ and $\tilde{l}_{2}(j)$ are the lower and upper limits on the time levels used at the point $i+j$. These limits are computed according the latest time level available and the number of time levels chosen at that point in the stencil. As an example, in Fig. 4 we identify a stencil to solve a system with four equations. At the two interior points the latest available time level is $n$, which has a zero delay. Thus, the limits are $\tilde{l}_{1}(j)=\tilde{l}_{2}(j)=0$, for $j\in\{-1,0\}$. As specified before, the latest available time level at the buffer point is given by ${\tilde{n}}=n-{\tilde{k}}_{i+1}$, and we use two successive time levels at this point. The limits on the time level at this point are then $\tilde{l}_{1}(1)={\tilde{k}}_{i+1}$ and $\tilde{l}_{2}(1)={\tilde{k}}_{i+1}+1$.

A choice of stencil will lead to a scheme only when there exists a solution to the resulting linear system $\boldsymbol{A}\boldsymbol{\tilde{c}}=\boldsymbol{b}$. Since $\boldsymbol{b}\neq\boldsymbol{0}$ due to the first condition in Eq. (7), the system is non-homogeneous and has a unique solution only when the matrix $\boldsymbol{A}$ is non-singular or has a full rank. If $N_{A}$ is the size of the linear system, then the matrix has full rank when $rank(\boldsymbol{A})=N_{A}$. We, then, obtain the scheme by solving the system and substituting the coefficients into Eq. (8). On the other hand, when $rank(\boldsymbol{A})<N_{A}$, the matrix is singular and the linear system possesses either no solution or infinite solutions. We can distinguish these two cases by computing the rank of the augmented matrix $\boldsymbol{A}|\boldsymbol{b}$. If $rank(\boldsymbol{A})\neq rank(\boldsymbol{A}|\boldsymbol{b})$, then the system is inconsistent and the choice of stencil does not result in a scheme. In the case where $rank(\boldsymbol{A})=rank(\boldsymbol{A}|\boldsymbol{b})$, the linear system is consistent, but has infinite solutions. The linear system, then, contains two or more equations that are linearly dependent. This means that for the choice of stencil, conditions on at least two of the terms are mathematically equivalent or a condition on at least one of the terms can be obtained from linear combination of others. In such a situation, we can get a scheme with greater accuracy with the same stencil. In some cases, it is possible to construct a smaller linear system (with a corresponding smaller stencil) comprised of linearly independent equations, which does in fact have a unique solution. The greater the number of linearly dependent equations, the smaller the linear system with linearly independent equations will be. This suggests that a judicious selection of grid points and time levels can be used to increase the number of linearly dependent equations in the resulting system, and thus reduce the stencil size which in turn reduces computations as well as the size of communication messages. This will be of interest in deriving AT schemes, which demand larger stencil due the presence of terms due to asynchrony.

Let us recall the second condition from Eq. (7), imposed to eliminate the terms due to asynchrony in deriving a scheme. After cancelling out the term $\Delta x^{\eta}\Delta t^{\zeta}/\eta!\zeta!$, which is constant across the equation for a given $(\eta,\zeta)$, we get

[TABLE]

It is evident from the above equation that the existence of linearly dependent equations rests on the values of $j$ and $l$ which are defined by the stencil, as well as $\eta$ and $\zeta$ which represent the order of the derivative corresponding to the equation. As mentioned earlier, the function in the stencil can belong to multiple time levels. The time level of the function at the interior points has a zero delay, that is $l=0$, and hence, will not appear in equations corresponding to asynchrony terms. If we choose a single uniform time level with a delay ${\tilde{k}}_{i+j}=\tilde{k}$ for all $i+j\in B$ in the stencil, then $\tilde{l}_{1}(j)=\tilde{l}_{2}(j)=\tilde{k}$ which leads to a uniform value of $l^{\zeta}$ in Eq. (9). The equation then reduces to

[TABLE]

which is independent of $\zeta$, and shows that for a given $\eta$, equations corresponding to $\zeta\neq 0$ are linearly dependent. With reference to Fig. 3, when we eliminate a term in the red or green groups, all the other terms in the corresponding row that are in the same group are also eliminated. This illustration shows that it is possible to choose a stencil which results in linearly dependent equations in a system.

We conclude this section by summarizing the steps to derive AT schemes.

List the terms on which conditions have to be imposed for a given $d$, $a$ and $r$ 2. 2.

Identify an appropriate stencil ($J_{1},J_{2},\tilde{l}_{1}(j),\tilde{l}_{2}(j)$) according to the terms in the list 3. 3.

Compute the rank of the matrix $\boldsymbol{A}$

• (a)

$rank(\boldsymbol{A})=N_{A}$: unique solution

• (b)

$rank(\boldsymbol{A})<N_{A}$ and $rank(\boldsymbol{A})\neq rank(\boldsymbol{A}|\boldsymbol{b})$: no solution, identify a new stencil

• (c)

$rank(\boldsymbol{A})<N_{A}$ and $rank(\boldsymbol{A})=rank(\boldsymbol{A}|\boldsymbol{b})$: infinite solutions, add more conditions to get greater accuracy or reduce the system size with only linearly independent equations and adjust the stencil size 4. 4.

Solve for $\boldsymbol{\tilde{c}}$ and substitute the coefficients into general scheme

3.3 Alternative approach

It is often necessary to use schemes with a specific structure in terms of stencil and the corresponding coefficients to either improve computational performance or satisfy numerical properties. In the context of AT schemes, it is desirable to use schemes at PE boundary points that are similar in nature to those at interior points. Such an implementation may improve the overall stability of a numerical method and relieve the natural tendency of concentrated errors in the spatial distribution near PE boundaries (e.g. see Fig. (3) in [11]). Though the method described earlier gives the flexibility to choose a particular structure for the stencil, there is not much control over the nature of the resulting coefficients in the scheme. This is because the necessary conditions imposed on the terms on the right hand side of Eq. (6) are all solved in a single linear system. And, explicit conditions on the coefficients have to be added to the linear system to address this issue.

In an alternative approach to derive schemes, we propose to impose necessary conditions (similar to Eq. (7)) on the set of terms arising from Taylor series in a step-by-step process. In each step, a subset of lower order terms are eliminated, while retaining the derivative order term using a particular stencil. This process is repeated until the desired accuracy is achieved. Linear systems of smaller size can be constructed in each step to enforces the conditions and obtain the coefficients. The procedure described in bullet $3$ in the summary of section 3.1 should be used in computing the solution of these systems.

We now proceed to outline the procedure to derive schemes similar to central differences using this approach, and will later provide a detailed illustration in example 3 of section 3.4.

Central difference schemes are widely used in solving parabolic and elliptic PDEs, and are shown to have low numerical dissipation, necessary to resolve all scales in multi-scale phenomena [14]. If we consider the structure of central difference schemes, it can be characterized by a symmetric stencil about the point of computation, and symmetric coefficients (in absolute value). A general synchronous central difference scheme can be expressed as

[TABLE]

where $J$ determines the size of stencil and $\phi_{j}$ are the appropriate coefficients. Let us consider this stencil in the presence of asynchrony. In practical simulations each PE, typically, is assigned a large number of grid points. When asynchrony is allowed in such cases, delays are experienced only on one side of the stencil, that is, either on the left or right about the point of computation $i$ for $i\in I_{B}$. If we assume the delay on the left side, which implies $i-j$ is a buffer point, then the terms in the sum in Eq. (11) take the form $\left({u_{i+j}^{n}}+(-1)^{d}{u_{i-j}^{n-l}}\right)$. To maintain the above mentioned symmetries in AT schemes, we use this sum to eliminate some of the lower order terms and retain the derivative order term in the Taylor series in the first step. As delay is present only at $i-j$, none of the terms due to asynchrony in the expansion of ${u_{i-j}^{n-l}}$ are cancelled out in the sum of the function at the two points. However, some of the terms, which are not a function of $\Delta t$, cancel out depending on the order of derivative $d$. If $d$ is odd, then terms that correspond to even power of $\Delta x$ cancel out, as shown next.

Consider the difference ${u_{i+j}^{n}}-{u_{i-j}^{n}}$ where we choose $l=0$ to simplify the analysis. The conclusions, though, are valid for arbitrary delays $l>0$. A Taylor series expansion can then be written as

[TABLE]

Similarly, if we consider the sum, ${u_{i+j}^{n}}+{u_{i-j}^{n-l}}$, terms with odd powers of $\Delta x$ will vanish. This reduces some of the terms on which conditions need to be imposed, as we move on to the next step. A further decrease in the number of conditions can be achieved by artificially imposing the same delay of $l$ levels on the other side of the stencil, that is, $\left({u_{i+j}^{n-l}}+(-1)^{d}{u_{i-j}^{n-l}}\right)$. The Taylor series expansion of this difference, for odd $d$, is

[TABLE]

which shows that all the terms with even powers of $\Delta x$, regardless of the power of $\Delta t$, are absent. Indeed, imposing this artificial delay on the function at the interior point, though demands additional storage of more time levels at each grid point, will lead to a smaller number of constraints. Thus, schemes with delay on both sides will need a smaller stencil to compute derivatives.

In [12], this approach was used to recover the drop in accuracy due to delay in communication in a particular application using central difference schemes. The authors further suggested that imposing delay on both sides of the stencil in central differences would suffice to maintain the accuracy under asynchronous conditions. However, it can be shown from a Taylor series expansion that the schemes cannot be accurate beyond second order under the conditions they presented. It is essential to increase the stencil size to achieve higher order accuracy when asynchrony is present, as shown in this work.

The remaining lower order terms can be eliminated by expanding the stencil with additional terms of the form $\left({u_{i+j}^{n}}+(-1)^{d}{u_{i-j}^{n-l}}\right)$ for different values of $j$ or $l$. This can be done either in a single or multiple steps, and both of them will ensure symmetry in the coefficients. Assuming delay on the left of the stencil, the resultant AT scheme takes the form:

[TABLE]

It is often useful to derive AT schemes that reduce to central difference schemes when all delays are zero, that is ${\tilde{k}}_{j}=0$ for $j\in B$. Such schemes can be derived by expanding the stencil in space, using the sum at different $j$, to eliminate terms that are not a function of $\Delta t$. And, use the sum at different levels in time to cancel out the terms due to asynchrony. This approach, which contains the essence of the alternative procedure presented in this section, will be illustrated in detail as Example 3 below.

3.4 Classification of schemes

In arriving at AT schemes there are several choices available in terms of choosing points and time levels in a stencil and on the nature of coefficients. We first provide a simple classification of AT schemes based on these choices, and then we present some examples.

Let us consider the stencil of the general AT scheme in Eq. (8), which is given by the limits $J_{1}$ and $J_{2}$ in space and $\tilde{l}_{1}(j)$ and $\tilde{l}_{2}(j)$ in time. If $J_{1}=J_{2}$, the number of points are equal on either sides of the point of computation $i$. We refer to this as a symmetric stencil in space. Else, $J_{1}\neq J_{2}$ and the stencil is asymmetric. Regarding the nature of the delays, schemes can potentially have different delay values at different points in a stencil. However, enforcing a uniform delay across all the buffer points in a scheme, that is ${\tilde{k}}_{i+j}={\tilde{k}}$ for all $i+j\in B$, may lead to linearly dependent conditions and a simpler implementation of schemes. We can, thus, classify schemes according to the presence or absence of uniform delay in schemes. In addition to the uniformity of delays, schemes can also be classified with respect to the time levels chosen at interior points. The function at these points can be either at the synchronous time level or at artificially imposed levels which, as we have shown, provide some numerical advantages.

Schemes can also be classified on the basis of the nature of coefficients. When asynchrony is present, a stencil with symmetric points may not necessarily give rise to symmetry in its coefficients. This is due to the non-uniform time levels in the stencil at these points. To obtain symmetry in coefficients, as in standard central difference schemes, we have earlier proposed to use a sum or difference of the function at symmetric grid points. In this regard, we classify schemes with symmetric coefficients as the ones which have $|{\tilde{c}}_{-j}^{0}|=|{\tilde{c}}_{j}^{0}|$ (i.e. when ${\tilde{k}}_{i+j}=0$). A summary of these classifications is given in Table 1.

From the discussions in previous sections, it is clear that each of these classifications will have consequences in terms of numerical properties and computational performance of schemes. We will now proceed to derive three AT schemes and demonstrate how the conditions corresponding to the classification can be implemented in arriving at them.

Example 1: first derivative - second order accurate ($d=1$, $a=2$)

Using $r=2$, conditions in Eq. (7) are imposed on terms that satisfy the inequality $\eta+2\zeta<3$. This gives rise to a linear system with four equations corresponding to the terms with $(\eta,\zeta)=\{(0,0),(1,0),(2,0),(0,1)\}$ in the Taylor series. The next step is to select a stencil with the function defined at four different combinations of points and time levels. Let, as before, $n-{\tilde{k}}_{i+j}$ be the latest available time level at a point $i+j\in B$ with $j>0$. Further, let us choose the function set $\{{u_{i-1}^{n}},{u_{i}^{n}},{u_{i+1}^{n-{\tilde{k}}_{i+1}}},{u_{i+2}^{n-{\tilde{k}}_{i+2}}}\}$ to construct the linear system $\boldsymbol{A}\boldsymbol{\tilde{c}}=\boldsymbol{b}$. This results in

[TABLE]

The rank of the coefficient matrix in the above equation is $4$, which is equal to the size of the system. Thus, the choice of stencil results in a scheme without any further adjustments. After solving for the coefficients, we obtain the scheme as

[TABLE]

Note that the coefficients are a function of the random delay at buffer points. It is easy to see by inspection that this scheme has to be complemented in two specific circumstances. First, when $3{\tilde{k}}_{i+1}-{\tilde{k}}_{i+2}=0$ the approximation has an infinite value. This can be avoided by artificially altering the delays such that $3{\tilde{k}}_{i+1}-{\tilde{k}}_{i+2}\neq 0$. Second, when the function at both buffer points is at a synchronous time level, i.e., no delay, the above scheme will result in an indeterminate form. In that case one can use

[TABLE]

which is obtained by substituting ${\tilde{k}}_{i+1}={\tilde{k}}_{i+2}={\tilde{k}}$ and simplifying Eq. (16). It is interesting to see that the coefficients in the above scheme, with a uniform delay across the buffer points, are independent of the delay value, which eliminates the limitations of the scheme in Eq. (16). Similar schemes can be derived by considering delays on the left of the stencil.

Example 2: second derivative - fourth order accurate ($d=2$, $a=4$)

The relationship between the time step and grid spacing is assumed as $\Delta t\sim\Delta x$, that is $r=1$. For this set of parameters, Eq. (7) enforces conditions on 21 terms in the Taylor series, which are highlighted in red in Fig. 5.

The resulting linear system has 21 equations which, in principle, will need the function at 21 combinations of points and time levels. However, from Eq. (9) we can see that a stencil with only two time levels, $n$ for interior points and $n-{\tilde{k}}$ for buffer points, will lead to linearly dependent equations in the system. We use this choice of stencil to reduce the size of the linear system. Choosing the limits $\{J_{1},J_{2}\}=\{5,6\}$ in space and assuming the buffer points are on the right side of the stencil, leads to a smaller linear system with 11 equations that has a unique solution. Upon solving the system, the resulting AT scheme is

[TABLE]

Note that a single linear system has been used to obtain the above scheme.

Example 3: second derivative - fourth order accurate ($d=2$, $a=4$)

In this example, we will use the alternative step-by-step approach described in section 3.3 to derive an AT scheme that reduces to a standard central difference scheme in the absence of delays. If we consider the Taylor series of $u$ at a generic point and time level in a stencil, and assume $\Delta t\sim\Delta x^{2}$, conditions have to be imposed on terms with

[TABLE]

In order to maintain a symmetry in the stencil points and coefficients, we use the sum $\left({u_{i+j}^{n}}+{u_{i-j}^{n-l}}\right)$ in the first step, which eliminates the terms with $(\eta,\zeta)=\{(1,0),(3,0)\}$ upon Taylor series expansion. In the second step, conditions are enforced on the terms that are only a function of $\Delta x$ or the terms with $\zeta=0$ by expanding the stencil in space. These are the three terms corresponding to $(\eta,\zeta)=\{(0,0),(2,0),(4,0)\}$, which result in three equations using the function $\left({u_{i+j}^{n}}+{u_{i-j}^{n-l}}\right)$ for $j=0,1,2$:

[TABLE]

After solving the equations we obtain a linear combination of the function that is free from lower order synchronous terms. We can express the linear combination as

[TABLE]

When $l=0$, the terms in the truncation error due to asynchrony disappear from the above expression, and we recover the standard fourth order central difference scheme. In the next step, we eliminate the remaining lower order terms which appear due to asynchrony. If we expand the stencil further in space, which corresponds to $j>2$, then the scheme would possess the required symmetries, but will not reduce to the fourth order central difference in the absence of delay. On the other hand, when the stencil size is increased in time, i.e., $l\in\{{\tilde{k}},{\tilde{k}}+1,{\tilde{k}}+2,\dots\}$, we get a scheme that does resemble a standard central difference. The conditions on the six asynchrony terms from the set in Eq. (19), need Eq. (21) at six time levels. However, with the use of multiple time levels in the stencil, the resulting linear system has three linearly dependent conditions. We find that conditions on terms with the same $\zeta$ are all mathematically equivalent. This reduces the size of the linear system that uses the linear combination in Eq. (21) at $l\in\{{\tilde{k}},{\tilde{k}}+1,{\tilde{k}}+2\}$ to three equations:

[TABLE]

The solution to this linear system results in the scheme

[TABLE]

Clearly, in the absence of delay, ${\tilde{k}}=0$, the scheme reduces to a standard fourth order central difference scheme. Note that, like the scheme in Example 1, the coefficients in the above scheme are a function of the random delay ${\tilde{k}}$. However, unlike Eq. (16) in Example 1, this scheme can take any delay value in the range $[0,L-1]$.

Some other useful examples with their leading order term in the truncation error are collected in Tables LABEL:tab:atschemes-left and LABEL:tab:atschemes-right. These will be used later on when we assess the numerical performance of different AT schemes.

4 Error analysis

In previous sections, we presented a method to derive AT schemes of arbitrary accuracy. As explained earlier, these schemes are, typically, used at PE boundaries ($i\in{I_{B}}$) where asynchrony is experienced. The number of computations that are carried out asynchronously in a domain depend on the number of PEs used to solve the problem, the stencil size of schemes used at interior points and statistics of the random delays, which in turn depend on the characteristics of communications in a computing system. These dependencies bring new challenges when trying to understand the overall accuracy of these AT schemes. First, due to the random nature of the delay, the associated truncation error is also random in nature. Second, schemes to compute spatial derivatives at interior points are not the same as AT schemes at PE boundary points and have, thus, different truncation errors. These issues result in a non-homogeneity of error in the domain, both in space as well as time.

In our previous work [11], we have proposed a statistical description to analyze the overall error and determine the accuracy of the numerical solution. We follow a similar procedure in this work. Before we develop the error analysis, we present some necessary definitions that will be used. First, let us define the probability of having a time level ${\tilde{n}}=n-{\tilde{k}}_{i}$ at a grid point $i$ as ${p_{k[i]}}$. The sum of probabilities of all levels at point $i$ is obviously

[TABLE]

To obtain the statistics of the error, we define two types of averages for a variable $f$: a space average and an ensemble average. The space average can be performed over all points in the set $I$ or the subsets ${I_{I}}$ and ${I_{B}}$. If the average is over the entire domain, that is $i\in I$, it is denoted by angular brackets and given by $\langle f\rangle=\sum_{i=1,N}f_{i}/N$. On the other hand, the average over the points in the subsets ${I_{I}}$ and ${I_{B}}$ are given by $\langle{f}\rangle_{{}_{B}}=\sum_{i\in{I_{B}}}f_{i}/{N_{B}}$ and $\langle{f}\rangle_{{}_{I}}=\sum_{i\in{I_{I}}}f_{i}/{N_{I}}$, respectively. The random nature of delays is taken into account by ensemble averages, which is denoted by an overline $\overline{f}$.

A common measure of the error incurred by using a finite difference representation of the original PDE is given by the so called truncation error. Formally, it is given by the difference between the PDE and the approximate finite difference equation (FDE), that is $E=PDE-FDE$. As introduced in [11], the assessment of the error of asynchrony schemes which are random in nature and heterogeneous in space, can be done by applying the two averages described above. That is,

[TABLE]

where ${E_{i}^{n}}$ is the truncation error at the point $i$ and time level $n$. Due to the non-uniform expression for the truncation error at interior and PE boundary points, it is convenient to split the error according to the two sets of points:

[TABLE]

Note that the error due to interior points does not possess randomness due to delays and are, hence, unaffected by the ensemble average. On the other hand, errors at PE boundary points have both random asynchronous and deterministic synchronous components. This allows us to further split the error in the set ${I_{B}}$ as

[TABLE]

where the subscripts $s$ and $a$ denote the synchronous and asynchronous components, respectively. It is clear that in the absence of delays $\overline{{{\tilde{E}}_{i}^{n}}}|_{a}=0$.

The order of accuracy of a scheme will depend on the leading order term in each of the error terms in the above equation. These terms comprise the sum of the truncation error due to all the derivatives in the original PDE, including the time derivative. Thus, it is important to choose the accuracy of time integration to match the order of accuracy of space derivatives. We will discuss this topic next and then present an example to illustrate the effect of asynchrony on the error. We will end this section with a generalization of the results on accuracy of asynchronous schemes.

4.1 Time integration

To understand the effect of time discretization on the overall order of accuracy, let us consider the equation $\partial u/\partial t=f$, where $f$ depends on spatial derivatives of $u$, integrated using Euler scheme. The scheme is first order in time with the leading order term being $-u^{(0,2)}\Delta t/2$. As mentioned in section 3.1, if we assume a relation of the form $\Delta t\sim\Delta x^{r}$, then the leading order term is equivalent to $\mathcal{O}(\Delta x^{r})$ in space. When the accuracy of the space derivatives is greater than $r$, the total error will, very likely, be dominated by the temporal term and will dictate the order of accuracy of the solution.

Thus, if a certain order is desired for space derivatives, it is important to select a time discretization with the same (or greater) order to keep the overall order unchanged. We will follow this practice as we demonstrate the accuracy of the proposed AT schemes next. For this, we choose a linear multi-step method to compute the time derivative. A general expression with $T$ time steps is given by

[TABLE]

where the coefficients $\beta_{m}$ determine the particular temporal scheme [15].

The advantage of using a temporal scheme of the form Eq. (28) is that the terms $f_{i}^{n-m}$ can be computed using AT schemes and are thus, free of asynchrony errors to the desired order of accuracy. Thus, so will be the linear combination of $f_{i}$ at different time steps. For example, if one uses an AT scheme that is fourth order accurate, with $r=2$ (i.e. $\Delta t\sim\Delta x^{2}$), then one needs a temporal scheme with second order accuracy to maintain fourth order accuracy globally. This can be accomplished by a two-step Adams-Bashforth method

[TABLE]

which is readily shown to be second order in time [15]. The generalization to higher orders is straightforward.

4.2 Example: heat equation with fourth-order accurate AT schemes

Let us consider the 1D heat equation,

[TABLE]

where $u(x,t)$ is the temperature and $\alpha$ is the thermal diffusivity of the medium. The above equation is solved on a uniform grid shown in Fig. 1 with periodic boundary conditions. The equation is approximated with the second order Adams-Bashforth scheme shown in Eq. (29) and standard fourth order central difference for the space derivative at interior points. At the PE boundary points, the space derivative is computed with a fourth order AT scheme which, with delay in the left boundary, is given by Eq. (23) derived in Example 3.

Using Taylor series, the truncation error at interior points is

[TABLE]

As mentioned above, at PE boundary points, the truncation error can be split into the synchronous and asynchronous components,

[TABLE]

Because by construction, the AT scheme in Eq. (23) reduces to the standard central difference in the absence of delay, the synchronous component of the error, ${{E_{i}^{n}}}|_{s}$, is the same as Eq. (31). Note that since this scheme has a uniform delay at buffer points, we drop the subscript for ${\tilde{k}}$ in the above expression for simplicity. The asynchronous component of the error, considering delays only on the left side of the stencil, can be readily shown to be

[TABLE]

The leading order term in the error remains the same when the delays are experienced on the right of the stencil. Clearly, when ${\tilde{k}}=0$, we have ${{\tilde{E}}_{i}^{n}}|_{a,{\tilde{k}}=k}=0$ and thus also its ensemble average. On the other hand, if ${\tilde{k}}\geq 0$, then the ensemble average is

[TABLE]

where moments are given by $\overline{{\tilde{k}}^{n}}=\sum_{k=0,L-1}p_{k}k^{n}$. It is interesting that the average error under the presence of asynchrony, depends not just on the mean of the delay as in [11], but also on its higher order moments. The implication of this result is that in assessing the performance of asynchronous numerical schemes a certain degree of details about the architecture of the computing system would be needed, such as the probability density function of the delays ${\tilde{k}}$. Conversely, one can quantitatively compare the performance of different computing systems by comparing moments of ${\tilde{k}}$.

We now substitute the leading order terms in Eqs. (31) and (34) into Eq. (32). Assuming the statistics of the delays are homogeneous in space, the average error is

[TABLE]

The first two sums on the right hand side are due to synchronous computations and can be conveniently combined by noting that $I=I_{I}\cup{I_{B}}$. To determine the spatial accuracy of the solution, we use the stability parameter $r_{\alpha}=\alpha\Delta t/\Delta x^{2}$ to substitute the time step $\Delta t$ in terms of $\Delta x$. This corresponds to $r=2$ in the formulation presented in section 3.1. The above equation then reduces to

[TABLE]

which can be rewritten as

[TABLE]

The average error is clearly seen to possess components due to synchronous and asynchronous computations. Either of the terms can dominate the average error depending on physical parameters ($\alpha$, initial conditions, etc.), numerical parameters ($\Delta x$, $r_{\alpha}$, etc.), and simulation parameters ($P$, network performance, etc.). If the synchronous part dominates the overall error, then the resulting scheme is fourth order accurate, that is $\langle{\overline{E}}\rangle\sim\mathcal{O}(\Delta x^{4})$.

If, on the other hand, the asynchronous component dominates, the error is given by

[TABLE]

where we have used $N_{B}=(J_{1}+J_{2})P$, with $J_{1}$ and $J_{2}$ being the stencil size in space at interior points. Using $N=\mathcal{L}/\Delta x$, where $\mathcal{L}$ is the length of the domain, and for all other parameters kept constant, the average error is found to scale as

[TABLE]

Interestingly, the order of accuracy of the numerical method now depends on how the problem is scaled on a parallel machine. In the case of weak scaling, where the computational effort per PE is kept constant, that is $P/N=constant$, the error varies as $\Delta x^{4}$ and the method is fourth order accurate in space. On the other hand, when the total computational effort is kept constant ($N=constant$) and the simulations are carried out on increasingly large number of PEs, the average error is $\langle{\overline{E}}\rangle\sim\mathcal{O}(\Delta x^{5})$ and the method is fifth order accurate. We also observe that the error scales linearly with $P$.

In some situations the error due the synchronous and asynchronous components may be comparable. In such cases, the overall error will depend on the sign of each contribution. If the synchronous and asynchronous components have opposite signs, then it is possible to expect some error cancellation. The order of accuracy though would remain unaltered.

4.3 Generalization

We now proceed to generalize the expressions for average error ($\langle{\overline{E}}\rangle$) presented in Eq. (39). For this, we restate the conditions and assumptions that lead to Eq. (39). First, we assumed that the asynchronous component dominates the overall error. Second, the AT scheme used in the analysis is fourth order accurate, which lead to an $\mathcal{O}(\Delta x^{4})$ leading term due to asynchrony. We have also assumed a uniform random delay in the stencil, that is ${\tilde{k}}_{i+j}={\tilde{k}}$ for all $i+j\in B$. Also, the scheme uses three successive asynchronous time levels (with delays ${\tilde{k}}$, ${\tilde{k}}+1$, ${\tilde{k}}+2$), which results in a cubic polynomial in $k$ in the leading order error term.

With the above observations, we can arrive at a general case which uses AT schemes with $\mathcal{T}$ number of successive asynchrony time levels and is accurate to an order $a$. If the asynchronous component of the error dominates the average error, then it is easy to generalize Eq. (39) as:

[TABLE]

Note that the average error still scales linearly with the number of PEs. However, higher order moments of the delay are necessary to characterize the error when the stencil size of AT schemes is expanded in time. A minimum accuracy of order $a$ is then assured, regardless of how simulations are scaled up.

5 Numerical Simulations

In this section, we verify the numerical performance of AT schemes. Let us consider the general PDE:

[TABLE]

where $\cal{D}$ is the highest derivative and the coefficient $\beta_{d}$ determines the characteristics of the physical process associated with the $d$-th derivative. Of particular interest are the heat equation (${\cal D}=2$ with $\beta_{1}=0$ and $\beta_{2}=\alpha$), and the advection-diffusion equation (${\cal D}=2$ with $\beta_{1}=c$ and $\beta_{2}=\alpha$), where $\alpha$ is the thermal or viscous diffusivity and $c$ is the advection speed. When the coefficients $\beta_{d}$ are constant, Eq. (41) is linear and usually possesses an analytical solution, which will be used here to evaluate the error in numerically computed solutions. The so-called nonlinear viscous Burgers’ equation, which is widely used in understanding physical properties of fluid flows, is obtained with ${\cal D}=2$, $\beta_{1}=u(x,t)$ and $\beta_{2}=\alpha$. We also perform simulations of this equation to demonstrate the feasibility of AT schemes in solving multi-scale phenomena with non-linear couplings.

5.1 Simulation details

The equations described in the above section are solved in a periodic domain of length $2\pi$. For initial conditions, we use a multi-scale spectrum given by superimposed sinusoidal waves:

[TABLE]

where $\kappa$ denotes the wavenumber. $A(\kappa)$ and $\phi_{\kappa}$ are the amplitude and phase angle corresponding to each wavenumber $\kappa$. The phase $\phi_{\kappa}$ are included in order to avoid circumstances like the coincidence of PE boundaries with zero-gradients in the function, which may result in very special cancellations of some of the error terms due to asynchrony. The results presented below are in fact ensemble averages of multiple simulations with different phases. In addition to avoiding special cases in terms of accuracy as mentioned above, this procedure provides a probability space over which ensemble averages can be obtained.

Simulations are carried out using several configuration cases with different governing equations to study the behavior of AT schemes in different regimes. The synchronous computations of spatial derivatives are carried out using standard central difference schemes. Close to PE boundary points, the AT schemes summarized in Tables LABEL:tab:atschemes-left and LABEL:tab:atschemes-right are used. The time derivatives are discretized according to the procedure described in section 4.1. The details of each numerical experiment are tabulated in Table 2.

In numerical simulations, we use a random number generator to simulate communication delays (${\tilde{k}}_{j}$) at PE boundaries. This provides a complete control over the statistics of the delays, thus, allowing us to compare the results against the theoretical predictions in different parameter regimes. At each time advancement, the delay at a buffer point is computed from a random number drawn with a given initial seed from a uniform distribution in the interval $[0,1]$. This interval is divided into $L$ bins according to the probabilities $\{{p_{0[j]}},{p_{1[j]}},\dots,{p_{L-1[j]}}\}$ corresponding to delays ${\tilde{k}}_{j}=0,1,2,\dots,L-1$, respectively. When a random number is drawn, it is matched with the corresponding bin which determines the delay. As we use i.i.d. random sequences at different PE boundaries in the simulations, there is no dependence on the location and hence, we drop the subscript $j$ in probabilities for simplicity and write $\{p_{0},p_{1},\dots,p_{L-1}\}$. As an example, if we choose $L=3$ and the set $\{p_{0},p_{1},p_{2}\}=\{0.6,0.3,0.1\}$, then the probability of having ${\tilde{k}}=0$, ${\tilde{k}}=1$ and ${\tilde{k}}=2$ is $0.6$, $0.3$ and $0.1$, respectively. In the case of schemes which use a uniform delay in their stencil (like in Eq. (23) of Example 3), a single random number is drawn at each PE boundary to obtain the uniform delay at all the buffer points at that PE boundary.

The error is computed by comparing the numerical solution against the analytical solution. With periodic boundary conditions and an initial condition given in Eq. (42), the analytical solution (denoted by subscript $a$) for the linear advection-diffusion equation is

[TABLE]

For the heat equation, the analytical solution is given by the above expression with $c=0$. In the case of nonlinear Burgers equation, the error is evaluated against the solution from a highly resolved simulation. The error at a point $i$ and time level $n$ is computed as ${E_{i}^{n}}={u_{i}^{n}}-u_{a}(x_{i},t_{n})$. The overall error in the domain is obtained using the different averages presented in section 4.

5.2 Results

5.2.1 Linear equations

Fig. 6 shows results from simulations with three different schemes: synchronous (solid black lines), asynchronous-standard (dashed red lines) and AT (green lines) schemes. Note that by asynchronous-standard schemes we mean standard (synchronous) schemes used in an asynchronous fashion. The governing equation and schemes are those corresponding to Case 1 in Table 2. The simulation parameters are $N=128$, $\kappa=\{2,3,5\}$ (the vector $\kappa$ here contains the wavenumbers used in the initial condition defined by Eq. (43)), $P=4$, and three allowable time levels for asynchronous computations according to $\{p_{0},p_{1},p_{2}\}=\{0.5,0.3,0.2\}$. In part (a) of the figure, we show the time evolution of function $u$. We observe that the initial condition which is a combination of sine waves is convected with a wave speed $c$ and simultaneously damped due to diffusive action, as expected. To highlight the differences between these cases, we show the evolution of the error in part (b) of the figure. The error in the case of asynchronous standard schemes is an order of magnitude greater near the PE boundaries (indicated by vertical dash-dotted lines). As discussed in [11], this is due to the asynchrony in the data available at buffer points. This error, which is initially localized near PE boundaries, propagates into the interior with time. In the case of AT schemes, which are designed to mitigate the affect of asynchrony, the error at PE boundaries is of similar magnitude as the synchronous schemes.

To verify the formal order of accuracy of AT schemes and the effect of simulation parameters ($N$, $P$, $\overline{{\tilde{k}}}$, etc.) on the overall error, we now proceed to the statistical description of the error. An example of the effect of asynchrony on the overall error for standard central difference schemes (Case 4 in Table 2) is shown in Fig. 7. For the fourth-order scheme used in these simulations, the error for $p_{0}=1.0$ decreases with a slope of $-4$, as expected in synchronous computing. In the presence of asynchrony ($p_{0}<1.0$) the slope reduces to $-1$, depicting a first order accurate solution [11]. Also, the absolute error for a given grid resolution increases when asynchrony is increased ($p_{0}$ is reduced). This drastic decrease in accuracy is mitigated when AT schemes are used, as we show next.

Fig. 8 shows the effect on asynchrony in different configurations when AT schemes are used. The parameters used in these numerical experiments are: $\kappa=\{3,4,5\}$, $A(\kappa)=\{2.0,0.5,1.5\}$, $P=16$, $L=3$. Different colors in the graphs represent results from different sets of $p_{k}$ ($k=0,1,2$) and their values are given in the caption of the figure. Part (a) shows results from Case 2 in Table 2. A second-order central difference scheme for synchronous computations and a first-order asymmetric stencil AT scheme for asynchronous computations are used. The error in absence of delays (red line) decreases with a slope of $-2$ as expected. We also observe that, asymptotically, this is also the case for the asynchronous cases ($p_{0}<1$). The reason for second-order accuracy even with a first-order AT scheme can be explained with the strong scaling argument presented in Eq. (38). The effect of an increase in the amount of asynchrony, is seen to increase the magnitude of error leaving the asymptotic rate of convergence unchanged.

In part (b), we show results for Case 3 in Table 2. A second-order accurate asymmetric stencil AT scheme is used for asynchronous computations. In this case too, we see an asymptotic convergence rate of order 2. Also, the magnitude of error increases with the amount of asynchrony. Note that in both parts (a) and (b), the AT schemes are constructed by expanding the stencil in space to improve the accuracy in the presence of asynchrony. These schemes do not reduce to or have the same form as the central difference schemes when ${\tilde{k}}=0$.

In parts (c) and (d) of Fig. 8, results are shown for AT schemes derived by expanding the stencil in time (instead of space) to maintain accuracy in the presence of asynchrony. These are cases 4 and 5 in Table 2. These schemes have symmetric stencils and coefficients, and they reduce to central difference schemes when ${\tilde{k}}=0$. As expected from the theory, the error in part (c) converges with an accuracy of order 4. The effect of asynchrony is hardly noticeable at higher resolutions. This can be attributed to the fact that the AT schemes, in this case, reduce to the synchronous central schemes in absence of asynchrony and result in a homogeneous synchronous truncation error terms across the domain. A similar observation is also found in part (d), which uses a sixth order accurate scheme for space derivative.

Let us now recall Eq. (40), which describes the scaling of the average error $(\langle{\overline{E}}\rangle)$ with $P$, $N$ and ${\tilde{k}}$:

[TABLE]

Note that this scaling holds only when the error due to asynchrony dominates the overall error. Otherwise, the error in the leading order may have a synchronous component and may show a different dependence on simulation parameters.

In Fig. 9 we show numerical data from Case 1. According to Eq. (44), the order of accuracy is one more than the order of the AT scheme when $P$ is fixed. Part (a) of the figure shows the convergence of the error for different $P$. For low $P$, the leading order terms of the error contain both synchronous and asynchronous contributions, and thus shows a convergence slope between $-3$ and $-2$. However, when $P$ increases, the error due to asynchrony dominates and shows a convergence of $-3$, as predicted by the theory. The linear scaling of the error with $P$ is verified in the inset of part (a). Results for weak scaling, that is when both $P$ and $N$ in increase such that $P/N$ is kept constant, are shown in part (b). As expected, the error for different $P/N$ asymptotically converges to second order accuracy. In the inset of the figure, an inverse dependence of the error on $N/P$ is observed as predicted also by the theory.

Unlike standard synchronous schemes when asynchrony is allowed, for which the error due to asynchrony depends only on the average of delays ($\overline{{\tilde{k}}}$) [11], the error for AT schemes can also depend on higher order moments of ${\tilde{k}}$, as shown in section 4. This dependence is indeed confirmed in Fig. 10. Results in parts (a) and (b) of the figure are for Cases 3 and 4, respectively and can be understood as follows. For these schemes, which use two and three delayed time levels in the stencil ($\mathcal{T}=2$ and $3$), respectively, the average error scales as:

[TABLE]

To simplify the dependence of error to a single variable, the probability of occurence of a level $k$ for a given $L$ is chosen as $p_{k}=1/L$. For example, if $L=3$, then $\{p_{0},p_{1},p_{2}\}=\{1/3,1/3,1/3\}$. This reduces the scaling of the average error to

[TABLE]

As mentioned earlier, this scaling is only valid for the asynchrony component of the error or for the total error when the former is dominant to leading order. Thus, to compare to Eq. (46), we compute the total error and subtract it from a completely synchronous, but otherwise identical, simulation. In Fig. 10 we show the thus obtained asynchronous part of the error, that is $(\langle{\overline{E}}\rangle-\langle{\overline{E}}\rangle_{s})/\langle{\overline{E}}\rangle_{s}$. The dashed curves in the graphs correspond to Eq. (46). There is very good agreement between the theoretical prediction and the data from simulations.

5.2.2 Nonlinear equations

In the above numerical experiments, we have verified the performance of AT schemes for linear equations. However, a number of natural and engineering systems are governed by highly nonlinear processes like fluid turbulence phenomena. The viscous Burgers’ equation is often used as a proxy to understand these nonlinear effects in fluid flows with negligible pressure effects. Thus, we use this equation to assess AT schemes in a more realistic setup. Fig. 11 shows the convergence of the error for the fourth-order schemes described in Table 2 as Case 6. Clearly, even with an increase in the degree of asynchrony the error converges with fourth-order accuracy.

In nonlinear problems like turbulence, one is interested not only in statistical moments of the velocity field but also in its gradients as they exhibit very strong but localized fluctuations, a phenomenon known as intermittency [16]. Thus, we investigated the variation of central moments of velocity ($u$) and velocity gradients ($\partial u/\partial x$) with the resolution, $N$. An example is shown in Fig. 12. For this problem, most of the contribution to the velocity field comes from low wavenumbers, while most of the contribution for its gradients comes from high wavenumbers. Thus, it is not surprising that asynchrony effects are more evident for velocity gradients at low grid resolution. Nevertheless, both synchronous and asynchronous cases seem to converge at the same grid resolution ($N=256$). This is consistent with the results for lower-order statistics in [11].

While our numerical experiments show accurate results even for high order statistics of velocity gradients, this result is not expected to be general for other equations or in higher-dimensional spaces. This is indeed an area that needs further investigation.

6 Conclusions

A number of natural and engineering systems are governed by PDEs that present solutions with wide range of scales which can only be captured by high-fidelity simulations (using high-order numerical schemes) on massive computational systems. At extreme scales, global communications and synchronizations will likely become an obstacle to sustained performance for number of scientific codes. In this work we have presented a general methodology to analyze and derive schemes that remove these two main obstacles by allowing some tunable level of asynchrony.

The concept relies on finite differences to approximate derivatives of general order using values of the function from neighboring points. Close to PEs boundaries current computational methodologies are stalled until communications between PEs is completed. In previous work we have shown that one can relax this forced synchronization at the mathematical level such that computations can proceed using values from past time levels. Here we generalize the concept, established conditions under which schemes can be obtained, classified the resulting schemes in terms of their properties, and provided a general framework in which schemes of arbitrary order can be obtained for any derivative of a function. These schemes are referred as asynchrony-tolerant or AT schemes.

By analyzing in detail the truncation error of general finite differences when asynchrony is allowed, we described the mathematical conditions needed to obtain a scheme of arbitrary order under asynchronous conditions. In particular, we showed that asynchrony errors can be eliminated either by extending the stencil in space as well as in time. These two alternatives lead to schemes with different properties and limitations. However, depending on the order of the scheme and the type of asynchrony allowed (e.g. on both sides of the stencil, uniform across the stencil, etc.) not all expansions of the stencil size will result in an AT scheme. The kind of stencils that do lead to a scheme, has also been presented and requires identifying the nature of terms present in the truncation error. The coefficients are obtained by solving a linear system of equations. An alternative method was also presented where successive terms in the truncation errors are eliminated by a step-by-step method. The process ends when the desired accuracy is achieved.

The resulting schemes can be classified on the nature of their coefficients. We presented four conditions for the classification: (i) symmetric layout of grid points, (ii) unconstrained or uniform delay on boundary points, (iii) artificial delay at interior points, and (iv) symmetry on the coefficients. Each have different numerical and performance properties. Actual examples of these different schemes were also put forth.

The truncation error was analyzed in a statistical framework that takes into account both the stochasticity of delays as well as the non-uniformity of delays in space. We have further shown that multi-step time-integration methods can be used successfully to obtain solvers of arbitrary order using AT schemes. The general form of the error, given by Eq. (40), shows that the average error depends on the number of processors as well as moments of the distribution of delays, which in turn depend on the characteristics of the computing system simulations are run on.

Theoretical predictions on the accuracy of schemes were compared to numerical experiments for linear as well as non-linear equations. Good agreement was found across different parameter space. The work presented here provides a strong foundation for mathematically asynchronous computing methods for PDEs at extreme scales. Application of this method to more complex phenomena and realistic conditions is a part of our ongoing research.

7 Acknowledgments

The authors gratefully acknowledge NSF (Grant OCI-1054966 and CCF-1439145) for financial support. The authors also thank NERSC and XSEDE for computer time on their systems. The authors benefited from discussions with Lawrence Rauchwerger, Raktim Bhattacharya and Jacqueline H. Chen.