An efficient trajectory tracking algorithm for the backward semi-Lagrangian method of solving the guiding center problems
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 be the computational costs of solving the Cauchy problem and finding the required departure points. Costs captured by 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 . 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 , we introduce an indicator defined by where and are cubic and quadratic interpolation polynomials, respectively, and denotes the discrete -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 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 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.
Section snippets
Semi-Lagrangian method for the guiding center model
We consider the guiding center model described by where , is a bounded domain, is a given smooth initial function, and the advection field is defined by for the self-consistent electrostatic potential function φ. The density and the potential are coupled to each other in relations (2.1) and (2.2). Hence, the problem we have to solve
Explicit formula for
We devote this section to designing an explicit formula for the solution Ψ for the system (2.16), for which the properties of the Jacobian matrix defined in (2.8) are used. We begin the section by introducing a useful property of the Jacobian matrix . Lemma 2 For the Jacobian matrix defined in (2.8), we have where Proof One can easily complete the proof from the definition of the Jacobian matrix . □
An interpolation method for the departure points
The aim of this section is to introduce a numerical procedure to overcome the computational complexity caused by the discretization of the Cauchy problem discussed in the end of Section 2. We first explain our motivation as follows. We use the difference between a cubic polynomial interpolating four given departure points and a quadratic polynomial interpolating three departure points among the same four points as shown in Fig. 1. The figures make clear that the bigger the magnitude of the
Conservation of mass
The aim of this section is to develop a method for updating the approximate solution by improving mass conservation to compensate for an accuracy that a solution created from interpolation technique discussed in Section 4, rather than a solution based on the departures of all Cauchy problems, can be lost. For convenience of discussion, we begin this section with some notation. Let be the set of total indices of the grid points. For the sets A and B of the departure points
Numerical examples
In this section, we demonstrate the numerical accuracy and efficiency of the proposed numerical algorithm for tracking the trajectory of the characteristic curve. For the accuracy and the convergence, we calculate a -norm error defined by where are the spatial grid sizes, is the numerical approximation of representing the reference solution at time and grid point . When the problem we solve does not have
Conclusion
We developed a third-order multi-step collocation method for tracking the trajectory of the characteristic curve in the backward semi-Lagrangian method. To reduce computational costs, a new completely explicit formula for the solution of the Cauchy problem was designed using the divergence-free condition for the velocity. Also, to overcome the heavy computational cost in discretizing the Cauchy problem, we proposed an interpolation scheme instead of the direct solver for the Cauchy problem;
Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Acknowledgements
This research was supported by the basic science research program through the National Research Foundation of Korea (NRF) funded by the Korea government (MSIT) (2020R1A2C1A01008506). Also, the first author Piao was supported by the National Research Foundation of Korea (NRF) grand funded by the Korea government (MSIT) (NRF-2020R1H1A1006982) and the second author Kim was partially supported by the National Research Foundation of Korea (NRF) grand funded by the Korea government (MSIT) (
References (33)
- et al.
A high order characteristic method for the incompressible Navier-Stokes equations
Comput. Methods Appl. Mech. Eng.
(1994) - et al.
High order time discretization for backward semi-Lagrangian methods
J. Comput. Appl. Math.
(2016) - et al.
A second-order time-accurate ALE Lagrange-Galerkin method applied to wind engineering and control of bridge profiles
Comput. Methods Appl. Mech. Eng.
(2004) - et al.
A completely explicit scheme of Cauchy problem in BSLM for solving the Navier–Stokes equations
J. Comput. Phys.
(2020) - et al.
An iteration free backward semi-Lagrangian scheme for solving incompressible Navier-Stokes equations
J. Comput. Phys.
(2015) - et al.
One-step -stable temporal integration for the backward semi-Lagrangian scheme and its application in Guiding center problems
J. Comput. Phys.
(2018) - et al.
A conservative high order semi-Lagrangian WENO method for the Vlasov-equation
J. Comput. Phys.
(2010) - et al.
Conservative high order semi-Lagrangian finite difference WENO methods for advection in incompressible flow
J. Comput. Phys.
(2011) - et al.
The semi-Lagrangian method for the numerical resolution of the Vlasov equation
J. Comput. Phys.
(1999) - et al.
Multi-dimensional conservative semi-Lagrangian method of characteristics CIP for the shallow water equations
J. Comput. Phys.
(2009)
A semi-Lagrangian high-order method for Navier-Stokes equations
J. Comput. Phys.
Numerical analysis of convection-diffusion-reaction problems with higher order characteristics/finite elements. Part I: time discretization
SIAM J. Numer. Anal.
Numerical analysis of convection-diffusion-reaction problems with higher order characteristics/finite elements. Part II: fully discretized scheme and quadrature formulas
SIAM J. Numer. Anal.
A high-order characteristics/finite element method for the incompressible Navier-Stokes equations
Int. J. Numer. Methods Fluids
Comparison of semi-Lagrangian discontinuous Galerkin schemes form linear and nonlinear transport simulations
A high order semi-Lagrangian discontinuous Galerkin method for the two-dimensional incompressible Euler equations and the Guding center Vlasov model without operator splitting
J. Sci. Comput.
Cited by (4)
Development of a parallel CUDA algorithm for solving 3D guiding center problems
2022, Computer Physics CommunicationsCitation Excerpt :In other words, register memory and local memory must be allocated to perform the kernel function used in each method, and the kernel function may fail to execute if the block size is too large. It is well known that the BSLM (backward semi-Lagrangian method), including the recently introduced ECM3 [19], is suitable for parallel computing since it is based on the departure points of the characteristic curve where their trajectories are solved independently. This paper developed a parallel CUDA algorithm for the ECM3.
One- and two-step semi-Lagrangian integrators for arbitrary Lagrangian–Eulerian-finite element two-phase flow simulations
2022, International Journal for Numerical Methods in FluidsTennis Serve Trajectory Capture Algorithm Based on Wavelet Multiscale Decomposition
2022, Mathematical Problems in EngineeringTrajectory Tracking Method of Volleyball Player's Arm Hitting Image Based on D-P Algorithm
2021, Scientific Programming