Fourier Pseudospectral Method Research Articles

On the viscoelastic-electromagnetic-gravitational analogy

The analogy between electromagnetism and gravitation was achieved by linearizing the tensorial gravitational equations of general relativity and converting them into a vector form corresponding to Maxwell’s electromagnetic equations. On this basis, we use the equivalence with viscoelasticity and propose a theory of gravitational waves. We add a damping term to the differential equations, which is equivalent to Ohm’s law in electromagnetism and Maxwell’s viscosity in viscoelasticity, to describe the attenuation of the waves. The differential equations in viscoelasticity are those of cross-plane shear waves, commonly referred to as SH waves. A plane-wave analysis gives the phase velocity, the energy velocity, the quality factor and the attenuation factor of the field as well as the energy balance. To obtain these properties, we use the analogy with viscoelasticity; the properties of electromagnetic and gravitational waves are similar to those of shear waves. The presence of attenuation means that the transient field is generally a composition of inhomogeneous (non-uniform) plane waves, where the propagation and attenuation vectors do not point in the same direction and the phase velocity vector and the energy flux (energy velocity) are not collinear. The polarization of cross-plane field is linear and perpendicular to the propagation-attenuation plane, while the polarization of the field within the plane is elliptical. Transient wave fields in the space–time domain are analyzed with the Green function (in homogeneous media) and with a grid method (in heterogeneous media) based on the Fourier pseudospectral method for calculating the spatial derivatives and a Runge–Kutta scheme of order 4 for the time stepping. In the examples, wave propagation at the Sun–Earth and Earth–Moon distances using quadrupole sources is considered in comparison to viscoelastic waves. The Green and grid solutions are compared to test the latter algorithm. Finally, an example of propagation in heterogeneous media is presented.

GPU-enabled extreme-scale turbulence simulations: Fourier pseudo-spectral algorithms at the exascale using OpenMP offloading

Fourier pseudo-spectral methods for nonlinear partial differential equations are of wide interest in many areas of advanced computational science, including direct numerical simulation of three-dimensional (3-D) turbulence governed by the Navier-Stokes equations in fluid dynamics. This paper presents a new capability for simulating turbulence at a new record resolution up to 35 trillion grid points, on the world's first exascale computer, Frontier, comprising AMD MI250x GPUs with HPE's Slingshot interconnect and operated by the US Department of Energy's Oak Ridge Leadership Computing Facility (OLCF). Key programming strategies designed to take maximum advantage of the machine architecture involve performing almost all computations on the GPU which has the same memory capacity as the CPU, performing all-to-all communication among sets of parallel processes directly on the GPU, and targeting GPUs efficiently using OpenMP offloading for intensive number-crunching including 1-D Fast Fourier Transforms (FFT) performed using AMD ROCm library calls. With 99% of computing power on Frontier being on the GPU, leaving the CPU idle leads to a net performance gain via avoiding the overhead of data movement between host and device except when needed for some I/O purposes. Memory footprint including the size of communication buffers for MPI_ALLTOALL is managed carefully to maximize the largest problem size possible for a given node count.Detailed performance data including separate contributions from different categories of operations to the elapsed wall time per step are reported for five grid resolutions, from 20483 on a single node to 327683 on 4096 or 8192 nodes out of 9408 on the system. Both 1D and 2D domain decompositions which divide a 3D periodic domain into slabs and pencils respectively are implemented. The present code suite (labeled by the acronym GESTS, GPUs for Extreme Scale Turbulence Simulations) achieves a figure of merit (in grid points per second) exceeding goals set in the Center for Accelerated Application Readiness (CAAR) program for Frontier. The performance attained is highly favorable in both weak scaling and strong scaling, with notable departures only for 20483 where communication is entirely intra-node, and for 327683, where a challenge due to small message sizes does arise. Communication performance is addressed further using a lightweight test code that performs all-to-all communication in a manner matching the full turbulence simulation code. Performance at large problem sizes is affected by both small message size due to high node counts as well as dragonfly network topology features on the machine, but is consistent with official expectations of sustained performance on Frontier. Overall, although not perfect, the scalability achieved at the extreme problem size of 327683 (and up to 8192 nodes — which corresponds to hardware rated at just under 1 exaflop/sec of theoretical peak computational performance) is arguably better than the scalability observed using prior state-of-the-art algorithms on Frontier's predecessor machine (Summit) at OLCF. New science results for the study of intermittency in turbulence enabled by this code and its extensions are to be reported separately in the near future.

Energy Maximization Absorption of Wave Energy Converter Based on Fourier Pseudo-Spectral Method and Adaptive Dynamic Programming

In this paper, a novel noncausal control framework to address the energy maximization problem of wave energy converters (WECs) with physical constraints is proposed. The energy maximization problem of WECs is a constrained optimal control problem. The proposed control framework converts this problem into a reference trajectory-tracking problem through the Fourier pseudo-spectral method (FPSM) and utilizes the online tracking adaptive dynamic programming (OTADP) algorithm to realize real-time trajectory tracking for practical use in the ocean environment. Using the wave prediction technique, the optimal trajectory is generated online through a receding horizon (RH) implementation. A critic neural network (NN) is applied to approximate the optimal cost value function and calculate the error-tracking control by solving the associated Hamilton-Jacobi-Bellman (HJB) equation. The proposed WEC control framework improves computational efficiency and makes it feasible for the online control of the WEC in practice. Simulation results show the effects of the receding-horizon implementation of FPSM with different window lengths and window functions, while verifying the performances of tracking control and energy absorption of WECs in two different sea conditions.

Numerical study on the polarity change process during internal wave shoaling

Polarity change is an important mechanism for internal waves shoaling. In this study, a numerical model for simulating the real-scale internal wave passing over slope-shelf topography is established based on the Fourier pseudo-spectral method and weakly nonlinear theory. By numerical simulation, the effects of shelf height, initial wave amplitude, and inclination angle on the waveform characteristics and energy properties of the internal wave during its shoaling are investigated. In the polarity change process, the initial internal wave converts into a depression wave and a generated elevation wave behind it. The distance between the peak of the elevation wave and the trough of the depression wave is a key feature to describe the polarity change. In terms of energy properties, the energy ratio of depression and generated elevation waves compared with the initial wave as well as their relative magnitude is mainly determined by the shelf height. In addition, the initial wave amplitude also affects the generation of the elevation wave and the attenuation of the depression wave to a certain extent. The increase in the inclination angle hinders the formation of the elevation wave but has little effect on the depression wave energy.

Early-time resonances in the three-dimensional wall-bounded axisymmetric Euler and related equations

We investigate the complex-time analytic structure of solutions of the three-dimensional (3D)-axisymmetric, wall-bounded, incompressible Euler equations, by starting with the initial data proposed in Luo and Hou [“Potentially singular solutions of the 3D axisymmetric Euler equations,” Proc. Natl. Acad. Sci. U. S. A. 111(36), 12968–12973 (2014)], to study a possible finite-time singularity. We use our pseudospectral Fourier–Chebyshev method [Kolluru et al., “Insights from a pseudospectral study of a potentially singular solution of the three-dimensional axisymmetric incompressible Euler equation,” Phys. Rev. E 105(6), 065107 (2022)], with quadruple-precision arithmetic, to compute the time-Taylor-series coefficients of the flow fields, up to a high order. We show that the resulting approximations display early-time resonances; the initial spatial location of these structures is different from that for the tygers, which we have obtained in Kolluru et al. [“Insights from a pseudospectral study of a potentially singular solution of the three-dimensional axisymmetric incompressible Euler equation,” Phys. Rev. E 105(6), 065107 (2022)]. We then perform asymptotic analysis of the Taylor-series coefficients, by using generalized ratio methods, to extract the location and nature of the convergence-limiting singularities and demonstrate that these singularities are distributed around the origin, in the complex-t2 plane, along two curves that resemble the shape of an eye. We obtain similar results for the one-dimensional (1D) wall-approximation (of the full 3D-axisymmetric Euler equation) called the 1D Hou–Luo model, for which we use Fourier-pseudospectral methods to compute the time-Taylor-series coefficients of the flow fields. Our work examines the link between tygers, in Galerkin-truncated pseudospectral studies, and early-time resonances, in truncated time-Taylor expansions of solutions of partial differential equations, such as those we consider.

Numerical study of the impacts of stochastic forcing on the vortex in fluid flow

This paper focuses on a numerical study about the stochastic Navier–Stokes equations. Unlike previous studies, this paper focuses on studying these equations from the perspective of vortices. The vorticity–stream function method was proposed to deal with incompressible fluid flow. And a Crank–Nicolson Fourier pseudo-spectral method was put forward to solve the formulation of stream function equation. In addition, we have conducted some numerical experiments to observe the effects of random forcing on vortices in the fluid flow by utilizing stochastic solution.

Error estimates of exponential wave integrators for the Dirac equation in the massless and nonrelativistic regime

AbstractWe present exponential wave integrator Fourier pseudospectral (EWI‐FP) methods and establish their error estimates of the fully discrete schemes for the Dirac equation in the massless and nonrelativistic regime. This regime involves a small dimensionless parameter where , and is inversely proportional to the speed of light. The solution exhibits highly oscillatory behavior in time and rapid wave propagation in space in this regime. Specifically, the time oscillations have a wavelength of , while the spatial oscillations have a wavelength of , with a wave speed of . We employ (symmetric) exponential wave integrators for temporal derivatives and Fourier spectral discretization for spatial derivatives. We rigorously derive the error bounds which explicitly depend on the mesh size , the time step and the small dimensionless parameter . The error estimates for the EWI‐FP methods demonstrate that their meshing strategy requirement (‐scalability) necessitates setting and when . Finally, some numerical examples are provided to validate the error bounds.

Multi-Symplectic Method for the Two-Component Camassa–Holm (2CH) System

In this paper, the multi-symplectic formulations of the two-component Camassa–Holm system are presented. Both the multi-symplectic structure and two local conservation laws of the generalized two-component Camassa–Holm model are proposed for its first-order canonical form. Then, combining the Fourier pseudo-spectral method in the spatial domain with the midpoint method in the time dimension, the multi-symplectic Fourier pseudo-spectral scheme is constructed for the first-order canonical form. Meanwhile, the discrete scheme of the residuals of the multi-symplectic structure and two local conservation laws are also provided. By using the multi-symplectic Fourier pseudo-spectral scheme, the evolution of one- and two-soliton solutions for the generalized two-component Camassa–Holm model is regained. The structure-preserving properties and the reliability of the numerical scheme are illustrated by the tiny numerical residuals (less than 3.5 × 10−8) of the conservation laws as well as the tiny numerical variations (less than 1 × 10−9) of the amplitudes and the propagating velocities of the solitons.

High order energy-preserving method for the space fractional Klein–Gordon-Zakharov equations

The space fractional Klein–Gordon-Zakharov equations are transformed into the multi-symplectic structure system by introducing new auxiliary variables. The multi-symplectic system, which satisfies the multi-symplectic conservation, local energy and momentum conservation, is discretizated into the semi-discrete multi-symplectic system by the Fourier pseudo-spectral method. The second order multi-symplectic average vector field method is applied to the semi-discrete system. The fully discrete energy preserving scheme of the space fractional Klein–Gordon-Zakharov equation is obtained. Based on the composition method, a fourth order energy preserving scheme of the Riesz space fractional Klein–Gordon-Zakharov equations is also obtained. Numerical experiments confirm that these new schemes can have computing ability for a long time and can well preserve the discrete energy conservation property of the equations.

Conformal structure-preserving SVM methods for the nonlinear Schrödinger equation with weakly linear damping term

In this paper, by applying the supplementary variable method (SVM), some high-order, conformal structure-preserving, linearized algorithms are developed for the damped nonlinear Schrödinger equation. We derive the well-determined SVM systems with the conformal properties and they are then equivalent to nonlinear equality constrained optimization problems for computation. The deduced optimization models are discretized by using the Gauss type Runge-Kutta method and the prediction-correction technique in time as well as the Fourier pseudo-spectral method in space. Numerical results and some comparisons between this method and other reported methods are given to favor the suggested method in the overall performance. It is worthwhile to emphasize that the numerical strategy in this work could be extended to other conservative or dissipative system for designing high-order structure-preserving algorithms.

An adaptive time-stepping Fourier pseudo-spectral method for the Zakharov-Rubenchik equation

An adaptive time-stepping Fourier pseudo-spectral method for the Zakharov-Rubenchik equation

Determination of an unknown coefficient in the Korteweg–de Vries equation

Abstract In this paper, a space-time spectral method for solving an inverse problem in the Korteweg–de Vries equation is considered. Optimal order of convergence of the semi-discrete method is obtained in L 2 {L^{2}} -norm. The discrete schemes of the method are based on the modified Fourier pseudospectral method in spatial direction and the Legendre-tau method in temporal direction. The nonlinear term is computed via the fast Fourier transform and fast Legendre transform. The method is implemented by the explicit-implicit iterative method. Numerical results are given to show the accuracy and capability of this space-time spectral method.

Second-order Sobolev gradient flows for computing ground state of ultracold Fermi gases

Based on density functional theory, the ground state for ultracold Fermi gases with dipole–dipole interaction is a functional minimization problem with orthonormality constraints. Many methods such as direct minimization, self-consistent iteration method and first-order gradient flow had been proposed to solve such problem. In this paper, we introduce the second-order Sobolev gradient flows, solve them to the steady state and get the ground state of ultracold Fermi gases numerically. Two contributions have been made: firstly we extend the usual first-order gradient flows to the second-order gradient flows; secondly, we extend the usual L2 gradient to the Sobolev gradient and get the Sobolev gradient flows. We prove that the second-order Sobolev gradient flows have the properties of orthonormality preserving and ‘modified’ energy diminishing, which are desirable for the computation of the ground state solution. Several numerical techniques have been used to discretize the time-dependent second-order Sobolev gradient flows. These include explicit leap-frog method and stabilized semi-implicit method for temporal discretization and Fourier pseudo-spectral method for spatial variables. Extensive numerical examples for the ground state are reported to demonstrate the power of the numerical methods.

A novel and simple spectral method for nonlocal PDEs with the fractional Laplacian

We propose a novel and simple spectral method based on the semi-discrete Fourier transforms to discretize the fractional Laplacian (−Δ)α2. Numerical analysis and experiments are provided to study its performance. Our method has the same symbol |ξ|α as the fractional Laplacian (−Δ)α2 at the discrete level, and thus it can be viewed as the exact discrete analogue of the fractional Laplacian. This unique feature distinguishes our method from other existing methods for the fractional Laplacian. Note that our method is different from the Fourier pseudospectral methods in the literature which are usually limited to periodic boundary conditions (see Remark 1.1). Numerical analysis shows that our method can achieve a spectral accuracy. The stability and convergence of our method in solving the fractional Poisson equations were analyzed. Our scheme yields a multilevel Toeplitz stiffness matrix, and thus fast algorithms can be developed for efficient matrix-vector multiplications. The computational complexity is O(2Nlog⁡(2N)), and the memory storage is O(N) with N the total number of points. Extensive numerical experiments verify our analytical results and demonstrate the effectiveness of our method in solving various problems.

Improved uniform error bounds of a Lawson-type exponential wave integrator method for the Klein-Gordon-Dirac equation

For the Klein-Gordon-Dirac equation (KGDE) with small coupling constant ε∈(0,1], we propose a Lawson-type exponential wave integrator Fourier pseudo-spectral (LEWIFP) method and establish the improved uniform error bounds in the time domain at O(1/ε). We first convert the KGDE to a coupled system and then consider LEWIFP method for the coupled system. The LEWIFP method is proved to be time symmetric which is an important structure in numerical geometric integration. Through careful and rigorous convergence analysis, we establish the error bounds O(hm+ετ2+τ0m) for the full-discretization, where m is determined by the regularity conditions. If the solution is sufficiently smooth with m sufficiently large, we obtain the errors with improved uniform bounds at O(hm+ετ2) in the long-time domain up to O(1/ε). These error bounds are much better than the classical bounds O(hm+τ2) provided by the traditional analysis for the non-Lawson-type exponential wave integrators equipped with Fourier pseudo-spectral method. Combined with the classical analysis tools such as mathematical induction and energy method, we complete our error analysis by adopting the regularity compensation oscillation (RCO) technique which controls the high frequency modes by the regularity of the solution and low frequency modes by phase cancellation. By applying the LEWIFP method to some problems, we show the numerical results to support our error bounds. In addition, the numerical results also show that the discrete mass and energy are stable in the time domain which is long enough. Finally we extend our method to the oscillatory problem.

Direct/split invariant-preserving Fourier pseudo-spectral methods for the rotation-two-component Camassa–Holm system with peakon solitons

The Fourier pseudo-spectral method is well suited to solve PDEs under the periodic boundary condition due to its high-order accuracy and easy-to-implement feature. In this paper, we explore as well as comparatively study four classes of Fourier pseudo-spectral schemes for solving the rotation-two-component Camassa–Holm system which possibly owns peakon solitons. Via exploiting inherent structural properties of the system, we reformulate it into two kinds of different equivalent forms and then apply the Fourier pseudo-spectral method to derive two spatial semi-discrete systems, both of which are proved to preserve the corresponding invariants including mass, momentum and energy. Subsequently, we construct two linearly implicit schemes based on Strang splitting technique and two nonlinear schemes, respectively, for both semi-discrete systems. Owing to the different equivalent forms in the structure, one of the nonlinear schemes preserves discrete mass and momentum, while the other one is shown to preserve all three invariants. Numerical results under the situation of smooth/nonsmooth initial values are provided for distinct types of solutions to test the accuracy in long time simulation and to verify the capacity of predicting water wave propagation, as well as advantages in preserving these invariants. For instance, the present schemes are shown to be at least 14 significant digits, improving upon 10 from ones in previous references.

Improved uniform error bounds on an exponential wave integrator method for the nonlinear Schrödinger equation with wave operator and weak nonlinearity

In this paper, we establish the improved uniform error bounds for the Lawson-type exponential wave integrator Fourier pseudo-spectral (LEWIFP) method for the nonlinear Schrödinger equation with wave operator (NLSW) and weak nonlinearity characterized by a small constant ε∈(0,1]. We first convert the NLSW to a coupled system and then consider the LEWIFP method for the coupled system. The LEWIFP method is proved to be time symmetric and mass-conservative which is very important from the point of view of numerical geometric integration. Through careful and rigorous convergence analysis, we establish improved uniform error bounds for the full-discretization at O(hm+ε2pτ2) in the long-time domain up to O(1/ε2p) where m is determined by the regularity conditions. These error bounds are much better than the classical error bounds O(hm+τ2) provided by the traditional error analysis. In error analysis, in addition to the classical tools such as mathematical induction and energy method, we adopt the regularity compensation oscillation (RCO) technique to analyze the accumulation of errors carefully. The numerical experiments support our error estimates and demonstrate the discrete mass-conservation. In addition, the numerical results show that the discrete energy is stable in the time domain up to O(1/ε2p).

Numerical Simulation of Flows Using the Fourier Pseudospectral Method and the Immersed Boundary Method

The present work proposes the application of a computational methodology based on the coupling of the Fourier Pseudospectral Method (FPSM) and the Immersed Boundary Method (IBM) for conducting flow simulations over slender airfoils. This methodology, termed IMERSPEC, leverages the benefits of both high accuracy and low computational cost inherent in pseudospectral methods, thanks to the utilization of the Fast Fourier Transform algorithm. IBM is employed to impose non-periodic boundary conditions in the Navier–Stokes equations, addressing the requirement of periodicity at boundaries for FPSM convergence and to accurately represent the immersed slender airfoil in the flow. The aerodynamic behavior of the analyzed profiles was assessed by calculating lift and drag coefficients, which were then compared with existing literature results. Consistently favorable outcomes were observed, particularly in flows at lower Reynolds numbers, demonstrating the effectiveness of the IMERSPEC methodology for simulating complex flows computationally. Additionally, weight functions, fundamental to IBM, are employed flexibly for aerodynamic force calculations. Specifically, within the same simulation, a Cubic function is utilized for drag calculation while a Hat function is employed for lift calculation, yielding results more closely aligned with the literature’s findings. This approach offers an alternative to previously proposed methods for IBM implementation.

An exact in time Fourier pseudospectral method with multiple conservation laws for three-dimensional Maxwell’s equations

Maxwell’s equations describe the propagation of electromagnetic waves and are therefore fundamental to understanding many problems encountered in the study of antennas and electromagnetics. The aim of this paper is to propose and analyse an efficient fully discrete scheme for solving three-dimensional Maxwell’s equations. This is accomplished by combining Fourier pseudospectral methods in space and exact formulation in time. Fast computation is efficiently implemented in the scheme by using the matrix diagonalisation method and fast Fourier transform algorithm which are well known in scientific computations. An optimal error estimate which is not encumbered by the CFL condition is established and the resulting scheme is proved to be of spectral accuracy in space and exact in time. Furthermore, the scheme is shown to have multiple conservation laws including discrete energy, helicity, momentum, symplecticity, and divergence-free field conservations. All the theoretical results of the accuracy and conservations are numerically illustrated by two numerical tests.

Coupling of the immersed boundary and Fourier pseudo-spectral methods applied to solve fluid–structure interaction problems

The current manuscript addresses fluid–structure interaction (FSI) problems to the analysis of vortex-induced vibration (VIV), employing the methodology designated as IMERSPEC2D. The approach seamlessly integrates the Fourier pseudo-spectral method (FPSM) and the immersed boundary method (IBM) to solve the Navier–Stokes (NS) and continuity equations. Simultaneously, the effects of cylindrical structure, involving 1 and 2 degrees of freedom (d.o.f.) are simulated, by incorporating structural motion equations into NS via the IBM source term and implementing a two-way FSI algorithm. In the paper are presented the comparison of low and high-order temporal integration methods applied to solve the dynamic behavior of the structure and the dimensionless of the cylinder’s structural motion equations, resulting in an increase in computational time step from 1E(-9\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$-9$$\\end{document}) to 1E(-3\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$-3$$\\end{document}). Validation results to refined mesh simulations include the accurate representation of the lock-in phenomenon to simulations of 1 d.o.f. obtaining the range between reduced velocity 5.456 and 5.568 and the maximal amplitude is 0.352. Notably, the eight-shape pattern of the 2 d.o.f. cylinder displacement is obtained and the formation of C2S wake patterns is achieved, as presented by other authors. In conclusion, the dimensionless approach reduces computational time, while the high-order both of IMERSPEC methodology and time advancement integration to FSI improve the results of the cylinder’s motion and alters distinct vortex shedding patterns in fluid dynamics.