The use of numerical simulation for prediction of characteristics of chaotic dynamical systems inherently involves unpredictable processes. In this work, we develop a model for the expected error in the simulation of ergodic, chaotic ordinary differential equation (ODE) systems, which allows for discretization and statistical effects due to unpredictability. Using this model, we then generate a framework for understanding the relationship between the sampling cost of a simulation and the expected error in the result and explore the implications of the various parameters of simulations. Finally, we generalize the framework to consider the total cost—including unsampled spin-up timesteps—of simulations and consider the implications of parallel computational environments to give a realistic model of the relationship between wall-clock time and the expected error in simulation of a chaotic ODE system.

Progress in computational processing capabilities and numerical methods has led to increasing use of unsteady, scale resolving simulation methods for turbulent flows for applied research.1–3 In spite of advances, the cost of these simulations, particularly large-eddy simulation (LES), remains an obstacle to industrial use. Cost reduction remains a key research area for these systems, including more aggressive modeling approaches.1,4 On the other hand, advancements in numerical methods have demonstrated that a stable, high-order solution of turbulent systems is possible5 with ongoing developments enabling new efficiencies.6 As these advances approach maturity, a key outstanding question is how to use them to extract estimates as efficiently as possible when simulations are unsteady and chaotic and discretization methods are not significantly constrained by explicit stability limits.

The resolution of the turbulent—and therefore chaotic—nature of the underlying flows is a primary complication for estimation using these more complex simulations. For chaotic systems, estimation of long-time behavior is challenging, because the ordinary differential equations (ODEs) to which they can be reduced have limited predictability.7 Of the general class of chaotic systems, a subset is ergodic systems, whose long-term states are drawn from a stationary distribution, independent of initial conditions (ICs).8 For ergodic chaotic problems, we frequently want to quantify the unique infinite-time average of some instantaneous quantity of interest of the system

J=limT1T0Tg(u(t))dt,
(1)

where g is the instantaneous output functional and the state u(t) is governed by a dynamical system of the form

dudt=f(u),
(2)

with a given initial condition (IC), u(0)=uIC.

Often, the complexity of chaotic systems of interest is high, and accordingly, the cost of an accurate computational estimate of J becomes formidable.9–12 As the scale, cost, and complexity of a computational simulation increase, efficient discretization methods become critical for accurately estimating quantities of interest.

Understanding the error in approximations of J is nontrivial because statistical errors (errors due to a finite-time approximation) and discretization error (error due to a numerical approximation of solutions) are always simultaneously present. In the largest direct numerical simulation (DNS) and large eddy simulation (LES) cases, for example, it is typical to fix the period of simulation time over which an average is performed at some large number of characteristic times and validate that discretization error converges as expected, assuming negligible sampling errors.13–16 A recent work has sought to quantify the effect of statistical error more robustly, using a turbulent flow theory,17 advanced spatiotemporal statistical post-processing methods,18 statistical windowing techniques,19 or by extending the concept of Richardson extrapolation to chaotic flows using auto-regressive models and Bayesian methods.20 The latter work is notable for its use to estimate the statistical errors in the DNS of a high- turbulent channel flow.21 

The objective of this paper is to investigate the behavior of statistical and discretization errors as a function of the computational cost for ergodic systems. Following a similar approach of Oliver et al.,20 we propose a simple error model for finite-time, discrete approximations of infinite-time averages on attractors. Using the Lorenz system as an example, we demonstrate that the discretization error converges as the time step size decreases. However, if we denote by “sampling time” the amount of simulation time on which the integral is taken (after discarding spin-up), we show that error does not increase exponentially with sampling time as might be expected from classical numerical analysis but rather asymptotes to a constant, time step-dependent value with respect to the sampling time. Furthermore, for a given computational cost (e.g., the number of timsteps), an optimal choice of discretization (i.e., time step) exists that minimizes the expected error in simulation, when accounting for both the effects of discretization error and sampling error. We show that this optimal choice results in a convergence rate with respect to computational cost that is bounded by the sampling convergence rate with a minor impact from the discretization order of accuracy. Finally, we consider the implications of spin-up time (i.e., unsampled time needed to arrive at the stationary distribution) and parallelism on the optimal error. We develop a method for estimating transient-related errors and then evaluate optimal choices incorporating the results.

To approximate J, we compute finite-time, discrete estimates of the outputs of interest of the true system as

JT,hp=1TsIt0t0+Ts(ghp(uhp(t))),
(3)

where the notation Iab(·) here represents the quadrature approximation of the integral ab(·)dt of a quantity (·) between a and b. Here, we have made a discrete approximation of the state using an order-p discretization with a temporal grid with characteristic size h=Δt, where an order-p discretization is one for which the discretization error behaves as

maxt[0,t0+Ts]|ghp(uhp(t))g(u(t))|=O(hp),
(4)

when the discretization is applied to a well-posed (non-chaotic) system. Then we sample a discrete state over a finite sampling period, Ts, starting at some initial time t0. We can define the error that is incurred as

eT,hp=JT,hpJ.
(5)

By introducing a third value,

JT=1Tst0t0+Tsg(u(t))dt,
(6)

we can re-write the error using an identity

eT,hp=(JT,hpJT)+(JTJ)=ehp+eT.
(7)

Here, we define the “discretization error” and “sampling error,” respectively, as

ehpJT,hpJT,
(8)
eTJTJ.
(9)

We can take an absolute value of both sides of (7), followed by a manipulation using the triangle inequality

|eT,hp|=|ehp+eT||ehp|+|eT|.
(10)

Thus, the total error incurred by the approximation is bounded by the sum of the absolute discretization and sampling errors. Next, we define the attractor of the operator f, A, as the set of long-term states toward which all trajectories converge independently of the initial condition.22 We can define the expectation EA[ϕ(u0)] for a generic function ϕ as the expectation taken over all the trajectories that can result from starting from points on the attractor, A,

EA[ϕ]=1|A|u0Aϕ(u0)du0.
(11)

For the case in question, we will be considering either

ϕ(u0)=1Ts|t0t0+Tsg(u(t))dtJ|

or

ϕ(u0)=1Ts|It0t0+Ts(ghp(uhp(t)))t0t0+Tsg(u(t))dt|,

with, for these examples, u(t0)=u0A. Given these definitions, we can now take the expectation of (10), giving

EA[|eT,hp|]EA[|ehp|]+EA[|eT|],
(12)

by linearity.

From here, we propose asymptotic forms for the two right-hand side terms in (12). Consider the definition of eT in (9)

eT=1Tst0t0+Tsg(u(t))dtJ.
(13)

Assuming that we choose t0 such that each u0 is effectively an independent sample from the attractor's stationary distribution, then the quantity g(u(t)) is a random variable drawn from a stationary distribution. The states of ergodic systems, in general, are not independent in time, but as long as the system has satisfactorily strong mixing properties, the central limit theorem (CLT) can be applied to finite time averages of its outputs. This is the case whenever the condition of α-mixing is met,23,24 which has been shown for the Lorenz system.25 Thus, we can write eT as

eTN(0,(π2A0Ts1/2)2),
(14)

with a constant coefficient A0, where N(μ,σ2) gives the normal distribution with mean μ and variance σ2. If we take the absolute value of this random variable, the result is a halfnormal distribution

|eT|H((π2A0Ts1/2)2),
(15)

where H(σ2) gives a halfnormal distribution such that |X|H(σ2) when XN(0,σ2). The expectation of the half-normal distribution is well defined, allowing

EA[|eT|]A0Ts1/2,
(16)

as Ts goes to infinity.

Now consider the use of a time-stepping method to give a discrete approximation uhp(tn) of u(tn) for each tn=n(Δt). Following the classical analysis,26 we might expect that the discretization error should take a form

|ehp|Cp(exp(ΛTs)1Λ)(Δt)p.
(17)

This analysis is based on bounding the growth of local truncation error at each time step by the Lipschitz constant, Λ, of the underlying system, with Cp being a constant parameter that depends on the choice of method. However, Viswanath showed27 that, the global error could be modeled by a form

|ehp|E(Ts;p)(Δt)p,
(18)

where E(Ts;p) could be bounded by a constant for some nonlinear but non-chaotic systems that are exponentially stable. While this result has not been extended to an ergodic system, the expected convergence onto the attracting set suggests a bound of the form

EA[|ehp|]Cp(Δt)p.
(19)

As our results in Sec. III will show, (19) is a good description of the expected discretization error.

Thus, taking (12), (16), and (19), we assume a bound of the form

EA[|eT,hp|]emodel=Cq(Δt)q+A0Tsr=Cq(Δt)q+A0Nsr(Δt)r,
(20)

that bounds EA[|eT,hp|] when Δt is small enough and Ts is large enough to satisfy the asymptotic assumptions. Here, we have used Ts=NsΔt, where Ns is the number of timesteps computed for sampling; additionally, q refers to the observed discretization convergence behavior, which, in practice, may differ from p due to numerical cancelations or if the solutions of the system are insufficiently regular. Similarly, r is an observed sampling convergence rate coefficient, which we expect to be 1/2 asymptotically under the CLT.

In this section (Sec. III), we will fit numerical results for the Lorenz system to determine q, r, Cq, and A0 and show that this model is representative of the observed behavior. The Lorenz system is given by:28 

dudt=f(u;α)=(α0(u1u0)u0(α1u2)u1u0u1α2u2),
(21)

where u=[u0,u1,u2] and α=[α0,α1,α2]. The Lorenz system is known to be chaotic for the classic Lorenz parametrization:29,α=[10,28,8/3], which is used everywhere in this text. For the output, we choose g(u)=u2. We consider a set of explicit methods: forward Euler (FE, p = 1), third-order Runge–Kutta (RK3, p = 3), and fourth-order Runge–Kutta (RK4, p = 4). In all of these methods, we expect asymptotic convergence of JT,hp to JT to be at least O(Δtp) for non-chaotic systems.30 The stability limit for each of these methods is Δt101, which is above the limits of asymptotic convergence for all of the methods here; tests with implicit methods indicated similar error trends at higher cost and are, therefore, omitted.

For any given discrete instance, we will start the simulation at an initial state at t = 0 that is sampled randomly from a normal distribution

uinit(N(1.0,5.02)N(1.0,5.02)N(1.0,5.02)).
(22)

To guarantee that the initial sampling state u0 at t0 is on the attractor (as well as further guaranteeing the independence from the other Monte Carlo instances), we evolve the state of any given Lorenz system discretization from its starting state uinit for t0=100 before proceeding to the sample; we refer to the process of evolving the solution until it is on the attractor as “spin-up.” Then, we evolve the state over the next Ts, during which we integrate and compute (3) using the same numerical integration scheme that was used for the state itself.

To approximate eT,hp, we must first estimate J by a reference value Jref. Jref is calculated using an ensemble mean of JT,hp over Mens=5122 instances of the Lorenz system. Each instance is started from a different uinit as given in (22) and simulated using RK4 with Δt=17.7×106 and Ts=6646.9. The resulting Jref is

Jref=23.549916±0.000074,
(23)

with a 95% confidence estimate based on the ensemble mean estimator.

The computation of Jref allows us to estimate errors eT,hpJT,hpJref. For a given Δt, Ts pair, we then approximate EA[|eT,hp|] using a Monte Carlo method over M=10000 independent instances of the discrete system, each started from initial states drawn from (22) and spun-up to independent sampling starting points on the attractor u0(m),

E[|eT,hp|]1Mm=1M|JT,hp(u0(m))Jref|.
(24)

In Figs. 1–3, we compare the results of simulations with the FE, RK3, and RK4 discretizations with different values of Ns. In these figures, Ts scales with Δt for a given Ns, so the Ts values on the x-axis will vary between lines on the plot. The fits shown are computed with truncated data in order to attempt to eliminate non-convergent data at small Ts or large Δt; the limits used for truncation are found in Table I. The results of the nonlinear least squares fits for Ns=104, 105, and 106 are given in Table II. In the table, we observe that r1/2 as the discretization error is reduced, either by increasing Ns or by pushing p higher. These figures demonstrate that (19) has explanatory value, as the errors in the discretization-dominated region collapse independently of Ts. It is also worth noting that Table II demonstrates higher-than-expected discretization error convergence rates for FE and RK4. In Fig. 4, we can examine the sampling error behavior between discretization methods for a single shared choice of Ns. Here, we can see that the sampling error effects on the left-hand side of the plot collapse independently of the discretization method. This indicates that the statistical effects are properties of the dynamical system, not artifacts of the discretization, as we might expect in the limit as Δt0. We note that Table II shows some Ns-dependence in the fit values of r; we expect that this is an artifact of the fit process, in which discretization effects and statistical effects are intermingled. As Ns or q become small, the region in which the convergent CLT behavior is the dominant effect becomes smaller, and for these cases, the value of r detected is subject to error and deviation from the expected rate r = 1/2.

FIG. 1.

Expected relative error as a function of Δt for forward Euler discretization of the Lorenz equations. Nonlinear least squares fit based on Ns=106 data.

FIG. 1.

Expected relative error as a function of Δt for forward Euler discretization of the Lorenz equations. Nonlinear least squares fit based on Ns=106 data.

Close modal
FIG. 2.

Expected relative error as a function of Δt for RK3 discretization of the Lorenz equations. Nonlinear least squares fit based on Ns=106 data.

FIG. 2.

Expected relative error as a function of Δt for RK3 discretization of the Lorenz equations. Nonlinear least squares fit based on Ns=106 data.

Close modal
FIG. 3.

Expected relative error as a function of Δt for RK4 discretization of the Lorenz equations. Nonlinear least squares fit based on Ns=106 data.

FIG. 3.

Expected relative error as a function of Δt for RK4 discretization of the Lorenz equations. Nonlinear least squares fit based on Ns=106 data.

Close modal
TABLE I.

Fit boundaries for nonlinear least squares fits.

MethodΔtmaxTs,min
FE 5.0×103 1.0 
RK3 5.0×102 1.0 
RK4 9.0×102 1.0 
MethodΔtmaxTs,min
FE 5.0×103 1.0 
RK3 5.0×102 1.0 
RK4 9.0×102 1.0 
TABLE II.

Values of error model coefficients computed from nonlinear least squares fits to Monte Carlo study data.

FERK3RK4
(a) Ns=104 
A0 2.19 1.74 1.63 
r 0.975 0.721 0.683 
Cq 4995 942 85 900 
q 1.65 2.70 4.83 
(b) Ns=105 
A0 1.94 1.50 1.41 
r 0.820 0.648 0.620 
Cq 1410 1310 96 100 
q 1.40 2.76 4.84 
(c) Ns=106 
A0 1.52 0.978 0.918 
r 0.693 0.553 0.538 
Cq 714.6 2740 165 000 
q 1.273 2.96 5.02 
FERK3RK4
(a) Ns=104 
A0 2.19 1.74 1.63 
r 0.975 0.721 0.683 
Cq 4995 942 85 900 
q 1.65 2.70 4.83 
(b) Ns=105 
A0 1.94 1.50 1.41 
r 0.820 0.648 0.620 
Cq 1410 1310 96 100 
q 1.40 2.76 4.84 
(c) Ns=106 
A0 1.52 0.978 0.918 
r 0.693 0.553 0.538 
Cq 714.6 2740 165 000 
q 1.273 2.96 5.02 
FIG. 4.

Expected relative error as a function of Δt for discretizations of the Lorenz equations.

FIG. 4.

Expected relative error as a function of Δt for discretizations of the Lorenz equations.

Close modal

Finally, we attempt to compare the computational costs across the various discretizations. In this case, the number of timesteps Ns is not a good proxy for fixed cost, since the computation time for a time step will vary between methods. Instead, we now fix Us, the total number of evaluations of the right-hand side f used in sampling timesteps. For the explicit schemes used in this work, we will have p right-hand side evaluations (e.g., forward Euler has p = 1 right-hand side evaluations) and, thus, Us=pNs. In Fig. 5, we can see the effect of changing Δt at fixed sampling cost Us across discretizations. The error that can be achieved with the Runge–Kutta methods is lower than that of the forward Euler scheme, a factor of 4.8 improvement in the error from FE to RK4. However, the best-case improvement for going from third-order to fourth-order Runge–Kutta schemes is the only factor of about 1.4. Moreover, the results show that to achieve the lowest possible error, the optimal time step will be discretization dependent. We investigate this further in Sec. IV.

FIG. 5.

Expected percent error as a function of Δt for discretizations of the Lorenz equations at a number of sampling residual evaluations. All fits evaluated at Us=1.2×106.

FIG. 5.

Expected percent error as a function of Δt for discretizations of the Lorenz equations at a number of sampling residual evaluations. All fits evaluated at Us=1.2×106.

Close modal

We now study the implications of the error model (20), specifically seeking to understand the convergence of the error with respect to computational effort. In this analysis, we will assume that r = 1/2.

Consider a non-dimensional form of the error model in which the error is normalized by the standard deviation of the instantaneous output σg, and the timescales Δt and Ts are normalized by decorrelation time Td. The decorrelation time relates the amount of variance from independent draws from the distribution on the attractor and the amount of variance in the finite-time mean estimators based on the correlated output signal, given by the relation:31 

Var[JT]=TdTsσg2.
(25)

Furthermore, combining (15) and (25) allows us to write

A0=2πσgTd1/2.
(26)

In general, Td is hard to estimate accurately; this is a crux of the work of Oliver et al.20 In our formulation of the error model, we identify A0, which avoids outright estimation of Td. However, for the purpose of understanding the behavior of the error, Td is an intrinsic timescale, which can be used to normalize Δt and Ts.

The resulting non-dimensional form of the error model is

emodelσg=CqTdqσg(ΔtTd)q+2π(TsTd)12.
(27)

We can also write the optimizers and optimal value of (27) in terms of the non-dimensional variables. These are given by

(ΔtTd)opt=(12π)12q+1(qCqTdqσg)22q+1Ns12q+1,(TsTd)opt=(12π)12q+1(qCqTdqσg)22q+1Ns2q2q+1,(emodelσg)opt=(12π)q2q+1(2+1q)(qCqTdqσg)12q+1Nsq2q+1.
(28)

In terms of convergence with respect to sampling costs, the error model will scale at best as

(emodelσg)optNsq2q+1.

In the limit as q, the rate q/(2q+1)1/2: the CLT limits the convergence rate. Table III gives the rates of convergence (28) for various values of q.

TABLE III.

Convergence rates for combined error with respect to sampling timesteps implied in (28) at common high-order discretization error convergence rates.

q12345
q2q+1 1/3 2/5 3/7 4/9 5/11  1/2 
q12345
q2q+1 1/3 2/5 3/7 4/9 5/11  1/2 

Using the reference simulation, we can also find

Var[JT]Var[JT,hp]=1.1692×104,σg2σ̂g2=74.34804±0.00018,
(29)

where σ̂g is an estimate of the standard deviation of g. Together, these allow us to estimate

Td1.0170×102,σg8.6225.
(30)

With these values, we can plot the non-dimensional error model with fixed r = 1/2, which is given for Ns=105 in Fig. 6.

FIG. 6.

Expected non-dimensional error as a function of the non-dimensional time step for discretizations of the Lorenz equations. r = 1/2 assumed.

FIG. 6.

Expected non-dimensional error as a function of the non-dimensional time step for discretizations of the Lorenz equations. r = 1/2 assumed.

Close modal

We now consider the implications of these results for increasing Ns. To focus solely on control of the discretization error, increases in Ns can be used to refine Δt=Ts/Ns, with Ts fixed. On the other hand, to focus solely on controlling the sampling error, Ts=NsΔt can be increased, holding Δt fixed. In Fig. 7, the two approaches are compared with the optimal use of resources. In orange is the discretization error control strategy. In this approach, the simulations converge at a high-order rate in Ns toward the optimal error behavior; once the error reaches this optimum, however, it asymptotes to a constant: statistical errors limit the estimation of JT,hp. On the other hand, the sampling error control approach is shown in blue. In this approach, the central limit convergence rate of 1/2 is initially achieved until the error asymptotes to a constant: discretization errors limit the estimation of JT,hp. In other words, this latter case represents having increasingly many significant figures of certainty about an answer that is significantly offset from the true quantity of interest, J. In the literature for large simulations, discussed in the introduction, simulations tend to be planned using either the discretization or statistical error control approach. What (20) implies and Fig. 7 demonstrates is that, in fact, there is a particular optimal scheme in which Δt and Ts are simultaneously varied that will extract the most accurate estimate of J as Ns increases.

FIG. 7.

Refinement study comparison for fixed Δt, fixed Ts, and optimized Δt and Ts using RK3 discretization to compute the expectation of the Lorenz system output g = u2.

FIG. 7.

Refinement study comparison for fixed Δt, fixed Ts, and optimized Δt and Ts using RK3 discretization to compute the expectation of the Lorenz system output g = u2.

Close modal

In this section, we show that our simulations of chaotic, ergodic ODEs are consistent with a bounded relationship between the local and global discretization errors. Consider an estimate of the global error based on Ns timesteps

ehp1Nsn=0Nsη=nNsG(tη,tn)°eLT,p(n),
(31)

where

eLT,p(n)uhp(tn+1)u(tn+1),
(32)

and u(tn+1) is the exact solution integrated from uhp(tn) through Δt

u(tn+1)=uhp(tn)+tntn+1f(u(t))dt.
(33)

In (31), we have assumed that the error from any given local state perturbation is propagated forward in time by the dynamics, before being transformed into an error in the output; this process is captured by an operator G. Because the effect of local error propagates forward and not backward in time, G(t,tn)=0 for t<tn, and moreover, we assume that due to ergodicity G(t,tn)=0 when ttnTd, where Td is the decorrelation time associated with the attractor. This allows us to write

ehp1Nsn=0Nsη=nn+Td/ΔtG(tη,tn)°eLT,p(n).
(34)

Now, we assume that a constant Gmax exists such that

|G(tη,tn)°v|Gmax||v||,
(35)

for all tn,tη, and vB(u(tn))d, where B(u) is the set of states possible by perturbation of u that remain in the basin of attraction of the attractor A of f. When this is the case, we can create a bound on the magnitude of ehp

|ehp|TdΔtGmax1Nsn=0Ns||eLT,p(n)||TdΔtGmaxmaxn||eLT,p(n)||.
(36)

We now attempt to bound the value of Gmax for the Lorenz system by approximating the local truncation error. To make an estimate, we compute both the solution at the next time step as well as a surrogate for the true solution at each time step: uhp(tn+1) and ũ(tn+1), where the former is computed with one time step of the method of interest and the latter is always computed with the highest available accuracy method, RK4, and subdividing t[tn,tn+1] into ten consecutive timesteps rather than one. Both uhp(tn+1) and ũ(tn+1) are always advanced from uhp(tn). This allows us to estimate eLT,p(n+1) locally

eLT,p(n+1)ẽLT,p(n+1)=uhp(tn+1)ũ(tn+1).
(37)

In Fig. 8, we characterize the convergence of local error estimates. Computations are run with Ts = 100 and t0=100 fixed, varying Δt. At each time step, the local truncation error is estimated by computing (37). The figure shows the computed maxn||ẽLT,p(n)|| and demonstrates that the expected rate of (p+1) is nearly exactly achieved.

FIG. 8.

Convergence of estimated local truncation error with respect to Δt. Fits to cpΔtp+1 shown (with offset for presentation).

FIG. 8.

Convergence of estimated local truncation error with respect to Δt. Fits to cpΔtp+1 shown (with offset for presentation).

Close modal

Using (36), we can estimate a bounding value for Gmax by

GmaxE[|ehp|]maxn||eLT,p(n)||ΔtTd=CqΔtqcpΔtp+1ΔtTd,
(38)

where cp is the leading truncation error coefficient fit in Fig. 8 and Cq and q are taken from Table II. Of course when q>qtheory=p, there will be Δt dependence. However, as (38) requires that the discretization error has an asymptotic behavior, we will only consider Δt in the asymptotic convergence regions given in Table I to compute Gmax. In Fig. 9, we show the values of the right-hand side quantity in (38), which allow us to make an estimate

Gmax4.3.
(39)
FIG. 9.

Estimation of bounding value Gmax.

FIG. 9.

Estimation of bounding value Gmax.

Close modal

Next, we use classical truncation error estimates26 to relate the discretization error to properties of the solution. We will assume that the local truncation error is bounded by a form

maxn||eLT,p(n)||CLT(p+1)!dp+1udtp+1Δtp+1,
(40)

where CLT is a local truncation constant term dependent on the numerical method and the ||·|| in this context refers to the maximum value in time of the inf-norm of a vector-valued, time-dependent quantity (·). The derivatives of u(t) can be computed by evaluating f(u) and its derivatives using solutions from a reference RK4 solution of the Lorenz system with Ts = 1000, t0=100, and Δt=104. Norms of the derivatives are shown in Fig. 10. The resulting values of CLT that can now be derived by fitting the asymptotic behavior in Fig. 8 can be found in Table IV. The result of these estimates is that we can reliably bound the global error of a dynamical system as an accumulation of the local errors over a region of correlation.

FIG. 10.

Norm of analytic derivatives of u computed on the attractor of f. State u computed with RK4 at Δt=104 and Ts = 1000 after discarding t0=100.

FIG. 10.

Norm of analytic derivatives of u computed on the attractor of f. State u computed with RK4 at Δt=104 and Ts = 1000 after discarding t0=100.

Close modal
TABLE IV.

Rate and coefficient fit for convergence of local truncation error of the discrete Lorenz system. CLTdp+1udtp+1 estimated by cp(p+1)! using cp fit from Fig. 8.

prate(observed)CLTdp+1udtp+1CLT
2.00 7.61×103 7.33 
4.02 3.28×106 156 
4.93 4.50×107 76.5 
prate(observed)CLTdp+1udtp+1CLT
2.00 7.61×103 7.33 
4.02 3.28×106 156 
4.93 4.50×107 76.5 

We now want to consider how the global error behavior demonstrated here might extrapolate to more complicated systems by evaluating the spectral behavior of the Lorenz system. Using a discrete Fourier transform (DFT) with a Hann window function,32 we perform a spectral analysis on the states of the Lorenz system with sampling time Ts = 1000, t0=100, and Δt=103. The resulting spectrum can be found in Fig. 11. The Lorenz system tends to have the most content in the frequencies with f101 with a region of exponential decay in the range 1f300. On scales with f300, machine precision plateaus are observed and omitted here.

FIG. 11.

Fourier spectrum of u(t). Computed with DFT using the Hann window function on data from RK4 discretization of the Lorenz system with Ts = 1000, t0=100, and Δt=103. Gray dashed line: fit assuming |û(f)|exp(af+b) with a = 0.872 and b = 2.58.

FIG. 11.

Fourier spectrum of u(t). Computed with DFT using the Hann window function on data from RK4 discretization of the Lorenz system with Ts = 1000, t0=100, and Δt=103. Gray dashed line: fit assuming |û(f)|exp(af+b) with a = 0.872 and b = 2.58.

Close modal

The fact that the Lorenz spectrum is an exponentially decreasing function of frequency f makes the use of high-order methods theoretically appealing for the spectral convergence of hp-refinement strategies.33 Unfortunately, the effect of statistical error in (28) limits the impact of this exponential decay, such that the benefits of higher-order discretization methods are limited compared to their steady-state and non-chaotic application. The convergence to the central limit rates can be seen in Fig. 12, which shows the convergence of (28) with the total sampling cost. The effect of increasing order improves the convergence rate in (28) toward the CLT-implied asymptotic rate of 1/2, as well as decreasing the value of the leading constant and the error never achieves the spectral rates possible with hp-refinement in the steady case. Nevertheless, the cost to achieve a given amount of error in expectation—in terms of function evaluations—is significantly less with higher-order methods. Managing to achieve 1% non-dimensional error in expectation is possible with RK4 at a cost ten times less than would be possible using FE; this factor grows larger than 100 when the tolerance is tightened to 104.

FIG. 12.

Convergence of optimal error with sampling costs for FE, RK3, and RK4 discretizations of the Lorenz output g = u2. Asymptotic 1/2 rate implied by the central limit theorem is shown.

FIG. 12.

Convergence of optimal error with sampling costs for FE, RK3, and RK4 discretizations of the Lorenz output g = u2. Asymptotic 1/2 rate implied by the central limit theorem is shown.

Close modal

In this section, we will consider how the error behaves when ensemble averaging (over multiple parallel instances) and when spin-up effects are present.

Sampling error can be reduced at a fixed wall clock time by ensemble averaging across multiple parallel processes.34 Consider a Monte Carlo approach to approximate J with a set of Mens independent realizations

JMC=1Mensm=1MensJT,hp(m).
(41)

We can write a modified version of (20) to approximate the error that we expect in the Monte Carlo estimator in (41)

E[|JMCJ|]emodel,MC=Cq(Δt)MCq+A0MensTs,MCr,
(42)

with an equivalent non-dimensional version, assuming r1/2

(emodelσg)MC=CqTdqσg(ΔtTd)q+2πMens12(TsTd)12,
(43)

and an optimum given by

(emodelσg)MC,opt=(12π)q2q+1(2+1q)(qCqTdqσg)12q+1Mensq2q+1Nsq2q+1,
(44)

at

(ΔtTd)opt=(12π)12q+1(qCqTdqσg)22q+1Mens12q+1Ns12q+1
(45)

and

(TsTd)opt=(12π)12q+1(qCqTdqσg)22q+1Mens12q+1Ns2q2q+1.
(46)

Equation (44) shows that, for finite values of q, the Monte Carlo method will have a mitigated return compared to its purely stochastic application as in Ref. 34; the optimal error scales as Mensq/(2q+1) as opposed to Mens1/2. However, parallelization can achieve perfect scaling in the expected error, in the sense that the effect of running Mens ensembles with Ns sampling timesteps each will have an equivalent error in expectation to simulate MensNs timesteps in series. As Mens is varied on the set of optimal solutions, (45) and (46) indicate that the time step and sampling time should be adjusted with the same factor Mens1/(2q+1) to achieve perfect scaling.

So far, we have considered the error and cost on the attractor, neglecting the impact of spin-up from t = 0 to t = t0. This spin-up is necessary because simulations of ergodic systems invariably need some time for the state to proceed onto the attractor from the initial condition.

Consider u(t), a solution of the ergodic chaotic system f from an arbitrary initial condition u(0)=uIC in the basin of attraction of an attractor, A. The existence of the attractor implies the nonlinear stability of the system, such that all uIC will converge to trajectories on the attractor A. Denote by uA(t) a trajectory that is on the attractor for all t and to which u(t) collapses as t. The perturbation δuA(t)u(t)uA(t) that describes the IC, therefore, exists in a stable subspace of perturbations to uA and can be associated with the negative Lyapunov exponents of the system. Thus, we can assume that such perturbations are governed asymptotically by

||δuA(t)||exp(tTλ),
(47)

with Tλ being a characteristic time associated with the stable Lyapunov modes. In practice, we are interested in averages of quantities on the attractor g(uA(t)), but we can only calculate quantities g(u(t)) that will include some effect—if small—of the spin-up transient.

Next, we seek to quantify the effect of this gap on estimates JTJ. Consider the computation of JT. In (7), we have effectively found an estimate of

JTA=1Tst0t0+Tsg(uA(t))dt,
(48)

by choosing t0 sufficiently large. We now want to consider an error model of the form

eT,hp=(JT,hpJT)ehp+(JTJTA)eλ+(JTAJ)eT,
(49)

where a new error eλ is introduced, associated with the spin-up transient. The model for eT in (9) will apply without modification, while the model for ehp will be subject to slightly different assumptions, where in (8), Cq was bounded by the value on the attractor, A, here we must assume that Cq is bounded from t = 0 to t=t0+Ts, including both the attractor and the transient part of the trajectory. We only require that the transient part be in the basin of attraction of A,B(A). We assume that a model of the form used in (8) applies in expectation when the transient component is included.

Next, we concentrate on eλ,

JTJTA=t0t0+Ts(g(u(t))g(uA(t)))dt.
(50)

We now assume that, like u, g will decay exponentially in t as (47), such that

g(u(t))g(uA(t))δgA(t)Aλexp(tTλ)
(51)

will apply for t[0,), with Aλ being a constant that can be related to the deviation between g(u(0)) and g(uA(0)).

From this assumption,

eλ=1Tst0t0+Tsg(u(t))g(uA(t))dt1Tst0t0+TsAλexp(tTλ)dt=AλTλTsexp(t0Tλ)(1exp(TsTλ)).
(52)

Taking the absolute value, we can find a bounding model

|eλ|=|Aλ|TλTsexp(t0Tλ).
(53)

As before, manipulation of (49) allows

|eT,hp|=|ehp+eλ+eT|
(54)
|ehp|+|eλ|+|eT|.
(55)

Now, we take an expectation of the absolute value of eT,hp,

E[|eT,hp|]EB(A)[|ehp|]+EIC[|eλ|]+EA[|eT|],
(56)

where EB(A) gives the expectation on the basin of attraction of A. Here, the expectation of |eT,hp| does not reduce to an expectation on the attractor. The statistical term is handled on the attractor as before, and we have assumed that the discretization error is bounded by the same form in expectation on B(A) as on A. Finally, the expectation of |eλ| is taken on the set of initial conditions used. This allows us to take the expectation of (53) to complete (56). Because we anticipate Tλ will be bounded by a constant for a given system, this is given by

EIC[|eλ|]=EIC[|Aλ|]TλTsexp(t0Tλ).
(57)

If a Aλ and Tλ can be identified by observation of g(u(t)) given an initial condition uIC,|eλ| is no longer stochastic and the E[|eT,hp|]|eT,hp| as in (53).

Putting all the pieces together, we can now give an error model that incorporates the effects of spin-up and ensemble estimation

emodel,MC=ÃλTλTs,MCexp(t0Tλ)+Cq(Δt)MCq+A0MensTs,MCr,
(58)

where Ãλ can be either estimated on an instance-by-instance basis or by estimating the expectation on the family of initial conditions. Under this model, eλ will scale with the exponent of a large negative value when t0Tλ. Even when t0̸Tλ, (53) suggests that the decay-induced error term will still scale with Ts1, faster than the expected CLT rate of Ts1/2, and thus, it will be dominated it as Ts1. This implies two “paths” to control spin-up errors: either choosing t0 long enough to shrink the mean offset error from t0 or choosing Ts long enough so that the mean offset contribution to the simulation error is small in spite of the error at t = t0.

We will now develop a method to fit the error model. In order to do so, consider observations gng(tn) and gnAgA(tn) for tn in {t0,t0+NskipΔt,,t0+Ts}. We will assume that Nskip is large enough that the solution at each tn is effectively independent. If this is the case, then we can assume that each gnA will be an independent and identically distributed (i.i.d.) draw from a bounded, stationary distribution with mean J. The distributions of gA(u(t)) and g(u(t)), in general, are not known. In order to facilitate an estimate of the mean behavior, we will assume gnA are i.i.d. drawn from a normal distribution with mean value J. Then, we have

gnN(J+δgA(tn),σg2),
(59)

where the relationship between gn and gnA is taken from (51).

In order to understand the implications of this model, we can use a set of reference RK4 simulations of the Lorenz system with Nt=105 timesteps sampled without spin-up over a period T = 100 from initial conditions similar to those given in (22) with a scaled-up standard deviation of 100 in all three variables to highlight the initial transient. In order to treat each of J, σg, Aλ, and Tλ in (59) as unknowns, we use Hamiltonian Monte Carlo with the likelihood function implied in (59). We discard from t = 0 to t = 5 and then take 10 000 equispaced samples from t = 5 to t = 100. For prior models, we start by computing naïve estimators of the mean and standard deviation of the trace, J̃ and σ̃ using the downsampled trace signal {gn} and then use

JN(J̃,σ̃2),σgΓ(ασ,βσ),AλN(0,max(ghp)min(ghp)),TλΓ(αT,βT),
(60)

where

(ασ,βσ)(μσ=σ̃,σσ=σ̃10),(αT,βT)(μT=10.0,σT=10.0).

It should be noted that in this specification, the Bayesian fit only requires a user-supplied prior for the decay time and for the uncertainty in the standard deviation, assumptions upon which the fitting method only requires to be reasonable.

A sample fit and trace are found in Fig. 13, for which the maximum a posteriori (MAP) estimate gives Tλ=0.312 and Aλ=0.925. For the Lorenz system, the initial transient onto the attractor is very rapid, almost negligible. Applying the Bayesian fit procedure to an ensemble of 1000 runs generated in the same way as Fig. 13, we can find maximum a posteriori (MAP) estimates of the variables Tλ and |Aλ| in the decay model. In Figs. 14 and 15, histograms of these variables are shown, which are needed to determine (53). We can see that the fit procedure identifies values

Tλ<4.03,|Aλ|<38.7,
(61)

for greater than 97% of initial conditions, up to two standard deviations above the mean. Using these values as a conservative estimate for the mean offset, we can now model the effect of the transient behavior.

FIG. 13.

g=u2(t) trace in the transient region with Bayesian method fit.

FIG. 13.

g=u2(t) trace in the transient region with Bayesian method fit.

Close modal
FIG. 14.

MAP estimate Tλ for Lorenz system transient. Collected over 1000 Lorenz trajectories with Δt=102, Ts = 100, and randomized uIC. Outliers truncated, greater than 97% of data in the pictured range.

FIG. 14.

MAP estimate Tλ for Lorenz system transient. Collected over 1000 Lorenz trajectories with Δt=102, Ts = 100, and randomized uIC. Outliers truncated, greater than 97% of data in the pictured range.

Close modal
FIG. 15.

MAP estimate |Aλ| for Lorenz system transient. Collected over 1000 Lorenz trajectories with Δt=102, Ts = 100, and randomized uIC. Outliers truncated, greater than 97% of data in the pictured range.

FIG. 15.

MAP estimate |Aλ| for Lorenz system transient. Collected over 1000 Lorenz trajectories with Δt=102, Ts = 100, and randomized uIC. Outliers truncated, greater than 97% of data in the pictured range.

Close modal

Now we can consider how the cost and error impact of spin-up is incorporated into the model for error at a fixed cost. The spin-up time requires the use of N0 timesteps

N0=t0(Δt)MCt0(Δt)MC,
(62)

with N being the total number of timesteps used, given by

N=N0+NMC=t0(Δt)MC+NMC,
(63)

where NMC is the number of timesteps during sampling for t0 to t0+Ts on a given instance.

By normalizing (58) then substituting (63), we arrive at a transient-inclusive non-dimensional model for the error

(emodelσg)MC=ÃλσgTλTd(N(ΔtTd)MCt0Td)1exp(t0/TdTλ/Td)+CqTdqσg(ΔtTd)MCq+2πMens12(N(ΔtTd)MCt0Td)12.
(64)

Using this result, we can solve numerically for (Δt)MC,opt and eMC,opt via (64).

Consider a Lorenz simulation on which a budget of U=pN=1.2×106 right-hand side evaluations are available on each of Mens parallel processors. We start by studying the error under (64) as Δt and t0 vary with a conservative estimate for the transient behavior using the bounding values in (61). In Fig. 16, we show emodel for forward Euler at a fixed cost of U=1.2×106. (The optimum is denoted by a red star.) Moving to the right, discretization error becomes the dominant factor as ΔtTd. The diagonal boundary gives the region of feasibility at which, under the cost constraint, sampling no longer occurs (Ts = 0). Moving from the optimum toward the bottom left, t00,Ts0, and ΔtTd; thus the transient error and sampling error become dominant. Similar plots for RK3 and RK4 are found in Figs. 17 and 18. The optimal errors and optimizing simulations are described in Table V. We can see from these results that, at a fixed budget with U=1.2×106, the effect of increasing the discretization order makes a smaller error possible with a larger time step, which means fewer timesteps to traverse the spin-up time. These two effects combine to allow for an increase in the sampling time available Ts, allowing significantly less sampling error for RK3 compared to FE, and an additional—albeit smaller—benefit moving from RK3 to RK4, holding cost fixed.

FIG. 16.

Dependence of normalized error expectation emodel,MC/σg on normalized time step Δt/Td and normalized spin-up time t0/Td with total cost set at U=1.2×106 for forward Euler. Red star denotes optimum, and dashed line indicates optimal t0 given Δt.

FIG. 16.

Dependence of normalized error expectation emodel,MC/σg on normalized time step Δt/Td and normalized spin-up time t0/Td with total cost set at U=1.2×106 for forward Euler. Red star denotes optimum, and dashed line indicates optimal t0 given Δt.

Close modal
FIG. 17.

Dependence of normalized error expectation emodel,MC/σg on normalized time step Δt/Td and normalized spin-up time t0/Td with total cost set at U=1.2×106 for third-order Runge–Kutta. Red star denotes optimum, and dashed line indicates optimal t0 given Δt.

FIG. 17.

Dependence of normalized error expectation emodel,MC/σg on normalized time step Δt/Td and normalized spin-up time t0/Td with total cost set at U=1.2×106 for third-order Runge–Kutta. Red star denotes optimum, and dashed line indicates optimal t0 given Δt.

Close modal
FIG. 18.

Dependence of normalized error expectation emodel,MC/σg on normalized time step Δt/Td and normalized spin-up time t0/Td with total cost set at U=1.2×106 for fourth-order Runge–Kutta. Red star denotes optimum, and dashed line indicates optimal t0 given Δt.

FIG. 18.

Dependence of normalized error expectation emodel,MC/σg on normalized time step Δt/Td and normalized spin-up time t0/Td with total cost set at U=1.2×106 for fourth-order Runge–Kutta. Red star denotes optimum, and dashed line indicates optimal t0 given Δt.

Close modal
TABLE V.

Optimal Lorenz simulations for output g = u2 under budget of U=1.2×106 right-hand side evaluations using Mens=1.

MethodpemodelΔtt0Ts
FE 0.0502 2.54×104 30.2 275 
RK3 0.0130 8.52×103 35.5 3370 
RK4 8.89×103 0.0224 36.7 6670 
MethodpemodelΔtt0Ts
FE 0.0502 2.54×104 30.2 275 
RK3 0.0130 8.52×103 35.5 3370 
RK4 8.89×103 0.0224 36.7 6670 

In Fig. 19, we take another perspective on these results for RK3 by varying Δt and plotting the optimal t0, Ts, and emodel. As Δt gets large, the optimal choice of t0 has logarithmic growth, and when Δt/Td1, the optimal choice of t0 rapidly falls to zero. Parallelization has a small but non-zero effect on the optimal choice of sample time. The sampling time also has a small effect from parallelization, in this case constrained to a small region. Outside that Δt region, Ts scales with Δt both as Δt0 and as Δt.

FIG. 19.

Dependence of normalized spin-up time t0/Td, sampling time Ts/Td, and model error emodel/σg on normalized time step Δt/Td with total cost set at U=1.2×106 for third-order Runge–Kutta.

FIG. 19.

Dependence of normalized spin-up time t0/Td, sampling time Ts/Td, and model error emodel/σg on normalized time step Δt/Td with total cost set at U=1.2×106 for third-order Runge–Kutta.

Close modal

The bottom plot of Fig. 19 shows the variation of error with Δt. In this plot, we can see three distinct regions. For Δt/Td102, discretization error is the dominating error, and the convergence goes with the discretization error rate. Approaching the optimum, sampling error becomes the dominant error contribution, starting at Δt2×102 until Δt103. In this region, the convergence is around the CLT-implied 1/2 rate, and the effect of parallelization is clearly seen. For Δt103, however, the spin-up error becomes the dominant error contribution. The optimal choice of t0 begins to fall rapidly, as the sampling and spin-up must compete for computational resources under the budget. Once the spin-up error dominates, the paradigm by which (53) is controlled shifts from the exp(t0) term to the Ts1 term as Δt/Td0, since resolving Ts delivers both spin-up and sampling error control.

This interdependence will evidently have an effect on the overall scaling between cost and error, which we now seek to understand. Here, we study the variation of emodel,MC with U under the optimal choices and evaluate how well emodel,MC approximates experimental data for E[|JMCJ|]. In Fig. 20, the variation of emodel,MC computed via (64) as a function of Mens and U is shown. From this figure, we can see that, in the limit of small error, the sampling costs dominate and the best possible rate is given by the estimate in (44), limited by the CLT. On the other hand, when the cost is more moderate, scaling of the error is close to the discretization error convergence rate in (19). In this region, the spin-up costs are significant, and high-order discretization brings the state more efficiently to the start of sampling. In the spin-up dominated region, the effect of the parallel ensemble approach is minimal since spin-up must be overcome on each processor.

FIG. 20.

Optimal non-dimensional error under model as a function of total cost U for RK3. Theory totem on the left-hand side: discrete convergence rate, 1/q; on right-hand side: 2(q+r)q rate from (44).

FIG. 20.

Optimal non-dimensional error under model as a function of total cost U for RK3. Theory totem on the left-hand side: discrete convergence rate, 1/q; on right-hand side: 2(q+r)q rate from (44).

Close modal

Now, we validate the total error model for the Lorenz system by a final numerical experiment. At each choice of Mens and U, we generate 1000 individual realizations of JMC at the computed (Δt)MC,opt and NMC,opt and using the model fit given in Table II. In Figs. 21–23, we show the predictions and the results of Monte Carlo estimates of E[|JMCJ|] for our three discretizations. These results validate the model with significant discrepancies only when the asymptotic assumptions—Δt small and Ts large—do not hold, due to budget limitations in the limit of small U.

FIG. 21.

Total cost model and Monte Carlo validation as a function of total cost U for FE.

FIG. 21.

Total cost model and Monte Carlo validation as a function of total cost U for FE.

Close modal
FIG. 22.

Total cost model and Monte Carlo validation as a function of total cost U for RK3.

FIG. 22.

Total cost model and Monte Carlo validation as a function of total cost U for RK3.

Close modal
FIG. 23.

Total cost model and Monte Carlo validation as a function of total cost U for RK4.

FIG. 23.

Total cost model and Monte Carlo validation as a function of total cost U for RK4.

Close modal

In this manuscript, we have developed a theoretical framework for the total error incurred by the discrete sampling of mean outputs of ergodic ODEs. These findings are validated by Monte Carlo studies of the Lorenz system using Runge–Kutta methods. We incorporate effects of parallelization and spin-up and validate that the models match observed results in experiments. Using these models, we are able to develop a comprehensive understanding of the relationship between the wall-clock cost of a simulation and the amount of error in expectation that it might achieve.

The first result of these studies is to demonstrate the central limit theorem that gives a lower limit of 1/2 on the convergence rate of error with respect to sampling cost in a simulation. The effect of discretization error is to penalize the observed convergence rate away from this limiting convergence rate, causing a higher convergence rate than the 1/2 rate for the central limit theorem. In order to minimize the error that might be expected in a given simulation, there exists a specific choice of Ts=NsΔt that achieves the minimum in expectation given a number of timesteps Ns available for computation.

In addition to this finding, we show that ensemble simulations reduce the sampling error, albeit with a similar penalty for discrete simulations. Moreover, the effect of spin-up costs—assuming optimal choices—is to increase the error at a given amount of the total simulation cost, which cannot be reduced by parallelization. The effect of spin-up, however, is negligible in the limit of infinite costs, but the results suggest it should be considered in cost-constrained contexts.

A key problem with the applicability of this research presented in this paper is the expense of identifying the parameters of the error model. In order to overcome this, we believe that leveraging a Bayesian approach as in Oliver et al.20 can allow us to approximate the model in (20) at relatively small cost, and then exploit the result to conduct a high-fidelity simulation at (approximately) optimal discretizations. Furthermore, the framework must be extended to handle chaotic partial differential equation systems as opposed to ordinary differential equation systems. Though many discrete partial differential equation (PDE) systems are discretized in a form that reduces to an ODE system, a rigorous model for the error and cost of a PDE system should account for the contributions of both temporal discretization and spatial discretization. These will be the primary concerns of our forthcoming work. Finally, for the most expensive systems, it is likely that cost constraints will make simulations that reside comfortably inside the asymptotic convergence region prohibitively expensive. Further investigation of the convergence behavior as Δt leaves the asymptotic region or at the convergence limit is left to future research.

The authors would like to acknowledge the support of The Boeing Company (technical monitor Dr. Andrew Cary).

The authors have no conflicts to disclose.

Cory V. Frontin: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Software (equal); Writing – original draft (equal). David L. Darmofal: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1.
S. T.
Bose
and
G. I.
Park
, “
Wall-modeled large-eddy simulation for complex turbulent flows
,”
Annu. Rev. Fluid Mech.
50
,
535
(
2018
).
2.
A.
Garai
,
L. T.
Diosady
,
S. M.
Murman
, and
N. K.
Madavan
, “
Scale-resolving simulations of bypass transition in a high-pressure turbine cascade using a spectral element discontinuous Galerkin method
,”
J. Turbomach.
140
,
031004
(
2018
).
3.
R.
Stoll
,
J. A.
Gibbs
,
S. T.
Salesky
,
W.
Anderson
, and
M.
Calaf
, “
Large-eddy simulation of the atmospheric boundary layer
,”
Boundary-Layer Meteorol.
177
,
541
581
(
2020
).
4.
K.
Duraisamy
, “
Perspectives on machine learning-augmented Reynolds-averaged and large eddy simulation models of turbulence
,”
Phys. Rev. Fluids
6
,
050504
(
2021
).
5.
L. T.
Diosady
and
S. M.
Murman
, “
Higher-order methods for compressible turbulent flows using entropy variables
,” in
53rd AIAA Aerospace Sciences Meeting
(
AIAA
,
2015
), p.
0294
.
6.
M.
Franciolini
and
S. M.
Murman
, “
Multigrid preconditioning for a space-time spectral-element discontinuous-Galerkin solver
,” in
AIAA Scitech 2020 Forum
(
AIAA
,
2020
), p.
1314
.
7.
M. J.
Lighthill
, “
The recently recognized failure of predictability in Newtonian dynamics
,”
Proc. R. Soc. London, Ser. A
407
,
35
50
(
1986
).
8.
J.-P.
Eckmann
and
D.
Ruelle
, “
Ergodic theory of chaos and strange attractors
,” in
The Theory of Chaotic Attractors
(
Springer
,
1985
), pp.
273
312
.
9.
D. R.
Chapman
, “
Computational aerodynamics development and outlook
,”
AIAA J.
17
,
1293
1313
(
1979
).
10.
P.
Spalart
,
W.
Jou
,
M.
Strelets
,
S.
Allmaras
 et al, “
Comments on the feasibility of LES for wings, and on a hybrid RANS/LES approach
,” in
Advances in DNS/LES
,
1997
.
11.
H.
Choi
and
P.
Moin
, “
Grid-point requirements for large eddy simulation: Chapman's estimates revisited
,”
Phys. Fluids
24
,
011702
(
2012
).
12.
X. I.
Yang
and
K. P.
Griffin
, “
Grid-point and time-step requirements for direct numerical simulation and large-eddy simulation
,”
Phys. Fluids
33
,
015108
(
2021
).
13.
J.
Kim
,
P.
Moin
, and
R.
Moser
, “
Turbulence statistics in fully developed channel flow at low Reynolds number
,”
J. Fluid Mech.
177
,
133
166
(
1987
).
14.
A.
Lozano-Durán
and
J.
Jiménez
, “
Effect of the computational domain on direct simulations of turbulent channels up to
Reτ=4200,”
Phys. Fluids
26
,
011702
(
2014
).
15.
J. C.
Del Álamo
,
J.
Jiménez
,
P.
Zandonade
, and
R. D.
Moser
, “
Scaling of the energy spectra of turbulent channels
,”
J. Fluid Mech.
500
,
135
144
(
2004
).
16.
K. A.
Goc
,
O.
Lehmkuhl
,
G. I.
Park
,
S. T.
Bose
, and
P.
Moin
, “
Large eddy simulation of aircraft at affordable cost: A milestone in computational fluid dynamics
,”
Flow
1
,
E14
(
2021
).
17.
R. L.
Thompson
,
L. E. B.
Sampaio
,
F. A.
de Bragança Alves
,
L.
Thais
, and
G.
Mompean
, “
A methodology to evaluate statistical errors in DNS data of plane channel flows
,”
Comput. Fluids
130
,
1
7
(
2016
).
18.
S.
Russo
and
P.
Luchini
, “
A fast algorithm for the estimation of statistical error in DNS (or experimental) time averages
,”
J. Comput. Phys.
347
,
328
340
(
2017
).
19.
C.
Mockett
,
T.
Knacke
, and
F.
Thiele
, “
Detection of initial transient and estimation of statistical error in time-resolved turbulent flow data
,” in
Proceedings of the 8th International Symposium on Engineering Turbulence Modelling and Measurements
(
European Research Collaboration on Flow Turbulence and Combustion
,
2010
), pp.
9
11
.
20.
T. A.
Oliver
,
N.
Malaya
,
R.
Ulerich
, and
R. D.
Moser
, “
Estimating uncertainties in statistics computed from direct numerical simulation
,”
Phys. Fluids
26
,
035101
(
2014
).
21.
M.
Lee
and
R. D.
Moser
, “
Direct numerical simulation of turbulent channel flow up to
Reτ5200,”
J. Fluid Mech.
774
,
395
415
(
2015
).
22.
A. M.
Stuart
, “
Numerical analysis of dynamical systems
,”
Acta Numer.
3
,
467
572
(
1994
).
23.
M.
Denker
, “
The central limit theorem for dynamical systems
,”
Banach Cent. Publ.
23
,
33
62
(
1989
).
24.
R. C.
Bradley
, “
Basic properties of strong mixing conditions. a survey and some open questions
,”
Probab. Surv.
2
,
107
144
(
2005
).
25.
V.
Araújo
,
I.
Melbourne
, and
P.
Varandas
, “
Rapid mixing for the Lorenz attractor and statistical limit laws for their time-1 maps
,”
Commun. Math. Phys.
340
,
901
938
(
2015
).
26.
E.
Hairer
,
Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems
, 2nd ed., Springer Series in Computational Mathematics Vol. 14 (
Springer
,
Berlin, Germany
,
1993
).
27.
D.
Viswanath
, “
Global errors of numerical ODE solvers and Lyapunov's theory of stability
,”
IMA J. Numer. Anal.
21
,
387
406
(
2001
).
28.
E. N.
Lorenz
, “
Deterministic nonperiodic flow
,”
J. Atmos. Sci.
20
,
130
141
(
1963
).
29.
C.
Sparrow
,
The Lorenz Equations: Bifurcations, Chaos, and Strange Attractors
(
Springer Science & Business Media
,
1982
).
30.
J.
Dormand
,
R.
Duckers
, and
P.
Prince
, “
Global error estimation with Runge-Kutta methods
,”
IMA J. Numer. Anal.
4
,
169
184
(
1984
).
31.
K. E.
Trenberth
, “
Some effects of finite sample size and persistence on meteorological statistics. Part I: Autocorrelations
,”
Mon. Weather Rev.
112
,
2359
2368
(
1984
).
32.
F.
Harris
, “
On then use of windows for harmonic analysis with the discrete Fourier transform
,”
Proc. IEEE
66
,
51
(
1978
).
33.
G.
Karniadakis
and
S.
Sherwin
,
Spectral/hp-Element Methods for Computational Fluid Dynamics
(
Oxford University Press
,
2005
).
34.
V.
Makarashvili
,
E.
Merzari
,
A.
Obabko
,
A.
Siegel
, and
P.
Fischer
, “
A performance analysis of ensemble averaging for high fidelity turbulence simulations at the strong scaling limit
,”
Comput. Phys. Commun.
219
,
236
245
(
2017
).