Issue 
A&A
Volume 578, June 2015



Article Number  A106  
Number of page(s)  15  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/201425125  
Published online  12 June 2015 
Online material
Appendix A: Numerical methods
Appendix A.1: Integration scheme
We have adapted the standard MacCormack method to suit our specific problem. In the simplest case of a onedimensional vector q of conserved quantities being advected on and equidistant grid with a spacing of Δx, MacCormack’s method can be written as where is the value of q at the kth grid point and the nth time step, Δt the time step, f(q) the flux function, and we use the convention that any parenthesised upper index refers to a substep of the method instead of a timestep index. The method is linearly stable provided that the CFL condition Δt ≤ Δx/ρ(A) is met, where A is the Jacobian matrix of the flux vector and ρ(A) is the largest characteristic value of A. Nonlinear stability typically requires the addition of some form of artificial viscosity. MacCormack’s method is secondorder accurate both in space and time.
We discretise Eqs. (4), (6), and (37) on a collocated, twodimensional grid of M × N cells with constant cell spacing (Δx, Δz). The two spatial dimensions and the presence of source terms in the equations forces us to significantly extend the basic MacCormack scheme. We begin by advecting the vector of variables q = (u,w,ϑ) in both spatial directions using Strang splitting, where we have written out the explicit form of the flux terms. The indices k and l refer to the position along the x and z axes, respectively. We proceed by adding the source terms to the momentum equations, where we use secondorderaccurate central differences to keep up with the order of accuracy of the advection scheme, ν_{l} = Pr_{x}(z_{l}) and μ_{l} = Pr_{z}(z_{l}) are the coefficients of our anisotropic artificialviscosity prescription (see Sect. 2.3), and μ_{l + 1/2} = (μ_{l} + μ_{l + 1})/2. The new velocity field is, in general, slightly divergent. We correct for this divergence by subtracting the gradient of a pressurecorrection field, u^{(1)} = u^{(1d)} − Δt∇(Δp)^{(1)}. The condition ∇·u^{(1)} = 0 leads to a Poisson equation for the pressure correction, (A.9)Since we use central differences to compute partial derivatives, the discrete form of the Laplace operator in Eq. (A.9) should be derived by applying the central differences twice. That would, however, lead to a sparse operator and cause oddevendecoupling problems on our collocated grid. Therefore we use the standard compact Laplacian and solve the approximate pressurecorrection equation (A.10)Equation (A.10) is solved by a spectral solver, see Sect. A.3. Having computed the pressure correction, we apply it to the velocity field, The approximate nature of the pressurecorrection equation (Eq. (A.10)) causes ∇·u^{(1)} to be small, but nonzero. Practical experience has shown that the residual divergence is negligibly small in the flows analysed in this paper provided that the boundary conditions are treated properly, see Sect. A.2. We should also write at this point, but our numerical tests have shown that the residual divergence in the velocity field becomes much smaller if we set , so we use the latter form. The next step is to integrate the remaining two terms in the energy equation. We begin by adding the −w term, (A.13)where its latest available value, −w^{(1)}, has been used. The diffusion substep is given by the implicit equation (A.14)which is also solved by a spectral solver, see Sect. A.3. We have thus completed the first step of the MacCormack scheme, analogous to Eq. (A.1), and obtained the new variables u^{(1)}, w^{(1)}, p^{(1)}, and ϑ^{(1)}. The second step, which we do not do not go into detail on, differs from the first one at two points. First, advection is done using backwardspace flux differencing, as in Eq. (A.2) (compare with Eq. (A.1)). Second, the pressure field is updated in this step, i.e. . The final step of the MacCormack’s scheme, Eq. (A.3), is used in the same form, with q = (u,w,ϑ). We also update the pressure field in the same way, , so that we obtain an estimate of the pressure field for the next time step.
Fig. A.1
Effect of three different methods of treating the pressure at the solid top and bottom boundaries. In all three panels, the vertical velocity component, w, is plotted on a split logarithmic colour scale. We use Θ = 10^{3}, L = 10^{1}, constant kinematic viscosity and set the resolution to only 16 × 64 to make the spurious oscillations visible. We obtain the result plotted in the left panel using the simple symmetry conditions for pressure (Eqs. (A.19) and (A.20)). Preceding the pressuregradient computation by thirdorder pressure extrapolation to the ghost cells reduces the oscillations’ amplitude by a factor of ~100 (middle panel). Increasing the extrapolation order to six brings about another decrease by a factor of ~30 in the oscillations’ amplitude (right panel). The pressure gradient is in all three cases computed by the secondorder central differences in the whole computational domain. 

Open with DEXTER 
Finally, there is a simple way of increasing the accuracy of the scheme at a given grid resolution, which we use. The MacCormack method contains a builtin asymmetry: Eqs. (A.1) and (A.2) show that it always starts with forwardspace flux differencing and continues with backwardspace flux differencing. The two fluxdifferencing methods can be reversed, obtaining a “reverse” MacCormack method, without decreasing the order of accuracy of the overall scheme. We compute every time step using both the “direct” and the “reverse” methods and use the arithmetic average of the estimates given by the two methods.
Appendix A.2: Boundary conditions
The treatment of boundaries is restricted by our decision to use spectral solvers, which do not allow changing the differentiation operators anywhere in the computation domain. We use the ghostcell technique for this reason. The boundary conditions we impose on the differentialheating flow are summarised in Sect. 2.3. The periodic boundaries in the horizontal direction are trivial to implement. The solid boundaries on the top and bottom of the computational domain, however, require much more care. We implement them using reflective boundary conditions for the velocity vector, so that the imaginary walls are located at l = −1/2 and at l = N − 1/2. The conditions imposed on u also eliminate any shear on the boundary. The pressure field is required to be symmetric with respect to the solid boundaries, The conditions imposed by Eqs. (A.15)–(A.20) can easily be shown to be consistent with the pressurecorrection equation (Eq. (A.10); sum both sides over k = 0, 1,...,M and l = 0, 1,...,N). They typically do, however, bring about a cusp in the pressure field along the normal to the walls. The resulting discontinuity in the vertical pressure gradient then propagates to the rest of the domain and can be seen as a lowamplitude oscillatory field superimposed on the true pressure field (see the left panel of Fig. A.1). We tried to cure this problem by changing the discretisation of the verticalgradient operator at the walls, so that the ghost cells would not be used when computing the pressure gradient. This solution has met with very little success, most likely because the abrupt change in the operator brings about an abrupt change in the discretisation error so the problem remains. Quite surprisingly, preceding the pressuregradient computation by highorder pressure extrapolation to the ghost cells has turned out to be an effective solution, able to eliminate nearly all of the spurious oscillations (see the middle and right panels of Fig. A.1). We therefore use sixthorder extrapolation in the simulations with constant artificial viscosity and increase the extrapolation order to ten when we let the artificial viscosity decrease with height. This technique cannot be viewed, however, as an allpurpose solution, because it is likely to be too unstable to be useful when computing highly turbulent flows.
We require the temperature fluctuation ϑ to have a fixed sinusoidal profile at the bottom boundary and to vanish at the upper boundary, which translates into
Appendix A.3: Spectral solvers
We use spectral methods to solve the two equations involving the Laplace operator, the Poisson equation for the pressurecorrection equation (Eq. (A.10)) and the implicit heatdiffusion equation (Eq. (A.14)). We express both the knowns and unknowns as linear combinations of the Laplacian’s eigenfunctions that comply with the desired boundary conditions. The solution procedure is then much simplified and effective, provided that the transform to the eigenfunction basis can be computed efficiently.
In case of the pressurecorrection equation (Eq. (A.10)), we use the linear transform (A.23)and its inverse (A.24)to transform any field f_{k,l} to an array of complex amplitudes and back. We can see that the basis functions in Eq. (A.24) are periodic in k and even around l = −1/2 and l = N − 1/2; i.e., they comply with our boundary conditions on the pressure field (see Sect. A.2). Upon using the spectral decomposition defined by Eq. (A.24) on both sides of the pressurecorrection equation (Eq. (A.10)), we readily obtain its solution in the wavenumber space, (A.25)where we have omitted the upper indices because the expression applies to both steps of the MacCormack scheme, Ŝ_{k,l} is the transformed righthand side of Eq. (A.10). The eigenvalues λ_{m,n} of the Laplacian are (A.26)and can be precomputed. We set λ_{0,0} to a large number to prevent division by zero and make the undetermined component vanish.
In case of the heatdiffusion equation (Eq. (A.14)), we use the linear transform (A.27)\newpage
and its inverse (A.28)to transform any field g_{k,l} to an array of complex amplitudes ĝ_{m,n} and back. We can see that the basis functions in Eq. (A.28) are periodic in k and odd around l = −1/2 and l = N − 1/2; i.e., they comply with our boundary conditions on the temperature field in case of a vanishing heating amplitude (see Sect. A.2). To allow for an arbitrary heating profile at the bottom boundary, we take out the known boundary term from the Laplacian on the righthand side of Eq. (A.14) and treat it as a source term. One can show that it is the same as replacing the diffusion equation ∂ϑ/∂t = ∇^{2}ϑ by the equivalent equation ∂(ϑ − ζ) /∂t = ∇^{2}(ϑ − ζ), where ζ is the static solution to the diffusion equation ∂ζ/∂t = ∇^{2}ζ with the desired boundary conditions (ζ can be precomputed for a fixed heating profile). The boundary conditions on the difference ϑ − ζ are then identically zero, and the spectral decomposition defined by Eq. (A.28) can be used. This way we obtain an explicit expression for the solution of the implicit Eq. (A.14) in the wavenumber space, (A.29)where the eigenvalues Λ_{m,n} of the Laplacian are (A.30)and can be precomputed. An equation analogous to Eq. (A.29) relates to .
In the practical implementation, we use the FFTW library (Frigo & Johnson 2005) to compute the transforms in Eqs. (A.23), (A.24), (A.27), and (A.28). We combine standard, onedimensional transforms of different kinds to obtain the nonstandard, twodimensional transforms that we need. Namely, Eq. (A.23) is implemented as a series of DCTII transforms over the rows of the input array, after which the columns of the resulting array are transformed by a series of DTF transforms. The backward transform (Eq. (A.24)) is then computed by a series of DFTs followed by a series of DCTIIIs. The transforms for the diffusion equation (Eqs. (A.27) and (A.28)) are implemented in the same way, but simply replacing the DCTIIs by DSTIIs and DCTIIIs by DSTIIIs. The transforms from the FFTW library do not include the normalisation factor (2MN)^{1}.
© ESO, 2015
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.