{'fraction_non_alphanumeric': 0.029077911413804763, 'fraction_numerical': 0.015118581848041346, 'mean_word_length': 4.604764482945317, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Development of a symplectic and phase error reducing perturbation finite-difference advection scheme

ABSTRACT The aim of this work is to develop a new scheme for solving the pure advection equation. This scheme formulated within the perturbation finite-difference context not only conserves symplecticity but also preserves the numerical dispersion relation equation. The employed symplectic integrator of second-order accuracy in time enables calculation of a long-time accurate solution in the sense that the Hamiltonian is conserved at all times. The generalized high-order spatially accurate perturbation difference scheme optimizes numerical phase accuracy through the minimization of the difference between the numerical and exact dispersion relation equations. Our proposed new class of phase error reducing perturbation difference schemes can in addition locally capture discontinuities underlying the concept of applying a shope/flux limiter. The high-order spatial accuracy can be recovered in a smooth region. Besides the Fourier analysis of the discretization errors, anisotropy and dispersion analyses are both conducted on the dispersion-relation and symplecticity-preserving pure advection scheme to shed light on the distinguished nature of the proposed scheme. Numerical tests are carried out and the results compare well with the exact solutions, demonstrating the efficiency, accuracy, and the discontinuity-resolving ability using the proposed class of high-resolution perturbation finite-difference schemes.


Introduction
In science and engineering areas, many transport processes involve various length scales. Development of a scheme of high-order spatial and temporal accuracy is essential to resolve highwave-number (or high-frequency) and small-amplitude wave components. The need for developing a high-order scheme is particularly necessary, for example, while performing a direct simulation of flow turbulence, in investigating transition in wall-bounded shear flow, and in modeling global weather. In addition to yielding a high order of truncation, it is also desirable to apply a high resolution numerical scheme, which is characterized by employing a large number of grid points per wavelength, to resolve the shortest-wave component. Good-quality solution can be therefore obtained at low computational cost by using this class of high-order numerical schemes for small length scales.
In the past three decades, a tremendous effort has been devoted to develop high-order schemes with high resolution. In these numerical schemes, the spectral and compact schemes are often referred to. In contrast to the finite-difference and finite-element methods developed normally under the comparative low-order representation of the functions in a small subdomain, spectral methods A numerically less explored high-order perturbation finite-difference (PFD) method was proposed by Gao in 2000 [21]. His underlying idea of developing this novel method [22,23] has been extended to the finite-volume method [24][25][26] in flow simulations. This class of flux resolution schemes modifies the classical convective and diffusive numerical fluxes by introducing a series of perturbation terms, which are functions of the power of grid spacing. The convection-diffusion differential equation under approximation is solved using the proposed numerical perturbation algorithm aiming at eliminating nonphysical oscillations resulting from the numerical convective instability and developing a high-order scheme without increasing the number of stencil points. The objective of this work is to refine the original PFD scheme to reconstruct a new set of coefficients which not only yields higher solution accuracy but also get improved dispersion characteristics. High resolution of the scheme is also addressed so that the number of grid points per wavelength required to resolve the shortest error component in the practical computation can be reduced as much as possible.
The rest of this article is organized as follows. In Section 2, the mathematical nature of the working equation is reviewed. The symplectic PFD scheme characterized by having the optimized dispersion relation equation (DRE) in a given grid stencil is developed in Section 3. Theoretical analysis of this optimized difference scheme is detailed in Section 4. The discontinuity-capturing perturbation finitedifference scheme accommodating the optimized numerical dispersion relation equation (DC-PFD2-DRE) is presented in Section 5. The efficiency, accuracy, and convergence characteristics of the proposed pure advection scheme are demonstrated in Section 6. Finally, concluding remarks are given in Section 7 based on the simulated results.

Working equation and solution algorithm
The linear pure advection equation considered as the modal equation to develop the proposed numerical scheme is given as In the above hyperbolic equation for /(x, y, t), the velocity components u, v in the x and y directions are assumed to be constant for the sake of simplicity to describe the proposed scheme. Given an initial condition /(x, y, t ¼ 0), the solution will be sought subject to its preserved symplectic property. This strong demand on the predicted solution is essential to get a long-term accurate solution since Eq. (1) has been well known to be one of the important members of the Hamiltonian equations.
While solving a multidimensional differential equation, we can adopt a class of time-spiltting methods to attain a good-quality solution at low computational cost. The motivation for the use of time-splitting methods in the numerical calculation of Eq. (1), for example, can avoid a less restrictive stability condition while using an explicit scheme. Provided that an implicit scheme is employed, it is our objective to retain the tridiagonal nature of the matrices so as to reduce the computational cost in solving the algebraic system of equations. The time-splitting method given below will be adopted in this study to solve the time marching solution iteratively from the following two fractional steps:

Numerical method
In light of Eqs. (2)-(3), the overall performance of applying the proposed time splitting method depends very much on the numerical method applied to solve the following model equation:

Symplectic temporal discretization
The hyperbolic equation (4) is endowed with the conserved Hamiltonians. This fundamental property embedded in the equation leads naturally to the adoption of a symplectic temporal scheme. Given the need for designing a scheme that can preserve symplecicity, in this study the following second-orderaccurate symplectic scheme for the model equation (4) is employed, leading to

Third-order-accurate perturbation finite-difference scheme (PFD3)
The main idea behind developing the perturbation finite-difference scheme is to eliminate the leading discretization error terms from the modified equation and therefore get high-order accuracy. If the first-order upwinding scheme is used to discretize the advection term in (5), the resulting discretized equation is written as follows for u > 0: When constructing the corresponding perturbation difference scheme, the convective velocity u in (6) is replaced with u p , thereby leading to In the above, u p (¼u þ A 1 Δx þ A 2 Δx 2 ) is approximated as the power series of the grid spacing Δx. The coefficients A 1 and A 2 will be determined below.
The terms / nþ1 i , / nþ1 iÀ 1 , and / n iÀ 1 are expanded first in Taylor series with respect to / n i , respectively, at node i. These expansion terms are then substituted into the discrete equation (7) to get the following modified equation after a lengthy derivation: Dt 2 þ ð:::ÞDt 3 þ ð:::ÞDt 4 þ ::: Two leading discretization error terms (Δx n terms, n ¼ 1, 2) shown in the above modified equation are eliminated to get A 1 ¼ (u/2)/ xx // x and A 2 ¼ À (u/6)/ xxx À (1/2)A 1 / xx // x . By virtue of the approximation of derivative terms as the perturbation finite-difference scheme yielding third-order accuracy can then be derived according to the resulting modified equation.

Second-order-accurate perturbation finite-difference scheme with the numerical dispersion relation equation (PFD2-DRE)
While approximating the propagation equation it is essential to reduce the cumulative type of numerical dispersion error. In this study we aim to minimize the difference between the exact and numerical dispersion relation equations for the linear equation / t þ u / x ¼ 0. Note that the dispersion relation is used to rigorously link the angular frequency ω of a wave with the wavenumber k.
To derive the numerical dispersion relation equation, the plane wave solution cast in the form of e i(ωtÀ kx) \ for / t þ u / x ¼ 0 is used to get the difference between the exact and numerical dispersion relation equations. To this end, derivation of the last algebraic equation for the linear equation can be transformed to its corresponding algebraic expression in (ω, k), where ω and k denote the angular frequency and wavenumber, respectively. Through the transformation between the (x, t) and (ω, k) planes, the exact dispersion relation (or Von Neumann stability relation) equation is derived as ω À uk ¼ 0.
To improve dispersive accuracy for the adopted PFD scheme (7), we refine the scheme derived in Section 3.2 by rewriting / xxx as ðc 1 =2Dx 3 Þ (/ iþ2 À 2/ iþ1 þ 2/ iÀ 1 À / iÀ 2 ). By substituting the discrete into Eq. (7), the numerical dispersion relation equation for the equation / t þ u / x ¼ 0 can be derived as In the above derived numerical dispersion relation equation, the numerical angular frequency ω num is no longer linearly proportional to the wavenumber k. This derivation clearly enlightens that numerical dispersion error has been introduced, and it can definitely affect the predicted solution.
One can directly reduce the dispersion error by decreasing the difference between the exact and numerical dispersion relation equations or indirectly by minimizing the difference between the numerical and exact phase velocities or group velocities. In this study the chosen error function (ω num À ω exact ) 2 is minimized in an integrated sense as follows within the proper interval of the scaled wavenumber: By applying the limiting condition ð@E=@c 1 Þ ¼ 0, the free parameter c 1 introduced into the approximation of / xxx can then be derived from the above error minimization procedure. The resulting derived coefficient c 1 will not be shown explicitly in the article because of the derived complex functions of Cr ð� uDt=DxÞ. As an example, the values of c 1 are 1.25, 1.45, and 1.95 are used for the cases of Cr ¼ 0.02, Cr ¼ 0.1, and Cr ¼ 0.5, respectively.
It is worth nothing here that the value of c 1 is not fixed throughout the computational domain but is rather varied locally according to the magnitude of the flow velocity. Given the users' prescribed Δt and Δx, the local Courant number Cr is determined from the local convective velocity u. Use of the resulting locally varying weighting coefficients enables us to obtain a solution that can theoretically optimize the numerical dispersion relation equation at every grid point. The numerical angular frequency ω num is plotted with respect to the wavenumber k in Figure 1 at different Courant numbers using the currently proposed PFD method. One can see from the figure that the PFD2-DRE scheme has a much better numerical dispersion relation than that of the PFD3 scheme at all the scaled Courant numbers.

Fourth-order-accurate perturbation finite-difference scheme with the numerical dispersion relation equation (PFD4-DRE)
To increase the prediction accuracy, the fourth-order-accurate perturbation finite-difference scheme with numerical dispersion relation equation is developed below. By performing Taylor series expansion on the terms / nþ1 i , / nþ1 iÀ 1 , and / n iÀ 1 with respect to / n i , and then substituting these expansion terms into the discrete equation (7), the modified equation then can be derived as The leading three discretization error terms are then eliminated from the above modified equation, thereby leading to Derivation of scheme PFD4-DRE is followed by introducing the free parameter c 2 for the approxi- 5 . We then minimize the resulting error function E shown in Eq. (10). Following the same procedure described in Section 3.3, the values of c 2 are 1.1, 1.3, 1.5, and 1.7 derivaed for the cases of Cr ¼ 0.02, Cr ¼ 0.1, Cr ¼ 0.3, and Cr ¼ 0.5, respectively. The resulting modified equation analysis enables us to know that the scheme has formal accuracy order of four.

A class of high-order perturbation finite-difference schemes
Following the same idea described in Sections 3.2-3.4, a generalized PFD scheme can be derived by perturbing u as follows: The resulting derived A i for the proposed upwinding PFD schemes with/without the optimization of the numerical dispersion relation equation can be derived below: The coefficients A i shown above are summarized in Table 1 for N ¼ 1À 4.

Dispersion analysis
Given that k 2 ¼ k 2 x þ k 2 y , both wavenumber components can be expressed as k x ¼ k cos θ and k y ¼ k sin θ. For the sake of comparison and discussion of the computed results, three extra parameters need to be defined. The first one is N λ ¼ λ/h, which denotes the number of points per wavelength λ ( ¼ 2 π/ k). The other two are the CFL numbers Cr x ¼ u Δ t/h ¼ 0. 2 and Cr y ¼ vDt h ¼ 0:2. Based on these definitions, the numerical phase velocity V p , which is defined as the ratio of the numerical angular frequency and the wavenumber k, can be derived. The ratio of the numerical phase velocity The numerical phase velocity V p1 for the PFD3 scheme can be also rewritten in terms of θ. Two ratios, V p /V e and V p 1 /V e are plotted with respect to the angle θ in Figure 4 for N λ ¼ 4. One can clearly  observe from Figure 4 that the PFD2-DRE scheme has better agreement between the numerical and exact dispersion relation equations.

Discontinuity-capturing perturbation finite-difference scheme (DC-PFD2-DRE)
When solving the convection-dominated flow equations, in order to obtain a high-order accuracy and at the same time to resolve discontinuity, the resulting discretized equation is written as follows for the case of u > 0: where The value of / ½/ iþ1 À / iÀ 1 =/ i À / iÀ 1 � shown above will be determined according to the monotonicity-maintenance criteria. The consequence is that the resulting nonlinear scheme can eliminate the nonphysical oscillatory solutions using the PFD2-DRE scheme in the region across discontinuities and retain high-order accuracy in smooth region. For u < 0, the proposed discontinuity-capturing PFD scheme can be similarly derived.

One-dimensional verification studies
The solution of the linear equation u t þ u x ¼ 0 is sought subject to the following initial condition: In this one-dimensional verification study, the periodic boundary condition is imposed for u. The time step is chosen to be Δt ¼ 0.02Δx. Figure 5 shows the exact waveform and the waveform obtained To validate that the proposed symplectic and dispersion error reducing perturbation finite-difference scheme is applicable to solve problems with nonsmooth solution, the solution of the same linear advection equation is sought subject to the initial waveform of the discontinuous type given below: The above equation is solved subject to the specified periodic boundary condition u(0, t) ¼ u(4, t). Figure 6 shows the predicted waveforms at t ¼ 1 for the case investigated at the constant grid and time spacings Δx ¼ 1/50 and Δt ¼ 0.02Δx, respectively. Since the proposed PFD2-DRE scheme is not classified to be monotone, the predicted oscillations near the root of the investigated square wave is under our expectation. Provided that the discontinuity-capturing slope limiter is used, one can see from Figure 6b that the proposed DC-PFD2-DRE scheme can capture the sharply varying solution quite well. One can also clearly see from Figure 6b that the predicted wave shape remains almost unchanged.  The linear advection equation u t þ u x ¼ 0 is also solved under the following initial condition [27][28][29][30]: In Eq. (21), G(x, z) ¼ e À β(x À z) 2 and Fðx; aÞ ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi Maxð1 À a 2 ðx À aÞ 2 ; 0Þ q . The constants are taken as a ¼ 0.5, z ¼ À 0.7, δ ¼ 0.005, α ¼ 10, and β ¼ log(2)/36Δ 2 . This investigated initial condition consists of a discontinuous square pulse and several continuous profiles. Figure 7 plots the solutions predicted from the PFD2-DRE and DC-PFD2-DRE schemes at time t ¼ 2 in 400 uniformly distributed grids. Under our expectation, the DC-PFD2-DRE scheme can completely resolve oscillations near the jumps.
The solution of the above two-dimensional pure advection equation is sought in a square of unit length, where the vortex flowfield is generated by the prescribed velocity vector (u, v) given by [32] In the above, T is the time required for the solution to return back to its initial solution. A circle of radius 0.15 centered at (0.5, 0.75) is prescribed initially. This problem, known as the vortex flow in a box, has been known to be computationally challenging in the sense that the solution profile will be stretched and torn by the prescribed vortex flow. A very thin filament with the width of a mesh size can be generated with time. This test problem is solved in 150 � 150 grids at t ¼ 0.5, 1.0, 1.5, and 3. The computed results match very well with the exact solutions and are plotted in Figure 8. Note that the maximal stretching occurs at t ¼ T ¼ 3.
The Zalesak's problem [33] is another well-known benchmark case for testing the ability of the developed multidimensional advection scheme. A slotted disk of radius 15.0 and width of 5.0 is initially located at the (50, 75) in a square domain [100,100]. The prescribed velocity field is given as The slotted disk under investigation rotates clockwise. The results predicted in a domain of 250 � 250 and 400 � 400 mesh points are plotted in Figure 9a and Figure 9b, respectively, at a time after completing one revolution. Our computed solution is seen to be in good agreement with the exact (or initial) solution, as shown in Figure 9b. Figure 10a shows the solution profile /(x, 70) along the x axis for 0 � x � 100 using the DC-PDF2-DRE scheme in 400 � 400 grids. Figure 1b shows again that the proposed scheme has excellent Hamiltonian preservation property.

Two-dimensional incompressible flow problem
The Taylor vortex flow problem regarding the transport of decaying Taylor vortices is frequently used as the benchmark test for the verification of numerical schemes developed to solve the convectiondominated incompressible viscous flow equations [34]. The initially periodic vortex flow is convected and decayed exponentially with time due to the nonzero fluid viscosity in the flow. The exact solution of this problem is varied with the Reynolds number Re as follows: uðx; y; tÞ ¼ À cosðxÞ sinðyÞ expð À 2 t Re Þ vðx; y; tÞ ¼ cosðyÞ sinðxÞ expð À 2 t Re Þ; pðx; y; tÞ ¼ À 0:25½cosð2xÞ þ cosð2yÞÞ expð À 4 t Re �

<
: The results computed in a mesh of 250 � 250 uniform grids are presented at Re ¼ 100. Figures 11a  and 11b show good comparison between the computed and exact u-and v-velocity components along the y-and x-axis lines at different times.

Concluding remarks
In this study a generalized perturbation finite-difference scheme has been developed for solving the pure advection equation in a single dimension. This scheme of considerable novelty has been easily extended to predict the two-dimensional pure advection equation with greater success. Besides the ability of predicting a very high spatial accuracy using the generalized pure advection PFD scheme of varying accuracy orders, we also aimed to compute the long-time accurate solution in a flow containing discontinuities. To this end, the second-order-accurate symplectic temporal scheme has been applied to conserve Hamiltanians in the pure advection scheme. To capture discontinuous solution, a slope limiter has been included within the framework of the perturbation finite-difference method. More important, the developed scheme theoretically minimizes the difference between the numerical and exact dispersion relation equations in a grid of fewer number of stencil points. The generalized scheme of varying high-order spatial accuracy has been demonstrated numerically and theoretically to be suitable to predict wave and Navier-Stokes equations.