Higher-order optimized hybrid Robert-Asselin type time filters

Journal of Computational Physics 399 (2019) 108941

Contents lists available at ScienceDirect

Journal of Computational Physics www.elsevier.com/locate/jcp

Higher-order optimized hybrid Robert-Asselin type time filters Praveen K. Maurya, Manoj K. Rajpoot ∗ , Vivek S. Yadav Department of Mathematics & Statistics, Rajiv Gandhi Institute of Petroleum Technology, Amethi 229304, UP, India

a r t i c l e

i n f o

Article history: Received 5 April 2019 Received in revised form 25 July 2019 Accepted 5 September 2019 Available online 11 September 2019 Keywords: Hybrid time filters Spurious dispersion Grid staggering Rotary shallow water equation Navier-Stokes equation Direct numerical simulation

a b s t r a c t In this paper, a class of optimized filters based on hybrid blending of implicit and explicit filters are formulated for the three time-level leapfrog time integration scheme. As the leapfrog method is the most recurrently employed in ocean and atmospheric simulations, however, it also admits a spurious mode in the numerical computations. To subdue the computational mode, the developed hybrid time filters are also optimized using constraints arising from error dynamics. The optimized filters are adept in suppressing the computational mode(s) more efficaciously without altering the physical mode. Finally, the developed hybrid filters coupled with centered and compact spatial discretization schemes are also authenticated for the numerical solutions to one-dimensional (1D) convection equation, two-dimensional (2D) dispersive rotating shallow water equation, and unsteady 2D incompressible Navier-Stokes equation at different Reynolds numbers. Present solutions are also compared with results reported in the literature. © 2019 Elsevier Inc. All rights reserved.

1. Introduction Sizable research continues to explore the efficient and accurate time integration methods for the direct numerical simulation (DNS) of flow and wave propagation problems. As, obtaining the high accuracy solutions for propagation of disturbances is a primary objective in fluid dynamics, and it is reported in [1,2] that to ascertain the enhanced accuracy, adopted methods should be analyzed for hyperbolic equations using 1D convection equation as the model system. Being non-dispersive, this equation provides proving ground to gauge the efficiency of a numerical method for its appropriateness for the DNS of the propagation of disturbances. Authors in [3–5] reported that in atmospheric and ocean models, results are highly vulnerable to numerous time advancement methods. Therefore, errors engendered by time advancement methods perform crucial role in the geophysical simulations [6]. The Robert-Asselin (RA) filtered leapfrog method [7–10] is largely used in geophysical simulations, and the dispersion analysis of RA and Robert-Asselin-Williams (RAW) filters by considering space-time discretization together is discussed in [11]. The key issue with the use of leapfrog time integration method is that besides physical mode it also introduces computational mode. The central postulation in using RA filtered leapfrog time integration method is that it attenuate the computational mode in discrete computing. However, such filters also attenuate the physical mode deteriorating the accuracy of the numerical solutions [9,10]. Williams [9,12] devised a variation to the RA filter that reduces the precluded attenuation of the physical mode. Author in [13] also discussed the functioning of RAW filter for implicit-explicit (IMEX) temporal advancement method by using explicit leapfrog method for the low time-frequency region and implicit CrankNicolson (CN) method for the high time-frequency region. Using spectral analysis, author in [14] presented the problematic

*

Corresponding author. E-mail addresses: [email protected], [email protected] (M.K. Rajpoot).

https://doi.org/10.1016/j.jcp.2019.108941 0021-9991/© 2019 Elsevier Inc. All rights reserved.

2

P.K. Maurya et al. / Journal of Computational Physics 399 (2019) 108941

aspects of IMEX method for the DNS of transitional and turbulent flows using explicit third-order Runge-Kutta method and implicit CN method coupled with spatial discretization schemes. Li & Trenchea [10] developed higher order RA type filter and noted that the performance of developed scheme is similar to the third-order Adams-Bashforth method having two computational modes. Marsaleix et al. [15] formulated alternatives to the RA filtered leapfrog time integration methods using Laplacian approach to introduce the artificial viscosity. Moustaoui et al. [16] presented the leapfrog method based on fourth-order implicit time filter of decentered type, which is centered around nth time-level, while RA/RAW filters are centered about (n − 1)th level. Easy implementation and broader stability properties of the method developed in [16] form a fair substitute for geophysical simulations, and are used for numerical solution to the Lorenz’s and shallow water systems. Moreover, method in [16] is also used for the time advancement of Maxwell equations for unstaggered temporal grids [17]. The main objective of this paper is to develop the hybrid filtered leapfrog methods to suppress the computational mode more efficiently without disturbing the neutral stability of the physical mode. The hybrid formulation is based on the weighted average [18] of explicit and implicit time filters. In the present study, a hybrid family of filtered leapfrog methods coupled with second-order centered (CD2 ) and compact spatial discretization schemes are analyzed by considering space and time discretization together. Fourier-spectral analysis is used to assess the numerical properties of the developed methods. Moreover, for improved accuracy, values of the filter parameters are fixed by performing the constrained optimization [19,20] based on error dynamics. The optimized filters provide the near-neutral stability for the physical mode, while further attenuating the dominating computational mode. Finally, to validate the robustness of the developed methods, propagation problems following model dispersive and non-dispersive dynamical systems are also considered on collocated and staggered meshes. Also, flow inside 2D lid-driven cavity (LDC) following Navier-Stokes equation at different Reynolds number is considered to demonstrate the efficiency of the present methods for the nonlinear case. The paper is organized in the following manner. In section 2, a brief introduction about the representative model systems is presented. In section 3, analysis and development of the hybrid time filtered leapfrog methods are presented. In section 4, Fourier-spectral analysis is used to discuss the spectral properties of the developed schemes with 1D convection equation. Furthermore, in section 5, constrained optimization based on error dynamics [21] is considered in the spectral plane to fix the optimal values of filter parameters. Validation of the developed methods is accomplished by solving the propagation problems governed by dispersive and non-dispersive systems in sections 6 & 7. Moreover, benchmark flow problem governed by 2D Navier-Stokes equation is reported in section 8, and finally, summary and conclusion complete the manuscript. 2. DNS of convection dominated flows The model 1D system for the linear hyperbolic transport problem is given as

∂u = f (u ) ∂t

(1)

for 1D linear advection equation f (u ) = −c ∂∂ux with c > 0, and the efficiency of the developed methods are demonstrated with the help of Eq. (1). Being the non-dissipative and non-dispersive, engaged numerical method for Eq. (1) needs to be neutrally stable and devoid of phase and dispersion errors. Furthermore, the governing equations for the unsteady incompressible flow problems are obtained by considering the conservation of mass and momentum. In the present case, we have considered the 2D incompressible Navier-Stokes equation and rotating shallow water equation as the representative systems for flow and wave propagation problems for the chosen domain topography. The Navier-Stokes equation in stream function-vorticity (ψ − ω) formulation is given as

ψ = −ω,

Dω Dt

=

1 Re

ω

(2)

D where,  is the Laplacian operator, Dt = ∂∂t + V .∇ represent the material derivative, ω is the component of the vorticity vector normal to the plane, and ψ is the stream function with (u , v ) as velocity components. The (ψ, ω) formulation is favored over primitive variables due to its computational capability in satisfying mass conservation exactly. Furthermore, single-layer approximation of rotating shallow water equation [22] given as

DV

+ f kˆ × V + g ∇ η˜ = 0,

D η˜

+ η˜ ∇. V = 0 (3) Dt Dt where, V = (u˜ , v˜ ) denotes the velocity components; η˜ represent the surface height; g is the acceleration due to gravity, and f is the Coriolis parameter. The developed optimized hybrid filters coupled with CD2 and near-spectral accurate compact scheme are used to solve the flow and wave propagation problems governed by systems (1)-(3) on collocated and staggered meshes. In the next section, detailed analysis and development of hybrid filters is discussed. 3. Formulations of hybrid time filters The hybrid filters are obtained from the weighted average [18] of explicit and implicit filters. The higher order explicit filters are based on stencils centered around (n − 1)th time-level, however the implicit filters are based on the stencils

P.K. Maurya et al. / Journal of Computational Physics 399 (2019) 108941

3

centered around nth time-level. Therefore, the explicit filter are of non-centered nature and implicit are of the centered nature. In general, formulation of time filters are based on system given by Eq. (1) and can be represented as

 t ∗ ∂ un  ∂ i u  n = f ( u ) + κ ∂ t L F ∂ti  where,

κ denotes the diffusion factor for the filter. For i = 2, Eq. (4) produces the classical RA filter [7,8] by choosing = γ (unF−2 − 2un−1 + un ), and is given by

n−1

κ ∂∂ t 2u  2

(4)

un+1 = unF−1 + 2t f (un ) unF−1 = un−1 + γ (unF−2 − 2un−1 + un )

(5)

unF

n

where, u and represent the unfiltered and filtered solution at the nth time-level, respectively, and γ denotes the filter parameter. Therefore, in the case of RA filter, the numerical diffusion of the computational mode is obtained via the second-order derivative centered at (n − 1)th time-level. Similar to RA filter, higher order diffusion is used by choosing higher order derivatives centered at (n − 1)th level. It is noted that such decentered filtered are of explicit type, because the filtered values at (n − 1)th level does not appear on the right hand side expression in Eq. (5), however, filtered values at lower time-level may appear. Moreover, as discussed in [15,16], the implicitness in the filters occur due to the presence of filtered values on the both sides at the same time-level. Similarly, implicit time filtering methods are derived by considering the approximation which is centered at nth timelevel having stencil points up to (n + 1)th level. In opting, second-derivative is given as

n−1

κ ∂∂ t 2u  2

= γ (−unF−1 + 2un − un+1 ), the implicit filter based on

un+1 = unF−1 + 2t f (un ) unF−1 = un−1 + γ (−unF−1 + 2un − un+1 )

(6)

The first hybrid filter (IFL) is obtained by taking the weighted average of explicit and implicit filters in Eqs. (5)-(6), given as, ω1 (6) + (1 − ω1 )(5), with ω1 representing the weight factor. Resultant IFL is obtained as

un+1 = unF−1 + 2t f (un )



unF−1 = un−1 + γ (1 − ω1 )unF−2 − ω1 unF−1 + (2ω1 − 2)un−1 + (1 + ω1 )un − ω1 un+1



Unlike to explicit filter, here un+1 and unF−1 both are unknown and depend on each other, and can be represented as



−1 1 + γ ω1

1

γ ω1



u n +1 unF−1





=

c1 c2

 



where, c 1 = 2t f (un ) and c 2 = un−1 + γ (1 − ω1 )unF−2 + (2ω1 − 2)un−1 + (1 + ω1 )un , and the coefficient matrix is easily invertible, as



u n +1 unF−1

 =

1



(1 + 2γ ω1 )

1 + γ ω1

−γ ω1

1 1



c1 c2

 (7)

which is almost equivalent to solving a explicit system of two equations. Similarly, the second hybrid filter (IIFL) based on third-order derivative can also be developed by choosing i = 3 in Eq. (4), as

un+1 = unF−1 + 2t f (un )



unF−1 = un−1 + γ (ω2 − 1)unF−3 + (3 − 2ω2 )unF−2 − 3ω2 unF−1 + (3ω2 − 3)un−1 + (1 + 2ω2 )un − ω2 un+1

 

which is obtained from the weighted average of explicit and implicit filters, respectively, given as, un+1 = un−1 + γ −unF−3 +

 n



3unF−2 − 3un−1 + u + 2t f (un ) and un+1 = un−1 + γ unF−2 − 3unF−1 + 3un − u filter (IIIFL) based on fourth-order derivative is given as

un+1 = unF−1 + 2t f (un )



 n+1

+ 2t f (un ). Similarly the third hybrid



unF−1 = un−1 + γ (1 − ω3 )unF−4 + (3ω3 − 4)unF−2 − 6ω3 unF−1 + (4ω3 − 4)un−1 + (1 + 3ω3 )un − ω3 un+1   which is based upon the explicit and implicit filters given as, un+1 = un−1 + γ unF−4 − 4unF−3 + 6unF−2 − 4un−1 + un +   2t f (un ) and un+1 = un−1 + γ −unF−3 + 4unF−2 − 6unF−1 + 4un − un+1 + 2t f (un ). Finally, the fourth hybrid filter (IVFL) based on fifth-order derivative is given as

4

P.K. Maurya et al. / Journal of Computational Physics 399 (2019) 108941

Table 1 Stability behavior of various order explicit and implicit filters. Filter

Explicit

Implicit

IFL IIFL IIIFL IVFL

stable unstable stable unstable

unstable stable unstable stable

un+1 = unF−1 + 2t f (un )



unF−1 = un−1 + γ (ω4 − 1)unF−5 + (5 − 4ω4 )unF−4 + (5ω4 − 10)unF−3 + 10unF−2 − 10ω4 unF−1

+ (5ω4 − 5)un−1 + (1 + 4ω4 )un − ω4 un+1



which is obtained by considering the weighted average of the following explicit and implicit filters, given as un+1 = un−1 +    γ −unF−5 +5unF−4 − 10unF−3 + 10unF−2 − 5un−1 + un + 2t f (un ) and un+1 = un−1 + γ unF−4 − 5unF−3 + 10unF−2 − 10unF−1 + 5un − un+1 + 2t f (un ), respectively. Here, ωi denotes the weight factor for the ith filter. It can be noted from Eq. (4) that the higher order time filtering to suppress the computational mode for leapfrog method is achieved by adding the numerical diffusion based on higher derivatives. In general, the added numerical diffusion/antidiffusion based on odd (upwinded) and even (centered) order derivatives is given as [23,24],

⎧  n−(i −1)/2 ⎪ i  i +1 n−(i −1)/2 ⎪ ⎪ + κ E ∂∂ t i+1u CD i = 1, 3, 5 . . . ⎨ ∂∂ tui  2

n ∂ i u  CD2 = ∂ t i UD1 ⎪ ⎪ i n−(i /2) ⎪ ⎩∂ u

(8)

i = 2, 4, 6 . . .

∂ t i CD2

where, UD1 denotes first-order upwinding and κ E represent the diffusion coefficient. Such diffusion terms are used to design the explicit filters having stencils up to nth time-level. In addition to suppressing the computational mode, the use of explicit time filtered leapfrog scheme also attenuates the physical mode. To circumvent this problem with explicit time filtering, implicit filters based on the odd and even order derivatives are obtained as

⎧ i  i +1 n−(i −3)/2 n−(i −3)/2 n+1 ⎪ − κ I ∂∂ t i+1u CD i = 1, 3, 5 . . . − ∂∂ tui CD ⎨  2 2 ∂u − i  = ∂ t UD1 ⎪ ⎩− ∂ i u n−(i −2)/2 i = 2, 4, 6 . . . ∂ t i CD i

(9)

2

with κ I representing the diffusion coefficient for the implicit filtering. The implicit filtering based on Eq. (9) with stencils up to (n + 1)th level amplifies the physical mode, and further attenuate the computational mode. Therefore, the hybrid filters based on weighted average of explicit and implicit time filtering strategy are capable of suppressing computational mode more effectively without any attenuation for the physical mode. Furthermore, for the cases of upwinding, the sign of the numerical diffusion or anti-diffusion is decided by feedback stabilization. As stated in [23], at the jth node Eq. (1) can be expressed as ∂∂ut = k a j u j ±k + b j , and in computing added numerical diffusion is indispensable to avoid the numerical instability. The added numerical diffusion can be defined in terms of error e j , as, e j = C e −αt , representing a negative feedback whenever α > 0. Therefore, for the first filter based on the second-order diffusion term, positive value of the filter parameter (γ ) introduces the numerical diffusion, while negative value corresponds to the anti-diffusion. Moreover, for the fourth-order diffusion term, negative value introduces the numerical diffusion while positive value is for anti-diffusion and so on for higher order diffusion terms. The stability behavior of the developed various order explicit and implicit filters based on the added numerical diffusion or anti-diffusion is given in Table 1. It can be noted from Table 1, that for each combination the stability behavior of the explicit and implicit filters are of opposite nature for the physical mode, and therefore the hybrid version of the filters provide near-neutral stability for the physical mode. In the following section, bilateral Fourier-Laplace transform [25] is used to appraise the spectral properties of the physical and computational modes in the spectral plane. 4. Fourier-spectral stability analysis In this section, Fourier-spectral analysis [26,23] is used to investigate the accuracy of the developed schemes in the spectral plane. Using Fourier transform, unknown, u j , and its derivative are represented as

∞ u (x j , t ) = −∞

Uˆ (k, t )e ikx j dk,

u  (x j , t ) =

∞

−∞

ik Uˆ (k, t )e ikx j dk

P.K. Maurya et al. / Journal of Computational Physics 399 (2019) 108941

5

where, Uˆ (k, t ) stands for the Fourier transform of the unknown and k denotes the wavenumber. As reported in [23,27], numerical approximation of the derivative is also represented as

km



u N (x j , t ) =

ikeq Uˆ (k, t )e ikx j dk

−km

where, keq denotes the equivalent wavenumber for first-order derivative, and ±km are construed through Nyquist limit. In discrete computing, the first-order derivative is written as, u j =

1 h

N

l=1

C jl ul , where N denotes the number of mesh-points,

h represent the mesh-spacing, and C is the corresponding discretization matrix. Utilizing nodal amplification factors (G j = Uˆ j (k, t + t )/Uˆ j (k, t )), the general numerical solution at nth time-level is written as



u (x j , nt ) =

Uˆ 0 (k) |G j |n e i (kx j −nβ j ) dk

where, Uˆ 0 (k), denotes the Fourier transform of the initial condition with β j = −Arg(G j ). Expressions for the nondimensional phase speed and group velocity at the jth mesh-point for Eq. (1) are given as,

 cN  c

β

=  j t , j



V gN c



j

=

1 dβ j N c d(kh)

[21], where  denotes the circular frequency and c N , V g N represent the numerical phase speed and group velocity, respectively. Using the Fourier transform, nodal amplification factors for the IFL are expressed by the following equation













1 + 2γ ω1 G 2 + 2 −γ + 2iN c φ(kh)(2γ ω1 + 1) G + 2γ (1 − ω1 ) + 2i γ N c φ(kh)(ω1 − 1) − 1 = 0

(10)

c t h

where, N c = is the Courant-Friedrichs-Lewy (CFL) number and φ(kh) is non-dimensional equivalent wavenumber for the first-order spatial derivative defined in section 6 for the used spatial discretization schemes. Roots (G 1,2 ) of Eq. (10) are given in Appendix A, where, G 1 and G 2 represent the physical and computational mode, respectively. Theoretically, Eq. (1) admits only one mode, which governs the propagation of the solution. For second filter (IIFL) nodal amplification factors are given by the following cubic equation











1 + 4γ ω2 G 3 + 2iN c φ(kh)(3γ ω2 + 1) − 4γ G 2 + 4γ (1 − ω2 ) + 2i γ N c φ(kh)

   (2ω2 − 3) − 1 G − 2i γ N c φ(kh)(ω2 − 1) = 0

(11)

and the roots are given in Appendix A. Out of three roots, one corresponds to the physical mode and other two portray computational modes. For IIIFL, expression for the nodal amplification factors using MATLAB symbolic toolbox, is given as









1 + 7γ ω3 G 4 + −7γ − γ ω3 + 2iN c φ(kh)(6γ ω3 + 1) G 3 +



2

N c φ(kh)(ω3 − 3) − 1 G +





γ (8 − 7ω3 ) + 4i γ 





γ (ω3 − 1) + 2i γ N c φ(kh)(4 − 3ω3 ) G + 2i γ N c φ(kh)(ω3 − 1) = 0 (12)

having four roots of which one represent the physical mode and other three depict the computational modes. Roots of Eq. (12) are given in Appendix A. Finally, nodal amplification factors for the IVFL are given by the following quintic equation, as











1 + 11γ ω4 G 5 + −11γ − 4γ ω4 + 2iN c φ(kh)(10γ ω4 + 1) G 4 + 5γ (3 − 2ω4 )

    − 20i γ N c φ(kh) − 1 G 3 + γ (4ω4 − 5) + 10i γ N c φ(kh)(2 − ω4 ) G 2 + γ (1 − ω4 )    + 2γ N c φ(kh)(4ω4 − 5) G − 2i γ N c φ(kh)(ω4 − 1) = 0

(13)

owing to intricacy of the Eq. (13), Taylor’s series expansion is used to analyze the numerical behavior of physical mode (G phys ) given as

G phys (γ , ω4 , N c φ(kh)) =

∞ 



n

cn (γ , ω4 ) −iN c φ(kh)

n =0

where, coefficients (cn ) are real and behavior of the physical mode would be dictated by the corresponding leading order terms in the expansion. By substituting the series expansion into the Eq. (13), coefficient corresponding to 0th order term is given as

(1 + 11γ ω4 )c 05 − γ (11 + 4ω4 )c 04 + (15γ − 10γ ω4 − 1)c 03 + γ (4ω4 − 5)c 02 + γ (1 − ω4 )c 0 = 0 admitting five roots of which one corresponds to physical mode (c 0 = 1), and other four values of c 0 belong to computational modes. Comparison of higher order coefficients, considering first few terms, are given as

6

P.K. Maurya et al. / Journal of Computational Physics 399 (2019) 108941

Table 2 Stability and temporal order of accuracy of the physical mode for the hybrid filters for ωi = 0.5, i = 1, 2, 3, 4.

c 0 = 1, c6 =

c 1 = 1,

Filter

Stability

IFL IIFL IIIFL IVFL

unconditionally unconditionally unconditionally unconditionally

c2 =

1 2

,

c 3 = 0,

unstable stable stable unstable

1 c4 = − , 8

c5 =

Amp error/ time step

Phase error/ time step

4 4 6 6

3 3 3 3

γ (1 − 2ω4 ) 2 − 8γ (1 − 2ω4 )

1 + 4γ (−5 + 8ω4 − 56γ ω4 + 14γ + 56γ ω42 ) 16(1 − 8γ + 16γ 2 + 64γ 2 ω42 + 16γ ω4 − 64γ ω42 )

Therefore, the absolute value of the physical mode is given as

|G phys | = 1 +

1 16

 6 (16c5 − 16c6 + 1) N c φ(kh) 

6

For the weight factor, ω4 = 0.5, absolute value is given by |G phys | = 1 + γ N c φ(kh) . Consequently, the amplitude and phase mismatch errors are given as

|G phys | − |G exact | =

1 16

 6 (16c5 − 16c6 + 1) N c φ(kh)

Arg(G phys ) − Arg(G exact ) = −

( N c φ(kh))3 6

where, G exact = exp (iN c kh), represent the exact amplification factor for Eq. (1). Here, the amplitude error/time-step for the physical mode is proportional to O (t )6 , and phase error/time-step varies as O (t )3 . Moreover, stability and temporal order of accuracy for the various developed filters are shown in Table 2. As, for ω4 = 0.5 it is perceived that the physical mode is unconditionally unstable for finite values of N c and kh, which is also the case for the IFL. Furthermore, it is also observed from Table 2 that the temporal order of accuracy of the developed methods is higher than the second- and third-order Adams-Bashforth methods [10]. Here, we would like to mention that though use of hybrid filters attenuate the computational mode without affecting the neutral stability of physical mode, even if the presence of highly attenuated computational mode(s) affects the atmospheric simulations, as reported in [27,23]. To accomplish the near neutrally-stable physical mode, constrained optimization problem based on error dynamics is formulated in next section to fix the values of filter parameters (γ , ωi ). 5. Mathematical formulation of the optimization problem To implement the optimization problem as in [19,20] for the improved phase, dispersion and stability properties, following objective function is considered

δ |G phys − G exact |2 d(kh)

F (γ , ωi , N c ) =

(14)

0

Here, G phys is the amplification factor corresponding to physical mode and G exact denotes the corresponding exact amplification factor, as defined in the previous section. The objective function given by Eq. (14) is also subjected to the following explicit constraints, given as

δ1   F 1 (γ , ωi , N c ) = |G phys | − 1d(kh) ≤ 1

(15)

0

δ1 F i (γ , ωi , N c ) =

|G i comp |d(kh) ≤ i , i = 2, 3, 4...

(16)

0

where, G i comp represent the ith computational mode(s) for the corresponding hybrid filters. Here, δi and i are constants based on the prescribed error tolerance bar in the spectral plane. Sufficiently small tolerance values (i ), as given in [19,20],

P.K. Maurya et al. / Journal of Computational Physics 399 (2019) 108941

7

Table 3 Standard values for filter parameters for explicit, implicit and hybrid filters of various orders. Filter

Filter parameter (explicit)

Filter parameter (implicit)

Hybrid (γ , ωi = 0.5)

IFL IIFL IIIFL IVFL

0.1 0.05 1/30 1/40

0.125 0.0625 1/33 1/39

1/9 1/18 2/63 2/79

Table 4 Optimized values for the hybrid filter parameters. Filter

Filter parameter (γ )

Weight factor (ω)

IFL IIFL IIIFL

0.19203521 0.20254816 0.20351572

0.40133920 0.52013561 0.51025443

are chosen to ensure the neutral stability for the physical mode and highest possible attenuation for the computational modes. The optimization problem is solved using grid-search technique [23] in which the feasible region in (γ , ωi )-plane is fixed by using the constraints for the arbitrary but fixed values of N c . For the present case, δ = π /2, δ1 = π for the objective function and constraints, respectively. The chosen error tolerances (i ) are of the order 10−4 . The standard and optimized values of parameters (γ , ωi ) are given in the following Tables 3 & 4. As reported in [9,11], γ (= ν /2 with ν = 0.2) provides the best result for the numerical solutions to the dispersive and non-dispersive linear systems for the RA/RAW filters. In the present case, standard values of the parameters for developed filters are obtained via comparing with the value of RA filter parameter (γ = 0.1). Numerical properties of the developed optimized methods coupled with spatial discretization schemes are discussed next. 6. Numerical properties and solution to the model 1D convection equation Numerical properties of the hybrid filters coupled with CD2 and compact spatial discretization schemes in the

( N c , kh)-plane are discussed next. For the collocated mesh, optimized upwind compact scheme (OUCS3) given in [26] for first-order spatial derivative is obtained by fixing the coefficients of an original sixth order compact scheme [28,29] by minimizing L 2 norm of its deviation from the derivative obtained by spectral method up to the Nyquist limit using different initial conditions for 1D convection equation. The OUCS3 scheme for the first-order spatial derivative for interior nodes is given as

pu i −1 + u i + pu i +1 =

q2 h

( u i +2 − u i −2 ) +

q1 h

( u i +1 − u i −1 )

(17)

where, p = 0.3793894912, q2 = 0.045801298 and q1 = 0.787786895. The explicit near-boundary stencils for the OUCS3 scheme to solve non-periodic problems are given by Eq. (23). Non-dimensional equivalent wavenumbers for interior nodes using OUCS3 is given as

 φ(kh) = 2

q2 sin(2kh) + q1 sin(kh)



1 + 2p cos(kh)

and non-dimensional equivalent wavenumber for the CD2 scheme is given as, φ(kh) = sin(kh). Absolute value of the nodal amplification factors for the IFL coupled with CD2 and OUCS3 schemes are shown in Fig. 1(a). Top four frames of Fig. 1(a) represent the numerical properties with CD2 scheme, and bottom four frames are with OUCS3. From Fig. 1(a) it is noted that for the optimized values of filter parameters computational mode is further damped, as compared to the standard values. Also, the physical mode for the optimized filter parameters is near neutrally stable, while there is mild instability for the standard values of filter parameters. Similar phenomena are also observed for OUCS3 spatial discretization. Moreover, phase error and group velocity contours for the IFL with CD2 and OUCS3 are displayed in Fig. 1(b). As OUCS3 provides the near spectral resolution, therefore phase and dispersion properties with the OUCS3 are superior to CD2 , as shown in Fig. 1(b). For the IIFL method, magnitudes of the amplification factors with OUCS3 are shown in Fig. 2. In this case, G 1 represent the physical mode while G 2,3 are the computational modes. Also, |G 3 | is almost negligible and the dominating computational mode (G 2 ) is more effectively attenuated for the optimized values of the parameters. Finally, for IIIFL-OUCS3, there are three computational modes G 2,3,4 of which G 2,3 are sufficiently small and the dominating computational mode (G 4 ) is much more damped for the optimized values of the parameters, as shown in Fig. 3. Numerical error in wave propagation introduced by discrete schemes can be gauged by considering the propagation of a wave-packet following Eq. (1), given as

8

P.K. Maurya et al. / Journal of Computational Physics 399 (2019) 108941

(a) Fig. 1. (a) Contour plots showing numerical properties in the spectral plane for physical and computational modes using indicated schemes for standard and optimized values of parameters. (b) Contour plots showing non-dimensional phase error and group velocity for physical and computational modes using indicated schemes for standard values of parameters.

u (x, t = 0) = e −32(x−5) cos(ko x) 2

To solve this propagation problem for interior nodes, a computational domain 0 ≤ x ≤ 30, with 3000 equispaced mesh-points is considered with N c = 0.5 and ko h = 0.5. Numerical solutions using IFL, IIFL and IIIFL with CD2 scheme are shown in the top six frames of Fig. 4 for the optimized values of the parameters. Being three-time level method, for IFL time advancement from t = 0 to 2t is performed using fourth-order Runge-Kutta (RK4 ) method, and for IIFL, IIIFL and IVFL, RK4 is used for first three, four and five time steps, respectively. It is observed from Fig. 4, that IFL method introduces large diffusion at higher time instants, while results using IIFL and IIIFL are indistinguishable. Numerical solutions displayed in the top frames of Fig. 4 are in conformity with the numerical properties shown in Figs. 1 & 2, and also compared with the exact solutions.

P.K. Maurya et al. / Journal of Computational Physics 399 (2019) 108941

9

(b) Fig. 1. (continued)

Numerical solutions at indicated time instants using IFL, IIFL and IIIFL along with OUCS3 are shown in bottom frames of Fig. 4. From the left-column frames, it can be noticed that IFL - OUCS3 introduces higher diffusion and in result packet gets attenuated at subsequent times, which is also evident from the property contours in Fig. 1. However, for IIFL & IIIFL, numerical solutions are identical with the exact results for the indicated time instances. Numerical solutions to the dispersive linearized rotating shallow water equation (LRSWE) are discussed next.

10

P.K. Maurya et al. / Journal of Computational Physics 399 (2019) 108941

Fig. 2. Contour plots showing numerical properties in the spectral plane for physical and computational modes for IIFL − OUCS3 scheme for standard and optimized values of parameters.

7. Solution to 2D dispersive LRSWE The rotating shallow water equation based on single-layer approximation [22,30,31] are used for the modeling of atmosphere and ocean dynamics. The shallow water equation given by Eq. (3) form a dispersive hyperbolic system of conservation laws and used to model free surface flows in layers using the hydrostatic pressure approximation. Moreover, shallow water equation describe the flow of a fluid layer of constant density, where horizontal scale is large enough as compared to fluid depth. In vector notations, 2D shallow water equation linearized about mean state are given as

zt + Azx + Bz y + C z = 0 where z = [u , v ,

η]T , and

(18)

P.K. Maurya et al. / Journal of Computational Physics 399 (2019) 108941

11

Fig. 3. Contour plots showing numerical properties in the spectral plane for physical and computational modes for IIIFL − OUCS3 scheme for standard and optimized values of parameters.

⎡

0 A=⎣ 0 H

0 0 0

⎤

⎡

g 0 0 ⎦, B = ⎣ 0 0 0

0 0 H

⎤

⎡

0 0 g ⎦ and C = ⎣ f 0 0

−f 0 0

⎤

0 g⎦ 0

here, (u , v ) are the perturbation velocity component, η is the perturbation surface elevation, and H represent the mean height. Using the bilateral Fourier-Laplace transform [32], the dispersion relation for Eq. (18) is given as

 ( 2 − c 2 | K |2 − f 2 ) = 0

(19)

12

P.K. Maurya et al. / Journal of Computational Physics 399 (2019) 108941

Fig. 4. Exact (dotted lines) and numerical (solid lines) solutions to one-dimensional advection equation using indicated schemes with N c = 0.5 and ko h = 0.5.

√

 denotes the circular frequency, K = (kx , k y ) and c =  g H , with kx and k y representing wavenumber components. Here, 1 = 0 represents the geostrophic mode, and 2,3 = ± f 2 + c 2 | K |2 symbolize the inertia-gravity modes. Group where,

velocity and phase speeds components for the inertia-gravity modes [32,11] are obtained from Eq. (19) as ( V gx, y )2,3 = ±(kx, y c 2 )/( f 2 + c 2 | K |2 )1/2 and (c x, y )2,3 = ±( f 2 + c 2 | K |2 )1/2 /kx, y , and the resultant group velocity for the inertia-gravity modes is given as, V g = c 2 | K |/( f 2 + c 2 | K |2 )1/2 . Mesinger & Arakawa [33] suggested various collocated/staggered meshes, and among them C -grid remains more up-andcoming [32,11]. In the literature, various forms of C -grid with rectangular, triangular and hexagonal [34–38] cell structures have been investigated. Moreover, for the Navier-Stokes equation C -grid is equivalent to the mesh used in Marker and Cell (MAC) [39] method. Using bilateral Fourier-Laplace transform, Eq. (18) with spatial discretization scheme can be represented as

P.K. Maurya et al. / Journal of Computational Physics 399 (2019) 108941

⎡ ∂ Zˆ /∂ t + D Zˆ = 0, with D = ⎣

0 f eq i Hk xeq

13

⎤

− f eq

igkxeq igk y eq ⎦ 0

0 i Hk y eq

ˆ ] T , and kxeq , k yeq representing equivalent wavenumbers. For two-level time advancement method, the where, Zˆ = [Uˆ , Vˆ , H

expression combining values of unknown Zˆ at (n + 1)th and nth levels is given as, Zˆ n+1 = Q Zˆ n . The stability/instability for geostrophic and inertia-gravity modes are given by eigenvalues of corresponding evolution matrix [32]. The evolution matrix ( Q ) for the IFL is obtained from roots given in Appendix A by replacing iN c φ(kh) with t D, and is given by

Q 1,2 IFL =

1

(1 + 2γ ω1 )

  ±σ + γ I − (1 + γ ω1 )t D

with



σ = (γ 2 ω12 + 2γ ω1 + 1)t 2 D 2 + 2γ ω1 (−2γ ω1 + γ − 1)t D + (4γ 2 ω12 − 4γ 2 ω1 + 4γ ω1 + γ 2 − 2γ + 1) I

 12

were, I denotes 3 × 3 identity matrix. Here, eigenvalues of Q 1IFL represent the nodal amplification factors for physical geostrophic and inertia-gravity modes, while eigenvalues of Q 2IFL represent spurious modes corresponding to geostrophic and inertia-gravity modes. Amplification factors for the higher order filtered schemes can be obtained in the similar manner. On staggered meshes, optimized staggered compact scheme (OSCS) with optimized staggered interpolation scheme for spatial discretization [32] derived from the staggered compact scheme given in [40], is used. Staggered compact scheme for evaluating the first-order spatial derivative at jth node is given as [40]

α1 u j−1 + u j + α1 u j+1 = b1



u j +3/2 − u j −3/2





+a1

3 x

u j +1/2 − u j −1/2



(20)

x

representing a single parameter (α1 ) family of fourth-order schemes with a1 = 38 (3 − 2α1 ) and b1 = 18 (−1 + 22α1 ). For the present case we have chosen α1 = 0.18, as in [32]. As, staggered meshes also require information at midpoint locations obtained by interpolation, given by [40]

α2 u Ij−1 + u Ij + α2 u Ij+1 =

b2 2





u j +3/2 + u j −3/2 +

a2



2



u j +1/2 + u j −1/2

(21)

where u Ij denotes interpolated value of u at jth node. This approximation is of fourth-order accuracy if a2 = 18 (9 + 10α2 ) and b2 = 18 (−1 + 6α2 ). In the present case, the optimized value α2 = 0.35 [32] is used for the computations. The resultant equivalent wavenumbers and transfer functions for interpolation using Eqs. (18)-(19), respectively, are given as

 φ(ki i ) = 2

a1 sin(ki i /2) +

 T F i (ki i ) =

b1 3

sin(3ki i /2)



1 + 2α1 cos(ki i ) a2 cos(ki i /2) + b2 cos(3ki i /2) 1 + 2α2 cos(ki i )

 , i = x, y

where x and  y denote the mesh-width in x− and y − directions, respectively. For C-grid, evaluation of the first-order derivatives does not require midpoint interpolation. However, momentum conservation equations having terms f v and f u in Eq. (18) require midpoint interpolations in x- and y-directions, which give rise to, f eq t = f t T F x (kx x) T F y (k y  y ), for C-grid, and for A-grid it is equal to f t. To access the dispersion properties of different numerical schemes, propagation of a 2D wave-packet governed by LRSWE is considered next. The 2D wave-packet given by

u (x, y , t = 0) = 0, v (x, y , t = 0) = 0

η(x, y , t = 0) = e−0.05

 (x−x0 )2 15

+( y − y 0 )2



sin(kx x + k y y )

(22)

For the present case, values of the required parameters as given in [32] are used. Here, computational domain of size (360 × 360) is considered with uniform mesh-width. As, in [32,41], stretched wave-packet is chosen to manifest the effects of q-waves with propagation angle, θ = 45o and aspect ratio of λ =  y /x = 1. Other parameters are chosen as, g = 10 m/s2 , H = 2.5 m, x =  y = 0.3 m, f t = 1.2 × 10−2 and Nc x = Nc y = 0.2. Initially, wave-packet is located at (k0 x, k0  y ) = (2.3365, 2.3365), which is close to the inception limit (k0 x = 2.3910) of q-waves for OUCS3. Thus, there exist wavenumbers in the proximity of (k0 x, k0  y), which will provide the source for the genesis of q-waves. Computations of this propagation problem are performed for collocated A- and staggered C-grids, and compared with the exact solutions obtained via 2D FFT of initial condition. Numerical solutions using IIFL, IIIFL & IVFL with OUCS3 on A-grid are

14

P.K. Maurya et al. / Journal of Computational Physics 399 (2019) 108941

Fig. 5. Exact and numerical solutions to two-dimensional dispersive shallow water equations on collocated A- and staggered C-grids using indicated schemes with N c = 0.2 and ko h = 2.39.

shown in the left column frames of Fig. 5. It is evident from the left frames of Fig. 5 that for chosen values of parameters, numerical solutions admit q-waves on collocated A-grid, which skewed the numerical solutions at higher time instant. To further, gauge the efficiency of the developed schemes, numerical solutions using IIFL, IIIFL & IVFL with OSCS on the staggered C-grid are shown in the right column frames in Fig. 5. From max/min values marked in each sub-frame, it is observed from Fig. 5 that IIFL - OSCS generates higher diffusion in comparison to IIIFL - OSCS, and IVFL - OSCS provides the best results showing the propagation of inertia-gravity modes with correct speed. This is in the conformity with the dispersion properties of the OSCS scheme on staggered C-grid, as reported in [32]. Next, DNS of the 2D unsteady Navier-Stokes equation for the flow inside LDC problem is solved at different Reynolds numbers.

P.K. Maurya et al. / Journal of Computational Physics 399 (2019) 108941

15

Table 5 Maximum and minimum values of stream-function (ψ) using indicated methods for flow inside LDC at Re = 1, 000. Methods, mesh-size & time-step (t )

ψmax

ψmin

Botella & Peyret [48] ( N = 160) Ghia et al. [49] (129 × 129) IFL - OUCS3 (300 × 300), t = 5 × 10−4 IFL - OUCS3 (400 × 400), t = 1 × 10−4 IFL - OUCS3 (512 × 512), t = 1 × 10−4

0.11893 0.11792 0.11863 0.11875 0.11882

−1.7297 × 10−3 −1.7510 × 10−3 −1.7582 × 10−3 −1.7456 × 10−3 −1.7394 × 10−3

8. Direct numerical simulation of 2D LDC problem DNS of the flow inside lid-driven cavity is carried out using stream function-vorticity (ψ − ω) formulation of the NavierStokes equation as given by Eq. (2). Numerical simulation exhibiting typical features of flow inside LDC are studied in [42,43] and the experimental evidences for the presence of polygonal/triangular core vortex are given in [44,45]. Being the non-periodic problem, DNS of flow inside LDC require a near-spectral accurate compact scheme with appropriate boundary closures. As pointed out in [23,43], compact scheme in [28] perform poorly for non-periodic problems and this is due to the lack of suitable boundary closures. Detailed grid sensitivity analysis and various sources of error for the LDC are presented in [46] at different Reynolds numbers. Here, in the vorticity transport equation given by Eq. (2), first- and second-order spatial derivatives are approximated using OUCS3 and CD2 schemes, respectively. Also, to deal with the non-periodic problem, explicit boundary and near-boundary closures [26] for the OUCS3 scheme are given as

u 1 =

1

1

(−3u 1 + 4u 2 − u 3 ), u N = (3u N − 4u N −1 + u N −2 ) 2h 2h        1 2β2 1 8 β2 1 8β2 1 2β2 u 2 = − u1 − + u 2 + (4β2 + 1)u 3 − + u4 + u5 h

u N −1 = −

3

1 h



3

2β N 3

−

1 3



3

uN −



2

8β N 3

+

1 2



3

6

u N −1 + (4β N + 1)u N −2 −



3

8β N 3

+

1 6

 u N −3 +

2β N 3

 u N −4

(23)

Equation (23) represent the boundary stencils and for global accuracy and numerical stability, where β2 = −0.025 and β N = 0.09. Furthermore, temporal advancement is performed using hybrid filtered leapfrog time integration method. A variant of conjugate gradient method known as, Bi-CGSTAB method [47], is used in solving the stream function equation and the residual convergence criteria is held as 10−6 . First, this problem is solved using IFL - OUCS3 for Re = 1, 000 with 512 × 512 mesh-points and the stream function and vorticity contours at marked time instants are shown in Fig. 6. Numerical results for Re = 1, 000 using different mesh-sizes are also compared with the benchmark results reported by Botella & Peyret [48] by comparing the maximum and minimum values of ψ inside the cavity, as given in Table 5. Additionally, computed results are also compared with those reported in [49] for Re = 1, 000. For all the reported results, solutions are advanced in time until ∂ ω/∂ t reduced to 10−12 . It is observed from Table 5 that there is good match with the benchmark results reported in [48] up to the three decimal points for optimized values of filter parameters. For the genuine non-steady case, LDC problem is also solved for Re = 10, 000 and vorticity contours are shown in Fig. 7 with 300 × 300 mesh-points using IFL and IIIFL with OUCS3 scheme. Results are shown till time t = 700 exhibiting triangular/polygonal vortical structure in the core of cavity starting at t = 400 surrounded by revolving secondary vortices, as in [50,42]. For the IFL - OUCS3, due to higher in-build numerical diffusion the central vortical structure disappears at higher time instants, which is not the case for IIIFL - OUCS3. The behavior shown in Fig. 7 is in consonance with the analysis proclaimed in [51] using proper orthogonal decomposition (POD) conforming the presence of various polygonal vortices as the eigensolutions. 9. Summary and conclusion In this paper, a family of hybrid filtered leapfrog methods are developed, and to achieve the better numerical properties a constrained optimization is considered in the spectral plane. Developed methods also have been used for the numerical solution to 1D non-dispersive, 2D dispersive linear model dynamical systems. Moreover to assess the robustness of the presented methods for the nonlinear case, flow inside lid-driven cavity governed by Navier-Stokes equation is also considered at various Reynolds number. It is observed that for standard values of the filter parameters, there is small instability in the physical mode and the dominating computational mode is also not adroitly damped. However, the instability issue is undertaken with the optimized values of the filter parameters, as the physical mode is near neutrally stable and computational modes are highly damped.

16

P.K. Maurya et al. / Journal of Computational Physics 399 (2019) 108941

Fig. 6. Stream function-vorticity contours using IFL - OUCS3 for the flow inside driven-cavity problem at Re = 1000 with (512 × 512) mesh-points.

As, the optimized hybrid filters strikingly attenuate the computational mode and also the spectral weight [23] associated with computational mode is sufficiently small, therefore, existence of computational modes would not influence the numerical solutions significantly. Furthermore, propagation problem following the dispersive rotating shallow water equation is also solved for collocated (A) and staggered (C) grids and the numerical solutions are also compared with exact solutions. Comparison of numerical solutions with the results reported in the literature display good match, which authenticate the veracity of the developed methods. Finally, optimized hybrid filters are employed for the direct numerical simulation of Navier-Stokes equation for flow inside 2D lid-driven cavity at different Reynolds numbers. Computed solutions display the good match with the benchmark results reported in the literature. Implementation of present methods to the physical modeling of ocean and atmospheric models will be reported in the near future. The computational cost of the present hybrid filtered leapfrog methods is comparable to the explicit two-stage Runge-Kutta method, and computational complexity of these schemes are of same type as of two-stage second-order Runge-Kutta methods- with no extra storage requirement. However, the developed methods allow the use of higher time-step in the computations as compared to Runge-Kutta method. 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 Authors gratefully acknowledge the financial support from the SERB, Government of India, grant ECR/2016/000381. First author also acknowledge the support from CSIR (grant 09/1088(0003)/2017-EMR-I), New Delhi. Appendix A. Roots of the hybrid filters using symbolic toolbox The roots of the quadratic expression given by Eq. (10) are computed as

G 1,2 =

1

(1 + 2γ ω1 )

  ±σ1 + γ − iN c φ(kh)(1 + γ ω1 )

(24)

P.K. Maurya et al. / Journal of Computational Physics 399 (2019) 108941

17

Fig. 7. Stream function-vorticity contours using marked schemes for the flow inside driven-cavity problem at Re = 10, 000 with (300 × 300) mesh-points and t = 1 × 10−3 .

18

P.K. Maurya et al. / Journal of Computational Physics 399 (2019) 108941

with



σ1 = N c2 φ 2 (kh)(−γ 2 ω12 − 2γ ω1 − 1) + N c φ(kh)(−4i γ 2 ω12 + 2i γ 2 ω1 − 2i γ ω1 )

1 + (4γ 2 ω12 − 4γ 2 ω1 + 4γ ω1 + γ 2 − 2γ + 1) 2

The three roots of the depicting the nodal amplification factors for second hybrid filter (IIFL) are computed as

√ 



−σ1 σ3 3 σ3 σ3 G 1 = σ1 − σ2 − , G 2,3 = ± σ1 − i − σ2 − σ1 2 2 σ1 2σ1 where,







σ1 = −σ5 + −σ33 + σ5 +

2  12 1 σ63 σ63 + σ − − σ4 3 4 3 3 27σ8 27σ8

σ2 σ6 σ7 σ6 σ7 i γ N c φ(kh)(1 − ω2 ) , σ3 = + 6 2 , σ4 = , σ5 = 3 σ8 3σ8 σ8 9σ8 6σ82 σ6 = −4γ + 2i γ N c φ(kh)(1 + 3γ ω2 ), σ7 = 4γ (ω2 − 1) + 2i γ N c φ(kh)(3 − 2ω2 ) + 1 σ2 =

σ8 = 1 + 4γ ω2 Finally, the roots of the quartic expression for the nodal amplification factors for the third hybrid filter (IIIFL) are given as

G 1,2 = −σ3 +

σ12 4σ15

± σ1 , G 3,4 = σ3 +

σ12 4σ15

± σ2

where,



σ1 =

2

−9σ63

√

√

√

1



σ2 =

1

√

 12

1

√

 12

σ5 − σ92 σ5 + 12σ8 σ5 + 12σ9 σ63 σ5 − σ4 1

6σ66 σ54 2

−9σ63

√

√

√

σ5 − σ92 σ5 + 12σ8 σ5 + 12σ9 σ63 σ5 + σ4 1

1

6σ66 σ54 1

σ3 =

σ52 1 6

√ √  1 2 2 , σ4 = 3 6σ10 −72σ9 σ8 + 3 3σ7 − 2σ93 + 27σ10

6 σ6

σ5 = − σ6 = −

12σ11

σ15 4σ9 σ8 3

1

2

+ 6σ9 σ63 + 9σ63 + σ92 − √ +



3σ7

18

−

σ93 27

+

2 σ10

2

4 9σ12 4 15

64σ

, σ10 =

−

3σ12 σ14 2 15

σ

−

2 3σ12 σ13 3 4σ15

σ3 σ14 σ12 σ13 + 123 + 2 σ15 8σ15 2σ15

2 4 2 σ7 = −4σ93 σ10 + 16σ94 σ8 + 256σ83 + 128σ92 σ82 + 27σ10 − 144σ9 σ10 σ8

 12

4 2 2 3σ12 σ13 σ11 σ12 σ14 σ12 σ13 3σ12 + + + , σ9 = + 4 2 3 2 σ15 256σ15 σ15 8σ15 4σ15 16σ15 σ11 = 2i γ N c φ(kh)(1 − ω3 ), σ12 = 7γ + γ ω3 − 2iN c φ(kh)(1 + 6γ ω3 )

σ8 =

σ13 = 1 − 8γ + 7γ ω3 + 4i γ N c φ(kh)(3 − ω3 ), σ14 = γ − γ ω3 − 2i γ N c φ(kh)(4 − 3ω3 ) σ15 = 1 + 7γ ω3 References [1] R. Vichnevetsky, J.B. Bowles, Fourier Analysis of Numerical Approximations of Hyperbolic Equations, SIAM, Philadelphia, 1982. [2] G. Strang, G.J. Fix, An Analysis of the Finite Element Method, second ed., Wellesley-Cambridge, 2008. [3] R.L. Pfeffer, I.M. Navon, X.L. Zou, A comparison of the impact of two time-differencing schemes on the NASA GLAS climate model, Mon. Weather Rev. 120 (1992) 1381–1393. [4] R.X. Huang, J. Pedlosky, On aliasing Rossby waves induced by asynchronous time stepping, Ocean Model. 5 (2003) 65–76. [5] B. Zhao, Q. Zhong, The dynamical and climate tests of an atmospheric general circulation model using the second order Adams-Bashforth method, Acta Meteorol. Sin. 23 (2009) 738–749.

P.K. Maurya et al. / Journal of Computational Physics 399 (2019) 108941

19

[6] J. Teixeira, C.A. Raynolds, K. Judd, Time stepping sensitivity of nonlinear atmospheric models: numerical convergence, truncation error growth, and ensemble design, J. Atmos. Sci. 64 (2007) 175–189. [7] A.J. Robert, The integration of a low order spectral form of the primitive meteorological equations, J. Meteorol. Soc. Jpn. 44 (1966) 237–244. [8] R. Asselin, Frequency filter for time integrations, Mon. Weather Rev. 100 (1972) 487–490. [9] P.D. Williams, A proposed modification to the Robert-Asselin time filter, Mon. Weather Rev. 137 (2009) 2538–2546. [10] Y. Li, C. Trenchea, A higher order Robert-Asselin type time filter, J. Comput. Phys. 259 (2014) 23–32. [11] M.K. Rajpoot, Dispersion analysis of Robert-Asselin type filters for linear non-dispersive and dispersive systems, Comput. Fluids 130 (2016) 49–83. [12] P.D. Williams, Achieving seventh-order amplitude accuracy in leapfrog integrations, Mon. Weather Rev. 141 (2013) 3037–3051. [13] P.D. Williams, The RAW filter: an improvement to the Robert-Asselin filter in semi-implicit integrations, Mon. Weather Rev. 139 (2011) 1996–2007. [14] T.K. Sengupta, A critical assessment of simulations for transitional and turbulent flows, in: T.K. Sengupta, S.K. Lele, K.R. Sreenivasan, P.A. Davidson (Eds.), Proc. of the IUTAM Symp. on “Advances in Computations, Modeling and Control of Transitional and Turbulent Flows”, 2015, pp. 491–532. [15] P. Marsaleix, F. Auclair, T. Duhaut, C. Estournel, C. Nguyen, Alternatives to the Robert-Asselin filter, Ocean Model. 41 (2012) 53–66. [16] M. Moustaoui, A. Mahalov, E.J. Kostelich, A numerical method based on leapfrog and a fourth-order implicit time filter, Mon. Weather Rev. 142 (2014) 2545–2560. [17] A. Mahalov, M. Moustaoui, Time-filtered leapfrog integration of Maxwell equations using unstaggered temporal grids, J. Comput. Phys. 325 (2016) 98–115. [18] N. Sharma, A. Sengupta, M.K. Rajpoot, R.J. Samuel, T.K. Sengupta, Hybrid sixth order spatial discretization scheme for non-uniform Cartesian grids, Comput. Fluids 157 (2017) 208–231. [19] M.K. Rajpoot, T.K. Sengupta, P.K. Dutt, Optimal time advancing dispersion relation preserving schemes, J. Comput. Phys. 229 (2010) 3623–3651. [20] T.K. Sengupta, M.K. Rajpoot, Y.G. Bhumkar, Space-time discretizing optimal DRP schemes for flow and wave propagation problems, Comput. Fluids 47 (2011) 144–154. [21] T.K. Sengupta, A. Dipankar, P. Sagaut, Error dynamics: Beyond von Neumann analysis, J. Comput. Phys. 226 (2007) 1211–1218. [22] J. Pedlosky, Geophysical Fluid Dynamics, Springer Verlag, 1979. [23] T.K. Sengupta, High Accuracy Computing Methods: Fluid Flows and Wave Phenomena, Cambridge Univ. Press, 2013. [24] T. Kawamura, H. Takami, K. Kuwahara, A new higher-order upwind scheme for incompressible Navier-Stokes equations, Fluid Dyn. Res. 1 (1985) 145–162. [25] B.V.D. Pol, H. Bremmer, Operational Calculus Based on the Two-Sided Laplace Transform, Cambridge Univ. Press, 2008. [26] T.K. Sengupta, G. Ganeriwal, S. De, Analysis of central and upwind compact schemes, J. Comput. Phys. 192 (2003) 677–694. [27] T.K. Sengupta, A. Dipankar, A comparative study of time advancement methods for solving Navier-Stokes equations, J. Sci. Comput. 21 (2004) 225–250. [28] S.K. Lele, Compact finite difference schemes with spectral-like resolution, J. Comput. Phys. 103 (1992) 16–42. [29] Z. Haras, S. Ta’asan, Finite difference schemes for long-time integration, J. Comput. Phys. 114 (2) (1994) 265–279. [30] A.E. Gill, Atmosphere-Ocean Dynamics, International Geophysics Series, Academic Press, New York, 1982. [31] G.K. Vallis, Atmosphere and Ocean Fluid Dynamics: Fundamentals and Large Scale Circulations, Cambridge Univ. Press, 2006. [32] M.K. Rajpoot, S. Bhaumik, T.K. Sengupta, Solution of linearized rotating shallow water equations by compact schemes with different grid-staggering strategies, J. Comput. Phys. 231 (2012) 2300–2327. [33] F. Mesinger, A. Arakawa, Numerical Methods Used in Atmospheric Models, vol. 1, GARP Publ. Ser., vol. 17, WMO, Geneva, 1976. [34] A.J. Adcroft, C.N. Hill, J.C. Marshall, A new treatment of the Coriolis terms in C-grid models at both high and low-resolutions, Mon. Weather Rev. 127 (1999) 1928–1936. [35] V. Casulli, R. Walters, An unstructured grid, three dimensional model based on shallow water equations, Int. J. Numer. Methods Fluids 32 (2000) 331–341. [36] R. Nicolaides, Direct discretization of planar div-curl problems, SIAM J. Numer. Anal. 29 (1992) 32–56. [37] R. Walters, V. Casulli, A robust finite element model for hydrostatic surface-water flows, Commun. Numer. Methods Eng. 14 (2001) 931–940. [38] H. Shirkhani, A. Mohammadian, O. Seidou, H. Qiblawey, Comparison of 2D triangular C-grid shallow water models, Comput. Fluids 161 (2018) 136–154. [39] F.H. Harlow, J.E. Welch, Numerical calculation of time dependent viscous incompressible flow of fluids with free surface, Phys. Fluids 8 (1965) 2182–2189. [40] S. Nagarajan, S.K. Lele, J.H. Ferziger, A robust high-order compact method for large eddy simulation, J. Comput. Phys. 19 (2003) 392–419. [41] T.K. Sengupta, Y.G. Bhumkar, M.K. Rajpoot, V.K. Suman, S. Saurabh, Spurious waves in discrete computation of wave phenomena and flow problems, Appl. Math. Comput. 218 (2012) 9035–9065. [42] T.K. Sengupta, V.V.S.N. Vijay, S. Bhaumik, Further improvement and analysis of CCD scheme: dissipation discretization and de-aliasing properties, J. Comput. Phys. 228 (2009) 6150–6168. [43] L. Lestandi, S. Bhaumik, G.R.K.C. Avatar, M. Azaiez, T.K. Sengupta, Multiple Hopf bifurcations and flow dynamics inside a 2D singular lid driven cavity, Comput. Fluids 166 (2018) 86–103. [44] M. Beckers, G.J.F.V. Heijst, The observation of a triangular vortex in a rotating fluid, Fluid Dyn. Res. 22 (1998) 265–279. [45] G.F. Carnevale, R.C. Koolsterziel, Emergence and evolution of triangular vortices, J. Fluid Mech. 259 (1994) 305–331. [46] V.K. Suman, S. Viknesh, M.K. Tekriwal, S. Bhaumik, T.K. Sengupta, Grid sensitivity and role of error in computing a lid-driven cavity problem, Phys. Rev. E 99 (2019) 013305. [47] H.A.V.D. Vorst, BI-CGSTAB: a fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear system, SIAM J. Sci. Stat. Comput. 13 (1992) 631–644. [48] O. Botella, R. Peyret, Benchmark spectral results on the lid-driven cavity flow, Comput. Fluids 27 (1998) 421–433. [49] U. Ghia, K.N. Ghia, C.T. Shin, High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method, J. Comput. Phys. 48 (1982) 387–411. [50] T.K. Sengupta, V. Lakshmanan, V.V.S.N. Vijay, A new combined stable and dispersion relation preserving compact scheme for non-periodic problems, J. Comput. Phys. 228 (2009) 3048–3071. [51] T.K. Sengupta, V.V.S.N. Vijay, N. Singh, Universal instability modes in internal and external flows, Comput. Fluids 40 (2011) 221–235.