Accelerated evolution of convective simulations

High Peclet number, turbulent convection is a classic system with a large timescale separation between flow speeds and the thermal relaxation time. In this paper, we present a method of fast-forwarding through the long thermal relaxation of convective simulations, and we test the validity of this method. This accelerated evolution (AE) method involves measuring the dynamics of convection early in a simulation and using its characteristics to adjust the mean thermodynamic profile within the domain towards its evolved state. We study Rayleigh-B\'enard convection as a test case for AE. Evolved flow properties of AE solutions are measured to be within a few percent of solutions which are reached through standard evolution (SE) over a full thermal timescale. At the highest values of the Rayleigh number at which we compare SE and AE, we find that AE solutions require roughly an order of magnitude fewer computing hours to evolve than SE solutions.


I. INTRODUCTION
Astrophysical convection occurs in the presence of disparate timescales. Studying realistic models of natural systems through direct numerical simulations is infeasible because of the large separation between various flow timescales and relaxation times. Stiffness in astrophysical systems can manifest in multiple ways, some of which can be handled by clever choices of numerical algorithms, and some which cannot. For example, flows in the convection zones of stars like the Sun are characteristically low Mach number (Ma) in the deep interior. Initial value problems solved using explicit timestepping methods are bound by the Courant-Friedrich-Lewy (CFL) timestep limit corresponding to the fastest motions (sound waves), resulting in timesteps which are prohibitively small for studies of the deep, low-Ma motions. These systems are numerically stiff, and the difference between the sound crossing time and the convective overturn time has made studies of low-Ma stellar convection difficult. This stiffness can be avoided using approximations such as the anelastic approximation, in which sound waves are explicitly filtered out [1,2]. Recently, advanced numerical techniques which use fully implicit [3][4][5] or mixed implicit-explicit [6][7][8] timestepping mechanisms have made it possible to study convection in the fully compressible equations at low Mach numbers, and careful studies of deep convection which would have been impossible a decade ago are now widely accessible.
Unfortunately, astrophysical convective systems are stiff in more than one timescale. Specifically, the Peclet number (Pe), the ratio of the thermal diffusion timescale to the convective velocity timescale, is large. In a high Pe system, many convective timescales must pass before the domain thermally relaxes into a steady state in which evolution of the thermal structure of the convective region has ceased. As the timestep size of implicit methods is bound to the fastest nonlinear flows (convection), fully implicit methods cannot be used to address this form of stiffness [3][4][5].
Resolving dynamics in atmospheres which are sufficiently thermally relaxed therefore remains a challenging problem. Solar convection is a prime example of this phenomenon, as dynamical timescales in the solar convective zone are relatively short (convection overturns every ∼ 5 min at the solar surface) compared to the Sun's thermal relaxation timescale, which is O(10 7 ) years [9]. In such a system, it is impossible to resolve the convective dynamics while also meaningfully evolving the thermal structure of the system using traditional timestepping techniques alone. As modern simulations aim to model natural, high-Pe convection in the high-Rayleigh-number (Ra) regime, the thermal diffusion timescale (t κ , defined in Sec. III) becomes intractably large compared to dynamical timescales such as the freefall time (t ff , defined in Sec. II) [7], Furthermore, as dynamical and thermal timescales separate, simulations become more turbulent. Capturing appropriately resolved turbulent motions requires finer grid meshes and smaller timesteps. Thus, the progression of simulations into the high-Ra regime of natural convection is slowed by two simultaneous effects: timestepping through a single convective overturn time becomes more computationally expensive and the number of overturn times required for systems to thermally relax grows. The vast difference between convective and thermal timescales has long plagued numericists studying convection, and an abundance of approaches has been employed to study thermally relaxed solutions. One popular method arXiv:1807.06687v2 [physics.flu-dyn] 22 Aug 2018 for accelerating the convergence of high-Ra solutions is by "bootstrapping" -the process of using the relaxed flow fields and thermal structure of a low Ra solution as initial conditions for a simulation at high Ra. This method has been used with great success [10,11], but it is not without its faults. In systems in which there are multiple stable solutions, such as a roll state and a shearing state of convection, the large-scale convective structures present in the low Ra solution imprint onto the dynamics of the new, high Ra solution. It is possible that this puts the high Ra solution into a different stable state than it would naturally reach from an initial, hydrostatically stable configuration. Another commonly-used tactic in moderate-Ra simulations is to first solve a simple model of the system in question, and then use the solution of that model as initial conditions for full nonlinear direct numerical simulations. For example, past studies have used initial conditions from a linear eigenvalue solve [12] in plane-parallel studies, or axisymmetric solutions in studies of convection in three-dimensional (3D) cylinders [11]. In other systems, particularly when convective zones are adjacent to stable regions, authors often choose initial conditions which are not in a classic, hydrostatic, conductive state. Rather, either through knowledge of low-Ra solutions [13] or broader convective theories such as Mixing Length Theory [14], initial conditions can be adjusted such that the stratification within the convective domain is closer to a relaxed state than the largely unstable hydrostatic state.
Despite the numerous methods that have been used, the most straightforward way to achieve a thermally relaxed solution is to evolve a convective simulation through a thermal timescale. Some modern studies do just that [2]. Such evolution is computationally expensive, and state-of-the-art simulations at the highest values of Ra can only reasonably be run for hundreds of freefall timescales [15], much less the thousands or millions of freefall times required for thermal relaxation.
In this work, we study a method of achieving accelerated evolution of convective simulations. Our technique is similar in essence to the approach used in asymptotically reduced models of rapidly rotating convection, in which the mean temperature profile evolves separately from the fast convective dynamics [16,17]. We couple measurements of the dynamics of unstable, evolving convective simulations with knowledge about energy balances in the relaxed solution to instantaneously and self-consistently adjust the mean vertical thermodynamic profile toward its relaxed state. While a technique of this kind has been used by many studies previously (e.g., [18]), the details of implementation, the convergence properties, and whether or not the thermally relaxed state achieved from accelerated evolution corresponds to the relaxed state in standard evolution are not documented. In Sec. II, we describe our convective simulations and numerical methods. In Sec. III, we describe the procedure for achieving accelerated evolution. In Sec. IV, we compare accelerated evolution solutions to solutions obtained from standard evolution through a full thermal diffusion timescale. In Sec. V, we compare the numerical cost of accelerated solutions to standard solutions. Finally, in Sec. VI, we offer concluding remarks and discuss extensions of the methods presented here.

II. EXPERIMENT
In this work we study a simple form of thermal convection: incompressible Rayleigh-Bénard convection under the Oberbeck-Boussinesq approximation, such that the fluid has a constant kinematic viscosity (ν), thermal diffusivity (κ), and coefficient of thermal expansion (α). The density of the fluid is a constant, ρ 0 , except in the buoyancy term, where it is ρ = ρ 0 (1 − αT 1 ). The gravitational acceleration, g = −gẑ, is constant. The equations of motion are [19]: where u = ux + vŷ + wẑ is the velocity, T = T 0 (z) + T 1 (x, y, z, t) are the initial and fluctuating components of temperature, and P is the kinematic pressure. The initial temperature profile, T 0 , decreases linearly with height. We non-dimensionalize these equations such that length is in units of the layer height (L z ), temperature is in units of the initial temperature jump across the layer (∆T 0 = L z ∂ z T 0 ), and velocity is in units of the freefall velocity (v ff = αgL 2 z ∂ z T 0 ). By these choices, one time unit is a freefall time (t ff = L z /v ff ). We introduce a reduced kinematic pressure, ≡ (P/ρ 0 + φ + |u| 2 /2)/v 2 ff , where the gravitational potential, φ, is defined such that g = −∇φ. In non-dimensional form, and substituting u · ∇u = ∇(|u| 2 /2) − u × (∇ × u) and ∇ 2 u = −∇ × (∇ × u), Eqs. (3) and (4) become where ω = ∇ × u is the vorticity. The dimensionless control parameters R and P are set by the Rayleigh (Ra) and Prandtl (Pr) numbers, We hold Pr = 1 constant throughout this work, such that P = R. P and R are related to the inverse Reynolds and Peclet numbers of the evolved flows. In Eqs. (2), (5), and (6), linear terms are grouped on the left-hand side of the equations, while nonlinear terms are found on the right-hand side. We timestep linear terms implicitly, and nonlinear terms explicitly. We utilize the Dedalus 1 pseudospectral framework [20] to evolve Eqs. (2), (5), and (6) forward in time using an implicit-explicit (IMEX), third-order, four-stage Runge-Kutta timestepping scheme RK443 [21]. The code used to run the simulations in this work is included in the supplemental materials as a zip file [22].
Variables are time-evolved on a dealiased Chebyshev (vertical) and Fourier (horizontal, periodic) domain in which the physical grid dimensions are 3/2 the size of the coefficient grid. We study two-(2D) and three-dimensional (3D) convection in which the domain is a cartesian box, whose dimensionless vertical extent is z ∈ [0, 1], and which is horizontally periodic with an extent of x, y ∈ [0, Γ], where Γ = 2 is the aspect ratio, as has been previously studied [10,23]. In 2D simulations, we set v = ∂ y = 0. We specify no-slip, impenetrable boundary conditions at both the top and bottom boundaries, The temperature is fixed at the top boundary, and the flux is fixed at the bottom boundary, such that For this choice of boundary conditions, the critical value of Ra at which the onset of convection occurs is Ra crit = 1295.78 [24], and the supercriticality of a run is defined as S ≡ Ra/Ra crit . Studies of convection which aim to model astrophysical systems such as stars often employ mixed thermal boundary conditions, as we do here [12,25,26]. Our choice of the thermal boundary conditions in Eqn. (9) was motivated by the fact that accelerated evolution is simpler when both the thermal profile and the flux through the domain are fixed at a boundary (see Sec. III). The initial temperature profile is linearly unstable, T 0 (z) = 0.5 − z. On top of this profile, we fill T 1 with random white noise whose magnitude is 10 −6 P, and which is vertically tapered so as to match the thermal boundary conditions. This ensures that the initial perturbations are much smaller than the evolved convective temperature perturbations, even at large Ra. We filter this noise spectrum in coefficient space, such that only the lower 25% of the coefficients have power; this low-pass filter is used to avoid populating the highest wavenumbers with noise in order to improve the stability of our spectral timestepping methods.

III. THE METHOD OF ACCELERATED EVOLUTION
Here we describe a method of Accelerated Evolution (AE), which we use to rapidly relax the thermodynamic state of convective simulations. We compare this AE method to Standard Evolution (SE), in which we evolve the atmosphere from noise initial conditions for one thermal diffusion time, t κ = P −1 . Both AE and SE simulations begin with identical initial conditions, as described in Sec. II. As Ra increases, and P decreases, SE solutions become intractable, while the timeframe of convergence for an AE solution remains nearly constant in simulation freefall time units (see table II in appendix A).
We study in depth a 2D simulation at S = 10 5 to demonstrate the power of AE. We compare kinetic energy (KE, black line) and mean temperature (blue line) traces from a SE run in Fig. 1(a) to an AE run in Fig. 1(c). In Fig. 1(a), the time evolution of the SE simulation is shown. The KE grows exponentially from white noise during the first ∼ 25 t ff . The solution then saturates and begins to slowly relax toward the saturated isothermal profile in the interior of the domain. This slow relaxation is evident in the behavior of the blue line, which measures T − T top , where T is the volume-average of T , and T top = −0.5 is the temperature at the upper boundary. The mean atmospheric temperature and kinetic energies are fully converged when t = 4000t ff = 0.35t κ . We show roughly the first thousand freefall time units of evolution, as well as the evolved thermodynamic state reached after a full thermal time of evolution. In Fig. 1(c), similar traces are shown for the corresponding AE solution at the same parameters. The same linear growth phase occurs, but shortly after the peak of convective transient we accelerate the convergence of the atmosphere through the process which we describe below. We adjust the 1D vertical profile of the atmosphere three times, as denoted by the three labeled arrows in the graph numbered 1-3. The third profile adjustment associated with arrow 3 is small enough (see appendix B) that we assume the atmosphere is sufficiently converged, and we begin to sample the evolved convective dynamics.
The horizontally averaged profiles of the vertical conductive flux, F κ = −κ∂ z T x,y , and the vertical convective enthalpy flux, F E = wT x,y , are the basis of the AE method. Here we use x,y to represent a horizontal average. We measure both of these fluxes early in a simulation, retrieving profiles such as those shown in Fig. 1(b). As the atmosphere relaxes towards the isothermal profile specified by the upper (cold) boundary condition, excess temperature throughout the atmosphere must leave the domain. This excess thermal energy leaves through the upper boundary, as seen in Fig. 1(b), where the amount of flux exiting at the top of the domain is nearly 20 times larger the flux entering the bottom of the domain. Once the atmospheric temperature profile reaches its evolved state, the flux entering the bottom boundary is equal to the flux exiting through the upper boundary. In general, this evolution is slow in SE [ Fig. 1(a)], but AE [ Fig. 1(c)] can rapidly advance a system whose fluxes are in a strongly disequilibrium state [ Fig. 1(b)], into a near-equilibrium state, as shown in Fig. 1(d). In this final state both boundaries conduct the same amount of flux. The converged fluxes achieved through AE are at most 5% different from the SE solution, as shown in Fig. 1(e). This is a very small difference considering the short timescales on which convergence is reached and the strongly disequilibrium state used to inform the AE process.
In order to adjust the temperature profile to achieve AE, we calculate the total flux, F tot = F E + F κ , and then derive the profiles which have the systematic asymmetries [ Fig. 1(b)] removed. These profiles describe which parts of the atmosphere depend on convection to carry flux (where f E (z) = 1 and f κ (z) = 0). We presume that the early convection occupies roughly the same volume as the evolved convection, and thus that the extent of the early thermal boundary layers (where f κ (z) = 1 and f E (z) = 0) will not change significantly over the course of the atmosphere's evolution. Under this assumption, in order to reach the converged state, the flux through the atmosphere must be decreased by some amount, where F B = P is the amount of flux that enters the atmosphere at the bottom, fixed-flux boundary. For example, in Fig. 1b, F tot ≈ 19P at the upper boundary, but in the relaxed state it should just be P, so ξ ≈ 1/19 at that depth.
To reduce the conductive flux by ξ, we examine the horizontally-and time-averaged Eqs. (5) and (6) in the timestationary state; after neglecting terms which vanish due to symmetry, these equations become Here, we construct F E, ev = ξF E and u × ω x,y, ev = ξ u × ω x,y . Using these profiles and Eqs. (12) and (13), we solve for x,y and T 1 x,y . Our choice of fixing the temperature at the top boundary and the flux at the bottom boundary ensures that there is a unique solution for T 1 x,y given the evolved fluxes. If flux were fixed at both boundaries, there would be infinitely many temperature solutions. If temperature were fixed at both boundaries, F B would not be precisely known a priori.
Solving the above equations adjusts the mean thermal state and F κ of the atmosphere while leaving F E unchanged. In the bulk of the domain, convective enthalpy flux dominates transport, and so there we assume that F tot ≈ F E . Upon removing terms which vanish due to symmetry, we note that the convective enthalpy flux is carried solely by the velocity field and perturbations in temperature away from the mean state, F E = w(T 1 − T 1 ) . In order to instantaneously decrease the convective enthalpy flux in a manner which is consistent with the conductive flux, we multiply both the velocity, u, and temperature perturbations about the mean state, T 1 − T 1 x,y , by √ ξ. This scaling of the velocity field is the reason that we multiply u × ω x,y , which is nondimensionally of order u 2 /L, by ξ while constructing u × ω x,y, ev .
In general, the AE method occurs in the following steps. First, we pause a convective simulation and construct ξ(z, t) from measured flux profiles in the convective domain. We solve a 1D boundary value problem (BVP) consisting of Eqs. (12) and (13) to obtain an evolved thermodynamic profile, and its corresponding conductive flux. We multiply both the temperature perturbations around the mean and the convective velocity flows by √ ξ. After adjusing the fields of a simulation in this manner, we continue timestepping forward. For specifics on the precise implementation of the AE method, we refer the reader to appendix B.

IV. RESULTS
We study evolved standard evolution (SE) solutions whose supercriticalities (S) are S ∈ (1, 10 5 ] in 2D and S ∈ (1, 10 4 ] in 3D. We compare their properties to accelerated evolution (AE) runs at S ∈ (1, 10 7 ] in 2D and S ∈ (1, 10 4 ] in 3D. We refer the reader to appendix A for a full list of simulations. The Nusselt number (Nu) quantifies the efficiency of convective heat transport and is defined as where the volume average of a quantity, η, is shown as η . The time evolution of Nu in AE and SE is compared to the mean temperature evolution in Fig. 2. In Fig. 2(a), we show the evolution of the SE run at S = 10 5 . The black trace shows the instantaneous value of Nu, and the overplotted gray line shows a moving time average of Nu. The time average is taken using a centered boxcar window whose width is 50 t ff . The mean value of Nu evolves towards its relaxed value (dashed horizontal gray line) rapidly compared to T − T top (blue), which approaches its evolved state (dashed horizontal blue line) slowly. As such, it is possible to measure the statistically stationary mean value of Nu after only a few hundred simulated freefall times, as has been done previously [27]. However, fluctuations in Nu around the mean value at early times have both a larger magnitude and higher frequency than the fluctuations in the relaxed state. In Fig. 2(b), we show the AE traces of Nu and T − T top at S = 10 5 . The times at which AE solves occur are marked by arrows. Nu reaches its evolved mean value slightly more rapidly than in the SE case, and the frequency and magnitude of fluctuations in Nu away from the mean resemble the final relaxed state of SE. When S < 10 3.67 in 2D and for all runs in 3D, the evolved system is characterized by a time-stationary value of Nu, and is thus in a state of constant convective heat transport. At larger S in 2D, the value of Nu varies significantly over time even in the relaxed state (as seen in Fig. 2). We examine the shaded region of Fig. 2(b) in more detail in the top left panel of Fig. 3, as well as a comparable time trace at S = 10 7 (bottom, left panel). We find that these systems exhibit large Nu during states in which temperature fluctuations travel in their natural buoyant direction (Fig. 3, Ia and IIa, where cold elements fall and hot elements rise). However, when wrongly-signed temperature perturbations are entrained in an upflow or downflow with oppositely signed fluid, Nu is suppressed (Fig. 3 plumes in these systems naturally oscillate horizontally over time, switching between transport being dominated by a counterclockwise cell, as pictured in Fig. 3, and a clockwise cell. Our choice of no-slip boundary conditions prevents the fluid from entering a full domain shearing state [23], and the oscillatory motions of the plumes are a long-lived, stable phenomenon. However, thanks to our choice of periodic boundary conditions and despite the no-slip conditions, the full system of the upflow and downflow plumes is free to slowly migrate to the left or right over time, and we observe this phenomenon in our solutions. The 2D SE simulations exhibit the same horizontally oscillatory behavior as the AE solutions for the same initial conditions. This time-dependent behavior of Nu is not seen strongly in our 3D solutions, however most 3D simulations we conducted were at low S compared to the runs in which this behavior was observed in 2D. The time-and volume-averaged values of Nu, the RMS Peclet number (Pe), and the mean temperature are shown for AE solutions in Figs. 4(a)-4(c). Mean values are shown by the symbols (purple circles and red stars), and the vertical lines represent the standard deviation of the measurement over time. Nu is shown as a function of Ra and S in Fig. 4a, while Pe = |u| /P is shown in Fig. 4(b), and T − T top (the mean temperature value minus its value at the upper, fixed temperature boundary) is shown in Fig. 4(c). We report Nu ∝ Ra 1/5 , Pe ∝ Ra 0.45 , and ( T − T top ) ∝ Ra −1/5 . The average temperature, T , is dominated by its value in the isothermal interior, so this measurement serves as a probe of the temperature jump across the boundary layers. Thus, the inverse scaling of average temperature and Nu that we find here is expected for thermally equilibrated solutions [28].
In Figs. 4(d)-4(f), we report the fractional difference between measurements in the AE and SE solutions. The mean values of Nu and T − T top measured in AE are accurate to SE values to within ∼ 1%. Pe measurements show marginally greater error, with AE measurements being ≤ 2% different from SE measurements.
For the select 3D runs conducted in this study, the scaling of Nu, Pe, and T − T top reported in Figs. 4(a)-4(c) is nearly identical to the 2D simulations. Errors between AE and SE solutions in 3D fall within the same range as errors in 2D in Figs. 4(d)-4(f). AE is therefore equally effective in both 2D and 3D, and we restrict much of our study to 2D here in order to more thoroughly sample parameter space.
The measurements presented in Fig. 4 demonstrate that AE can be powerfully employed in parameter space studies in which large numbers of simulations are compared in a volume-averaged sense. We now turn our examination to a more direct comparison of AE and SE for 2D convection at S = 10 5 , the time, flux, and Nu evolution of which are shown in Figs. 1, 2, and 3. All comparisons that follow for these two runs occur over the times shaded in green in Figs. 1(a) and 1(c). Measurements are sampled every 0.1 freefall time units for 500 total freefall time units. As AE is fundamentally a 1D adjustment to the thermodynamic structure of the solution, we compare the horizontally-and time-averaged temperature profiles attained by AE and SE in Fig. 5(a). The boundary layer width and structure are nearly identical between the two solutions, but the the mean temperature in the isothermal interior differs by roughly 0.5% [ Fig. 5(c)].
The probability distribution functions (PDFs) of point-by-point temperature measurements are compared for the two runs in Fig. 5(b). To construct these PDFs, we interpolate the full temperature field at each measurement time onto an evenly spaced grid, determine the frequency distribution of all T values over the duration of the 500 t ff measurement window, and then normalize the distribution such that its integral is unity. The two PDFs have noticeably different maxima, as is expected from Fig. 5(a). Over long timescales, the 0.5% difference between the two profiles would disappear, as the AE solution evolves to be exactly the SE solution; this is evident in the asymmetry of the AE PDF near the maxima in Fig. 5(b) and also the trend of the mean temperature over time in Fig. 1(c).
One method of comparing two PDFs to determine if they are drawn from the same underlying sample distribution is through the use of a Kolmogorov-Smirnov (D KS ) test [29]. We calculate the D KS statistic for a PDF of some value, q, as where F stands for cumulative distribution function (CDF), the integral of the PDF. A traditional Kolmogorov-Smirnov statistic is just a single value, D KS (q) = |D KS (q)| ∞ = max|D KS (q)|, and we use both the profile KS(q) and D KS (q) to gain insight into the likeness of two PDFs. We show D KS (T ) in Fig. 5(d), and the CDFs used to construct it overlay the PDFs in Fig. 5(b). Near the maxima of the temperature PDFs, D KS (T ) = 0.495, which is very large and implies that roughly half of all measurements in the AE case are at a lower T than those in the SE case. While this difference is significant, it is also expected from Fig. 5(a). Fortunately, D KS (T ) is very small away from the maxima, indicating that the temperature fluctuations off of the maxima, which are the primary drivers of convective transport, are nearly identical. In addition to adjusting the 1D thermal profile, the AE method also scales the simulation velocities and temperature fluctuations by √ ξ. In Fig. 6 we examine the velocities and heat transport found in the evolved states. Shown are the PDFs of vertical velocity [w, Fig. 6(a)], horizontal velocity [u, Fig. 6(b)], and the nonlinear vertical convective flux [w(T 1 − T 1 x,y ), Fig. 6(c)]. Each PDF here shows a strong peak near zero due to the no-slip, impenetrable velocity boundary conditions (Eq. (8)). The CDFs of each profile are overplotted, and corresponding KS profiles are shown in Figs. 6(d)-6(f). We report D KS (w) = 0.00615, D KS (u) = 0.0349, and D KS (w(T 1 − T 1 x,y )) = 0.0263. The difference in vertical velocity and heat transport between AE and SE is negligible, which is unsurprising in light of the Nu measurements of Figs. 4(a) and 4(d). This also confirms that the large D KS (T ) in Fig. 5(d) is not of concern, and that the AE run achieves the same relaxed convective solution as the SE run. We find that the difference in D KS (u), which consistently has more probability of flows moving left (in the -x direction), appears to be caused by a more prominent migration of the full roll system in the -x direction in the AE run than in the SE run. This migration does not appear to affect the vertical transport appreciably.

V. COMPUTATIONAL TIME-SAVINGS OF AE
Computational time-saving is the primary reason to use AE rather than evolving all solutions through SE. In table I, we compare cpu-hour cost for select 2D and 3D runs in which both AE and SE solutions were computed. Times reported for AE and SE runs only include the time required to reach a relaxed state, and do not include the time over which measurements were taken in that state (e.g., the green shaded regions of Figs. 1(a) and 1(c) are not included in the S = 10 5 times). All simulations were performed on Broadwell nodes on NASA's Pleiades supercomputer (Intel Xeon E5-2680v4 processors). The key metric which highlights the usefulness of AE is the number of cpu-hours used for the AE run divided by cpu-hours used for the SE run (t CPU,AE /t CPU,SE ). We see that at low resolution and low supercriticality, t CPU,AE /t CPU,SE > 1, and AE is not useful. However, as S grows, t CPU,AE /t CPU,SE shrinks. At the highest supercriticalities for which AE and SE were compared in this work, AE runs cost roughly an order of magnitude less computing time than SE runs.
Integrating information about the mean state in time (fluxes, etc.) decreases the rate at which our solver timesteps early in the AE cases. However, the first application of AE in a given simulation [e.g., Fig. 1(c), at the arrow labeled "1"] drastically increases the average timestep by fastforwarding the simulation into a more stable state with lower convective velocities. For the S = 10 5 case we examined in detail, the average time step grew by a factor of 2-3. At S = 10 7 , the AE solve immediately improved the timestep size by nearly a factor of 4. In this work we have studied a method of Accelerated Evolution (AE) which can be employed to achieve rapid thermal relaxation of convective simulations. We compared this technique to the Standard Evolution (SE) of convection through a full thermal diffusion timescale, and we showed that AE rapidly obtains solutions whose dynamics are statistically similar to SE solutions. The AE method is valid at low values of S, where SE solutions converge quickly due to the short thermal timescale, and AE remains applicable at high values of S, where SE solutions are intractable. As discussed, AE is equally applicable in 2D and 3D; here we have restricted most of our study to 2D to extend our parameter space coverage. At the largest values of S in which AE and SE are compared in this work, we find time savings of nearly an order of magnitude.
Here we studied the simplest possible case for the application of AE: Rayleigh-Bénard convection at low aspect ratio with mixed thermal boundary conditions. We anticipate that this technique will be powerful in its extensions to more complicated studies. To achieve AE in more complicated systems, one need only derive the steady-state, horizontallyaveraged equations governing the convective dynamics [e.g., the analogs to Eqs. (12) and (13)] and couple those equations with knowledge of the boundary conditions and current dynamics as described in Sec. III and appendix B. In general, AE should be useful in studies where there are two disparate timescales which must both be resolved and which cannot be overcome through clever timestepping techniques. Some avenues in which extensions of AE could be beneficial for expanding the available parameter space of exploration include studies of internally heated convection [24], convection with height-dependent conductivities [30], penetrative convection [13,14,18], or fully compressible, stratified convection [7]. As AE is fundamentally a horizontally uniform adjustment to the thermodynamic structure of the convective domain, it is unlikely that these techniques should be applied straightforwardly to nonperiodic convective domains.
We conclude by noting that AE should be extended to these more complicated studies with caution. While AE was extremely effective in this simple case studied here (where the aspect ratio was low, the bounds of the convective domain were pre-defined, and the solutions were simple rolls), this may not be the case for more complicated systems. For example, at higher aspect ratios, multiple stable solution branches may exist, and there is no guarantee that AE and SE will arrive at the same solution. Some assumptions which inform the AE solution, such as the assumption that the convection intially occupies the same space as the evolved convection, may not hold in studies of penetrative convection, despite the fact that similar methods have long been used in those studies [18]. Extensions to fully compressible convection in which there are two true thermodynamic variables [7] must contain a very careful treatment of AE pressure adjustments, so as to avoid wave-launching pressure mismatches. Our work here serves as a basis for determining if AE techniques are effective in more complex studies of convection. II. Simulation parameters. We report the supercriticality (S), Rayleigh number (Ra), and coefficient resolution (nz, nx, and ny are the number of coefficients in the z, x, and y directions respectively). Simulation run times required to reach convergence are reported for the SE solutions (tSE) and the AE solutions (tAE). The amount of time over which simulations measurements were taken in the evolved state is listed (tavg). All times are in freefall time units. The volume-averaged Nusselt number (Nu) of the AE and SE solutions are shown. In the upper part of the table, information pertaining to 2D runs is reported, while information pertaining to 3D runs is in the lower part of the 6. Multiply the velocity field, u = ux + vŷ + wẑ, and the temperature fluctuations, T − T x,y , by √ ξ in the DNS to properly reduce the convective flux.

Continue running the DNS.
We refer to this process as an "AE BVP solve." While the use of a single AE BVP solve rapidly advances the convecting state to one that is closer to the relaxed state, we find that repeating this method multiple times is the best way to ensure that the AE solution is truly converged. For all runs in 2D at S < 10 5 , we set t transient = 50, completed an AE BVP solve with t min = 30 and P = 0.1, and then repeated the procedure. For all 3D runs and 2D runs with S ∈ [10 5 , 10 6 ], we did a first AE BVP solve with t transient = 20, t min = 20, and P = 1 in order to quickly reach a near-converged state and vastly increase our timestep size. After this first solve, we completed two AE BVP solves, with t transient = 30, t min = 30, and P = 0.1 to get very close to the solution (as in Fig. 1c). At very high S = 10 7 , we ran two AE BVP solves with t min = 20 and P = 1. For the first solve, we set t transient = 20, and for the second we set t transient = 30. We used fewer solves at this high value of S in part to reduce the computational expense of the run, and in part because a third BVP generally did not greatly alter the solution (as in Fig. 1c, arrow 3). We wait 50 freefall times after the final AE BVP solve of each run before beginning to take measurements.
In general, to use AE, a threshold, f , should be chosen. When the fractional change of the mean temperature profile from an AE BVP solve becomes less than f , the solution can be considered converged on its chosen solution branch. In other words, once the solution is converged. In this work, we chose our number of AE iterations such that f ≈ 10 −2 . In general, the user of AE could set f smaller, and in doing so reduce the separation between the AE and SE solutions, which can be seen in e.g., Fig. 5a&b. However, smaller values of f require additional AE BVP solves, and likely require smaller values of P , resulting in longer wait times while the horizontal averages are computed.