High Order Semi-Lagrangian Methods for the Incompressible Navier–Stokes Equations | Journal of Scientific Computing Skip to main content
Log in

High Order Semi-Lagrangian Methods for the Incompressible Navier–Stokes Equations

  • Published:
Journal of Scientific Computing Aims and scope Submit manuscript

Abstract

We propose a class of semi-Lagrangian methods of high approximation order in space and time, based on spectral element space discretizations and exponential integrators of Runge–Kutta type. The methods were presented in Celledoni and Kometa (J Sci Comput 41(1):139–164, 2009) for simpler convection–diffusion equations. We discuss the extension of these methods to the Navier–Stokes equations, and their implementation using projections. Semi-Lagrangian methods up to order three are implemented and tested on various examples. The good performance of the methods for convection-dominated problems is demonstrated with numerical experiments.

This is a preview of subscription content, log in via an institution to check access.

Access this article

Subscribe and save

Springer+ Basic
¥17,985 /Month
  • Get 10 units per month
  • Download Article/Chapter or eBook
  • 1 Unit = 1 Article or 1 Chapter
  • Cancel anytime
Subscribe now

Buy Now

Price includes VAT (Japan)

Instant access to the full article PDF.

Fig. 1
Fig. 2
Fig. 3
Fig. 4
Fig. 5
Fig. 6
Fig. 7
Fig. 8

Similar content being viewed by others

Notes

  1. We also assume that the characteristic paths are integrated to high accuracy.

  2. Extension to higher order methods is straightforward (See [30] and references therein).

  3. See 6.2 for the definitions of these norms.

  4. The projection does not need to be orthogonal, but should be guaranteed not to compromise the order of the method, an orthogonal projection will have this property. The target of the orthogonal projection map on the discrete divergence-free subspace is the element of shortest distance to the point which is projected. Since \(y_{n+1}=\Pi Y_s\) and \(\Vert y(t_{n+1})-Y_s\Vert = {\mathcal {O}}(h^{r+1}),\) where \(r\) is the order of the integration method, then

    $$\begin{aligned} \Vert y(t_{n+1})-y_{n+1}\Vert \le \Vert y(t_{n+1})-Y_s\Vert +\Vert Y_s-y_{n+1}\Vert = {\mathcal {O}}(h^{r+1}), \end{aligned}$$

    because

    $$\begin{aligned} \Vert y_{n+1}-Y_s\Vert \le \Vert y(t_{n+1})-Y_s\Vert . \end{aligned}$$

References

  1. Ascher, U.M., Ruuth, S.J., Spiteri, R.J.: Implicit-explicit Runge–Kutta methods for time-dependent partial differential equations. Appl. Numer. Math. 25, 151–167 (1997)

    Article  MathSciNet  MATH  Google Scholar 

  2. Ascher, U.M., Ruuth, S.J., Wetton, B.T.R.: Implicit–explicit methods for time-dependent partial differential equations. SIAM J. Numer. Anal. 32(3), 797–823 (1995)

    Article  MathSciNet  MATH  Google Scholar 

  3. Brown, D.L., Minion, M.L.: Performance of under-resolved two-dimensional incompressible flow simulations. J. Comput. Phys. 122(1), 165–183 (1995)

    Article  MathSciNet  MATH  Google Scholar 

  4. Bruneau, C.-H., Saad, M.: The 2d lid-driven cavity problem revisited. Comput. Fluids 35(3), 326–348 (2006)

    Article  MATH  Google Scholar 

  5. Canuto, C., Hussaini, M.Y., Quarteroni, A., Zang, T.A.: Spectral Methods in Fluid Dynamics. Springer Series in Computational Physics. Springer, New York (1988)

    Book  Google Scholar 

  6. Celledoni, E.: Eulerian and semi-Lagrangian commutator-free exponential integrators. CRM Proc. Lect. Notes 39, 77–90 (2005)

    MathSciNet  Google Scholar 

  7. Celledoni, E., Kometa, B.K.: Semi-Lagrangian Runge–Kutta exponential integrators for convection dominated problems. J. Sci. Comput. 41(1), 139–164 (2009)

    Article  MathSciNet  MATH  Google Scholar 

  8. Celledoni, E., Kometa, B.K.: Order conditions for semi-Lagrangian Runge-Kutta exponential integrators, Preprint series: Numerics, 04/2009, Department of Mathematics, NTNU, Trondheim, Norway. http://www.math.ntnu.no/preprint/numerics/N4-2009.pdf (2009)

  9. Celledoni, E., Kometa, B.K.: Semi-Lagrangian multistep exponential integrators for index 2 differential-algebraic systems. J. Comput. Phys. 230(9), 3413–3429 (2011)

    Article  MathSciNet  MATH  Google Scholar 

  10. Childs, P.N., Morton, K.W.: Characteristic Galerkin methods for scalar conservation laws in one dimension. SIAM J. Numer. Anal. 27, 553–594 (1990)

    Article  MathSciNet  MATH  Google Scholar 

  11. Chorin, A.J.: Numerical solution of the Navier–Stokes equations. Math. Comput. 22, 745–762 (1968)

    Article  MathSciNet  MATH  Google Scholar 

  12. Chorin, A.J.: On the convergence of discrete approximations to the Navier–Stokes equations. Math. Comput. 23, 341–353 (1969)

    Article  MathSciNet  MATH  Google Scholar 

  13. Deville, M.O., Fischer, P.F., Mund, E.H.: High-Order Methods for Incompressible Fluid Flow. Cambridge Monographs on Applied and Computational Mathematics. Cambridge University Press, Cambridge (2002)

    Book  Google Scholar 

  14. Eich-Soellner, E., Führer, C.: Numerical Methods in Multibody Dynamics, European Consortium for Mathematics in Industry. B. G. Teubner, Stuttgart (1998)

    Book  Google Scholar 

  15. Falcone, M., Ferretti, R.: Convergence analysis for a class of high-order semi-Lagrangian advection schemes. SIAM J. Numer. Anal. 35(3), 909–940 (1998). (electronic)

    Article  MathSciNet  MATH  Google Scholar 

  16. Falcone, M., Ferretti, R.: Semi-Lagrangian approximation schemes for linear and Hamilton-Jacobi equations. SIAM (to appear)

  17. Fischer, P., Mullen, J.: Filter-based stabilization of spectral element methods. C. R. Acad. Sci. Paris Sér. I Math. 332(3), 265–270 (2001)

    Article  MathSciNet  MATH  Google Scholar 

  18. Fischer, P.F.: An overlapping Schwarz method for spectral element solution of the incompressible Navier–Stokes equations. J. Comput. Phys. 133(1), 84–101 (1997)

    Article  MathSciNet  MATH  Google Scholar 

  19. Fischer, P.F., Kruse, G.W., Loth, F.: Spectral element method for transitional flows in complex geometries. J. Sci. Comput. 17(1–4), 81–98 (2002)

    Article  MathSciNet  MATH  Google Scholar 

  20. Ghia, U., Ghia, K.N., Shin, C.T.: High re-solutions for incompressible flow using the Navier–stokes equations and a multigrid method. J. Comput. Phys. 48(3), 387–411 (1982)

    Article  MATH  Google Scholar 

  21. Giraldo, F.X.: The Lagrange–Galerkin spectral element method on unstructured quadrilateral grids. J. Comput. Phys. 147(1), 114–146 (1998)

    Article  MathSciNet  MATH  Google Scholar 

  22. Giraldo, F.X., Perot, J.B., Fischer, P.F.: A spectral element semi-Lagrangian (SESL) method for the spherical shallow water equations. J. Comput. Phys. 190(2), 623–650 (2003)

    Article  MathSciNet  MATH  Google Scholar 

  23. Guermond, J.L., Shen, J.: A new class of truly consistent splitting schemes for incompressible flows. J. Comput. Phys. 192(1), 262–276 (2003)

    Article  MathSciNet  MATH  Google Scholar 

  24. Guermond, J.L., Minev, P.D., Shen, J.: An overview of projection methods for incompressible flows. Comput. Methods Appl. Mech. Eng. 195, 6011–6045 (2006). doi:10.1016/j.cma.2005.10.010

    Article  MathSciNet  MATH  Google Scholar 

  25. Hairer, E., Lubich, Ch., Wanner, G.: Geometric Numerical Integration, Second ed., Springer Series in Computational Mathematics, vol. 31, Structure-Preserving Algorithms for Ordinary Differential Equations. Springer, Berlin (2006)

    Google Scholar 

  26. Hairer, E., Wanner, G.: Solving Ordinary Differential Equations. II, Second ed., Springer Series in Computational Mathematics. Springer, Berlin (1996)

    Book  Google Scholar 

  27. Kanevsky, A., Carpenter, M.H., Gottlieb, D., Hesthaven, J.S.: Application of implicit–explicit high order Runge–Kutta methods to discontinuous-Galerkin schemes. J. Comput. Phys. 225(2), 1753–1781 (2007)

    Article  MathSciNet  MATH  Google Scholar 

  28. Kennedy, C.A., Carpenter, M.H.: Additive Runge–Kutta schemes for convection–diffusion–reaction equations. Appl. Numer. Math. 44(1–2), 139–181 (2003)

    Article  MathSciNet  MATH  Google Scholar 

  29. Kometa, B.K.: Semi-Lagrangian methods and new integration schemes for convection-dominated problems, PhD thesis at NTNU, ISSN 1503–8181; 2011:270, Norwegian University of Science and Technology, Department of Mathematical Sciences, http://urn.kb.se/resolve?urn=urn:nbn:no:ntnu:diva-15729 (2011)

  30. Kvarving, A.M.: Splitting schemes for the unsteady Stokes equations: a comparison study, Preprint series: Numerics, 05/2010, Department of Mathematics, NTNU, Trondheim, Norway, http://www.math.ntnu.no/preprint/numerics/N5-2010.pdf (2010)

  31. Maday, Y., Patera, A.T., Rønquist, E.M.: An operator-integration-factor splitting method for time-dependent problems: application to incompressible fluid flow. J. Sci. Comput. 5(4), 263–292 (1990)

    Article  MathSciNet  MATH  Google Scholar 

  32. Morton, K.W.: Generalised Galerkin methods for hyperbolic problems. Computer Methods Appl Mech Eng 52, 847–871 (1985)

    Article  MathSciNet  MATH  Google Scholar 

  33. Rannacher, R.: On Chorin’s projection method for the incompressible Navier–Stokes equations, Lecture Notes in Mathematics, vol. 1530. Springer, Berlin (1991)

    Google Scholar 

  34. Restelli, M., Bonaventura, L., Sacco, R.: A semi-Lagrangian discontinuous Galerkin method for scalar advection by incompressible flows. J. Comput. Phys. 216(1), 195–215 (2006)

    Article  MathSciNet  MATH  Google Scholar 

  35. Shen, J.: On error estimates of projection methods for Navier–Stokes equations: first-order schemes. SIAM J. Numer. Anal. 29(1), 57–77 (1992)

    Article  MathSciNet  MATH  Google Scholar 

  36. Témam, R.: Sur l’approximation de la solution des équations de Navier–Stokes par la méthode des pas fractionnaires. II. Arch. Ration. Mech. Anal. 33, 377–385 (1969)

    Article  MATH  Google Scholar 

  37. Xiu, D., Karniadakis, G.E.: A semi-Lagrangian high-order method for Navier–Stokes equations. J. Comput. Phys. 172(2), 658–684 (2001)

    Article  MathSciNet  MATH  Google Scholar 

Download references

Acknowledgments

This work was supported in part by the GeNuIn project, Grant from the Research Council of Norway.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Olivier Verdier.

Appendix

Appendix

1.1 IMEX and DIRK-CF Methods for Convection–Diffusion Equations

A DIRK-CF method is a IMEX Runge–Kutta method of exponential type. It is defined by two Runge–Kutta tableaus one of implicit type to treat the linear diffusion and one of explicit type to treat the nonlinear convection by composition of exponentials. These methods have a Runge–Kutta like format with two sets of parameters:

$$\begin{aligned} {\mathcal {A}}=\{a_{i,j}\}_{i,j=1,\dots ,s},\quad {\mathbf {b}}=[b_1,\dots ,b_s],\quad {\mathbf {c}}=[c_1,\dots , c_s] \end{aligned}$$

and

$$\begin{aligned} \alpha _{i,l}^j,\quad \beta _l^{i},\quad i=1,\dots s, \,\, j=1,\dots ,s,\,\, l=1,\dots , J,\qquad \hat{{\mathbf {c}}}=[\hat{c}_1,\dots ,\hat{c}_s]. \end{aligned}$$

When applied to the convection-diffusion problem

$$\begin{aligned} \dot{U}(t)+C(U(t))\,U(t)=AU(t), \end{aligned}$$

with linear diffusion and nonlinear convection the DIRK-CF methods have the following format:

figure e

The methods are associated to the two Butcher tableaus,

(23)

where we have defined

$$\begin{aligned} \hat{a}_{i,j}:\!\!{=}\sum _{l=1}^J\alpha _{i,l}^j, \qquad \hat{b}_{j}:\!\!{=}\sum _{l=1}^J\beta _{l}^j, \end{aligned}$$
(24)

for \(i=1,\dots , s\), \(\hat{{\mathcal {A}}}=\{\hat{a}_{i,j}\}_{i,j=1,\dots ,s}\) and \(\hat{{\mathbf {b}}}=[\hat{b}_1,\dots ,\hat{b}_s]\). The coefficients of the first tableau are used for the linear vector field \(Ay\) while the coefficients of the second tableau, split up in the sums (24), are used for the nonlinear vector field \(C(y)y\). We choose the first tableau to be a DIRK (diagonally implicit Runge–Kutta) method, this means we are solving only one linear system per stage. The tableaus (23) are typically chosen so that they define a classical IMEX method, which we call the underlying IMEX method, see [1] and [7] for more details. This IMEX method has the format

figure f

The order theory for classical IMEX methods reduces to the theory of partitioned Runge–Kutta methods, [25]. Given an implicit and an explicit method of order \(\kappa \) they must satisfy extra compatibility conditions in order for the corresponding IMEX method to have order \(\kappa \). The extension of this theory to the DIRK-CF methods has been discussed in [7] and [8].

1.2 Definition of Norms

For a square-integrable (respectively \(H_1\)) function \({\mathbf {u}} : \Omega \rightarrow {\mathbb {R}}^n,\) where \(\Omega \subset {\mathbf {R}}^m\) is bounded and connected, the \(L_2\)-norm (\(\Vert \cdot \Vert _{L_2(\Omega )}\)) and the \(H_1\)-norm (\(\Vert \cdot \Vert _{H_1(\Omega )}\)) are defined by

$$\begin{aligned} \Vert {\mathbf {u}}\Vert _{L_2(\Omega )}&:\!\!{=}\left( \sum _{i=1}^n\int _{\Omega }u_i^2\,d\Omega \right) ^{1/2},\end{aligned}$$
(25)
$$\begin{aligned} \Vert {\mathbf {u}}\Vert _{H_1(\Omega )}&:\!\!{=}\left( \sum _{i=1}^n\int _{\Omega }(u_i^2+\nabla {u}_i\cdot \nabla {u}_i)\,d\Omega \right) ^{1/2}. \end{aligned}$$
(26)

In the spectral element approximations the continuous integrals of numerical solutions are accurately computed using Gauss quadrature rules.

1.3 Boundary Conditions and Discrete Stiffness Summation

For the sake of completeness, we illustrate the strategy for implementing the boundary conditions in the context of spectral element methods. We use the spectral element notion known as the direct-stiffness summation (DSS), see for instance [13].

Suppose we have to impose periodic or homogeneous Dirichlet boundary conditions, and that the variable \(\bar{y}\) represents the values of the numerical solution at all discretization nodes in the computational domain (including boundary nodes). The variable \(y\) represents the restriction of \(\bar{y}\) to the minimum degrees of freedom \(k\) required to define the numerical solution, while \(\bar{y}\) contains typically redundant components. Thus if the number of components of \(\bar{y}\) is \({\mathcal {N}},\) then \(k<{\mathcal {N}}.\) We denote by \(Q\) a prolongation or “scatter” operator such that \(\bar{y} = Qy.\) Associated to \(Q\) is a restriction or “gather” operator denoted by \(Q^T.\) The operator \(Q\) is a \({\mathcal {N}}\times k\) constant matrix of rank \(k\le {\mathcal {N}}\). The variable \(\bar{y}\) is referred to as the local variable, while \(y\) is the global variable. The DSS operator \(QQ^T\) ensures inter-element continuity and the fulfillment of the appropriate boundary conditions. So, for example, if the boundary conditions are periodic, given a vector \(\bar{y}\) in the solution space or in the space of vector fields, \(QQ^T\bar{y}\) is periodic.

The spectral element discretization of the Navier–Stokes equations yields the discrete system

$$\begin{aligned} \bar{B}\dot{\bar{y}}&= \bar{A}\bar{y}+\bar{C}(\bar{y})\bar{y}-\bar{D}^Tz \end{aligned}$$
(27)
$$\begin{aligned} \bar{D}\bar{y}&= 0, \end{aligned}$$
(28)

where \(\bar{y}\) is assumed to be in the range of \(Q\) (i.e. \(\bar{y} = Qy\)). The relation between the local and global operators is \(B=Q^T\,\bar{B}\,Q\), \(A=Q^T\,\bar{A}\,Q\), \(C(y)=Q^T\, \bar{C}(Qy)\,Q\) and \(D=\bar{D}\,Q\).

Applying on both sides of (27) the DSS operator \(QQ^T\) we obtain

$$\begin{aligned} \Sigma \dot{\bar{y}}&= QQ^T \bar{A}\bar{y}+QQ^T\bar{C}(\bar{y})\bar{y}-QQ^T\bar{D}^Tz, \end{aligned}$$
(29)
$$\begin{aligned} \bar{D}\bar{y}&= 0. \end{aligned}$$
(30)

The matrix \(\Sigma =QQ^T\bar{B}\) is \({\mathcal {N}}\times {\mathcal {N}}\) and invertible on the range of \(Q\). In practice the integration methods are reformulated for the local variable \(\bar{y}\) and the local operators.

Indeed, in the computations the full data for the local variable \(\bar{y}\) is stored, since all computations involving the operators \(B^{-1}, {\mathcal {H}}^{-1}\) and \((DB^{-1}D^T)^{-1}\) must be done within the range of \(Q.\) These operators are symmetric and positive-definite, and so they can be inverted using a fast iterative solver such as the conjugate gradient method. For example the problem \(y = {\mathcal {H}}^{-1}f\) is reformulated as follows: Find \(\bar{y}\) such that

$$\begin{aligned} \left\{ \begin{aligned} QQ^T{\mathcal {\bar{H}}}\bar{y}&= Qf\\ Qy&= \bar{y} \end{aligned}\right. , \end{aligned}$$
(31)

where \({\mathcal {\bar{H}}} = \frac{1}{a_{i,i}h}\bar{B} - \bar{A}.\)

We refer to [18] for further details on DSS and boundary conditions. In the experiments reported in this paper, no special treatment has been taken to enforce pressure boundary conditions, since the discrete pressure space is not explicitly defined on discretization nodes on the boundary.

1.4 Reformulation of the Integration Methods in Local Variables

In this section we briefly discuss the correct implementation of the methods of Sect. 2 as applied to (29). The purpose of reporting here these implementation details is to explicitly highlight when care has to be taken in the implementation. See also the remark below.

We first perform the elimination of the discrete pressure Lagrangian multiplier from (29) in analogy to (9). We use

$$\begin{aligned} \bar{H}=\bar{D}^T\left( DB^{-1}D \right) ^{-1}DB^{-1}Q^T, \end{aligned}$$

and we get a system of ODEs for the variable \(\bar{y}\):

$$\begin{aligned} \Sigma \dot{\bar{y}}&= QQ^T\bar{A}\,\bar{y}+QQ^T \bar{C}(\bar{y})\,\bar{y} - QQ^T\bar{H}(\bar{A}\,\bar{y}+\bar{C}(\bar{y})\,\bar{y}) , \end{aligned}$$
(32)

Introducing the projection \(\bar{\Pi }=I-\bar{H}\) allows us to write the following projected system of ODEs

$$\begin{aligned} \Sigma \dot{\bar{y}}= QQ^T\bar{\Pi } \bar{A}\bar{y}+ QQ^T\bar{\Pi } \bar{C}(\bar{y})\bar{y}. \end{aligned}$$
(33)

Denote \(\tilde{\Pi } = Q\Pi B^{-1}Q^T.\) Applying the method of Sect. 2.2 to the ordinary differential equation (33) split in its projected convection and projected diffusion terms, and then rewriting it as a method for the differential-algebraic equation (29) we obtain

figure g

under the assumption \(\bar{y}_n = Qy_n.\) Since \(\bar{y}_{n+1}=\bar{Y}_s,\) the approximation of the velocity satisfies the discrete incompressibility constraint, \(\bar{D}\bar{y}_{n+1}=0\).

Remark 2

We observe that this is not equivalent to what we obtain applying directly the IMEX method to (29), which written in ODE form is

$$\begin{aligned} \Sigma \dot{\bar{y}}=QQ^T \bar{A}\bar{y}+QQ^T\bar{C}(\bar{y})\bar{y}-QQ^T\bar{H} \left( \bar{A}y+\bar{C}\left( \bar{y}\right) \bar{y}\right) . \end{aligned}$$

In fact if we apply the IMEX method to (29) we need to treat the term \(QQ^T\bar{D}^Tz\) either with the implicit method or with the explicit method, while in the approach outlined in this section we treated \(QQ^T\bar{H}\bar{C}(\bar{y})\bar{y}\) explicitly and \(QQ^T\bar{H}\bar{A}\bar{y}\) implicitly, see also [29].

Analogously, the method of Sect. 2.1 applied to for (33) becomes

figure h

The method of Sect. 6.5 may be reformulated in a similar way.

1.5 Projected DIRK-CF Methods for Navier–Stokes Equations

In this section we present an alternative semi-Lagrangian approach compared to the previous Sect. 2.1. Also in this case the integration methods are a variant of the exponential integrators methods of [7].

We consider (11), and rearrange the terms in the form

$$\begin{aligned} \dot{y}= B^{-1}A\,y-HB^{-1}(A+C(y))\,y+ B^{-1}C(y)\,y. \end{aligned}$$
(34)

This is a differential equation on the subspace of discrete divergence-free vector fields, i.e. \(Dy=0\) for all \(t\). To approximate the solution of this equation, we consider a projection method of the type reviewed in [25, IV.4], see also [14, Sect. 5.3.3] and [26, Sect. VII.2]. The idea is to use a one-step integrator \(\phi _h\) for advancing the numerical solution of (34) by one step, and an orthogonal projection on the subspace of divergence-free vector fields applying \(\Pi \) at the end of each step. We choose \(\phi _h\) to be the following integration method, in which the coefficients of both the DIRK-CF method and the underlying IMEX method are used:

  • the term \((I-H)B^{-1}A\,y\) is treated implicitly with the DIRK coefficients,

  • the term \(HB^{-1}C(y)\,y\) is treated explicitly with the coefficients of the underlying explicit method,

  • the term \(B^{-1}C(y)\,y\) is treated with the coefficients of the corresponding CF method.

The projection \(\Pi \) is used to guarantee divergence-free numerical approximations, i.e. \(Dy_{n+1}=0\). We obtain

figure i

where

$$\begin{aligned} a_{s+1,j} = b_j,\,\,\hat{a}_{s+1,j} = \hat{b}_j,\,\,\alpha _{s+1,\gamma }^k = \beta _{\gamma }^k,\,\,j,k = 1,\ldots ,s,\,\, a_{s+1,s+1} = 0,\,\,\hat{a}_{s+1,s+1} = 0. \end{aligned}$$

Since \(\Pi \) is an orthogonal projection, the order of the method is not affected by the use of the projection in the update step.Footnote 4 This approach requires only exponentials of pure convection problems, this will ease the implementation of the method as a semi-Lagrangian method, because now we simply use the semi-Lagrangian approximation

$$\begin{aligned} \exp \left( hB^{-1} C\big (w\big )\right) \,g= {\mathcal {I}}(g)(\Phi _h^w(\Gamma )). \end{aligned}$$

Observe that at each stage \(Y_i \) does not necessarily satisfy \(DY_i =0\).

Rights and permissions

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Celledoni, E., Kometa, B.K. & Verdier, O. High Order Semi-Lagrangian Methods for the Incompressible Navier–Stokes Equations. J Sci Comput 66, 91–115 (2016). https://doi.org/10.1007/s10915-015-0015-6

Download citation

  • Received:

  • Revised:

  • Accepted:

  • Published:

  • Issue Date:

  • DOI: https://doi.org/10.1007/s10915-015-0015-6

Keywords

Mathematics Subject Classification

Navigation