A high-order, conservative integrator with local time-stepping^††thanks: Submitted to the editors 10 October 2019. \funding This work was supported in part by the Sherman Fairchild Foundation and NSF grants PHY-1606654 and ACI-1713678 at Cornell.

Suppose we wish to numerically solve a set of coupled first-order ordinary differential equations

with $\mathbf{y}=\mathbf{y}_{0}$ at some initial time $t_{0}$. Here $\mathbf{D}(\mathbf{y})$, the time-derivative operator, is the right-hand side evaluated when the system is in state $\mathbf{y}$. (Throughout this work we will restrict our attention to autonomous systems. As any non-autonomous system can be recast in an autonomous form (see, for example, Section II.2 of [18]) we lose no generality from this restriction.) A common method is to solve for the variables at a (monotonic) sequence of times $t_{0},t_{1},\ldots$ using a $k$th-order Adams-Bashforth method

with $\Delta\mathbf{y}_{n}=\mathbf{y}_{n+1}-\mathbf{y}_{n}$, $\Delta t_{n}=t_{n+1}-t_{n}$ and the coefficients corresponding to the step given by [29]

Here

are Lagrange polynomials. The values required for evaluating the first step can be obtained using a standard Adams-Bashforth GTS self-start procedure [23], after which the step sizes can be adjusted to the desired ratio.

If different degrees of freedom require different time steps for stability, it may be desirable to evaluate these variables at different frequencies, in order to avoid unnecessary computations for the more stable variables. Suppose we divide the components of $\mathbf{y}$ into $S$ sets $\mathbf{y}^{1},\ldots,\mathbf{y}^{S}$, each of which is a collection of degrees of freedom that are always evaluated at the same times, $t^{s}_{0},t^{s}_{1},\ldots$. We can then split (1) into an equation for each of these sets:

where $\mathbf{D}^{s}$ is the result of $\mathbf{D}$ restricted to the set $s$. Any attempt to use this equation to perform an LTS evolution immediately encounters the problem that evaluating its right-hand side requires knowing the entire state of the system, which conflicts with the goal of independent evaluation times for different degrees of freedom.

A linear conserved quantity is a quantity $C$ expressible as an inner product of a vector $\mathbf{c}$ with the evolved variables (treated as a vector)

with

for all values of $\mathbf{y}$. Such a quantity is constant under exact integration of the system and under integration using Euler’s method. An integrator is called (linearly) conservative if all such quantities remain constant when integrating a system using it [28]. It is desirable for an integrator to keep such quantities precisely constant (up to roundoff error) rather than merely constant up to the truncation error of the scheme. Such quantities often have an intuitive physical meaning, and frequently even a small rate of drift can cause qualitative changes in the evolution of the system.

When solving a PDE representing a physical system, the most common linear conserved quantities are integrals over the computational domain of fields representing densities. The vector $\mathbf{c}$ in these cases is the vector of coefficients necessary to perform a numerical integral. In a discontinuous Galerkin scheme these coefficients would combine quadrature weights on the elements and factors arising from coordinate mappings of the elements relative to their canonical shapes.

Let us first consider as an example the case of a second-order scheme on two sets, $A$ and $B$, with $B$ being evaluated twice as often as $A$. Call their step sizes $\Delta t^{A}$ and $\Delta t^{B}=\Delta t^{A}/2$. This step pattern is shown in Figure 1, where for simplicity we consider the steps starting from $t=0$ leading up to $t=\Delta t^{A}$. There are three types of steps to consider: the large step on set $A$, labeled $(a)$, and the first and second halves of that step on set $B$, labeled $(b)$ and $(c)$. This case is considered in Sandu and Constantinescu [25], but the method presented there only provides a second-order value when sets have stepped to the same time; intermediate values are only accurate to first order.

We will start with the small step $(b)$. For a GTS Adams-Bashforth method, this step would be given by

where $\tilde{\mathbf{D}}^{B}(t)$ is an approximation to the derivative at time $t$ that we will now construct.

At time $t=0$ we have values for the entire system, so we can obtain a derivative by applying the derivative operator to the initial data, so we take $\tilde{\mathbf{D}}^{B}(0)=\mathbf{D}^{B}(\mathbf{y}^{A}(0),\mathbf{y}^{B}(0))$. We cannot evaluate $\tilde{\mathbf{D}}^{B}(-\Delta t^{B})$ in this manner, however, because we do not have data for $\mathbf{y}^{A}$ at $t=-\Delta t^{B}$, so we must construct it from the values at $t=0$ and $t=-2\Delta t^{B}$. Up through this point, the derivation matches the results of Sandu and Constantinescu [25], but they now choose to use the known value of $\mathbf{y}^{A}(-2\Delta t^{B})$ to evaluate the derivative, while we will search for a more accurate approximation. There are two reasonable choices of how to do this: average the known values of $\mathbf{y}^{A}$ to get a value at the desired time and use that to apply the derivative operator, or apply the derivative operator at both times (using the value of $\mathbf{y}^{B}$ at $-\Delta t^{B}$ both times) and average the results. We choose the latter, taking $\tilde{\mathbf{D}}^{B}(-\Delta t^{B})=[\mathbf{D}^{B}(\mathbf{y}^{A}(0),\mathbf{y}^{B}(-\Delta t^{B}))+\mathbf{D}^{B}(\mathbf{y}^{A}(-\Delta t^{A}),\mathbf{y}^{B}(-\Delta t^{B}))]/2$, so (8) becomes

The error in averaging the derivatives is of order $(\Delta t^{B})^{2}$, so it introduces an error of order $(\Delta t^{B})^{3}$ in the value after the step, preserving the second-order quality of the base GTS method.

The second small step, $(c)$, proceeds similarly, except that we now use a derivative at $\Delta t^{B}$ instead of $-\Delta t^{B}$. Instead of averaging the derivatives at different $\mathbf{y}^{A}$ we must therefore perform a (linear) extrapolation to obtain our approximate derivative $\tilde{\mathbf{D}}^{B}(\Delta t^{B})=[3\,\mathbf{D}^{B}(\mathbf{y}^{A}(0),\mathbf{y}^{B}(\Delta t^{B}))-\mathbf{D}^{B}(\mathbf{y}^{A}(-\Delta t^{A}),\mathbf{y}^{B}(\Delta t^{B}))]/2$. Thus, we obtain the rule

We could use the same procedure to evaluate the large step $(a)$, but, as this would not take into account the value $\mathbf{y}^{B}(\Delta t^{B})$ used for taking the second small step, there is no way this procedure could be conservative. This, however, gives us a hint as to how to proceed: we treat the large step as having two internal steps, one for balancing each of the small steps. In fact, in order to remain conservative, we must take each of these internal steps using the same scheme as for the corresponding small step, except using the part of the derivative corresponding to set $A$. This can be seen by considering a generic pair of methods for a step: $\Delta\mathbf{y}^{A,B}=\sum_{i}k^{A,B}_{i}\mathbf{D}^{A,B}(\mathbf{q}_{i})$, where the $k^{A,B}_{i}$ are given coefficients and the $\mathbf{q}_{i}$ are vectors obtained in some manner from the known values of $\mathbf{y}$. The change in a linear conserved quantity during that step is

The first term vanishes by (7), so the only way for two sets to take equal-sized steps in a conservative manner is if they use the same step rule. The procedure for the large step can therefore be found by summing (9) and (10), giving

Note that the coefficients have changed by a factor of $2$ compared to the previous equations because of the change of the leading coefficient to $\Delta t^{A}$. As the two small steps were accurate to second order and this is effectively their concatenation, it is also accurate to second order. A detailed analysis of the consistency and convergence of the general method described below, of which this is a special case, is presented in Appendix A.

Let us turn now to the task of finding a general conservative, high-order LTS integrator. First, we will consider the implications of requiring an Adams-Bashforth-like LTS scheme to be conservative. For such a scheme it only makes sense to evaluate (6) at times at which all the degrees of freedom are evaluated. We therefore introduce a new quantity $\tilde{\mathbf{y}}_{n}$ that is defined for the entire set of degrees of freedom for each time $\tilde{t}_{n}$ at which any set is evaluated, and is equal to the (as yet unknown) complete numerical solution at all grid points wherever it exists. If we provide a rule to compute $\tilde{\mathbf{y}}_{n+1}$ from earlier values then, as long as portions of $\tilde{\mathbf{y}}_{n}$ that do not correspond to values in one of the $\mathbf{y}^{s}_{m}$ are never used, we can obtain an LTS method by summing the $\Delta\tilde{\mathbf{y}}_{n}$ between evaluations. Furthermore, if the step from $\tilde{\mathbf{y}}_{n}$ to $\tilde{\mathbf{y}}_{n+1}$ is conservative, then the implied full method will be as well.

The condition for this small step to be conservative is

This is satisfied if we evaluate $\Delta\tilde{\mathbf{y}}_{n}$ using a standard Adams-Bashforth method, but that would require values of $\tilde{\mathbf{y}}_{n}$ that are not included in any of the $\mathbf{y}^{s}_{m}$. Comparing (7) and (14), we see that we will obtain a conservative method if we take

for some coefficients $\beta_{n;i^{1}\cdots i^{S}}$. The choices of these coefficients are not unique, but there is a natural choice. We evaluate each step using a standard order-$k$ Adams-Bashforth scheme, but instead of using the derivatives of the function that we cannot evaluate, we use approximate derivatives $\tilde{\mathbf{D}}_{n}$. As long as these are accurate to order $k-1$, we will lose no formal accuracy for the step. We evaluate $\tilde{\mathbf{D}}_{n}$ by treating $\mathbf{D}(\mathbf{y}^{1}(t_{1}),\ldots,\mathbf{y}^{S}(t_{S}))$ as a function of the times $t_{1},\ldots,t_{S}$ independently, and then performing a multidimensional interpolation from the known values in the space of the evaluation times on each set. To obtain the required accuracy, we must use evaluations from at least $k$ times from each set, and it is natural to choose the most recent values. The known values of $\mathbf{D}$ then form a lattice in the multidimensional space. Multidimensional interpolation from such a lattice is not unique, but a natural choice is to perform it as a series of one-dimensional interpolations.¹¹1This freedom arises from the fact that the system of equations defining this interpolation is underdetermined for $S$ and $k$ greater than 1: we must find $k^{S}$ fitting coefficients but there are only $\binom{k+S-1}{S}$ monomial terms of degree less than $k$ (which are the ones relevant for an order $k$ fit). A general choice of interpolation coefficients will result in an interpolating polynomial containing all terms of degree less than $k$ in each of the $t^{s}$ individually. We therefore have the freedom to modify the interpolation coefficients as long as the modification alters only terms in the interpolating polynomial of total degree at least $k$. This freedom could be used, for example, to set certain coefficients to zero to reduce the number of computations required or to decrease the effect of terms where the times on different sets have large mismatches. In the case where the step size on each set is constant, the alternative sets of interpolation coefficients can be obtained by adding high-order products of discrete Chebyshev polynomials [7] to the coefficients in (19). In the general case we know of no simple method to calculate alternative coefficients. We have not investigated the use of such alternative coefficients in either of these cases. Combining all these ideas, we have

where $\tilde{\alpha}_{ni}$ are Adams-Bashforth coefficients corresponding to the sequence of times $\tilde{t}_{n}$ and $m^{s}(n)$ is defined by $t^{s}_{m^{s}(n)}\leq\tilde{t}_{n}<t^{s}_{m^{s}(n)+1}$, i.e., it is the index of the last evaluation on set $s$ that can influence the step $\Delta\tilde{t}_{n}$ (see Figure 2). The interpolation coefficients are given by

For computational purposes, it is useful to rewrite these steps as

where the coefficient is

The full change in the value of a given set of degrees of freedom over an entire step can then be obtained by summing the contributions of all these small steps. This will give for each set of degrees of freedom an equation of the form

for some coefficients $a^{s}_{m;q^{1}\cdots q^{S}}$. Setting (20) equal to this sum over (18) gives the coefficients

The equations (20) involve many more evaluations of the derivative than the standard GTS Adams-Bashforth method, so in this form the LTS method is unlikely to be more efficient. However, if the couplings between the sets of degrees of freedom are inexpensive to calculate compared to the interactions within each set, then the required number of evaluations can be reduced.. Let us suppose that the derivative on set $s$ is split into a “volume” portion depending only on set $s$ itself and a “boundary” portion encoding the coupling to other sets:

These names are motivated by finite volume and discontinuous Galerkin methods, where the terms from the interior and boundaries of elements split in this manner. Substituting (22) into (18) and summing over the small steps, the volume contribution to the full step on set $s$ is

where $n^{s}(m)$ is defined by $\tilde{t}_{n^{s}(m)}=t^{s}_{m}$ (see Figure 2). This is the same form as the GTS Adams-Bashforth method (2) using the bracketed expression as coefficients (absorbing the $\Delta t$ factor). The bracketed expression does not depend on the form of the derivative, so to evaluate it we can take the boundary coupling $\mathbf{B}^{s}$ to be zero, in which case this is the only contribution to the step. As this is then a $k$th-order GTS method and the Adams-Bashforth method is the unique $k$th-order method of this form, the bracketed quantity must be the standard Adams-Bashforth coefficient. Returning to the general case with a coupling, this shows that a set of degrees of freedom can be evolved using the standard Adams-Bashforth method for the volume portion with only the coupling terms evaluated using (18).

This simplification applies in intermediate cases as well: if the full derivative can be split into portions each of which depends on only some of the degree-of-freedom sets, each of those contributions to the step can be calculated independently using (18) ignoring non-contributing sets. In calculations where the sets are only coupled pairwise, this implies that only the $S=2$ case need be considered.

As mentioned previously, a major intended application of this time-stepping scheme is to evolution of PDEs using discontinuous Galerkin methods. In such a method, the sets of degrees of freedom are only coupled pairwise and the update method reduces to a collection of standard Adams-Bashforth methods and LTS methods with $S=2$. For the two-set case, we call the sets $A$ and $B$ and define the selection functions $\Theta^{A}_{n}$, $\Theta^{B}_{n}$, and $\Theta^{AB}_{n}$ to be one if $\tilde{t}_{n}$ is an evaluation time for only set $A$, only set $B$, or both sets, respectively. By construction, a time evaluated on neither set can never occur. These selection functions sum to one, so we can write $\beta_{n;q^{A}q^{B}}=\beta^{A}_{n;q^{A}q^{B}}+\beta^{B}_{n;q^{A}q^{B}}+\beta^{AB}_{n;q^{A}q^{B}}$ with

By the definition of $m(n)$, $\tilde{t}_{n}$ is not older than $t^{A,B}_{m^{A,B}(n)}$, so, from the construction of the $\tilde{t}_{n}$ we see that $\tilde{t}_{n-i}\geq t^{A,B}_{m^{A,B}(n)-(k-1)}$. This implies that if $\tilde{t}_{n-i}$ is an evaluation time for either set, it is one of the control points in the corresponding Lagrange polynomial. We can therefore collapse those polynomials to obtain

and $\beta^{B}_{n;q^{A}q^{B}}=\beta^{A}_{n;q^{B}q^{A}}$. The expressions should be taken to be zero when the conditionals on the right are not satisfied. Some example values are shown in Tables 1 and 2. The meaning of, for example, the first entry for (a) in Table 1 is that in (20) the coefficient $a^{A}_{0;0,1}=115/64$ (where we have chosen to number the steps starting from $t=0$ so the step on set $B$ at $\Delta t^{B}$ is step 1), so the equation for this step begins

Similarly, the lower-left entry for (c) in Table 2 indicates that one term in the second small step is

Additional tables of coefficients can be found in Appendix D.

For our first test problem, we wanted a simple nonlinear problem to test the convergence and conservation properties of the method. Accordingly, we evolved the Burgers equation:

The DG elements were coupled using the numerical flux of Harten, Lax, and van Leer [19]. In the common case where the solution has the same sign on both sides of the interface between two elements, this reduces to an upwind flux. For this PDE, the spatial integral of the field $u$ is a linear conserved quantity, and this carries over to the discretized system when using a DG scheme.

For our first test, we considered the exact solution

We used the analytic value of $u(t=-1/8,x)$ as an initial condition and evolved $u$ from $t=-1/8$ to $t=3/2$ on the interval $-9/8\leq x\leq 1/8$, with free boundary conditions. The interval was divided into 16 elements of equal size, each of which contained 10 nodes distributed as Legendre-Gauss-Lobatto points. These values reduce the spatial truncation error to approximately $10^{-15}$, smaller than the floating-point roundoff error, so it can be ignored for this analysis.

To simplify bookkeeping, the size of each time step was restricted to be a power of two, but each element was allowed to independently adjust its step size to conform to the CFL limit. We found that increasing the step size too rapidly led to a growth in the error, so each step size increase was limited to a factor of two and an increase was only allowed if the previous $k-1$ steps were all of the same size, where $k$ is the integration order. All elements were initialized with a very small time step of $2^{-27}$. Figure 3 shows a typical step pattern.

The sole characteristic speed in the Burgers system is $\lvert u\rvert$, so the CFL condition predicts that the stable step size will be proportional to $1/u$. The coefficient of proportionality can be estimated from the point distribution and properties of the time stepper, but in order to investigate the convergence of the method we kept it as a free parameter. The convergence as a function of this parameter is shown in Figure 4.

The previous test does not demonstrate conservation because of flow through the domain boundaries. To eliminate this effect, we perform an evolution of the Burgers equation on a periodic domain. We use the same domain decomposition as for the previous example, but identify the endpoints at $-9/8$ and $1/8$. We initialize the field at $t=0$ to