Fourier Continuation Method Research Articles

Vector potential-based MHD solver for non-periodic flows using Fourier continuation expansions

A high-order method to evolve in time electromagnetic and velocity fields in conducting fluids with non-periodic boundaries is presented. The method has a small overhead compared with fast FFT-based pseudospectral methods in periodic domains. It uses the magnetic vector potential formulation for accurately enforcing the null divergence of the magnetic field, and allowing for different boundary conditions including perfectly conducting walls or vacuum surroundings, two cases relevant for many astrophysical, geophysical, and industrial flows. A spectral Fourier continuation method is used to accurately represent all fields and their spatial derivatives, allowing also for efficient solution of Poisson equations with different boundaries. A study of conducting flows at different Reynolds and Hartmann numbers, and with different boundary conditions, is presented to study convergence of the method and the accuracy of the solenoidal and boundary conditions.

The smooth forcing extension method: A high-order technique for solving elliptic equations on complex domains

High-order numerical methods for solving elliptic equations over arbitrary domains typically require specialized machinery, such as high-quality conforming grids for finite elements method, and quadrature rules for boundary integral methods. These tools make it difficult to apply these techniques to higher dimensions. In contrast, fixed Cartesian grid methods, such as the immersed boundary (IB) method, are easy to apply and generalize, but typically are low-order accurate. In this study, we introduce the Smooth Forcing Extension (SFE) method, a fixed Cartesian grid technique that builds on the insights of the IB method, and allows one to obtain arbitrary orders of accuracy. Our approach relies on a novel Fourier continuation method to compute extensions of the inhomogeneous terms to any desired regularity. This is combined with the highly accurate Non-Uniform Fast Fourier Transform for interpolation operations to yield a fast and robust method. Numerical tests confirm that the technique performs precisely as expected on one-dimensional test problems. In higher dimensions, the performance is even better, in some cases yielding sub-geometric convergence. We also demonstrate how this technique can be applied to solving parabolic problems and for computing the eigenvalues of elliptic operators on general domains, in the process illustrating its stability and amenability to generalization.

Fourier continuation method for incompressible fluids with boundaries

We present a Fourier Continuation-based parallel pseudospectral method for incompressible fluids in cuboid non-periodic domains. The method produces dispersionless and dissipationless derivatives with fast spectral convergence inside the domain, and with very high order convergence at the boundaries. Incompressibility is imposed by solving a Poisson equation for the pressure. Being Fourier-based, the method allows for fast computation of spectral transforms. It is compatible with uniform grids (although refined or nested meshes can also be implemented), which in turn allows for explicit time integration at sufficiently high Reynolds numbers. Using a new parallel code named SPECTER we illustrate the method with two problems: channel flow, and plane Rayleigh–Bénard convection under the Boussinesq approximation. In both cases the method yields results compatible with previous studies using other high-order numerical methods, with mild requirements on the time step for stability.

The Fourier extension method and discrete orthogonal polynomials on an arc of the circle

The Fourier extension method, also known as the Fourier continuation method, is a method for approximating non-periodic functions on an interval using truncated Fourier series with period larger than the interval on which the function is defined. When the function being approximated is known at only finitely many points, the approximation is constructed as a projection based on this discrete set of points. In this paper we address the issue of estimating the absolute error in the approximation. The error can be expressed in terms of a system of discrete orthogonal polynomials on an arc of the unit circle, and these polynomials are then evaluated asymptotically using Riemann–Hilbert methods.

Two-dimensional local Fourier image reconstruction via domain decomposition Fourier continuation method.

The MRI image is obtained in the spatial domain from the given Fourier coefficients in the frequency domain. It is costly to obtain the high resolution image because it requires higher frequency Fourier data while the lower frequency Fourier data is less costly and effective if the image is smooth. However, the Gibbs ringing, if existent, prevails with the lower frequency Fourier data. We propose an efficient and accurate local reconstruction method with the lower frequency Fourier data that yields sharp image profile near the local edge. The proposed method utilizes only the small number of image data in the local area. Thus the method is efficient. Furthermore the method is accurate because it minimizes the global effects on the reconstruction near the weak edges shown in many other global methods for which all the image data is used for the reconstruction. To utilize the Fourier method locally based on the local non-periodic data, the proposed method is based on the Fourier continuation method. This work is an extension of our previous 1D Fourier domain decomposition method to 2D Fourier data. The proposed method first divides the MRI image in the spatial domain into many subdomains and applies the Fourier continuation method for the smooth periodic extension of the subdomain of interest. Then the proposed method reconstructs the local image based on L2 minimization regularized by the L1 norm of edge sparsity to sharpen the image near edges. Our numerical results suggest that the proposed method should be utilized in dimension-by-dimension manner instead of in a global manner for both the quality of the reconstruction and computational efficiency. The numerical results show that the proposed method is effective when the local reconstruction is sought and that the solution is free of Gibbs oscillations.

Regularized integral formulation of mixed Dirichlet-Neumann problems

This paper presents a theoretical discussion as well as novel solution algorithms for problems of scattering on smooth two-dimensional domains under Zaremba boundary conditions, for which Dirichlet and Neumann conditions are specified on various portions of the domain boundary. The theoretical basis of the proposed numerical methods, which is provided for the first time in the present contribution, concerns detailed information about the singularity structure of solutions of the Helmholtz operator under boundary conditions of Zaremba type. The new numerical method is based on the use of Green functions and integral equations, and it relies on the Fourier continuation method for regularization of all smooth-domain Zaremba singularities as well as newly derived quadrature rules which give rise to high-order convergence, even around Zaremba singular points. As demonstrated in this paper, the resulting algorithms enjoy high-order convergence, and they can be used to efficiently solve challenging Helmholtz boundary value problems and Laplace eigenvalue problems with high-order accuracy.

A Domain Decomposition Fourier Continuation Method for Enhanced $$L_1$$ L 1 Regularization Using Sparsity of Edges in Reconstructing Fourier Data

$$L_1$$ regularization is widely used in various applications for sparsifying transform. In Wasserman et al. (J Sci Comput 65(2):533–552, 2015) the reconstruction of Fourier data with $$L_1$$ minimization using sparsity of edges was proposed—the sparse PA method. With the sparse PA method, the given Fourier data are reconstructed on a uniform grid through the convex optimization based on the $$L_1$$ regularization of the jump function. In this paper, based on the method proposed by Wasserman et al. (J Sci Comput 65(2):533–552, 2015) we propose to use the domain decomposition method to further enhance the quality of the sparse PA method. The main motivation of this paper is to minimize the global effect of strong edges in $$L_1$$ regularization that the reconstructed function near weak edges does not benefit from the sparse PA method. For this, we split the given domain into several subdomains and apply $$L_1$$ regularization in each subdomain separately. The split function is not necessarily periodic, so we adopt the Fourier continuation method in each subdomain to find the Fourier coefficients defined in the subdomain that are consistent to the given global Fourier data. The numerical results show that the proposed domain decomposition method yields sharp reconstructions near both strong and weak edges. The proposed method is suitable when the reconstruction is required only locally.

A boundary integral algorithm for the Laplace Dirichlet–Neumann mixed eigenvalue problem

We present a novel integral-equation algorithm for evaluation of Zaremba eigenvalues and eigenfunctions, that is, eigenvalues and eigenfunctions of the Laplace operator with mixed Dirichlet–Neumann boundary conditions; of course, (slight modifications of) our algorithms are also applicable to the pure Dirichlet and Neumann eigenproblems. Expressing the eigenfunctions by means of an ansatz based on the single layer boundary operator, the Zaremba eigenproblem is transformed into a nonlinear equation for the eigenvalue μ. For smooth domains the singular structure at Dirichlet–Neumann junctions is incorporated as part of our corresponding numerical algorithm—which otherwise relies on use of the cosine change of variables, trigonometric polynomials and, to avoid the Gibbs phenomenon that would arise from the solution singularities, the Fourier Continuation method (FC). The resulting numerical algorithm converges with high order accuracy without recourse to use of meshes finer than those resulting from the cosine transformation. For non-smooth (Lipschitz) domains, in turn, an alternative algorithm is presented which achieves high-order accuracy on the basis of graded meshes. In either case, smooth or Lipschitz boundary, eigenvalues are evaluated by searching for zero minimal singular values of a suitably stabilized discrete version of the single layer operator mentioned above. (The stabilization technique is used to enable robust non-local zero searches.) The resulting methods, which are fast and highly accurate for high- and low-frequencies alike, can solve extremely challenging two-dimensional Dirichlet, Neumann and Zaremba eigenproblems with high accuracies in short computing times—enabling, in particular, evaluation of thousands of eigenvalues and corresponding eigenfunctions for a given smooth or non-smooth geometry with nearly full double-precision accuracy.

A Fourier Continuation Method for the Solution of Elliptic Eigenvalue Problems in General Domains

We present a new computational method for the solution of elliptic eigenvalue problems with variable coefficients in general two-dimensional domains. The proposed approach is based on use of the novel Fourier continuation method (which enables fast and highly accurate Fourier approximation of nonperiodic functions in equispaced grids without the limitations arising from the Gibbs phenomenon) in conjunction with an overlapping patch domain decomposition strategy and Arnoldi iteration. A variety of examples demonstrate the versatility, accuracy, and generality of the proposed methodology.

Hybrid Fourier-Continuation Method and Weighted Essentially Non-oscillatory Finite Difference Scheme for Hyperbolic Conservation Laws in a Single-Domain Framework

We investigate a hybrid Fourier-Continuation (FC) method (Bruno and Lyon, J Comput Phys 229:2009---2033, 2010) and fifth order characteristic-wise weighted essentially non-oscillatory (WENO) finite...

High-order numerical methods for nonlinear acoustics: A Fourier Continuation approach

Dispersion errors, which result from the use of low-order numerical methods in wave-propagation and transport problems, can have a devastating impact on the accuracy of acoustic simulations. These errors are especially problematic in settings containing nonlinear acoustic waves that propagate many times their fundamental wavelength. In these cases, the use of high-order numerical schemes is vital for the accurate and efficient evaluation of the acoustic field. We present a class of high-order time-domain solvers for the treatment of nonlinear acoustic propagation problems. These solvers are based on the Fourier Continuation method, which produces rapidly-converging Fourier series expansions of non-periodic functions (thereby avoiding the Gibbs phenomenon), and are capable of accurately and efficiently simulating nonlinear acoustic fields in large, complex domains.

A spectral FC solver for the compressible Navier–Stokes equations in general domains I: Explicit time-stepping

We present a Fourier continuation (FC) algorithm for the solution of the fully nonlinear compressible Navier–Stokes equations in general spatial domains. The new scheme is based on the recently introduced accelerated FC method, which enables use of highly accurate Fourier expansions as the main building block of general-domain PDE solvers. Previous FC-based PDE solvers are restricted to linear scalar equations with constant coefficients. The FC methodology presented in this text thus constitutes a significant generalization of the previous FC schemes, as it yields general-domain FC solvers for nonlinear systems of PDEs. While not restricted to periodic boundary conditions and therefore applicable to general boundary value problems on arbitrary domains, the proposed algorithm inherits many of the highly desirable properties arising from rapidly convergent Fourier expansions, including high-order convergence, essentially spectrally accurate dispersion relations, and much milder CFL constraints than those imposed by polynomial-based spectral methods—since, for example, the spectral radius of the FC first derivative grows linearly with the number of spatial discretization points. We demonstrate the accuracy and optimal parallel efficiency of the algorithm in a variety of scientific and engineering contexts relevant to fluid-dynamics and nonlinear acoustics.

Numerical simulation of focused nonlinear acoustic beams: A Fourier continuation direct solver and comparisons to a paraxial approximation.

A recently developed Fourier-continuation (FC) method is used to develop a numerical algorithm, which can accurately solve the fully nonlinear acoustic vector wave equation for the type of hundred-wavelength domains arising in the field of focused medical ultrasound. The FC nonlinear-ultrasound solver is used to outline the domain of applicability of the widely used Khokhlov–Zabolotskaya–Kuznetsov (KZK) approximation in focused ultrasound particularly in configurations arising in the field of high intensity focused ultrasound (HIFU). Good agreement was found for small-radius transducers; but as the source radius was increased to dimensions equivalent to practical HIFU sources the KZK approximation leads to errors in the location of the predicted focus. Examples for media with multiple scattering centers also show the utility of the FC approach in modeling complex media. In general the FC method produces solutions with significantly reduced discretization requirements over those associated with finite-difference and finite-element methods, and they lead to improvements in computing times by factors of hundreds and to thousands over competing approaches. [Work supported by NSF Grant No. 0835795.]

Comparison of the Khokhlov–Zabolotskaya–Kuznetsov equation and Fourier-continuation for modeling high-intensity focused ultrasound beams.

High-intensity focused ultrasound (HIFU) employs acoustic amplitudes that are high enough that nonlinear propagation effects are important in the evolution of the sound field. A common model for HIFU beams is the Khokhlov–Zabolotskaya–Kuznetsov (KZK) equation which accounts for nonlinearity and absorption of sound beams. The KZK equation models diffraction using the parabolic or paraxial approximation. For many HIFU sources the source aperture is similar to the focal length and the paraxial approximation may not be appropriate. Here we compare results obtained using the “Texas code” a time-domain numerical solution to the KZK equation to a new method of solving nonlinear differential equations: the Fourier-continuation (FC) method. Here the FC method solves the underlying fluid dynamics equations and does not invoke the paraxial approximation. Results were obtained for a 1.1-MHz focused HIFU source. The medium parameters were taken to either be hyperviscous water or to match human liver. For a low focusing gain transducer (focal length 50λ and radius 10λ) the KZK and FC models showed excellent agreement through the focal region. As the source radius was increased, discrepancies started to appear and for a radius of 30λ, the KZK model did not capture diffraction of the HIFU source accurately. [Work supported by NSF 0835795.]

High-order unconditionally stable FC-AD solvers for general smooth domains II. Elliptic, parabolic and hyperbolic PDEs; theoretical considerations

A new PDE solver was introduced recently, in Part I of this two-paper sequence, on the basis of two main concepts: the well-known Alternating Direction Implicit (ADI) approach, on one hand, and a certain ''Fourier Continuation'' (FC) method for the resolution of the Gibbs phenomenon, on the other. Unlike previous alternating direction methods of order higher than one, which only deliver unconditional stability for rectangular domains, the new high-order FC-AD (Fourier-Continuation Alternating-Direction) algorithm yields unconditional stability for general domains-at an O(Nlog(N)) cost per time-step for an N point spatial discretization grid. In the present contribution we provide an overall theoretical discussion of the FC-AD approach and we extend the FC-AD methodology to linear hyperbolic PDEs. In particular, we study the convergence properties of the newly introduced FC(Gram) Fourier Continuation method for both approximation of general functions and solution of the alternating-direction ODEs. We also present (for parabolic PDEs on general domains, and, thus, for our associated elliptic solvers) a stability criterion which, when satisfied, ensures unconditional stability of the FC-AD algorithm. Use of this criterion in conjunction with numerical evaluation of a series of singular values (of the alternating-direction discrete one-dimensional operators) suggests clearly that the fifth-order accurate class of parabolic and elliptic FC-AD solvers we propose is indeed unconditionally stable for all smooth spatial domains and for arbitrarily fine discretizations. To illustrate the FC-AD methodology in the hyperbolic PDE context, finally, we present an example concerning the Wave Equation-demonstrating sixth-order spatial and fourth-order temporal accuracy, as well as a complete absence of the debilitating ''dispersion error'', also known as ''pollution error'', that arises as finite-difference and finite-element solvers are applied to solution of wave propagation problems.

Hybrid Fourier-Continuation method and WENO-Z finite difference scheme for multi-dimensional detonation structure simulations

Hybrid Fourier-Continuation method and WENO-Z finite difference scheme for multi-dimensional detonation structure simulations