precision loss in non-trigonometric, periodic functions using FFTW and NaNs after marching forward in time (Fortran)

Question

I have developed a pseudospectral solver of the Navier-Stokes equations using FFTW. I tested my formulation of right hand sides (RHS) of the NS equations against standard trigonometric functions (sines, cosines and their combinations). For example, I set

density = sin 5x
x_velocity = 5cos 5y + 6sin7z
y_velocity = 4sin4y + cos x
z_velocity = 1
pressure = cos z

Supplying these values to the solver, it computed the RHS of the NS equations. I did the same by hand and compared the results with that obtained by the solver. Results were to good agreement. The maximum error between the exact answer and that computed by the solver was of the order of E-13 for a 128*128*128 grid.

Next I used a different function of the following form:

density = constant1+constant2*(tanh(x-constant3)-tanh(x-constant4))
x_velocity = 0
y_velocity = 0
z_velocity = 0
temperature = 1
pressure -> from ideal gas equation connecting density, temperature and pressure

The density was adjusted suitably based on the constants, to have a period of 2*pi. On calculating the RHS of the x-momentum Navier Stokes based on these values given and comparing it with my answer (calculated by hand), I obtained a maximum error of the order of E-03.

Further, using these values as initial values of the variables and moving forward in time by a Runge-Kutta 4 scheme, I get values of the density that seem to diverge very quickly. After about 30 time steps, I get NaNs.

1. Is there a specific reason why I notice a decrease in precision when non trigonometric periodic functions are used ?

2. Is 1. related to why my code seems to produce unstable results when marching forward in time ?

I wouldn't mind pasting the code here but it's pretty large.

I thought I would plot the initial density and its variation. But turns out I can't as I do not have enough reputation to do so.

The initial plot (@t = 0.0s) is a density plot that looks like a rectangular wave with the tanh functions used to smoothen the wave at the various corners.

At around t = 0.10s (the time step is 0.01s so, after 10 iterations), it develops spikes and becomes non-differentiable (still continuous).

Answer 1

There are three issues that are likely to cause such problems in pseudospectral methods:

1. Gibbs oscillations
2. Aliasing
3. Time step too large

In any case you likely develop oscillations in the solution until some point ends up with a negative density, resulting in a NaN when computing the pressure or sound speed or some other term. The solution to 3 is obvious, decrease the time step until the time integration is stable. The other two are more nuanced.

Gibbs oscillations

Gibbs oscillations arise when computing the Fourier series of discontinuous functions. Gibbs oscillations arise in the derivatives if the function is non-smooth. If you have large jumps in your initial conditions then the Fourier series will match the values at grid points exactly, but the derivative will have large oscillations, leading to loss of precision in the derivative (right hand side) computations. See the image below for a demonstration of this, the values match but the derivative does not. As a rule of thumb, jumps must be smoothed out over about 10 grid to prevent this behavior.

Even if your initial conditions are smooth on the scale of the grid, the state variables may quickly steepen. In compressible Navier-Stokes, the viscous terms act to prevent shocks from forming, but if your simulation is not sufficiently resolved you will still develop jumps in your simulation. Sufficiently resolved means having grid spacing small enough to capture the viscous dissipation, which can be estimated by looking at the Kolmogorov scale, see this PDF. This quickly leads to large, non-physical oscillations and a divergent solution.

Aliasing

Aliasing occurs in pseudospectral methods due to the presence of nonlinear terms (e.g. $u u_x$) in the evolution equation. Computation of the derivatives in spectral space assumes that you are resolving a certain number of wavelengths. However, nonlinear terms continuously generate higher and higher wavenumbers. In a discrete problem, these higher wavenumbers are "aliased" back to affect the lower wavenumbers that can actually be represented at the chosen resolution. This corrupts the lower wavenumber values and can quickly lead to oscillations, non-physical results, and the simulation blowing up.

A simple demonstration of how nonlinear terms generate higher wavenumbers, and how those wavenumbers are aliased to lower wavenumbers is shown in the following scenario (see image below):
Take a grid of 7 points on the interval [0,1). Let $a = \cos(6\pi x )$ and $b = \sin(6\pi x)$. These terms both have (angular) frequencies of $6\pi$. The term $$ab = \cos(6\pi x )\sin(6\pi x ) = \frac{1}{2}\sin(12\pi x)$$ has frequency $12\pi$. This frequency cannot be resolved on the chosen grid, but the values of $ab$ on the grid are equal to the values of $-\sin(2\pi x)/2$ on the grid. So instead of a $12\pi$ frequency component (which is not captured in the discretized space), the result of $ab$ appears as a $2\pi$ frequency component.

There are multiple options available to prevent aliasing from corrupting your results, commonly termed dealiasing. The most common methods are zero padding (3/2-rule, effectively increasing grid resolution before nonlinear multiplications and then discarding the higher frequencies), 2/3-rule truncation (zeroing out the highest wavenumbers before nonlinear multiplications, original paper), or filtering procedures (e.g. this method).

Additionally, there is evidence that for well resolved simulations dealiasing is not crucial since the highest wavenumbers components (which may produce aliasing) are small due to the viscous dissipation.

This presentation (PDF) also provides a good overview of dealiasing, including some of the history.