Adapting a Fourier pseudospectral method to Dirichlet boundary conditions for Rayleigh--B\'enard convection

We present the adaptation to non--free boundary conditions of a pseudospectral method based on the (complex) Fourier transform. The method is applied to the numerical integration of the Oberbeck--Boussinesq equations in a Rayleigh--B\'enard cell with no-slip boundary conditions for velocity and Dirichlet boundary conditions for temperature. We show the first results of a 2D numerical simulation of dry air convection at high Rayleigh number ($R\sim10^9$). These results are the basis for the later study, by the same method, of wet convection in a solar still.


I. Introduction
Experimental observations [1] show that the onset of a turbulent convective flux can significantly enhance the efficiency of a basin-type solar still, but until now a theoretical explanation is lacking. Any adequate hydrodynamical simulation must incorporate the effects of moisture and condensation. Recent works [2,3] show that this can be achieved through a Boussinesq-like approximation, which simplifies considerably the problem. However, realistic simulations are still demanding, given the need to resolve fine flux details and cope with Rayleigh numbers up to ∼ 10 9 [4].
Spectral methods [5] are well suited for this kind of tasks, and have many attractive features: they are simple to implement, show much better resolution and accuracy properties than finite difference or finite volume methods, and are highly efficient in large-scale simulations [6]. Fourier-based pseu-dospectral methods are the simplest and fastest, since the discretized spatial differential operators are local, nonlinear terms can be computed through Fast Fourier Transform (FFT) convolutions, and solving the Poisson equation originating from the incompressibility (divergence-free) condition is almost trivial. Nevertheless, they usually work only for free (in fact, periodic) boundary conditions (BCs).
The presence of non-free (e.g., Dirichlet or Neumann) BCs introduces additional complications. For example, two-dimensional Rayleigh-Bénard convection with laterally periodic BCs can be treated by using a spectral Galerkin-Fourier technique in the horizontal coordinate and a collocation-Chebyshev method in the vertical one [7], but vertical derivatives must then be computed by matrix multiplication. On a grid with N horizontal and M vertical points, this needs the solution of a linear system of dimension M for each of the N horizontal Fourier modes, at each timeintegration step.
Another complication with non-free BCs arises from the need to fulfill numerically the divergencefree condition, leading mainly to two different groups of methods (see Ref. [6] and references therein). In a first group the velocity field is written in terms of scalar potentials such that the divergence-free condition is satisfied by construction, e.g., in the 2D streamfunction-vorticity formulation or the 3D decomposition into toroidal and poloidal velocity potentials. In these methods, pressure is not present in the equations, but they lead to systems of higher-order partial differential equations with coupled BCs. In a second group, a primitive variable formulation of the equations is adopted and projection methods [8] are used to decouple velocity and pressure. These methods use a specific splitting of the equation system based on the chosen time-integration scheme, and determine pressure by projecting an appropriate velocity field onto a divergence-free space, leading to predictor-corrector algorithms. Besides the problem of correctly specifying the pressure BCs [9,10], these methods require solving a Poisson equation for the pressure at each time-integration step. On a N = N × M grid, the best Fourierbased Fast Poisson Solvers (FPS) have operation counts O(N log 2 N ) for the lowest (second) order discretization (and significantly worse for higher orders) [11,12], and those using GMRES are KO(N ) with K 100 [13][14][15].
In this work, we will show how a Fourier-based pseudospectral method can be adapted to simple non-free (but periodic) BCs without losing its more appealing features. This is a first step towards building a pseudospectral simulation of wet air convection inside a basin-type solar still, and must be considered just as a proof of concept.

II. System
We consider 2D dry air convection in a Rayleigh-Bénard cell of width L = 1 m and height H = 0.5 m, close to room temperature and with temperature differences ∆T = T h − T c up to ∼ 65K between the hot lower (T = T h ) and cold upper (T = T c ) plates (roughly the parameters of a real still [1]). Discarding thermal fluctuations and the heat generated by viscous dissipation, and assuming an incompressible fluid across which all thermodynamical parameters change little, the dynamics is given by [16] In these equations, the dynamical variables are the density ρ, the velocity u, the temperature T , and the pressure P . The parameters are the shear (or dynamic) viscosity η, the constant pressure specific heat capacity c p , the thermal conductivity K, and the gravitatioinal acceleration g. The Boussinesq approximation [17] consists in discarding the dependence of η, c p , and K on temperature and density, keeping in the buoyancy term (−ρgẑ in Eq. 1) but otherwise setting ρ =ρ elsewhere. Here,T = 1 2 (T h + T c ) is a reference temperature,ρ is a reference density (that of air at normal temperature and pressure), and α is the thermal expansion coefficient. Dropping the bars signifying reference quantities, absorbing some constants into the pressure gradient term, and defining the viscous diffusivity (or kinematic viscosity) ν = η/ρ and the thermal diffusivity κ = K/(ρc p ), we obtain the (dimensional) Oberbeck-Boussinesq equations [17] (∂ t + u · ∇)u = −∇(P/ρ) + ν∇ 2 u − αgTẑ, (5) Assuming perfect thermal contact with the lower (z = 0) and upper (z = H) plates, these equations admit the stationary conductive solution whereP is a reference (e.g. normal atmospheric) pressure, and z ′ = z − H/2. We now scale lengths with the cell height H, times with the characteristic vertical thermal diffusion time t c = H 2 /κ, and temperatures with the temperature difference ∆T . The nondimensional lengths, times, and velocities are then We also define the nondimensional temperature θ and pressure P ′ by Substituting into Eqs. (5)- (7), discarding the primes for simplicity, and absorbing all pure gradient terms into the pressure, we get the dimensionless Oberbeck-Boussinesq equations [17] σ where is the Rayleigh number and σ = ν/κ is Prandtl's number (≃ 0.7 for dry air). The BCs we adopt are periodic in the horizontal direction, and homogeneous Dirichlet for both velocity u (no-slip BCs) and temperature θ (perfect thermal contact) on the lower and upper plates. Note that Eq. (13) is not a differential equation but a constitutive relationship, expressing the incompressibility of the flux. In fact, the pressure term in Eq. (11) is computed by enforcing Eq.(13), which gives where i, j = x, z. In primitive variable integration schemes, this Poisson equation must be solved with adequate BCs at each time-step [6,9], to insure ∇ · (∂ t u) = 0.

III. Helmholtz decomposition
Given a vector field f , twice continuously differentiable, Helmholtz's Theorem [18] states that it can be decomposed as Here, f and f ⊥ are the longitudinal (or irrotational) and the transverse (or solenoidal) components of the field [8], respectively, with We are now going to rewrite Eq. (11) in the form Then, if ∇ · u = 0 initially, the incompressibility condition Eq. (13) requires This amounts to requiring the field f − ∇P to be purely transverse, that is where we used (∇P ) ≡ ∇P , since the pressure gradient is purely longitudinal. This shows that the only effect of the term ∇P in Eq. (18) is to cancel the longitudinal component of f . Equation (18) can then be set as with no explicit reference to a pressure field. In 2D and in free space (that is, disregarding surface terms), the longitudinal and transverse components of f can be computed as the projections of its Fourier transformf along the unit vectorŝ However, on a finite domain, the surface terms cannot be ignored, since they are essential for f ⊥ to have the correct BCs, which by Eq. (22) are the same as those for u in Eq. (11). Using Eq. (24), the fieldf ⊥ in Eq. (23) can be seen to be a particular solution (the free-space solution) of the Poisson equation where F −1 stands for the inverse Fourier transform.
To be able to impose the required BCs, we need the general solution of this equation, which we can get by adding the general solution of the corresponding homogeneous equation The required transverse component of f can then be redefined as and where the last equation is needed to insure the transversality of w and hence of f ⊥ . This requirement can be automatically fulfilled by writing w explicitly asw with the scalar field w satisfying where c is a constant. Noting that Eq. (27) can also be rewritten as we can also rewrite Eq. (26) as which shows explicitly the transversality of f ⊥ . The determination of the value of c, and the treatment of possible divergences inṽ at k = 0, are closely related and will be dealt with in the next section. The BCs for w at the lower and upper plates can be obtained from those for f ⊥ , and are w = −v; the BCs on the horizontal direction are periodic but otherwise free, and will be automatically fulfilled by the constructive procedure for w given in the next section.
IV. Ultra-fast Laplace solver with c p , a, and b (possibly complex) constants. For convenience, and without loss of generality, we will rewrite it in the form where a p and b p are constants, and z ′ = z − H/2. Here, we have used that cosh(0) = 1 to absorb all constant terms in a 0 , and used that sinh(0) = 0 to absorb the linear term in z ′ as the particular term of the sum for p = 0, leaving Eq.
give the discretization of w on the wavenumber grid. The matricesC pq andS pq are given explicitly bỹ Then from Eqs. (31) and (32), the discretization of f ⊥ on the wavenumber grid will be expressed asf The constant c can now be formally chosen to cancel the possible divergence at k = 0; in practice, we set c = 0 and redefineṽ 00 to be an arbitrary but finite constant, which without loss of generality can also be taken to vanish. Imposing the BCs w = −v, or equivalently f ⊥ = 0, at z = 0 and z = H, is achieved as follows: First we note that any field discretized on the wavenumber grid through a DFT of its discretization on the coordinate grid is automatically periodic in both the horizontal and vertical directions, so the BCs at z = 0 and z = H will give identical sets of equations. Next we use the fact that together with the completitude of the DFT, to rewrite the BC at z = 0 as qf ⊥,pq = 0 ∀p.
We then use Eq. (40) and the parity ofC pq (even in q) andS pq (odd in q) to obtain the conditions from which the coefficients a p and b p can be immediately retrieved. Finally, substituting into Eq. (40) leads tõ where it is understood that we takeṽ 00 = 0, and we have introduced the normalized matrices Equations (44) and (31) show that the transverse fieldf ⊥,pq on the wavenumber grid can be obtained directly in terms of the non-transverse field f pq without needing, at any point, to return to the coordinate grid. Moreover, Eq. (45) shows that c pq and s pq are given matrices that can be computed just once at the start of the simulation, as is the denominator k 2 x,p + k 2 z,q in Eq. (31). The only tasks to be performed at each time-step are then: first, computing the scalarsṽ pq fromf pq , requiring three multiplications per grid point; second, computing the sums in Eq. (44), which requires one multiplication per grid point; third, multiplying them by c pq and s pq , costing two multiplications per grid point; and fourth, multiplying by k z,q and k x,p at each grid point. The total cost of obtainingf ⊥,pq is then 8N multiplications on a N = N ×M grid, thus outperforming even the best FPS by a significant factor on large grids.
It must be noted that the method introduced here has some similarity to the streamfunctionvorticity formulation [6], in the sense that the scalar fieldsṽ andw play a role similar to these potentials. However, in our method they are not taken as dynamical variables, and the evolution equations are not formulated in terms of them but of primitive variables. The method presented here shows also a strong similarity with projection methods [8,9], but differently from them, pressure is not computed along the time evolution and, in fact, does no longer appear in the evolution equations.

V. Algorithm outline
We present now an outline of the numeric algorithm as we implemented it in the simulations. Initialization: • Take zero velocity and Gaussian white noise for the temperature (with amplitude of the order of thermal noise), discretized on the coordinate grid, and take their FFT to getũ x,pq , u z,pq ,θ pq on the wavenumber grid.
• Compute r.h.s. of evolution equations as Here, we have denoted by (f * g) pq the convolution product of fieldsf pq andg pq discretized on the wavenumber grid, which is performed by FFT.
It must be noted that the spatial discretization of the evolution equations has been performed in a closed form, independent of the time-stepping algorithm to be employed to solve the resulting set of ordinary differential equations, outlined above. We must also note that this system does not involve multiplication by matrices with dimension N , the only nonlocal parts being the convolution products (handled through FFT) and the elimination of the longitudinal component of the velocity.
The time-stepping can then be performed by any algorithm designed to solve systems of ordinary differential equations. In our case, we opted for an adaptive-stepsize fifth order Runge-Kutta-Cash-Karp algorithm [19], which in our previous experience we have found efficient, stable, and flexible.

VI. Test runs
Even for a simple system like the one we are studying here, the phenomenology found is rich; we  At low R, a stationary regime state, consisting in two counter-rotating rolls, is reached in times ∼ t c or less. This time falls rapidly with increasing R, to ∼ 0.1t c at R ∼ 1000R c (see Figs. 1, 2 and 3).
Around R ∼ 5000R c , these rolls develop lateral oscillations, and the first "secondary structures" (small whirlpools) appear near the base of the ascending and descending plumes (see Fig. 4).
Above R ∼ 5000R c , the regime state becomes disordered and aperiodic, consisting of intermittent plumes and whirlpools in a wide size range (see Fig.  5).
At R ∼ 5×10 5 R c , the temperature difference is ∼ 65K; the smallest whirlpools are ∼ 1cm wide, and the typical wind speeds are ∼ 1m/s, in agreement  with experimental observations [1] (see Fig. 6). The time to reach this regime state is rather short, ∼ 0.001t c .
Note that in all figures, except for Figs. 1 and 2, the velocity grid has been decimated to enhance clarity.

VII. Code performance
Over the full range of R tested here (more than five orders of magnitude), the code maintained the typical velocity divergence and the field values at the bottom and top boundaries at essentially machine-precision zero, showing that the implemented method is sound.
Also over this range of R, the grid spacing needed to achieve "smooth" fields (i.e., to capture all the physical detail down to the smallest present scales) Figure 5: Temperature and velocity fields for R = 5 × 10 4 R c at t = 0.01t c on a 384×192 grid (velocity decimated to a 64×32 grid).
is consistent with the width of the (thermal) boundary layer. However, for coarser grids the code still gives qualitatively sound results; typically a checkerboard-like instability develops, but the algorithm keeps it quenched, showing very good stability even in presence of a severe accuracy loss.
The algorithm is also fast: the simulation for R ∼ 10 9 (see Fig. 6), on a 512×256 grid, took less than one day per simulated minute on a single core of the 3GHz PentiumD CPU on which all our runs were performed, with no code optimization.
It is difficult to find in the literature a directly comparable simulation for more accurate benchmarking. However, from Ref. [6] we can see that, for example, the relaxation to a (steady) regime state in Rayleigh-Bénard convection at R 30000 ∼ 20R c on a ∼ 20000-node grid takes 46 minutes on a similar processor at similar speed (3.2GHz Pentium 4), by the method implemented there. In the case of our code, at R = 20R c on a 200 × 100 grid, and with the sole optimization of grouping the real FFTs in pairs (see the 2FFT algorithm in Ref. [19]), the equivalent relaxation took 8 minutes, which is shorter by a factor of ∼ 6. But it must be taken into account that with our initial conditions (zero velocities and thermal noise in temperatures) the convection onset is slow, and is followed by a transient stage with strong and disordered convective patterns that decay to the regime state very slowly. Figure 6: Temperature and velocity fields for R = 5×10 5 R c at t = 0.002t c on a 512×256 grid (velocity decimated to a 64×32 grid).

VIII. Conclusions and outlook
We have been able to show that a Fourier-based pseudospectral method can be adapted to a (admittedly simple) non-free BC setting, at the cost of moderate analytical work on the solutions of Poisson's and Laplace's equations. The method is formulated in primitive variables, but the pressure is not explicitly computed nor referenced, like in a streamfunction-vorticity formulation. It also shares some properties with projection methods, but it decouples the implementation of the incompressibility condition from the time-stepping scheme, allowing great flexibility in the selection of the last. The resulting code is fast and stable, and implements the BCs and the incompressibility condition essentially to machine-precision.
Work on the extension of this scheme to a fully closed Rayleigh-Bénard cell (i.e., with non-free BCs also on the lateral walls) is currently under course.