Elsevier

Journal of Computational Physics

Volume 229, Issue 4, 20 February 2010, Pages 1292-1310
Journal of Computational Physics

Global field interpolation for particle methods

https://doi.org/10.1016/j.jcp.2009.10.031Get rights and content

Abstract

A common problem in particle simulations is effective field interpolation. Field interpolation is any method for creating an accurate representation of a given continuous field by a linear combination of overlapping basis functions. This paper presents two techniques for field interpolation, based on a radial basis function (RBF) formulation using Gaussians. The application in mind is vortex methods, where one needs to determine the circulation (or strength) of individual vortex particles with known position and scale to represent a given vorticity field. This process is required both to initially discretize a given vorticity field, and to replace a configuration of particles with another for the purposes of maintaining spatial accuracy. The first technique presented is formulated as an RBF collocation problem, and obtains a solution accurately and with excellent algorithmic efficiency by means of a preconditioned iterative method. The preconditioner is a sparse approximation, based on localization, to the dense coefficient matrix of the RBF system. The second technique uses approximate solutions to the reverse heat equation, recognizing that the standard regularization used in vortex methods (estimating particle strengths using the local value of vorticity multiplied by particle area/volume) corresponds to a Gaussian blurring of the original field. A single time step is used, thus avoiding amplification of high frequencies, and accurate solutions are produced using explicit finite difference methods. Computational experiments were performed in two dimensions, demonstrating the accuracy and convergence of the proposed techniques. Application in three dimensions is straightforward, as radial basis function interpolation is neutral to dimension, but will require more computational effort.

Introduction

Particle and meshfree (or meshless) methods have received great interest in recent years and proven to be very successful for solving a variety of problems. They have a natural application in problems which are inherently discrete, such as the simulation of astrophysical systems or the dynamics of molecules or atoms out of equilibrium. In these applications there is no need for discretization, and the computational particles correspond to the physical objects under study. One method which originated in the astrophysical field and has found great success is the smoothed particle hydrodynamics (SPH) method. This method first appears in the literature in the late 1970s, and was reviewed in [23]. The SPH method uses particles as a discretization technique for a continuum formulation of the problem, and it has been extended to many other applications beyond astrophysics. It joins a variety of other meshfree methods where the particle representation is used as a means of numerical discretization. Another salient member of this group is the vortex particle method, where a particle representation of the vorticity field is formulated in a Lagrangian method used to solve the Navier–Stokes equations. The vortex method traces its origins earlier than SPH, with the regularized vortex particle approach being introduced in [12] but based on the singular vortex particle method first appearing in the classic work of Rosenhead [28]. An overview of the main computational and theoretical issues with the vortex method can be found in the book [13], while standard review papers of the early development of the method are [19], [20].

The term “particle method” does not refer to a single algorithm but rather a diverse collection of computational techniques for solving evolution equations. Researchers analyze and develop different techniques for resolving issues related to the solution of partial differential equations with particles, including capturing diffusion, satisfying boundary conditions, capturing vortex stretching in three-dimensional flows and effective time-marching schemes for the particles. In this paper, we examine the problem of representing a field using particles and/or replacing one set of particles with another for the purpose of maintaining accuracy and stability in the overall simulation. This procedure corresponds to the general problem of finding a set of basis functions with their weights, such that their linear combination is able to reconstruct the field variable accurately. We will refer to this procedure as particle field interpolation, or simply field interpolation.

When a particle method is used as a numerical discretization of a continuum problem, one represents the field variable as a linear combination of localized basis functions (the particles). Overlap of the particles is necessary to accurately reconstruct a continuous field variable, a fact which has been amply recognized in the vortex method literature where convergence results rely on the overlap condition [7]. In practice, it has been observed that the spatial accuracy in the numerical representation of the continuous vorticity field using the vortex particles increases exponentially with decreasing overlap ratio, defined as the distance between particles divided by the characteristic size [5]. In a Lagrangian approach (moving basis functions), the initial particle configuration can become severely strained or sheared, leading to a deterioration of spatial accuracy in the continuum representation. Thus, most if not all investigators face the need of remediating configurations of particles so that the basis functions continue to overlap in space throughout the simulation.

While various investigators have developed techniques for replacing a deformed configuration of particles with one that is uniform or nearly so, this aspect of particle method simulation is seldom addressed in the literature. Possibly, spatial and temporal errors have dominated particle calculations and errors induced by the particle redistribution has been of lesser concern. Currently, vortex methods involving elliptical, deforming basis functions are available, providing fourth-order accuracy in space [31]. The accurate replacement of one particle configuration with another acquires greater relevance when maintaining accuracy in combination with such a method.

A central issue in particle field interpolation is that the associated linear systems of equations are ill-conditioned. In one of the first attempts to avoid inverting the ill-conditioned system, Beale proposed applying a residual correction method to interpolation problems with fourth-order kernels, which relies on having square lattices of particles [6]. Computational experiments using Beale’s approach were reported in [11], here using Gaussian kernels on a fixed mesh with changing values of the core size. The importance of preconditioning in field interpolation was later recognized in [21]. However, tensor-product interpolation is a popular alternative that avoids solving the linear system altogether. In this method, a disordered set of particles is replaced with a set of uniformly distributed ones by projecting particle weights (circulation) onto a uniform mesh. The procedure is called remeshing or particle redistribution or re-initialization by different authors. One of its first successful uses was in the simulation of cylinder flow in [18], and the approach has since been at the heart of many important results with vortex particle methods, such as a recent calculation of aircraft wakes with billions of particles [10]. Alternative spatial adaptation techniques that have been introduced include the splitting and merging of particles [29], [30]. However, the use of tensor-product interpolations continues to be the state-of-the-art in vortex methods. Some evidence is available suggesting that remeshing could pose an accuracy limitation on the overall calculations (experiments in [4] show a jump in the errors upon the first remeshing process), and for that reason we are interested in the alternatives. One alternative that has been recently proposed and carefully studied [5] is to re-project the field variable, which in the case of vortex methods is the vorticity, using the particles as interpolation centers. The technique is based on a radial basis function (RBF) formulation of the problem, but the issue of ill-conditioning was not investigated.

The motivation of the present work is providing a spatial adaptation technique appropriate for use with fourth-order viscous vortex methods with elliptically deforming bases [31]. However, the results are general and applicable to any interpolation problem utilizing Gaussian radial basis functions. The paper begins with an overview of the vortex method, and an explanation of how the present work advances us toward our goal. In Section 3, we overview some relevant aspects of the theory of radial basis function interpolation and issues related to the solution of the ill-posed reverse heat equation. In Section 4, we present methods of solution and numerical experiments for the field interpolation problem applied to the vortex particle method. We use mainly two experimental setups: one is a passive scalar streak which has been deformed by a differentially rotating flow, and evolved using the fourth-order elliptical vortex method [32]; the second is the high-gradient but smooth vortex patch of [22]. These two experimental setups are used to demonstrate several approaches to field interpolation, including: solving the dense linear system associated with an RBF interpolation problem using a preconditioned iterative method; and, using the reverse heat equation to deconvolve an initial guess obtained using local vorticity and particle volumes. We end with conclusions and indications of the questions left for future work.

Section snippets

Vortex methods and representing fields with smooth particles

The vortex method solves the Navier–Stokes equation at constant density in vorticity formulation, namely, the vorticity transport equation:ωt+u·ω=ω·u+1ReΔω.In the case of a two-dimensional and inviscid flow the right-hand-side of (1) is zero and the governing equation reduces to the simple form DωDt=0, where DDt stands for the material derivative. One can immediately see how a Lagrangian approach is ideal for such a problem: to satisfy this equation it suffices to allow small computational

Solving global field interpolation for vortex particle representations

In this section, we will develop two techniques for particle field interpolation: preconditioned radial basis function collocation, and a deconvolution formulated with the reverse heat equation. We will compare these techniques with the common low-order practice that initializes the particle strengths as γi=ωih2. One common technique that has not been assessed in this paper is the use of tensor-product formulations with splines-based kernels. In the vortex method community, most workers use an

Two model problems

To study the efficacy of the two field interpolation methods suggested, we will use two test fields. The first test field consists of data obtained from a numerical calculation of a convection–diffusion problem, using the elliptical vortex method [32]. Thus, the field is the superposition of deformed elliptical Gaussian basis functions, with aspect ratio a. The elliptically deforming vortex particles have been evolved for a number of time steps, at the end of which the greatest aspect ratio is

Conclusions

We have studied the important process of field interpolation for Gaussian basis functions. Field interpolation is a general technique, but it plays a special role in the computational ecosystem surrounding vortex methods. Field interpolation is necessary for initialization of particles in flow simulations and for re-initializing distorted configurations of particles during long-time simulations. We have identified two techniques that provide results that are orders of magnitude more accurate

Acknowledgments

LAB acknowledges support from EPSRC (Engineering and Physical Sciences Research Council, UK) under Grant contract EP/E033083/1, and additional travel funds from Boston University College of Engineering. The authors thank the anonymous referees for their helpful comments and suggestions.

References (35)

  • L.A. Barba, A. Leonard, C.B. Allen, Numerical investigations on the accuracy of the vortex method with and without...
  • L.A. Barba et al.

    Advances in viscous vortex methods – meshless spatial adaption based on radial basis function interpolation

    Int. J. Numer. Meth. Fluid

    (2005)
  • J.T. Beale

    On the accuracy of vortex methods at large time

  • J.T. Beale et al.

    Vortex methods I: convergence in three dimensions

    Math. Comput.

    (1982)
  • M.D. Buhmann

    Radial Basis Functions. Theory and Implementations

    (2003)
  • R.H. Chan et al.

    Cosine transform-based preconditioners for total variation deblurring

    IEEE Trans. Image Proc.

    (1999)
  • J.P. Choquin et al.

    Accuracy of a deterministic particle method for Navier–Stokes equations

    Int. J. Numer. Meth. Fluid

    (1988)
  • Cited by (16)

    • RBF-vortex methods for the barotropic vorticity equation on a sphere

      2015, Journal of Computational Physics
      Citation Excerpt :

      An active research frontier is to accelerate RBF interpolation. Cost can also be greatly reduced by using fast multipole methods and treecodes such as [1,50,81] to sum the “far field” parts of RBF series. Unfortunately, treecodes are not accelerative for interactions between neighboring RBFs [15].

    • Periodized radial basis functions, part I: Theory

      2014, Applied Numerical Mathematics
    • FMM-based vortex method for simulation of isotropic turbulence on GPUs, compared with a spectral method

      2013, Computers and Fluids
      Citation Excerpt :

      Since the Gaussian function decays rapidly, we use the FMM neighbor list to calculate Eq. (3) between neighboring particles only. For an overlap ratio of h/σ = 1, where h is the average distance between particles, the RBF system is well conditioned and converges in 5–10 iterations [3,22]. The six FMM kernels and the matrix–vector multiplication in RBF interpolation are done on the GPU.

    • A multi-moment vortex method for 2D viscous fluids

      2012, Journal of Computational Physics
      Citation Excerpt :

      To continue to ensure accuracy for later times, we will need to incorporate spatial adaptation of the vortex particles, for example by using the highly accurate methods recently developed in [5] or the higher order interpolation methods of PIC [44].

    • PetRBF - A parallel O(N) algorithm for radial basis function interpolation with Gaussians

      2010, Computer Methods in Applied Mechanics and Engineering
      Citation Excerpt :

      The fast convergence of this small but non-local Gaussian basis function can be achieved only with the use of rasm with carefully truncated matrix-vector operations. Indeed, the estimated condition numbers in the present calculations are in the order of 104 for h/σ = 1.0 and 106 for h/σ = 0.8, consistent with those reported in [7]. There are two common measures for evaluating the performance of parallel codes: the strong scaling and the weak scaling.

    View all citing articles on Scopus
    View full text