1. Introduction
It is well known (e.g. Kim & Leonard Reference Kim and Leonard2024) that direct numerical simulations (DNS) are a powerful tool for studying the fundamental physics of turbulence in canonical settings. The first DNS were performed for a three-dimensional (3-D) periodic domain more than 50 years ago by Orszag & Patterson (Reference Orszag and Patterson1972), with
$32^3$
grid points at Taylor-scale Reynolds number (
$R_\lambda$
) 64. Since then, with computer power growing approximately a million times every 25 years, the drive towards higher Reynolds numbers has been a prime motivation for turbulence researchers (Yokokawa et al. Reference Yokokawa, Itakura, Uno, Ishihara and Kaneda2002; Lee, Malaya & Moser Reference Lee, Malaya and Moser2013; Ravikumar, Appelhans & Yeung Reference Ravikumar, Appelhans and Yeung2019) seeking to harness resources that have marked the arrival of the terascale, petascale and more recently exascale eras (Atchley et al. Reference Atchley2023). At the same time, it has become clear (Yakhot & Sreenivasan Reference Yakhot and Sreenivasan2005) that scales finer than the Kolmogorov scales appear more copiously as the Reynolds numbers increase, which implies a need to resolve the small scales better than was often practised in the past, and limits the Reynolds number achievable at the required resolution. This consideration also makes the pursuit of extreme-scale simulation capabilities ever more important.
An important science question requiring access to 3-D data at increasing Reynolds numbers and improved spatial and temporal resolution at the smallest scales of turbulence is whether scaling and other trends observed at moderate Reynolds numbers continue at ever-increasing Reynolds number, or whether qualitative transitions occur at some point. An important example of such trends concerns enstrophy having been found to be more important than dissipation in most known data sources. Past literature (Nelkin Reference Nelkin1999) included arguments that since dissipation and enstrophy are both small-scale quantities belonging to the same symmetry group, the tails of their probability density functions (PDFs) must have the same scaling characteristics (i.e. their moments could differ by no more than a Reynolds-number-independent constant). Despite recent progress (Yeung, Sreenivasan & Pope Reference Yeung, Sreenivasan and Pope2018) that addressed the impact of finite resolution, high-fidelity data at significantly higher Reynolds numbers are still required.
The largest Fourier pseudo-spectral DNS dataset that can fit into the memory of the world’s first exaflop leadership class supercomputer (named Frontier) is
$40\,960^3$
, but we focus on
$32\,768^3$
, which is more realistic for our present purposes. However, the viability of a specific simulation still depends on both the cost per step (at least proportional to
$N^3$
on an
$N^3$
grid) and the number of time steps for a desired physical time span (as
$N$
increases). As reported in Yeung et al. (Reference Yeung, Ravikumar, Nichols and Uma-Vaideswaran2025), the best timing for a
$ 32\,768^3$
grid using 44 % of the machine is 13 s per step. At
${R_\lambda }\approx 2500$
, each Kolmogorov time scale (
$ \tau _\eta$
) takes approximately 1000 time steps, while the ratio of large-eddy turnover time (
$T_E$
) to
$ \tau _\eta$
is approximately 300. At this scale, a long simulation spanning
$10T_E$
enabling long-time averaging would involve two years or more of virtually non-stop computing. It is clear that a different strategy is necessary.
This paper is focused on recently available DNS datasets at high resolution and high Reynolds numbers that may be used to address fundamental questions pertaining to properties of the small scales of turbulence, such as on comparisons between the asymptotic scaling of tails of dissipation and enstrophy, at increasing Reynolds numbers. The data are obtained using the multi-resolution independent simulations (MRIS) approach introduced by Yeung & Ravikumar (Reference Yeung and Ravikumar2020), here extended to facilitate attainment of higher Reynolds numbers. The technique will be described in more detail in § 2, but it suffices here to say that this method includes ensemble averaging over short segments (referred to as MRIS segments) with statistically equivalent initial conditions. The short duration of the simulations calls for a separate assessment of inertial range properties, including some tests of statistical convergence, but dissipative scales that are spatially and temporally finer appear to be adequately captured. We present these results in §§ 3 and 4, including fundamental properties of the energy spectrum, the scaling of the mean energy dissipation rate, and new insights into the PDFs of energy dissipation and squared vorticity. We conclude the paper in § 5, where we also note additional simulation capabilities developed on Frontier.
2. Technical approach
In this section, we first briefly review our basic DNS algorithm and its high-performing implementation on the leadership-class graphics processing unit (GPU) platform on Frontier. The use of the MRIS protocol for both improving resolution in time and space and for achieving progressively higher Reynolds number is then discussed in some detail.
2.1. The GPU-enabled extreme-scale turbulence simulations on Frontier
The GPU-enabled extreme-scale turbulence simulations code (Yeung et al. Reference Yeung, Ravikumar, Nichols and Uma-Vaideswaran2025) on Frontier performs DNS of incompressible isotropic turbulence on a 3-D periodic domain. Numerical integration in time is carried out in wavenumber space, using the Fourier pseudo-spectral formulation of Rogallo (Reference Rogallo1981), with solenoidal numerical forcing at the largest scales (in the lowest few wavenumber shells). The principal mathematical operation is a discrete 3-D fast Fourier transform (FFT), which is invoked multiple times between physical space and wavenumber space at each Runge–Kutta substep. Distributed-memory parallelism requires a solution domain partitioned in one or two directions simultaneously, while FFTs are taken one direction at a time. To complete a 3-D FFT, the domain must be re-divided in different directions, which is achieved by exchanging data among multiple parallel processes via message-passing communication.
Data movements in communication calls often add substantially to the overall simulation cost, especially at large problem sizes or large node counts, ultimately limiting the size of problem accessible on a given machine. However, the arrival of exascale computing defined by a theoretical peak of
$10^{18}$
floating point operations per second has led to improved prospects, by enabling much faster computation, while increased memory per node allows using fewer parallel processes, thus reducing the communication overhead. In addition, while early GPU platforms often had limited memory, parity in memory between central processing unit (CPU) and GPU on Frontier allows the GPUs to perform almost all operations, including computation and communication, thus reducing data movement between CPUs and GPUs to a minimum. Further algorithmic details and per-step timing data are given in Yeung et al. (Reference Yeung, Ravikumar, Nichols and Uma-Vaideswaran2025).
2.2. Successive grid refinement and increasing the Reynolds number
We next consider the physical and numerical issues that determine the number of time steps required for specific purposes, and how best to push the limits of available computing power to obtain high-Reynolds-number results for specific purposes.
The MRIS approach is built on the premise that while long simulations at high resolution may not be attainable, long simulations at low resolution are much more affordable, and that high-resolution results can be obtained from the latter by refining resolution in space and/or time, and then running at high resolution for only a short period of time. In pseudo-spectral methods, grid refinement results in an additional set of high-wavenumber modes, which soon start to receive energy by spectral transfer from wavenumbers close to the maximum represented on the coarser grid. Some small-scale statistics will adjust quickly because their own time scales are typically short. Indeed, there is evidence (Yeung et al. Reference Yeung, Sreenivasan and Pope2018; Yeung & Ravikumar Reference Yeung and Ravikumar2020) that when resolution is improved by halving the time step size or the grid spacing, it takes only approximately
$ 2\tau _\eta$
for some equilibration to occur of small scales on the order of
$\eta$
of the coarser grid. The core MRIS strategy thus consists of a three-phased sequence: (1) running a reasonably long simulation at temporal and spatial resolution adequate for scales larger than
$\eta$
; (2) continuing but halving the time step and running for a few
$ \tau _\eta$
; and finally (3) halving the grid spacing and running at the highest resolution also for a small number of
$ \tau _\eta$
. For smoother transitions in the numerics, increased resolution can be implemented as a series of small increments: e.g. by starting from
$N_1$
grid points in each direction, and first refining to
$3N_1/2$
before ultimately reaching
$N_2=2N_1$
. We are also able to improve statistical sampling by repeating this sequence using different instantaneous velocity field snapshots for phase 1, and taking ensemble averages over multiple MRIS segments. Generally, phase 2 in each MRIS segment involves reducing the Courant number
$C = [|u|+|v|+|w|]_{max}\, \Delta t/\Delta x$
(where the velocity maximum is taken over all grid points) from 0.6 to 0.3, while phase 3 involves increasing
${ k_{max}\eta }=(\sqrt {2} N/3) \eta \approx 2.96 /(\Delta x /\eta )$
from 1.4 to 2.8. The value of
$\eta$
corresponds to the highest resolution, but it does not change significantly with resolution.
Table 1 shows the main parameters in the computations reported in this paper. Although the number of phase 3 simulations available is not large, there is still appreciable benefit from ensemble averaging over the initial conditions (which is rarely discussed in the literature). We have used the forcing scheme of Donzis & Yeung (Reference Donzis and Yeung2010), where the total modal energy in each of the first three wavenumber shells (
$k\lt k_F$
, where
$k_F= 2.5$
on a
$(2\pi )^3$
domain) is kept constant at values derived from long-time averages from stochastic forcing (Eswaran & Pope Reference Eswaran and Pope1988). A test of a single segment of
$R_\lambda=1000$
at
$8192^3$
and spanning
$27 { \tau _\eta }$
suggests that the results presented are insensitive to the values of
$T_3/{ \tau _\eta }$
in the table, presumably because the energy of the large scales is kept frozen by the forcing scheme. For consistency, we have used MRIS results at all Reynolds numbers except the lowest, even if longer simulations at lower Reynolds numbers are available from prior work. The MRIS procedure also has some conceptual similarities with scale enrichment in large-eddy simulations (Ghate & Lele Reference Ghate and Lele2020).
Basic parameters in this work: for each approximate
$R_\lambda$
, highest grid resolution (
$N^3$
) used, kinematic viscosity (
$\nu$
), number of MRIS segments (
$M$
), and (see text in § 2.2) time spans
$T_1$
,
$T_2$
,
$T_3$
in
$ \tau _\eta$
for phases 1–3, with
$\{{ k_{max}\eta },C\}\approx \{1.4,0.6\}$
,
$\{1.4,0.3\}$
and
$\{2.8,0.3\}$
, respectively. Here,
$T_1$
,
$T_2$
and
$T_3$
are averaged over different segments since they have varied somewhat. Phases 1 and 2 were not needed for the
$R_\lambda=650$
case, which started from a
$C=0.15$
stationary state in Yeung et al. (Reference Yeung, Sreenivasan and Pope2018). Early portions of data from phase 3 are not used in the time averaging.
It should be noted that all of the simulations reported used single precision arithmetic. Numerical tests in Yeung et al. (Reference Yeung, Sreenivasan and Pope2018) focused on the peak values of dissipation and enstrophy showed that effects of machine precision were small compared to those associated with spatial and temporal resolution as well as statistical sampling. Recently, there has been increasing interest in DNS regimes where round-off errors may accumulate over time (Okamoto et al. Reference Okamoto, Ishihara, Yokokawa and Kaneda2025
a), while a rigorous computer science approach can provide insights into the effects of different levels of machine precision on quantities of different types (Karp et al. Reference Karp2025). While rigorous testing is not readily available, based on the above, and since our MRIS segments are only approximately
$ 2\tau _\eta$
long, it seems unlikely that machine precision would have strongly affected the conclusions drawn in the present work.
To reach higher Reynolds numbers, we start from developed velocity fields at, say,
${R_\lambda }_1$
on an
$N_1^3$
grid, and reduce the viscosity while keeping both the forcing details and the resolution parameter
$ k_{max}\eta$
unchanged. We refine the grid to
$N_2^3$
with
$\eta \propto 1/N$
. Since the mean dissipation rate is nearly independent of viscosity in the stationary state, the definition
$\eta = (\nu ^3/{ \langle \epsilon \rangle })^{1/4}$
implies
$\nu \propto (1/N)^{4/3}$
. In addition, with the root mean square velocity
$u'$
also unchanged, the standard local isotropy relation
$\langle \epsilon \rangle = 15 \nu (u'/\lambda )^2$
predicts
${R_\lambda }\propto N^{2/3}$
. This implies that at a higher resolution, the expected new Reynolds number can be estimated by
${R_\lambda }_2\approx {R_\lambda }_1 (N_2/N_1)^{2/3}$
, which amounts to a 58 % increase in the case
$N_2=2 N_1$
.
(a) Evolution of kinetic energy (black), mean dissipation (red),
$S_\epsilon$
(green) and
$R_\lambda$
(blue), all normalised by initial values, during transition from
$R_\lambda=2080$
on a
$12\,288^3$
grid to
$R_\lambda=2500$
on
$16\,384^3$
, by reducing the viscosity. Time
$t$
is normalised by the final value of
$\tau _\eta$
. (b) Plot of
$D(k)$
compared with data at
$R_\lambda=2080$
data (gold), at
$t/\tau _\eta \approx 0.01$
, 5.7, 12.8, 20.5, 28.8, 37.8 (red, green, blue, black, cyan, magenta, respectively).
Figure 1 shows the manner in which quantities dominated by different scale sizes evolve during transition to a higher Reynolds number, when subjected to the procedures described above. We have chosen an example where the velocity field evolves from
$R_\lambda=2080 $
, which served as an intermediate step between
$R_\lambda=1600 $
and
$R_\lambda=2500 $
. In figure 1(a), it can be seen that the kinetic energy varies little, while mean dissipation recovers from the sudden reduction in viscosity within approximately
$ 20\tau _\eta$
. For a diagnostic more strongly linked to the small scales, we consider the ‘dissipation skewness’ defined from (despite its name) a fourth-moment integral of the energy spectrum, i.e.
$S_\epsilon = (4/35)(15\nu /{ \langle \epsilon \rangle })^{3/2} \int _0^\infty \nu k^4\, E(k)\,{\rm d}k$
, which is equal to the negative of the velocity gradient skewness in stationary isotropic turbulence. As expected,
$S_\epsilon$
recovers sooner. When grid refinement begins,
$R_\lambda$
rises artificially from its stable value
${R_\lambda }_1$
on a
$12\,288^3$
grid to
${R_\lambda }(0)={R_\lambda }_1/0.6814$
, then decreases to approximately
$0.84{R_\lambda }_0=1.23{R_\lambda }_1$
, which agrees well with the expectation of a 21.1 % increase (corresponding to
$N_2/N_1=4/3$
).
Unlike the quest for resolution improvement alone in Yeung & Ravikumar (Reference Yeung and Ravikumar2020), here a reduction of viscosity directly affects almost all Fourier modes. Figure 1(b) shows, without normalisation, that the dissipation spectrum settles very quickly onto a new and approximately stationary form at higher Reynolds number. Results in this figure correspond to ‘phase 1’ in our MRIS protocol, for a case of grid refinement by ratio
$N_2/N_1=4/3$
. In most cases, we operate with
$N_2/N_1=2$
, which requires a longer transition period, with the Reynolds number dropping further before increasing towards the new value expected on the finer grid.
One of the three
$32\,768^3$
snapshots (3-D fields of three velocity components and pressure, encompassing half a petabyte of data) has been ingested into a public database system (the Johns Hopkins Turbulence Database) where the data are accessible using web-services-enabled virtual sensors.
3. Basic scaling behaviours versus the Reynolds number
Scaling behaviours with respect to
$R_\lambda$
of (a) integral-scale Reynolds number, (b) length scale ratio (triangles) and time-scale ratio (circles), and (c) the normalised energy dissipation rate. Dotted lines in (a) and (b) show expected power laws. Horizontal dashed lines in (c) indicate the range in various simulations.
As seen in table 1 the DNS data reported here cover more than six-fold spread of
$R_\lambda$
, with the values at the upper end exceeding those in most of the DNS literature. It is thus useful to compare the new data with some basic classical scaling relations, to confirm the fidelity of the new results and to check if classical scaling continues to apply at these higher Reynolds numbers without qualitative change. Figure 2(a) shows the large-scale Reynolds number (
$R_\ell$
, where the length scale
$\ell$
is taken to be the longitudinal integral length scale
$L_1$
of the velocity field) versus the Taylor-scale Reynolds number
$R_\lambda$
. The dashed line corresponds to
$C_\epsilon {R_\lambda }^2/15$
, by assuming that
$C_\epsilon ={ \langle \epsilon \rangle } L_1/u'^3$
is 0.45 (figure 2
c). Figure 2(b) shows the length scale ratio
$L_1/\eta$
and the time scale ratio
$(L_1/u')/{ \tau _\eta }$
. The lines show the expected
$3/2$
power and the linear relations, again taking
$C_\epsilon = 0.45$
. Figure 2(c) shows the variation of the normalised mean energy dissipation rate
$C_\epsilon ={ \langle \epsilon \rangle } L_1/u'^3$
, vis-à-vis the dissipative anomaly or the zeroth law of turbulence (Sreenivasan Reference Sreenivasan1984, Reference Sreenivasan1998). The data show a weak trend up to
${R_\lambda }\approx 1500$
, beyond which they stabilise to a constant, slightly larger than 0.4. Continued validity of the zeroth law of turbulence to
$R_\lambda$
over 2500 has important implications for turbulence theory and modelling practices.
(a) The 3-D energy spectrum
$E(k)$
at
${R_\lambda }\approx390 $
(red), 650 (green), 1000 (blue), 1600 (black), 2500 (magenta), after averaging in time and over available MRIS segments. The dotted line has slope
$-5/3$
. (b) Compensated 3-D spectrum
$\langle \epsilon \rangle ^{-2/3}k^{5/3}E(k)$
(upper curve) and the one-dimensional longitudinal spectrum
$E_{11}(k)$
in the Kolmogorov variable. The dashed horizontal lines are the corresponding (generally accepted) values 1.62 and
$(1.62)(18/55)=0.53$
.
A most investigated but still elusive topic in fundamental turbulence theory is the form of the energy spectrum in the inertial range. Figure 3(a) suggests a range
$E(k)\propto k^{-5/3}$
, which widens with an increase of
$R_\lambda$
. Yet a closer look at figure 3(b) shows that the compensated quantity
$E(k){ \langle \epsilon \rangle }^{-2/3}k^{5/3}$
does not display a plateau. Instead, the data support a power law for
$E(k)$
with a slope steeper than
$-5/3$
, and a bump, called the bottleneck, which is strongest at
$k\eta$
between 0.1 and 0.2. The first observation is consistent with the phenomenology of intermittency (Kolmogorov Reference Kolmogorov1962) that predicts a correction (steepening) of the slope and has motivated extensive research, including analyses of experimental data on one-dimensional spectra at higher Reynolds numbers (see Sreenivasan & Antonia Reference Sreenivasan and Antonia1997). The second issue has also spawned past papers (Donzis & Sreenivasan Reference Donzis and Sreenivasan2010) and proposed spectral fitting functions that include both intermittency correction and bottleneck (Meyers & Meneveau Reference Meyers and Meneveau2008). It is understood that because of the
$5/3$
exponent, the peak of the bottleneck is at a wavenumber slightly below that for the dissipation spectrum
$D(k)\propto k^2\, E(k)$
, hence viscous effects certainly play a role. Figure 3(b) also shows similar features in the one-dimensional spectrum
$E_{11}(k)$
. We have verified consistency with isotropy relations such as
$E(k)=k^3\,{\rm d}[(1/k){d}\,E_{11}(k)]/{\rm d}k$
, and noted that the bottleneck in
$E_{11}(k)$
is weaker than
$E(k)$
. Similarly to the work by Ishihara et al. (Reference Ishihara, Morishita, Yokokawa, Uno and Kaneda2016), the high Reynolds number in the current data is expected to be useful for further work on this topic.
4. The PDFs of dissipation and enstrophy density
Base-10 logarithms of PDFs of enstrophy versus the same for dissipation at (a)
$R_\lambda=390$
and (b)
$R_\lambda=2500$
. Symbols in red denote data points where
$\epsilon$
and
$\varOmega$
are both below their mean values; blue indicates samples above the mean values. Black dashed lines of slopes
$1/3$
(upper) or 1 (lower) have been added for some comparisons (in the text).
The PDFs of dissipation rate
$\epsilon$
(formed from the square of strain rate) and enstrophy density
$\varOmega \equiv \omega ^2$
(formed from the square of vorticity – henceforth enstrophy, for brevity) are prominent indicators of small-scale turbulence. Earlier measurements have suggested that the PDF tails of both quantities are well represented by stretched exponentials, with enstrophy being the more intermittent variable.
To firmly understand how the two variables behave, high Reynolds numbers, high resolution and large datasets are necessary. In this work, we have extracted the PDFs of dissipation and enstrophy from a large number of instants of time even though the MRIS segments themselves are short. Although one may question whether these PDFs have fully equilibrated from forcing over the short time interval covered by the MRIS segments, we can mitigate such concerns by plotting the two PDFs against each other, This is done in figure 4, where different colours represent the left tails (red, for small
$\epsilon$
or
$\varOmega$
) and right tails (blue, for large
$\epsilon$
or
$\varOmega$
). Deviations from the dashed line of slope 1 indicate that the two PDFs do not have identical form. We will qualify this statement below. Both very low and very high values of enstrophy are more probable than those of dissipation. The behaviour in the left tail can be explained by noting (Yeung, Donzis & Srenivasan Reference Yeung, Donzis and Srenivasan2012; Gotoh & Yeung Reference Gotoh and Yeung2025) that zones of low strain rate and low enstrophy are nearly Gaussian, hence the PDFs behave like a chi-square distribution of order
$n$
, which holds for the sum of squares of
$n$
independent and identically distributed Gaussian variables. For small enstrophy
$n=3$
, while for dissipation
$n=5$
, since there are (due to the incompressibility constraint) five independent strain-rate components. Consequently (for small
$x$
), the chi-square result is
$f(x)\propto x^{n/2-1}$
, which then gives
$f_\epsilon (x)\propto x^{3/2}$
, while
$f_\varOmega (x)\propto x^{1/2}$
. Combining these results, we obtain
\begin{equation} f_\varOmega (x) \propto [(f_\epsilon (x)^{2/3}]^{1/2} \propto [f_\epsilon (x)]^{1/3}, \end{equation}
which is confirmed by a dashed line of slope
$1/3$
in figure 4. This result is kinematic and holds equally well at both Reynolds numbers of this figure.
The right tails of
$f_\epsilon ({\cdot })$
and
$f_\varOmega ({\cdot })$
can be represented quite well by expressions of the form
\begin{equation} f_\epsilon (x)\propto \exp (-a_\epsilon x^{b_\epsilon }),\quad f_\varOmega (x)\propto \exp (-a_\varOmega x^{b_\varOmega }). \end{equation}
The data are reasonably represented by these expressions with
$b_\varOmega =b_\epsilon$
, in which case extremely high
$\varOmega$
, being more probable, implies
$a_\varOmega \lt a_\epsilon$
. The right tail portion in figure 4(b) appears close to linear, with slope close to 0.9. Thus we can conclude that differences in high-order moments of the two variables occur only via the coefficients
$a_\varOmega$
and
$a_\epsilon$
.
Test of the stretched exponential form of the right tails of the PDFs of
$\epsilon /{ \langle \epsilon \rangle }$
(red) and
$\varOmega /{ \langle \varOmega \rangle }$
at approximately (a) 390, (b) 1000, (c) 2500. Here,
$f_X({\cdot })$
denotes the PDF of either
$\epsilon /{ \langle \epsilon \rangle }$
or
$\varOmega /{ \langle \varOmega \rangle }$
. The dotted lines are at slope 0.518, which corresponds to 0.225 reported in Gotoh & Yeung (Reference Gotoh and Yeung2025) using the common log instead of the natural log as here.
The stretched exponential expressions above are relevant in the ranges where
$\epsilon /{ \langle \epsilon \rangle }$
and
$\varOmega /{ \langle \varOmega \rangle }$
far exceed unity. For a closer examination of the parameters involved, figure 5 shows, for three Reynolds numbers, the quantity
$\ln (-\ln [f_\epsilon (x)/f_\epsilon (1)])$
versus
$\epsilon /{ \langle \epsilon \rangle }$
for
$x\gt 1$
only, and likewise for the enstrophy. Perfect fits give straight lines of slope equal to the exponent and intercepts on the
$y$
-axis equal to the prefactor. As
$R_\lambda$
increases, non-zero samples are found further to the right, while the lines for dissipation (red) and enstrophy (blue) become more nearly parallel, with the latter having a smaller intercept on the
$y$
-axis. These observations are consistent with the discussion of figure 4 above, while also demonstrating the value of high Reynolds number in establishing trends more conclusively than possible before.
Reynolds number dependence of (a) normalised moments of
$\epsilon$
, orders
$n=2,3,4$
(bottom to top), (b) normalised moments of
$\varOmega$
, and (c) ratio of normalised moments of
$\varOmega$
to those of
$\epsilon$
. Red triangles are from integration of time-averaged PDFs. Blue circles are from averaging in physical space over several instantaneous snapshots.
The PDFs obtained here are smooth enough, extending to large values of the respective variable, that they can be used to compute high-order moments as measures of intermittency. Having checked for satisfactory convergence of the moment integrands associated with dissipation and enstrophy PDFs, in figure 6 we show the behaviour of second, third, and fourth moments of
$\epsilon /{ \langle \epsilon \rangle }$
and
$\varOmega /{ \langle \varOmega \rangle }$
at all five Reynolds numbers of table 1. As expected from the contrast in their PDFs, the moments of the normalised enstrophy for any given order
$n\gt 1$
are higher than those of normalised dissipation. Their ratios in figure 6(c) are nearly independent of the Reynolds number, which is consistent with arguments by Nelkin (Reference Nelkin1999) that while enstrophy is more intermittent than dissipation, their scaling with Reynolds number must be the same asymptotically.
We note in passing that the skewness and the flatness factor of the longitudinal velocity gradient agree well with the data compiled by Sreenivasan & Antonia (Reference Sreenivasan and Antonia1997).
5. Concluding remarks
Three-dimensional turbulence governed by the Navier–Stokes equations continues to be a grand challenge in computational science, even as leadership-class computing power has recently advanced past the exascale. Direct numerical simulations (DNS) at
$32\,768^3$
resolution are now possible on the world’s first exascale computer (named Frontier). However, as the range of time scales also grows, full-length simulations spanning multiple large-eddy time scales at this high resolution remain prohibitively expensive, especially if rigorous resolution requirements (Yakhot & Sreenivasan Reference Yakhot and Sreenivasan2005; Yeung et al. Reference Yeung, Sreenivasan and Pope2018) are to be met. In this paper we have not fully addressed the implications of these rigorous requirements for all aspects of DNS, but we have shown that a protocol named MRIS (for multi-resolution independent simulations, introduced in Yeung & Ravikumar Reference Yeung and Ravikumar2020), consisting of ensemble averaging over results from a number of short simulation segments at high resolution, can produce reliable results on small-scale turbulence.
This work is focused on statistically stationary isotropic turbulence with a forcing scheme, which produces very stable energetics at the large scales. Reynolds numbers are increased by reducing the viscosity at modest resolution, followed by short simulations that improve accuracy by increasing both temporal and spatial resolution. The highest Taylor-scale Reynolds number reached is approximately 2500. Even though the MRIS segments for resolution improvement are sometimes as short as approximately two Kolmogorov scales, many fundamental scaling properties, including the property of dissipative anomaly, and effects of intermittency and a bottleneck in the energy spectrum, are reproduced well. This study also successfully confirms and extends prior knowledge of the scaling properties of the probability densities and moments of the energy dissipation rate and enstrophy, where extreme events are of great interest. Our results show that these PDFs behave like chi-square distributions at low amplitudes, but as stretched exponentials at high amplitudes and asymptotics, consistent with the arguments provided by Nelkin (Reference Nelkin1999).
Turbulence simulations have been advancing steadily in size,with the most prominent drivers reaching higher Reynolds number and achieving better resolution at the small scales. As Reynolds numbers increase, new results may either lead to confirmation of prior trends, thus removing doubts that might otherwise blunt the path forwards, or suggest new physics not apparent from results at lower Reynolds number or insufficient resolution. Examples in this work include figures 2(c) and 6(a,b), where conclusions would differ if data points at
${R_\lambda }\approx 1600$
and 2500 were absent. Both types of contributions continue to provide motivation for the meticulous work (as documented in Yeung et al. Reference Yeung, Ravikumar, Nichols and Uma-Vaideswaran2025) required to extract the best performance from the best machines available. (The magnitude of the computational effort exerted in obtaining the present results can be inferred from the triangles in figure 2(b), where a 2.5-fold increase in
$R_\lambda$
from 1000 to 2500 (third data point to the last) requires a 4-fold increase in the number of grid points in each direction – which translates to at least a (
$4^3=64$
)-fold increase in elapsed wall time per step, and a (
$4^4=256$
)-fold increase for a given physical time span, on a machine with at least 64 times more memory.)
We note again that this paper provides primarily dissipation range results, while a discussion of inertial range quantities is postponeed to a later paper. The MRIS technique may not be as directly applicable when one wishes to study phenomena that are allowed to evolve on time scales associated with motions in the inertial or energy-containing range. Examples include non-stationary turbulence or flows subjected to forcing schemes that supply a highly variable rate of energy input. It would seem that longer MRIS segments may help us to understand the dominant phenomena of interest that occur over intervals of several Kolmogorov time scales. Continuing work includes extensions of MRIS to study the small scales in passive scalar mixing, and short-time Lagrangian properties such as statistics of the fluid particle acceleration or rapidly separating particle pairs in relative dispersion.
Note added in proof
The authors have become aware of a new manuscript (Okamoto et al. Reference Okamoto, Ishihara, Yokokawa and Kaneda2025 b) which provides some new insights relevant to figure 3(b) in this paper.
Acknowledgements
We thank especially Dr B. Messer, Director of Science at OLCF, for his constant support and encouragement. The first author also wishes to thank Professors T. Gotoh, Y. Kaneda and S.K. Lele for helpful recent discussions.
This manuscript has been authored in part by U T-Battelle, LLC, under contract DE-AC05-00OR22725 with the US Department of Energy (DOE). The publisher acknowledges the US government licence to provide public access under the DOE Public Access Plan (http://energy .gov/downloads/doe-public-access-plan).
Funding
This research used resources allocated through the INCITE program at the Oak Ridge Leadership Computing Facility, which is a US Department of Energy Office of Science User Facility supported under contract DE-AC05-00OR22725. The work at Georgia Tech is further supported by a subaward from NSF via CSSI grant 2103874 to Johns Hopkins University. The research of K.R.S. was funded by New York University.
Declaration of interests
The authors report no conflict of interest.
Data availability statement
A single instantaneous snapshot of the
$32\,768^3$
velocity field (including pressure fluctuations derived therefrom) obtained in this work is publicly available from the Johns Hopkins Turbulence Database. See https://turbulence.idies.jhu.edu/datasets/homogeneousTurbulence/isotropic32768.
Open access
