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.
Similar content being viewed by others
Notes
We also assume that the characteristic paths are integrated to high accuracy.
Extension to higher order methods is straightforward (See [30] and references therein).
See 6.2 for the definitions of these norms.
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
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)
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)
Brown, D.L., Minion, M.L.: Performance of under-resolved two-dimensional incompressible flow simulations. J. Comput. Phys. 122(1), 165–183 (1995)
Bruneau, C.-H., Saad, M.: The 2d lid-driven cavity problem revisited. Comput. Fluids 35(3), 326–348 (2006)
Canuto, C., Hussaini, M.Y., Quarteroni, A., Zang, T.A.: Spectral Methods in Fluid Dynamics. Springer Series in Computational Physics. Springer, New York (1988)
Celledoni, E.: Eulerian and semi-Lagrangian commutator-free exponential integrators. CRM Proc. Lect. Notes 39, 77–90 (2005)
Celledoni, E., Kometa, B.K.: Semi-Lagrangian Runge–Kutta exponential integrators for convection dominated problems. J. Sci. Comput. 41(1), 139–164 (2009)
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)
Celledoni, E., Kometa, B.K.: Semi-Lagrangian multistep exponential integrators for index 2 differential-algebraic systems. J. Comput. Phys. 230(9), 3413–3429 (2011)
Childs, P.N., Morton, K.W.: Characteristic Galerkin methods for scalar conservation laws in one dimension. SIAM J. Numer. Anal. 27, 553–594 (1990)
Chorin, A.J.: Numerical solution of the Navier–Stokes equations. Math. Comput. 22, 745–762 (1968)
Chorin, A.J.: On the convergence of discrete approximations to the Navier–Stokes equations. Math. Comput. 23, 341–353 (1969)
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)
Eich-Soellner, E., Führer, C.: Numerical Methods in Multibody Dynamics, European Consortium for Mathematics in Industry. B. G. Teubner, Stuttgart (1998)
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)
Falcone, M., Ferretti, R.: Semi-Lagrangian approximation schemes for linear and Hamilton-Jacobi equations. SIAM (to appear)
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)
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)
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)
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)
Giraldo, F.X.: The Lagrange–Galerkin spectral element method on unstructured quadrilateral grids. J. Comput. Phys. 147(1), 114–146 (1998)
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)
Guermond, J.L., Shen, J.: A new class of truly consistent splitting schemes for incompressible flows. J. Comput. Phys. 192(1), 262–276 (2003)
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
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)
Hairer, E., Wanner, G.: Solving Ordinary Differential Equations. II, Second ed., Springer Series in Computational Mathematics. Springer, Berlin (1996)
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)
Kennedy, C.A., Carpenter, M.H.: Additive Runge–Kutta schemes for convection–diffusion–reaction equations. Appl. Numer. Math. 44(1–2), 139–181 (2003)
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)
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)
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)
Morton, K.W.: Generalised Galerkin methods for hyperbolic problems. Computer Methods Appl Mech Eng 52, 847–871 (1985)
Rannacher, R.: On Chorin’s projection method for the incompressible Navier–Stokes equations, Lecture Notes in Mathematics, vol. 1530. Springer, Berlin (1991)
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)
Shen, J.: On error estimates of projection methods for Navier–Stokes equations: first-order schemes. SIAM J. Numer. Anal. 29(1), 57–77 (1992)
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)
Xiu, D., Karniadakis, G.E.: A semi-Lagrangian high-order method for Navier–Stokes equations. J. Comput. Phys. 172(2), 658–684 (2001)
Acknowledgments
This work was supported in part by the GeNuIn project, Grant from the Research Council of Norway.
Author information
Authors and Affiliations
Corresponding author
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:
and
When applied to the convection-diffusion problem
with linear diffusion and nonlinear convection the DIRK-CF methods have the following format:
The methods are associated to the two Butcher tableaus,
where we have defined
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
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
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
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
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
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
and we get a system of ODEs for the variable \(\bar{y}\):
Introducing the projection \(\bar{\Pi }=I-\bar{H}\) allows us to write the following projected system of ODEs
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
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
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
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
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
where
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
Observe that at each stage \(Y_i \) does not necessarily satisfy \(DY_i =0\).
Rights and permissions
About this article
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
Received:
Revised:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1007/s10915-015-0015-6