Properties of High-Order Finite Difference Schemes and Idealized Numerical Testing

Daosheng XU,Dehui CHEN,and Kaixin WU

1Guangzhou Institute of Tropical and Marine Meteorology/Key Laboratory of Regional Numerical Weather Prediction, Guangzhou 510640, China

2National Meteorological Center, Beijing 100081, China

ABSTRACT Construction of high-order difference schemes based on Taylor series expansion has long been a hot topic in computational mathematics,while its application in comprehensive weather models is still very rare.Here,the properties of high-order finite difference schemes are studied based on idealized numerical testing,for the purpose of their application in the Global/Regional Assimilation and Prediction System (GRAPES) model.It is found that the pros and cons due to grid staggering choices diminish with higher-order schemes based on linearized analysis of the one-dimensional gravity wave equation.The improvement of higher-order difference schemes is still obvious for the mesh with smooth varied grid distance.The results of discontinuous square wave testing also exhibits the superiority of high-order schemes.For a model grid with severe non-uniformity and non-orthogonality,the advantage of high-order difference schemes is inapparent,as shown by the results of two-dimensional idealized advection tests under a terrain-following coordinate.In addition,the increase in computational expense caused by high-order schemes can be avoided by the precondition technique used in the GRAPES model.In general,a high-order finite difference scheme is a preferable choice for the tropical regional GRAPES model with a quasi-uniform and quasi-orthogonal grid mesh.

Key words:high-order difference scheme,dispersion,uniform,orthogonal,computational efficiency

1.Introduction

Numerical differentiation is an elementary issue in numerical analysis.Theoretical research in computational mathematics has demonstrated the superior performance of higherorder difference schemes in reducing the truncation error and improving the ability of resolving high-frequency waves (Morinishi et al.,1998;Li,2005;Lin and Zhan,2008).Construction of high-order finite difference schemes falls into two categories:(1) The traditional finite scheme,which can be directly expressed by values on the computational grid.However,the number of grid points used for the computation stencil will be increased for higher-order schemes,and the computational efficiency may be reduced.(2) The compact difference scheme,which can limit the number of grid points used in the computation stencil to less than three,regardless of the order of difference schemes.Its weakness lies in the difference scheme cannot be expressed by grid values explicitly.As a result,a set of linear tridiagonal equations need to be solved (Sun et al.,2014).

High-order difference schemes have been widely applied in ocean models.McCalpin (1994) compared 2ndorder and 4th-order pressure gradient algorithms in a sigma coordinate ocean model and found that the higher-order scheme was helpful to improve the accuracy and stability of the model.Chu and Fan (1997,2000) further implemented a 6th-order difference scheme in a sigma ocean model to reduce the truncation error in calculating the horizontal pressure gradient term.A control-volume ocean model with a 4th-and 5th-order finite difference scheme was presented by Sanderson and Brassington (2002).

The accuracy of discretization for dynamic equations in numerical weather prediction (NWP) is an important issue for its performance,so how to design accurate high-order schemes for NWP models has long been a research focus(Cullen and Davies,1991;Lin and Rood,1996;Yang et al.,2015).Schär et al.(2002) demonstrated the beneficial impact on reducing small-scale noise by higher-order difference schemes with an idealized advection test.Feng and Li(2007) developed a high-order upwind-biased difference scheme,and tested its computational performance with the one-dimensional advection equation and inviscid Burgers’equation.A high-order positive-definite conservative multimoment center constrained finite volume transport model was presented by Shu et al.(2020) to improve the moist transport in NWP models.

In recent years,construction of dynamic cores with high accuracy has become one of the main international trends in developing new-generation NWP models.The truncation error of the GFDL finite volume cubed-sphere dynamical core (FV3,Lin,2004 2008) reached 3rd order.Li et al.(2015) built a 3rd-and 4th-order multi-moment constrained finite-volume method and tested it with a global shallowwater model on the Yin—Yang grid.At present,a central difference scheme with 2nd-order accuracy is still applied for the horizontal discretization of the Global-Regional Assimilation and Prediction System (GRAPES) model,which severely limits its computational performance.It will be meaningful to introduce a higher-order difference scheme for the horizontal discretization in the GRAPES model.

Considering the difficulty of coupling the compact difference scheme with a semi-implicit and semi-Lagrange time discretization scheme in the GRAPES model,we adopted the traditional finite difference scheme to construct a higher-order dynamic core.As the application of high-order difference schemes in operational NWP models is still very rare,their properties need to be further investigated before being implemented in the GRAPES model,which is the main purpose of this paper.The equations of high-order schemes based on Taylor series is displayed in section 2.Dispersion properties of high-order schemes are presented in section 3,and the performance under non-uniform and non-orthogonal grid conditions is presented in sections 4 and 5.Section 6 discusses the impact of high-order finite difference schemes on computational efficiency when implemented in the GRAPES model.Finally,a summary and discussion are provided in section 7.

2.Equations of high-order finite difference schemes

The equations for finite difference schemes with 2nd-,4th-and 6th-order accuracy are displayed in Table 1 (for details of derivation,see Appendix A).The number of grid points used in the computational stencil increases from three to seven with the accuracy raised from 2nd to 6th order,so it is an important issue to deal with the computational efficiency problem for application of high-order finite difference schemes.

The Fourier analysis of error for difference schemes for first- and second-derivative approximation under an unstaggered grid is shown in Fig.1 (for the detailed process of derivation,see Appendix B).The improvement of the higher-order schemes in resolving short waves is significant,especially from 2nd-to 4th-order accuracy.

3.Dispersion properties based on the onedimensional gravity wave equation

Dispersion of ageostrophic atmospheric energy with gravity waves is the physical mechanism of geostrophic adjustment.The dispersion property is an important criterion to evaluate the performance of the model dynamic core.In this paper,a one-dimensional gravity wave equation (Arakawa and Lamb,1977) is adopted to investigate the dispersion properties of high-order difference schemes:

where

u

and

v

are the horizontal components of wind in the

x

and

y

directions,

h

is the depth of fluid,f is a constant Coriolis parameter,g is gravity,and

H

is the mean of

h

.If the solution is assumed proportional to exp[i(kx−νt)](

k

is the wavenumber),then the angular frequency ν for the inertia gravity is given by

Table 1.Difference schemes with 2nd-,4th-and 6th-order accuracy for the first derivative under unstaggered and staggered grids.

Fig.1.Fourier analysis of error for the first derivative approximation with different schemes under an unstaggered grid.

where

L

is the wavelength in the

x

direction.

A grid+2nd order:

A grid+4th order:

A grid+6th order:

C grid+2nd order:

C grid+4th order:

C grid+6th order:

As shown by Fig.2,the frequencies get closer to the true solution when higher-order difference schemes are applied,both in the A grid and C grid for λ/d=2.The dispersion error of the 2nd-order scheme in the A grid is significant for low wavenumber area (wavelength longer than four grid-interval,as indicated by the blue rectangle in Fig.2),while the 4th-and 6th-order schemes in the A grid are as good as the C grid.For the grid scale wave (L/d ＜4),the dispersion error of the C grid is much better than that of the A grid.Considering high-frequency waves cannot be resolved well by the model grid,they need to be filtered out to avoid nonlinear instability.Therefore,poor dispersion properties of the A grid at the grid scale are not problematic as long as such waves can be selectively removed.This conclusion is consistent with the research of Chen et al.(2018).

4.One-dimensional advection equation test

Fig.2.The dependency of (nondimensional) frequency υ/f on the wavelength for the case λ/d=2.

The equations displayed in Table 1 are derived under the hypothesis of a uniform grid,but this hypothesis is not always strictly satisfied in NWP models.For example,the grid interval of the latitude—longitude grid applied in the GRAPES model will get smaller at higher latitude,and the distribution of vertical model levels are also designed to be non-uniform (Xue and Chen,2008).The accuracy of higher-order difference schemes will be degenerated under a non-uniform grid (Xu and Chen,2020).

The performance of higher-order schemes under a nonuniform grid is investigated with the following one-dimensional advection equation based on an unstaggered grid:

Here,u is an arbitrary scalar,and α is defined as a constant.

4.1.Test with a smooth wave solution

When the initial value is set as u(x,0)=Asech(kx),the true solution of Eq.(11) can be written asu(x,t)=Asech[k(x+αt)],where

A

is a parameter defined as a constant,and sech means the hyperbolic secant function (Li,2012).This equation describes the translational motion of a solitary wave being advected to the left,and the choice of the inverse hyperbolic secant function was based on the consideration of testing the high-order scheme with a modal,smooth solution.The parameters are set as follows:

A

=1.0,

k

=0.5,α=2.0,x ∈[−40,20],and an Euler scheme is adopted for time integration with the timestep ∆t=1×10.Three sets of grid configuration are designed as follows to assess the performance of high-order schemes under the non-uniform grid:

Uniform grid:

Nonuniform grid − 1:

Nonuniform − grid − 2:

The grid spacing is gradually reduced in Non-uniformgrid-1 from left to right,while the opposite is the case in Non-uniform-grid-2.As shown in Fig.3,the dissipation and phase errors of the nonuniform and uniform grid are very close to each other,which implies that the dissipation and dispersion properties of higher-order difference schemes will not be degenerated by the weak uniformity of grids.The improvement from the 2nd-order scheme to the 4th-order scheme is most obvious,while the difference between the 4th-order scheme and the 6th-order scheme is very small.This is consistent with the result of Fourier analysis in section 2.The results of the tests for nonuniform-grid-1 and nonuniform-grid-2 are very similar to each other,which also supports the conclusion that weak non-uniformity will not influence the performance of high-order difference schemes.

Fig.3.Numerical test of the one-dimensional advection equation at t=10 s in (a) a uniform grid,(b) non-uniform grid-1,and (c) non-uniform grid-2.

In general,the higher-order difference scheme are still effective for the non-uniform grid if the variation of grid distance is smooth,so they are applicable for the longitude—latitude grid of the GRAPES model and the non-uniform vertical difference process.

4.2.Test with a square wave solution

To study the performance of high-order difference schemes under more discontinuous conditions,a square wave test was further performed.In the square wave test,the initial value of the one-dimensional advection equation[Eq.(11)]was set as

The parameters were set the same as in section 4.1,and the advection test was performed under a uniform grid.Figure 4 shows the simulated square signals at

t

=1 s.The overshoot and undershoot noise were alleviated by the higher-order (4th-and 6th-order) schemes.This implies that the advantages of high-order schemes can be extended to the non-modal solutions,although the theoretical analysis is difficult.

5.Advection test under a terrain-following coordinate

The performance of high-order difference schemes under a vertical non-orthogonal grid with the terrain is further investigated with an idealized advection test (Schär et al.,2002) in this section.It provides a somewhat severe way of testing the effect of non-uniformity on a numerical scheme,as the non-orthogonality and non-uniformity of a terrain-following coordinate is much more severe than it will be on any horizontal grid.

The two-dimensional advection equation under a terrain-following vertical coordinate can be written as (Li et al.,2012)

Fig.4.The result of the square test at t=1 s.

where

and where zdenotes the maximum height of the obstacle.In the following test we use z=3 km,r=25 km,and L=8 km.As shown in Fig.5,the influence of the terrain on vertical levels is significant in the low level,and gradually disappears near the model top,so the computational mesh is nonorthogonal above the terrain.

In the idealized advection test,the topographic obstacle is submerged within a stagnant air mass,but aloft there is a uniform and purely horizontal flow directed from left to right (shown in Fig.6).The anomaly is assigned the shape

Fig.5.Schematic diagram of the non-orthogonal grid mesh for the two-dimensional advection test (black shaded area denotes the topographic obstacle).

where u=10 m s,z=4 km,z=5 km.The time step is 25 s,with grid increments of ∆x=1 km and ∆z=500 m.The total number of vertical levels is 51 and the top height is 25 km.

Fig.6.Vertical cross section for the idealized two-dimension advection test:(a) vertical distribution of velocity u;(b)analytical solution of the density as advected for three instances.

The error fields for the 2nd-,4th-,and 6th-order difference schemes are displayed in Fig.7.As the anomaly is located over the obstacle (at

t

=5000 s),the reduction of spurious small-scale distortion errors caused by higher-order schemes is obvious at the right-upper side of the obstacle(see Figs.7a,c and e),although the orthogonality of the grid is partly influenced by the terrain.The improvement is much more apparent when the computational accuracy is increased from 2nd to 4th order.When the non-orthogonality gets much more severe,the forecast error (see the error above the terrain at

t

=10000 s in Figs.7b,d and f) is almost the same for different schemes.

In general,the performance of high-order difference schemes highly depends upon the orthogonality of the model grid,and the advantage of high-order difference schemes can still be retained when the grid is nearly orthogonal.

6.Influence on the computational efficiency based on the forecast equation of the GRAPES model

The numbers of grids used in the calculation stencil will be increased if higher-order finite difference schemes are applied,so it is an important issue to consider the influence of high-order difference schemes on the computational efficiency for large-scale computation in operational NWP models.In this section,the stencils for the 2nd-and 4th-order difference scheme are compared based on the forecast equations of the GRAPES model,and the influence on the actual computational efficiency is discussed.

The dynamic equation of the GRAPES model under a spherical coordinate system can be written as (Chen and Shen,2006)

Fig.7.Numerical error of difference schemes for the two-dimensional advection equation:(a,b) 2nd-order scheme;(c,d) 4th-order scheme;(e,f) 6th-order scheme;(a,c,e) t=5000 s;(b,d,f) t=10000 s.

where B— Bare the coefficients of 19 grids of the stencil.If the 4th-order scheme is used for horizontal discretization,the Helmholtz equation contains 47 grids (the detail derivation process is omitted):

The generalized conjugate residual (GCR) method is used to solve the Helmholtz equation in the GRAPES model.To accelerate the speed of convergence,a preconditioned algorithm is applied to simplify the coefficient matrix.According to the fact that the coefficients at points(i,j,k),(i,j,k −1) and (i,j,k+1) are an order of magnitude larger than other points,the preconditioned matrix only retains these three points’ coefficients and then it becomes a segmented tridiagonal matrix (e.g.,Xue and Chen,2008,pp.132—136).The magnitude of coefficients at 47 points under the 4th-order scheme are displayed in Table 2,and we also find that the coefficients at points (i,j,k),(i,j,k −1) and(i,j,k+1) are an order of magnitude larger than other points.As a result,the preconditioned matrix is the same as the 2nd-order scheme,which implies that the time used to solving the Helmholtz equation will not be increased,and so a higher-order difference scheme is a feasible choice for the GRAPES model.

Table 2.Magnitude of coefficients in Eq.(27).

7.Conclusion

To go beyond the well-worn traditional linear shallowwater uniform-grid analyses into somewhat more realistic testing frameworks that consider situations of use in comprehensive weather and climate models,the properties of higherorder difference schemes were investigated for application in the GRAPES model in this paper.The one-dimensional gravity wave equation was used to study the dispersion properties of high-order schemes.Obvious reduction in the dispersion error was found in the higher-order schemes with an Arakawa-A grid for those well-resolved waves,especially from 2nd-to 4th-order accuracy,while the improvement for the Arakawa-C grid was not so apparent.Based on the discovery that the dispersion errors are almost equivalent for the A-and C-grid under 4th-and 6th-order schemes,the possibility of an unstaggered GRAPES dynamical core,which is more convenient to couple with physical processes (Purser and Leslie,1988),will be further investigated in the future.The preferable properties of high-order difference schemes are still obvious under a non-uniform grid when the variation of grid distance is smooth.The square wave test also indicated the superiority of high-order schemes.For a grid mesh with severe non-orthogonality and non-uniformity,the advantages of high-order schemes are not so obvious,as demonstrated by the results of idealized two-dimensional advection tests.The number of grid points used in the computational stencil will be increased for higher-order schemes,but it will not reduce the computational efficiency of the GRAPES model when the precondition technology being applied to reduce the magnitude of condition number and improving the convergence rate of the GCR.

In general,a high-order finite difference scheme is a feasible choice for a high-resolution,longitude—latitude grid model with a nearly uniform and orthogonal regional domain—for example,the Tropical Regional Area Model System (TRAMS;Xu et al.,2015) developed by the Guangzhou Institute of Tropical and Marine Meteorology,based on the regional version of GRAPES.

Although the result of the two-dimensional advection test under a terrain-following coordinate is in a similar spirit to the analysis of Reinecke and Durran (2009),it can still be misleading since reflections in the one-way advection equation are nearly always into grid-scale wavelengths that are easily diffused for a relatively long-term simulation (Harris and Durran,2010).Extending a similar test to a more dynamically active case supporting upstream physical wave propagation (Vichnevetsky,1987) will be more meaningful.A complete study of the performance by implementing high-order schemes in TRAMS will be the next stage in our work.A general finite high-order difference scheme applicable in a nonuniform grid will also be investigated in the future.

Acknowledgements

.This work was supported by the National Natural Science Foundation of China (Grant No.U1811464) .We thank the two anonymous reviewers and the editors for their comments,which greatly improved this paper.

APPENDIX A Derivation of High-Order Difference Schemes based on the Taylor Series Expansion Method

The general finite difference equations with arbitrary order accuracy are deduced first,followed by the specific forms of the 2nd-,4th-and 6th-order schemes.

1.Equations of traditional finite difference schemes with arbitrary order accuracy

The method to deduce high-order difference schemes at one point can be summarized as follows:the number of grid points needed for the required accuracy is determined first;then the expression of derivatives including indefinite coefficients are deduced;and finally these indefinite coefficients will be obtained by solving the set of simultaneous equations based on Taylor series expansion.

The following 2

n

+1 grid will be used to contract 2

n

order accuracy schemes for the first derivative:

The expansion of derivations can be written as (with coefficients undetermined):

The function f is then expanded at point i with a 2norder Taylor series:

Let n=1 in Eqs.(A1)—(A2),and the expressions of the second-order scheme can be obtained.The first derivative equation in Eq.(A1) can be rewritten as follows:

Substituting Eq.(A3) and Eq.(A6) into Eq.(A9),and omitting the remainder terms,we can obtain:

The following coefficients can be derived by equaling the corresponding terms in Eq.(A10):

The 2nd-order finite difference scheme for the first derivative can now be expressed as:

3.Fourth-order difference scheme

Let

n

=2 in Eqs.(A1)—(A2),and the expression of the 4th-order scheme can be obtained.The expansion of Eq.(A1) is

2.Second-order difference scheme

Substituting Eqs.(A3)—(A4) and Eqs.(A6)—(A7) into Eq.(A13),

where

The 4th-order finite difference scheme for the first derivative can now be expressed as:

4.Sixth-order difference scheme

Let

n

=3 in Eqs.(A1)—(A2),and the expression of the 6th-order scheme can be obtained.The expansion of Eq.(A1) is

Substituting Eqs.(A3)—(A5) and Eqs.(A6)—(A8) into Eq.(A17),and omitting the remainder terms,we can obtain:

where

The 6th-order finite difference scheme for the first derivative can now be expressed as:

For the staggered grid,the following 2nd-,4th-and 6th-order finite difference schemes can be obtained with a similar method:

APPENDIX B Fourier Analysis of Error for Finite Difference Schemes

For an arbitrary periodic function f(x) with its wavenumber defined as α,the Fourier transform (represented by

F

{})of f(x) and ∂f(x)/∂x is

According to the properties of Fourier transition,we can get

where ∆x is the grid distance.The differencing of the 2ndorder approximation for the first derivative under an unstaggered grid can be expressed as:

and its Fourier analysis error can be obtained by performing the Fourier transform on both sides of Eq.(B5):

Hence,

With the wavelength defined as L=2π/α,Eq.(B7) can be rewritten as

The following Fourier analysis error of the 4th-and 6th-order approximations can be derived with a similar method:

4th-order scheme: