Fourier Pseudo-Spectral Spatial Discretization

• Fourier pseudo-spectral spatial discretization is a high-order numerical technique that expands functions using truncated Fourier series and FFT to accurately approximate derivatives in periodic PDEs.
• It computes derivatives in Fourier space by exploiting the property that differentiation corresponds to multiplication by ik, while managing nonlinearities via physical space evaluation and dealiasing.
• Widely used in direct numerical simulations of turbulence, this method offers exponential convergence and computational speedups, making it essential for CFD applications.

Fourier pseudo-spectral spatial discretization is a numerical technique for approximating derivatives and solving partial differential equations (PDEs), particularly in domains with periodic boundary conditions and smooth solutions. The method leverages the fast Fourier transform (FFT) to achieve spectral accuracy in spatial discretization, making it highly attractive for direct numerical simulations (DNS) of incompressible turbulence and other PDE-dominated physical systems (Ketcheson et al., 2018). This approach is extensively utilized in computational fluid dynamics and turbulence research due to its combination of high order accuracy and computational efficiency.

1. Core Principles of Fourier Pseudo-Spectral Discretization

The Fourier pseudo-spectral method approximates spatial derivatives by expanding the unknown field (e.g., velocity, vorticity) in a truncated Fourier series,

$$u(x) = \sum_{k=-N/2}^{N/2-1} \hat{u}_k e^{ikx}, \qquad x \in [0, L],$$

where $\hat{u}_k$ are the discrete Fourier coefficients, $N$ is the number of collocation points, and $k$ runs over the discrete wavenumbers compatible with the periodic box. The method achieves “spectral accuracy”: when the underlying solution is smooth, the error in the spectral approximation decays faster than any algebraic order in $N$.

To compute derivatives, the method exploits the exact property that differentiation in the physical domain corresponds to multiplication by $ik$ in Fourier space:

$\widehat{u^{(n)}}_k = (ik)^n \hat{u}_k.$

In practical implementations, nonlinear terms are typically computed in physical space by transforming the solution from Fourier to physical space (via inverse FFT), evaluating products, and then transforming back (via forward FFT). This “pseudo-spectral” approach avoids the high computational complexity of evaluating convolutions in Fourier space.

2. Algorithmic Workflow and Implementation

A typical pseudo-spectral DNS workflow consists of the following steps (Ketcheson et al., 2018):

1. Grid and Initialization: Define an equispaced grid of $N$ points in space and initialize the physical field $u(x_j)$.
2. Forward FFT: Transform $u(x_j)$ to Fourier space to obtain $\hat{u}_k$.
3. Spectral Derivatives: In Fourier space, compute derivatives by multiplying each $\hat{u}_k$ by the appropriate factor of $ik$ (for first derivative) or $(ik)^n$ for higher derivatives.
4. Nonlinearity Evaluation: Transform field variables back to physical space, compute nonlinear products pointwise, then transform back by FFT.
5. Dealiasing: Apply dealiasing schemes, typically the 2/3-rule, by zeroing out upper third of the wavenumber spectrum after evaluating nonlinear terms, to eliminate aliasing errors.
6. Integration in Time: Advance in time using an appropriate ODE solver (commonly explicit or embedded Runge-Kutta methods tailored for stiff or non-stiff regimes).

Such an algorithm is efficiently implemented using FFT libraries, and the arithmetic complexity is $O(N\log N)$ per time step.

3. Application to Direct Numerical Simulation of Turbulence

In incompressible turbulence DNS (e.g., Navier–Stokes equations on a periodic box), Fourier pseudo-spectral methods are the prevailing standard. The method retains all resolvable scales accurately up to $k_{\max} = N/2$ and provides a natural framework for enforcing incompressibility via projection in Fourier space. Aliasing is addressed using spectral dealiasing techniques, essential for energy conservation and accuracy in turbulent cascade simulations.

The time integration of the resulting ODEs is typically performed using explicit Runge–Kutta methods, including high-order embedded pairs (e.g., Bogacki–Shampine BS5(4), Dormand–Prince DP5(4)) to achieve higher order and enable adaptive step size selection. Such embedded pairs estimate local error efficiently by comparing the difference between the high-order and embedded lower-order solutions at each time step (Ketcheson et al., 2018).

Speedups of $2\times$–$10\times$ over classical fourth-order Runge–Kutta with fixed step sizes have been demonstrated by deploying efficient adaptivity within pseudo-spectral DNS (Ketcheson et al., 2018).

4. Accuracy, Stability, and Error Control

Spectral methods deliver exponential convergence for smooth solutions, contingent on the absence of physical or numerical discontinuities. However, precision hinges on sufficient grid resolution and accurate dealiasing:

• Spectral Accuracy: Error decays faster than algebraic in $N$ for smooth, periodic functions.
• Aliasing Error: Nonlinear term evaluation can introduce spectral aliasing errors, mitigated by dealiasing schemes (e.g., 2/3 rule).
• CFL Stability Constraint: The stability of explicit time integration is subject to the CFL condition,

$$\Delta t \lesssim C_{\max}\frac{\Delta x}{\|u\|_\infty},$$

where $C_{\max}$ is determined by the time-stepping scheme. Adaptive step size controllers must always enforce this upper bound (Ketcheson et al., 2018).

Adaptive embedded methods compute the local error, $E_n = y_{n+1} - \hat{y}_{n+1}$, at each time step and adjust time steps according to

$$h_{\text{new}} = h \times \delta (1/\text{err})^{1/(q+1)},$$

where $q$ is the order of the embedded solver and $\delta$ a safety factor. This permits robust error control and optimizes computational cost (Ketcheson et al., 2018).

5. Representative Performance and Benchmark Problems

Fourier pseudo-spectral methods have been tested on canonical turbulence problems:

• Taylor–Green Vortex: Simulations across $\text{Re}=280$–$1600$ and $64^3$–$512^3$ grid points exhibit exponential convergence and highlight computational savings when adaptive embedded schemes are used.
• Rayleigh–Taylor Instability: For both 2D ($512 \times 2048$) and 3D ($256 \times 256 \times 1024$) setups, adaptive pseudo-spectral solvers dynamically adjust the step size through laminar and turbulent phases, efficiently tracking solution features while maintaining accuracy.
• Homogeneous Isotropic Turbulence: Adaptive pseudo-spectral DNS using BS5(4) or DP5(4) achieves maximal step sizes limited only by the CFL constraint during the developed turbulent cascade.

In all cases, adaptive pseudo-spectral DNS achieves equivalent (or better) accuracy at a fraction of the computational cost of traditional fixed-step methods, with the measured speedup consistently in the $2\times$–$10\times$ range (Ketcheson et al., 2018).

Test Case	Grid Size	Reference Integrator	Adaptive Integrator	Speedup
Taylor–Green Vortex	$64^3$ to $512^3$	RK4 (fixed step)	BS5(4) (adaptive/fixed step)	$2\times$–$10\times$
Rayleigh–Taylor Instability	$512\times2048$, $256^3\times1024$	RK4 (fixed step)	DP5(4), BS5(4) (adaptive)	$4\times$–$8\times$
Homogeneous Isotropic Turb.	$256^3$ to $1024^3$	RK4 (fixed step)	BS5(4) (adaptive/fixed step)	$4\times$–$10\times$

Performance metrics are measured as the number of right-hand-side function evaluations (i.e., FFT calls), with error measured in $L^2$ and max-norms against reference solutions (Ketcheson et al., 2018).

6. Implementation Considerations and Best Practices

Efficient pseudo-spectral solvers require:

• Use of high-performance FFT libraries for both forward and inverse transforms.
• Enforcement of dealiasing (e.g., 2/3-rule) after each nonlinear term computation to ensure spectral accuracy.
• Careful selection and implementation of time-adaptive, embedded explicit integrators (e.g., BS5(4), DP5(4), KCL5(4)8), including correct error normalization over all spectral modes and variables.
• Step-size controllers must enforce both error tolerance and the absolute CFL stability bound, with appropriate limiting and safety factors.

It is documented that with these practices, spectral DNS achieves substantially improved cost-accuracy balance without compromising resolution or dynamical fidelity, especially when combined with modern adaptive time-integration methods (Ketcheson et al., 2018).

7. Significance and Limitations

The Fourier pseudo-spectral method provides a compelling tool for the high-fidelity simulation of periodic, smooth-flow PDEs, particularly in incompressible turbulence research. Its spectral convergence and computational efficiency, especially when coupled with adaptive, embedded high-order time integration, yield significant improvements in performance and accuracy compared to classical approaches.

Limitations include poor handling of non-periodic or discontinuous solutions and the requirement for careful dealiasing and error control. The method is not generally applicable to complex geometries or boundary conditions, where alternative spectral element, finite difference, or finite volume approaches may be warranted.

References:

• "More efficient time integration for Fourier pseudo-spectral DNS of incompressible turbulence" (Ketcheson et al., 2018)