Oscillations on energy spectrum tail (high frequencies) for 2D Rayleigh-Benard convection

[tlunet]

Hi,

I'm currently using Dedalus to generate RBC data for a project, using different Rayleigh numbers, based on the example script on the website.

Problem : I have currently some weird results when extracting the energy spectrum from the generated data.

Just to drop it now, your code is great, I really enjoy working with it and it has inspired me a lot for many other projects ...

Also : this issue is linked to #291

Context : I have to be sure that for a given Rayleigh number, the space resolution is sufficient. So I use a similar approach as one I used in a previous work on DNS Navier-Stokes simulation (see Section 3.5). Idea is to set a constant mesh resolution, look at the energy spectrum for increasing Rayleigh numbers, up to the point where the tail of the spectrum (high frequency domain) goes "up", indicating an energy accumulation in the small scales, hence too coarse space resolution.

$\Rightarrow$ defines a "critical" Rayleigh number, corresponding to the maximum one we can simulate on this mesh size, considering a DNS-like resolution, hence a "correct" physical solution.

Methodology for energy spectrum

Based on the example script, I simulate up to 150 seconds, using a constant time-step size depending on the space resolution, ensuring no numerical instability from time integration. Base space resolution is $N_x, N_z = 256, 64$, and for this space resolution I use $\Delta t = 5e^{-3} sec$. Velocity fields is written into files using a file handler every 0.1 sec.

Then the energy spectrum is computed using an average approach using velocity fields from $t=50 sec$ to $t=150 sec$, dropping the initial transition phase.

Spectrum computation use RFFT for the $u_x$ and $u_z$ fields (np.fft.rfft) on the $x$ direction and multiply by the conjugate, average all those spectrum over the $z$ direction, then average over time. The final spectrum is computed by averaging spectrum of $u_x$ and spectrum of $u_z$, and normalizing by dividing by $N_x$.

Here is one spectrum example for the base resolution (in orange) and a finer resolution (in blue, 512 x 128) :

Looks consistent and physical, except for the high frequency tail, that is my focus : for the base space resolution (orange), I get this oscillation on the tail, that does not look physical even if the simulation seems to be correctly resolved, as the spectrum is similar for the finer resolution (so no energy accumulation in the small scales) :

Also, the finer resolution spectrum shows this high frequency oscillation too :

and furthermore :

• the last (high frequency) point of the spectrum is zero actually, so not show on the loglog plot
• it does not change if I reduce the time-step size, or change the time-integration method (from RK222 to RK443)
• the oscillation increases slightly with the Rayleigh number ... for instance, when reaching the critical Rayleigh number :

My problem (finally ...) : I want to script the detection of critical Rayleigh number by looking at when the tail goes up for the spectrum, which would be easy if the oscillation was not there. Because of the oscillation, I found a a workaround that uses the last part of the spectrum (1/4) without the last 3 points, fits it to a 2nd order polynomial, and check if the leading coefficient is negative (concave, DNS) or not (convex, under-resolved). This corresponds to the black dashed lines.

But this is not really satisfactory, as selecting eventually more or less point to throw changes a lot when the under-resolution is detected. Ideally, I would like to work with a spectrum that does not have those high frequency oscillations. So I wonder :

• is it due to the 2D nature of the problem (I haven't extended to 3D simulations yet, because of the computation cost).
• is it something from the computation of my spectrum ?
• is it something form the space discretization numerics in dedalus ?

Any hint or help on this would be awesome 😃

[kburns]

Hello, the mailing list is the right place for discussions like this, since we use the issue tracker to keep track of bugs in Dedalus. One issue may be how the spectrum is being computed in the first place -- the outputs of the Rayleigh Benard script are on a Chebyshev grid, so the signal is not periodic and the points are not evenly spaced, so putting this data into a 2D FFT does not correctly compute a 2D power spectrum. But please post on the mailing list if you'd like to continue discussing.