An efficient trajectory tracking algorithm for the backward semi-Lagrangian method of solving the guiding center problems

Highlights

• Construct a completely exlicit formula for the solution of each Cauchy problem in the backward semi-Lagrangian method.

• Construct an effective fast algorithm for the departure points.

• Propose a method to modify the solution to improve the physical quantities such as conservation of mass.

Abstract

In this paper, we develop an effective numerical algorithm that tracks the trajectory needed to solve guiding center models using the backward semi-Lagrangian method. In terms of numerical calculations, two appreciably fast algorithms for the departure points are designed. One is a completely explicit formula for numerical solutions of the discrete system for each Cauchy problem. This formula is characterized by numerical factors that are less than half the multiplication number used in the usual Gaussian elimination. The other finds the required departure points with an interpolation method that is at least 30% less expensive for heavy Cauchy problems. Last, we propose a method to modify the solution to improve estimation of physical quantities, such as conservation of mass, which can be lost in interpolation solutions calculated at the departure points. It turns out that the proposed method not only saves a great deal of computation time, but also preserves physical quantities such as mass and total kinetic energy much better than conventional methods. To demonstrate the numerical evidence, we use the proposed method to simulate several problems such as the incompressible Euler equation, Kelvin-Helmholtz instability, Diocotron instability and a three-dimensional guiding center model.

Introduction

The model problem we are interested in solving is the guiding center model, which was developed to efficiently describe low-frequency turbulence and resulting transport phenomena in strongly magnetized plasma. Assuming that the uniform external magnetic field is perpendicular to a plane, then the density of the guiding centers of charged particles can be modeled by a two dimensional hyperbolic equation, whose advection field is defined by an electrostatic potential function whose self-consistency satisfies the Poisson equation [17], [19]. The forcing term of the Poisson equation is characterized by the density of the guiding centers and results in a very strongly nonlinear problem having so-called self-consistency. The incompressible Euler equation as characterized by the vorticity-stream function, is also equivalent to the model problem if the sign of the Poisson equation is ignored. Moreover, the nonlinear structure is generic in descriptions of geophysical problems [5], [26], [33], [30] and the dynamics of plasma [27], [17], [19]. Hence, the model problem with which we are concerned may suggest the development of widely applicable numerical schemes. In addition, it serves as an ideal test problem for our newly suggested method.

Among many numerical methods for hyperbolic problems, the semi-Lagrangian method has long been used in a wide range of fields such as meteorology for weather prediction, oceanography, and more general problems in fluid dynamics [1], [2], [3], [4], [10], [21], [24], [25], [26], [28], [32], [33]. This approach was introduced at the end of the 1950s [31] and uses a fixed computational Eulerian grid. At each integration step, a discrete set of particles arriving at the grid points is tracked backwards over a single time step along its characteristic curve up to its set of departure points. The solution value at the departure point is then calculated by interpolating the solutions of the previous time step. The backward semi-Lagrangian method (BSLM), which implicitly treats the characteristic curve, is well known to use a large time step size, since there is no limitation due to the Courant-Friedrichs-Lewy (CFL) conditions by contrast with the explicit Eulerian scheme. Also, the BSLM has no mesh deformation as in pure Lagrangian methods, owing to the conformity of the arrival points with the grid points. According to the convergence and stability analysis found in [8], [17], the BSLM's accuracy depends on two decisive issues. One issue is the interpolation scheme for the calculation of the density and the advection field at departure points. There are many standard interpolation methods, such as cubic spline and Hermite interpolation [7] which includes cubic Lagrange interpolation. In addition, efforts have been made to develop a non-oscillatory interpolation method such as the WENO (weighted essentially non-oscillatory) method [22], [23].

The second critical issue is backward integration for the nonlinear Cauchy problem, the main feature of the BSLM. The Cauchy problem has the same self-consistency property as the model problem because the advection field described by the potential function is given as the velocity of the characteristic curve. Hence, it is difficult to use an explicit Runge-Kutta (RK) method directly since it gives poor stability and is not completely free from the CFL condition. In recent years, there has been a sustained effort to seek ever more efficient and stable numerical methods for the Cauchy problem. Backward integration is commonly solved using implicit schemes, as can be seen in the literature [17], [19], [32], [33]. Among such implicit schemes, there are two types of methods with third-order temporal accuracy; the multistep methods based on the Adams-Moulton method (AM) [9] and the three-stage one-step method recently developed by the authors [19]. Both methods were developed to solve transport equations for plasma physics. The multistep method [9] turns out to have poor stability and imposes constraints on the use of larger time step sizes compared to the one-step method. In contrast to the multistep method, the one-step method [19] with third-order accuracy and good stability involves an algorithm that is completed by finding the departure points at two intra-stages. Hence, it has a weak point in that the Poisson equation must be solved two more times at these intra-stages in each integration step, which requires a high computational cost to obtain the solution of the Poisson equation and the advection field.

As in the above mentioned research, various aspects of approximating the departure points require improvement. Let $cpu{t}_C$ be the computational costs of solving the Cauchy problem and finding the required departure points. Costs captured by $cpu{t}_C$ have two components and the main contribution of this paper is to develop a numerical approach to these components that reduces the total value of $cpu{t}_C$. The first component comes from solving the asymptotic linear system induced by the error correction technique [15], [17], [19]. We present a multistep method, inspired by our recent results [19], for backward integration of the Cauchy problem, that saves the computational cost associated with the advection fields. Also, we design a new explicit formula for the solution of the linear system that is quite flexible to implement, compared to that recently developed by the authors [14]. This explicit formula turns out to reduce the total number of multiplication by more than half. The second component is the discretization of the Cauchy problem, which requires heavy and expensive interpolation calculations. To resolve this problem, we develop a new method based on interpolation theory using some already calculated departure points to find the remaining ones, rather than using the solution of the Cauchy problem. For this purpose, we propose a method of comparing interpolation polynomials of different degree to decide whether it is reasonable to use an interpolation polynomial instead of the solution of the Cauchy problem. More precisely, for a given grid point ${\mathit{x}}_{i,j}$, we introduce an indicator $I_{i,j}$ defined by

$$I_{i,j}=3\Vert P_3({\mathit{x}}_{i+{\vartheta }_i(1),j+{\vartheta }_j(1)})-P_2({\mathit{x}}_{i+{\vartheta }_i(1),j+{\vartheta }_j(1)})\Vert ,$$

where $P_3$ and $P_2$ are cubic and quadratic interpolation polynomials, respectively, and $\Vert \cdot \Vert$ denotes the discrete $l_2$-norm. According to the magnitude of the indicator, we split the set of the departure points we have to find into two sets, say A and B. The departure points in set A are obtained by solving the Cauchy problem, while the departure points in set B are determined by evaluating the cubic interpolation polynomial $P_3$ that is constructed from set A. This overcomes the heavy computational cost of the interpolation procedure in the Cauchy problem. It turns out that the total number of Cauchy problems we have to solve at each integration step are reduced to less than 69% of their former number. Summarizing, the suggested procedures reduce the total the computational cost; we provide evidence for this in several numerical simulations.

Finally, we discuss how to improve the numerical solution in order to eliminate the possibility that the solution based on the departure points obtained by simple interpolation, as opposed to the solution of the Cauchy problem, could lose fidelity to physical properties such as the conservation of mass. In order to modify the solution, we partially correct the errors of the approximated solutions whose departure points are calculated using interpolation. For the error correction, we estimate the mean of the errors using the difference between the initial mass and that calculated from the approximated solutions. This error correction method is at least numerically verified to significantly improve the conservation of mass with only a slight difference in both accuracy and computational cost compared to the original unimproved method. In addition, it can be shown that the proposed method is superior to the existing Adams-Moulton (AM3) and three-stage one-step collocation (CM3) methods in terms of conservation of mass and kinetic energy.

The rest of this paper is organized as follows. In Section 2, we introduce a third-order multi-step method for the Cauchy problem by modifying the one-step collocation method [19]. Section 3 is devoted to developing a completely explicit and easily implemented formula for the departure points. In Section 4, we discuss the interpolation method for efficiently implementing the departure points without using solutions of Cauchy problems, for which we develop and analyze the indicator $I_{i,j}^{⁎}$ introduced above. Section 5 discusses how to improve the numerical solution so that modeling of physical properties such as conservation of mass is improved. Numerical experiments of the proposed scheme together with several comparisons are performed in Section 6. Finally, Section 7 concludes the paper.