unsteady flow problems
Tomas Lundquist · Jan Nordstr¨om
January 22, 2016
Abstract We make an initial investigation into the temporal efficiency of a fully dis-crete summation-by-parts approach for unsteady flows. As a model problem for the Navier-Stokes equations we consider a two-dimensional advection-diffusion problem with a boundary layer. The problem is discretized in space using finite difference ap-proximations on summation-by-parts form together with weak boundary conditions, leading to optimal stability estimates. For the time integration part we consider vari-ous forms of high order summation-by-parts operators and compare with an existing popular fourth order diagonally implicit Runge-Kutta method. To solve the result-ing fully discrete equation system, we employ a multi-grid scheme with dual time stepping.
Keywords Summation-by-parts in time· unsteady flow calcluations · temporal efficiency
Mathematics Subject Classification (2000) 65M06· 65M55
1 Introduction
Based on finite difference operators on summation-by-parts (SBP) form and the simultaneous-approximation-term (SAT) technique for imposing boundary conditions, the SBP-SAT technique constitutes a very robust framework for implementing high
T. Lundquist
Department of Mathematics, Computational Mathematics, Link¨oping University, SE-581 83 Link¨oping, Sweden
Tel.: +46 13 281456
E-mail: tomas.lundquist@liu.se J. Nordstr¨om
Department of Mathematics, Computational Mathematics, Link¨oping University, SE-581 83 Link¨oping, Sweden
Tel.: +46 13 281459 E-mail: jan.nordstrom@liu.se
order finite difference schemes on complex geometries. By construction, it leads to discrete energy estimates that perfectly imitates the corresponding continuous energy estimates, and hence to stability. See [17] for an overview of the development of the SBP-SAT technique until present, with focus on a wide range of applications in fluid dynamics.
The SBP-SAT technique was recently extended to initial value problems [14, 12], making it possible to formulate fully discrete approximations with optimal energy es-timates. The purpose of this work is to make an initial efficiency study of the temporal part of these schemes for a stiff model problem.
The numerical treatment of unsteady flow problems has gained increased atten-tion in later years due to increased computer resources making realistic calculaatten-tions more viable. However, the construction of efficient algorithms still remains a signifi-cant computational challenge. The most commonly used methods are based on either Newton iteration or dual time stepping. While Newton iterations typically achieves better convergence rates for small error tolerances, dual time stepping is more re-liable, at least for the initial iterations, and it can also be used for preconditioning purposes [11]. Some studies indicate that a combination of both these techniques gives the best results [2, 11, 3].
The organization of this report is as follows. In section 2 we review the basic properties of the SBP-SAT technique for time integration. In section 3 we analyze a linear model of the Navier-Stokes equations with a boundary layer. In section 4 a multi-grid dual time stepping scheme is formulated to solve the model problem iteratively. Numerical results for the model problem is presented in section 5, and finally in section 6 we draw conclusions.
2 Summation-by-parts in time
We introduce the SBP-SAT technique for time integration by considering the test equation ut+λu= 0, with initial condition u(0) = f . Let t = (0,∆t, . . . , N∆t= T ) be a uniform grid vector with N+ 1 grid points, and let (ej)Nj=0denote the unit vectors in the standard basis. An SBP-SAT approximation of the test problem is given by:
P−1Qu+λu= P−1σ(u0− f )e0. (2.1) The SAT penalty treatment on the right hand side of (2.1) forces the initial solution towards the initial data, and the SBP properties of the first derivative operator P−1Q
are the following:
Q+ QT= eNeTN− e0eT0 = Diag(−1,0,...,0,1), P = PT> 0.
These properties lead in an automatic way to a clean, optimally sharp energy esti-mate. Especially, forσ= −1 we get after multiplying (2.1) with u∗P and adding the
conjugate transpose:
where the norm is defined using the quadrature matrix P as||u||2P= u∗Pu. This
mim-ics the continuous energy estimate |u(T )|2+ 2Re(λ)kuk2dt = | f |2, wherekuk2= RT
0 |u|2dt (obtained by multiplying the test equation with u∗and then integrating).
As an alternative to the global formulation (2.1), we may also consider a multi-stage version in which the size of the equation systems to be solved remain limited, even for long time calculations. This can have computational advantages, such as reduced memory usage. With this approach we instead consider a sequence of points in time t0,t1, . . . ,tn,tn+1, . . ., and solve the problem by discretizing each subinterval using the SBP-SAT technique:
(P−1Q+λI)un+1= P−1σ(un0+1− unr)e0, (2.2) where un+1= (un+1
0 , un1+1, . . . , unr+1) for n = 0, 1, . . ., and the operator P−1Q is de-fined on the uniform subgrid tn+1= (tn,tn+∆t
n/r, . . . ,tn+1), where∆tn= tn+1−
tn. Starting with u0r = f , the scheme (2.2) then produces the sequence of solutions
u0r, u1
r, . . . , unr, unr+1, . . . approximating the solution in t0,t1, . . . ,tn,tn+1, . . ..
The size of the matrix P−1Q in the multi-stage formulation (2.2) is given by the
number of intermediate stages r+ 1 used for each subinterval, and thus remains
con-stant also for long time calculations. Conversely, the minimum number of stages stages depends on how small P−1Q can be made. The classical SBP operators are
based on a repeated central finite difference stencil together with boundary closures. An example is the second order operator given by
P=△t r 1 2 1 . .. 1 1 2 , Q= −1 2 1 2 −1 2 0 1 2 . .. . .. . .. −1 2 0 1 2 −1 2 1 2 .
Higher order operators need far more extensive boundary closures than in the exam-ple above, increasing the minimum number of stages required. Other types of opera-tors on SBP form do however exists, e.g. based on spectral collocation [6], and may alternatively be used to further decrease the necessary number of stages, as demon-strated in [4]. See also [8] for a detailed account on the construction of non-classical SBP operators based on any numerical quadrature.
The SBP-SAT technique presented above can be applied in an analogous way to any initial value problem, and the formal properties of the schemes obtained are summarized below (all proofs can be found in [12]):
– The SBP-SAT schemes (2.2) are always A-stable and L-stable. If P is diagonal they are also B-stable and preserve energy stability. Moreover, they lead to opti-mally sharp fully discrete energy estimates.
– The order of convergence follows the order of accuracy of the quadrature matrix
P. For classical SBP operators this coincides with the order of the interior scheme.
– The order of stiff convergence is given by the smallest local order of consistency in P−1Q , the so-called stage order. Using classical operators, the stiff
Note that this means that the non-stiff order of convergence is not reduced by the lower order at the boundaries of classical operators. This is explained by the fact that the SBP-SAT discretizations are dual consistent approximations of the continuous problem [12].
For the efficiency studies described later in this report, we have limited the dis-cussion to two types of SBP operators. First the classical operators with diagonal norms, denoted by SBP(2s,s), and secondly the spectral element operators based on Gauss-Lobatto quadrature, denoted by GL(2s,s). In both cases, 2s is the order of the quadrature P, and thus of the scheme itself, while s is the stage order. The GL(2s,s) schemes have recently been shown to be equivalent to the Lobatto IIIC type of clas-sical Runge-Kutta methods, also based on the same type of quadrature [5].
3 A stiff flow model in two dimensions
As a model of the Navier-Stokes equation, we study a fluid undergoing advective flow past a plate with fixed value:
ut+ ux=ε(uxx+ uyy) +ψ 0≤ x,y ≤ 1 t ≥ 0 u(0, x, y) = f (x, y) t= 0 u(t, 0, y) −εux(t, 0, y) = g1(t, y) ∂Ω1= {(x,y) : x = 0} εux(t, 1, y) = g2(t, y) ∂Ω2= {(x,y) : x = 1} u(t, x, 0) = 0 ∂Ω3= {(x,y) : y = 0} uy(t, x, 1) = 0 ∂Ω4= {(x,y) : y = 1}, (3.1)
whereε = 0.01. The solid boundary∂Ω3 is associated with a boundary layer of
width√εat y= 0, the inflow and outflow boundaries are∂Ω1and∂Ω2respectively,
and∂Ω4is a far-field boundary. An exact manufactured solution can be imposed by
appropriately specifying the forcing functionψ. We define the continuous norm bykuk2=R
Ωu2dΩ, and apply the energy method by multiplying the equation (3.1) with the solution u for the caseψ= 0. Integrating
the result over the physical domain and applying the divergence theorem leads to
kukt2+ 2 Z ∂Ω( 1 2u 2 −εuux,εuuy) · ndS = −2ε(kuxk2+ kuyk2).
We decompose the boundary integral into four parts corresponding to the four bound-aries ofΩ: Z ∂Ω1 (1 2u 2−εuu x,εuuy) · n1dS= Z ∂Ω1 (−1 2u 2+εuu x)dS = Z ∂Ω1 1 2((u − g1) 2− g2 1)dS Z ∂Ω2 (1 2u 2−εuu x,εuuy) · n2dS= Z ∂Ω2 (1 2u 2−εuu x)dS = Z ∂Ω2 1 2((u − g2) 2− g2 2)dS Z ∂Ω3 (1 2u 2−εuu x,εuuy) · n3dS= Z ∂Ω3 εuuydS= 0 Z ∂Ω4 (1 2u 2 −εuux,εuuy) · n4dS= Z ∂Ω4 −εuuydS= 0.
Thus the estimate becomes kukt2+ 2ε(kuxk2+ kuyk2) = Z ∂Ω1 (g21− (u − g1)2)dS + Z ∂Ω2 (g2 2− (u − g2)2)dS. (3.2)
The energy estimate (3.2) shows that the problem (3.1) is well-posed, since the growth of the solution is bounded in terms of the boundary data.
In order to properly resolve the boundary layer at y= 0 in (3.1), we introduce a
stretching functionηof the vertical coordinate, given by
y= 1 +tanh(B(η− 1)) tanhB ,
where B= 9/4, leading toηy(0) = 1/√ε. This stretching function is shown in Figure 3.1. See also Appendix A for a detailed discussion of this choice. After this change
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 η y
Fig. 3.1 Stretching of vertical coordinate to resolve the boundary layer.
of coordinate, the model problem (3.1) becomes
ut+ ux=ε(uxx+ηy(ηyuη)η) +ψ 0≤ x,η≤ 1 t ≥ 0 u(0, x,η) = f (x,η) t= 0 u(t, 0,η) −εux(t, 0,η) = g1(t,η) ∂Ω1= {(x,η) : x = 0} εux(t, 1,η) = g2(t,η) ∂Ω2= {(x,η) : x = 1} u(t, x, 0) = 0 ∂Ω3= {(x,η) :η= 0} ηyuη(t, x, 1) = 0 ∂Ω4= {(x,η) :η= 1}. (3.3)
3.1 Semi-discrete approximation
We follow the techniques, presented in [13, 15, 16] for the linearized Navier-Stokes equations to discretize the model problem (3.3) in space. Let Dx= Px−1Qxand Dη=
Pη−1Qη be SBP operators approximating the spatial derivatives in both coordinate directions on a cartesian grid. Organizing the solution vector U in lexicographical order, i.e. we define U= (UT
0, . . . ,UNTy)
T, where U
i= (U0,i, . . . ,UNx,i), all numerical
derivatives can be constructed using these one-dimensional operators through appli-cation of the Kronecker product. We thus define Ux= (Dx⊗ Iη)U, Uη= (Ix⊗ Dη)U and Uy= (Ix⊗ HyDη)U, where Hy= Diag(ηy) represents the stretching of the ver-tical coordinate. We also let(exi)Ni=0x and(eηj)
Nη
j=0 denote the standard unit vectors in the spaces corresponding to x andη, and use these vectors to restrict the solution to fixed coordinates in x andη. For example we define U|x=xi = (e
T
xi⊗ Iη)U and
U|η=ηj = (Ix⊗ e
T
ηj)U, restricting the solution to all points where x = xiandη=ηj respectively. We will use analogous definitions also for the derivatives.
In the analysis below, the vertical quadrature matrix Pηis assumed to be diagonal. This is necessary in order to mimic the estimate (3.2) in the discrete setting on the stretched coordinate grid [13]. An SBP-SAT approximation to (3.3) is obtained by replacing the derivatives with numerical ones, as well as imposing all four boundary conditions through weak penalty terms in the right-hand-side. We get
Ut+ (Dx⊗ Iη)U =ε((D2x⊗ Iη)U + (Ix⊗ (HyDη)2)U) + (Px−1⊗ Pη−1Hy)(Σx(t) +Ση(t)) +Ψ(t) U(0) = F, (3.4) where Ψ(t) =ψ(t, (x ⊗ 1η), (1x⊗η)) Σx(t) =σ0x(ex0⊗ Hy−1Pη(U|x=0−εUx|x=0− g1(t))) +σ1x(exNx⊗ Hy−1Pη(εUx|x=1− g2(t))) Ση(t) =σ0η(PxU|η=0⊗ eη0) +σ1η(PxUη|η=1⊗ eηNη) g1(t) = g1(t, (ex0⊗η)) g2(t) = g2(t, (exNx⊗η)) F= f ((x ⊗ 1η), (1x⊗η)).
The constant vectors 1xand 1ηcontain the value 1 on each position, and are of length
Nx+ 1 and Nη+ 1 respectively
Note that the matrix premultiplying the weak penalty terms in the right-hand-side of (3.4) defines a discrete normkUk2
Px⊗(Hy−1Pη)= U
T(P
x⊗ Hy−1Pη)U, approximat-ing the correspondapproximat-ing continuous normkuk2=R
Ω|u|2. We apply the discrete energy method by multiplying (3.4) with UT(Px⊗ Hy−1Pη), and then adding the conjugate transpose. A stable scheme is obtained with the following choice of penalty parame-ters:σ0x=σ1x= −1,σ1η= −1/2 andσ0η= −εηy(0)2/Pη00. The motivation for this
16] for details on the imposition of the various types of boundary conditions present here. Ignoring the forcing data, i.e. lettingΨ(t) = 0, the estimate becomes
(kUk2 Px⊗(Hy−1Pη))t+ 2ε(kUxk 2 Px⊗(Hy−1Pη)+ kUyk 2 Px⊗(Hy−1Pη)) + 2ε U|η=0 Uy|η=0 T (F ⊗ Px) U|U η=0 y|η=0 = kg1(t)k2 Hy−1Pη− kg(t)1−U|x=0k 2 Hy−1Pη + kg2(t)k2 Hy−1Pη− kg(t)2−U|x=1k 2 Hy−1Pη, (3.5) where Pη= Pη− Pη00e0ηeT0η, F= ηy(0) Pη00 1 2 1 2 Pη00 ηy(0) .
This estimate closely mimics the corresponding continuous estimate (3.2). Note that the value Pη00is “borrowed” from Pη into F in order to make the latter a positive
definite matrix. This method for imposing wall boundary conditions was used in [7, 16, 1].
3.2 Fully discrete approximation with SBP-SAT in time
Consider the semi-discrete energy estimate (3.5) again. Integrating over a finite time interval t∈ [a,b] yields
kU(b)k2 Px⊗(Hy−1Pη)+ 2ε Z b a (kUxk 2 Px⊗(Hy−1Pη)+ kUyk 2 Px⊗(Hy−1Pη))dt + 2ε Z b a U|η=0 Uy|η=0 T (F ⊗ Px) U|U η=0 y|η=0 dt = Z b a (kg1(t)k 2 Hy−1Pη− kg1(t) −U|x=1k 2 Hy−1Pη + kg2(t)k2 Hy−1Pη− kg2(t) −U|x=0k 2 Hy−1Pη)dt + kU(a)k2 Px⊗(Hy−1Pη). (3.6)
This form of the estimate is what we need to preserve in the fully discrete formulation. Before continuing, we write the semi-discrete problem (3.4) on a more compact form as the following linear system:
Ut+ BU = R(t) 0 < t ≤ T
U(0) = F,
where the matrix and the right-hand-side vector are given by
B= Px−1(Qx−σ0xex0eTx0−ε(Qx+σ0xex0eTx0−σ1xexNxe
T
xNx)Dx) ⊗ Iη
+ Ix⊗ Pη−1(−σ0ηeη0eTη0−ε(HyQηHy+σ1ηeηNηeTηNηHy)Dη)
−10 −8 −6 −4 −2 0 x 104 −1 −0.5 0 0.5 1x 10 4 Re(λ) Im( λ ) 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x y
Fig. 3.2 Left: spectrum of the semi-discrete spectrum. Right: The computational grid with stretched
ver-tical coordinate.
The semi-dicrete spectrum of a fifth order discretization (i.e. SBP(8,4) in space) with
Nx= Nη= 95 is shown in Figure 3.2. The spectral radius of almost 105indicates that an explicit time marching scheme would be an inefficient way to solve this problem. Instead we employ the SBP-SAT technique also for the time discretization. This au-tomatically leads to stability since the semi-discrete problem has an energy estimate (see Proposition 9 in [12]).
An implicit SBP-SAT time integration scheme with r+ 1 stages for (3.2) is given
by
(P−1Q⊗ IB)Un+1+ (It⊗ B)Un+1= (P−1σe0) ⊗ (U0n+1−Urn) + Rn+1 (3.7) where
Un+1= (U0n+1,U1n+1, . . . ,Urn+1)
Rn+1= (R(tn), R(tn+∆t/r), . . ., R(tn+∆t= tn+1)),
and as before we letσ= −1. The fully discrete norm is now given by kUn+1k2
P⊗Px⊗(Hy−1Pη)=
(Un+1)T(P⊗P
x⊗Hy−1Pη)Un+1, and we apply the energy method by multiplying (3.7) with(Un+1)T(P ⊗ P
x⊗ Hy−1Pη) and then adding the conjugate transpose. In this way we automatically preserve the discrete energy estimate, only adding a small damping term from the initial condition (see [14] for more examples and details). We get
kUn+1 r k2P x⊗(Hy−1Pη)+ 2ε(kU n+1 x k2P⊗P x⊗(Hy−1Pη)+ kU n+1 y k2P⊗P x⊗(Hy−1Pη)) + 2ε U n+1|η =0 Uny+1|η=0 T (P ⊗ F ⊗ Px) Un+1|η =0 Uny+1|η=0 = kGn+1 1 k2P⊗(Hy−1Pη)− kG n+1 1 − Un+1|x=1k2P⊗(H−1 y Pη) + kGn+1 2 k2P⊗(Hy−1Pη)− kG n+1 2 − Un+1|x=1k2P⊗(H−1 y Pη) + kUrnk2Px⊗(H−1 y Pη)− kU n+1 0 −U n rk2Px⊗(H−1 y Pη), (3.8)
where
Gn1+1= (g1(tn), g1(tn+∆t/r), . . ., g1(tn+1)) Gn2+1= (g2(tn), g2(tn+∆t/r), . . ., g2(tn+1)).
The estimate (3.8) perfectly imitates the semi-discrete estimate (3.6), only adding a small additional damping term. It of course also mimics the continuous estimate (3.2).
4 The dual time stepping scheme
Consider the compact formulation of (3.7): ˜ BUn+1= ˜R (4.1) where ˜ B= P−1(Q −σe0eT0) ⊗ IB+ It⊗ B ˜ R= −(P−1σe0) ⊗Un+ Rn+1.
We employ a multi-grid cycle for solving (4.1), where the smoothing step consists of a pseudo time stepping scheme. Thus, we add a pseudo time derivative to (4.1):
dUn+1 dτ + ˜BU
n+1= ˜R. (4.2)
The solution to (4.1) is now given by the steady-state solution to (4.2). To march forward in pseudo time, we use an explicit s-stage low storage Runge-Kutta smoother:
Wn0+1,m+1= Vn+1,m
Wnp+1,m+1= Vn+1,m+∆ταp( ˜R− ˜BWnp+1,m+1−1 ), p = 1, . . . , s Vn+1,m+1= Wns+1,m+1
(4.3)
The stability function of this scheme is given by S(z) = (1 +αsz(1 +αs−1z(...(1 +
α1z)...))), where z =λ∆τ. In order to match the specific shape of the semi-discrete
spectrum shown in Figure 3.2, we use the 4-stage smoother given byα= (0.0178571,
0.0568106 , 0.174513 , 1). It was proposed for viscous problems with similar
spec-trums in [10]. The stability region of this scheme is shown, and compared with the semi-discrete spectrum, in Figure 4. The reason that we focus on the semi-discrete spectrum is that adding the temporal part of the discretization results in very small perturbations. The stiff components of the semi-discrete problem remain as the dom-inating part of the fully discrete spectrum, making the stability condition independent of the physical time stepping scheme.
A time step restriction on the pseudo time stepping scheme is obtained by com-paring the largest (in magnitude) negative eigenvalue in the spectrum with the left boundary of the stability region in Figure 4. The same method is used on each of the coarse grids introduced in the multi-grid scheme. We have used the condition
Re(z) = −25 to determine the pseudo time step size, which leaves some safety
−10 −8 −6 −4 −2 0 x 104 −1 −0.5 0 0.5 1x 10 4 Re(λ) Im( λ ) Re(z) Im(z) −30 −20 −10 0 −15 −10 −5 0 5 10 15 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Fig. 4.1 Left: the semi-discrete spectrum of (3.4). Right: Stability region of Runge-Kutta smoother.
The full multi-grid scheme works as follows. Let Vn+1 be an approximation to the fully discrete solution Un+1. Marchingν steps in pseudo time with the ex-plicit scheme (4.3), starting from Vn+1, is defined by the update function Vn+1 =
Sν(Vn+1, ˜B, ˜R). A multi-grid cycle for (4.1) with L levels of grid coarsening can now be written as the fixed point iteration scheme Vn+1= MG(Vn+1, ˜R, L), where the multi-grid step function is defined as the recursive function x= MG(x, b, l), given
by:
– x= Sν(x, ˜Bl, b) (presmoothing) – if l> 0, do:
– r= Restrict(v − ˜Blx) (restriction)
– e= MG(0, r, l − 1) (coarse grid correction) – e= Prolong(e) (prolongation)
– x= x + e (fine grid correction) – x= Sν(x, ˜Bl, b). (postsmoothing)
5 Numerical results
We employ the manufactured solution u= sin(2π(x − t))e
1√−y
ε to (3.1), and compare
the numerical results for a selection of high order temporal schemes. The SBP opera-tors used are the following: SBP(4,2) with 8 implicit stages, SBP(8,4) with 16 implicit stages, GL(4,2) with 3 implicit stages, and GL(8,4) with 5 implicit stages. Remember that the two numbers in parenthesis indicate the order of convergence and the stage order, respectively. For the classical operators we use the minimum number of stages possible, i.e. only the boundary closures are included. As an additional comparison, we consider a fourth order diagonally implicit Runge-Kutta scheme ESDIRK4, with a stage order of 2, and 5 implicit stages. Schemes of this type is a popular choice for
implicit time integration of unsteady flow problems, and the tableau we have used can be found in [9]). 101 102 10−4 10−3 10−2 10−1 implicit stages ||e|| ∞ ESDIRK4 SBP(4,2) GL(4,2) SBP(8,4) GL(8,4) p=5 p=4
Fig. 5.1 Convergence in max norm at t= 1 of the high order temporal schemes, as a function of temporal
resolution. 101 102 101 102 103 implicit stages MG cycles/stage ESDIRK4 GL(4,2) SBP(4,2) GL(8,4) SBP(8,4)
Fig. 5.2 The amount of work in terms of multi-grid cycles required to converge each implicit stage in Fig
5.1, as a function of temporal resolution.
In order to minimize the spatial error component we use the fifth order discretiza-tion with Nx= Nη= 95 with spectrum shown in Figure 3.2. Note that the spectral
102 103 104 10−4 10−3 10−2 10−1 MG cycles ||e|| ∞ ESDIRK4 GL(4,2) SBP(4,2) GL(8,4) SBP(8,4)
Fig. 5.3 Convergence as a function of total work, defined by sum of multi-grid cycles over all implicit
stages.
radius of almost 105is a result both of the boundary layer and of this deliberate
over-resolution in space. The time integration is carried out using a multi-grid V-cycle on three grid levels, with refinement in the vertical coordinate only. On each grid, 10 steps of the explicit Runge-Kutta scheme is used as smoother, and the pseudo time step size is selected on each grid according to Re(z) = −25. The number of pseudo time iterations is set to make the iteration error small compared with the error from the physical time discretization.
In Figure 5.1 we measure the accuracy of the temporal schemes as a function of the total number of implicit stages in time. The accuracy is measured at t= 1, and in
all cases we observe some level of order reduction, with convergence rates slightly less than the order of the scheme (but higher than the stage order).
The amount of work that was required to converge each implicit stage on average is shown in Figure 5.2, again as a function of the total number of implicit stages. This can be interpreted as the efficiency of the multi-grid scheme as a function of the level of grid refinement. We see that the efficiency is lower for the eighth order schemes compared with the fourth order ones. This result is most likely explained by the fact that the lower order schemes need finer discretizations in physical time for a given error tolerance, leading to faster convergence in pseudo time. We also see that the multi-grid efficiency is worse for the classical SBP operators that use more stages in each implicit solve.
In Figure 5.3, we combine the results of the two previous figures to show the efficiency of the temporal schemes. Work is here defined as the sum of multi-grid cycles over all implicit stages required for each level of temporal accuracy. The re-duced efficiency of the classical SBP operators is apparent also here, but current work indicate that this might be possible to correct by modifying the multi-grid scheme in an appropriate way. We also note that there is no advantage of the diagonally
im-plicit ESDIRK4 method over the Gauss-Lobatto SBP schemes. Finally, the increased accuracy of the eighth order schemes is counterbalanced by the reduced multi-grid efficiency, resulting in very similar results as for the fourth order schemes.
6 Conclusions and further work
We have investigated the temporal efficiency of a fully discrete SBP-SAT approach to unsteady flow calculations. As a model problem we have used the advection-diffusion equation in two dimensions with a stiff boundary layer. A basic multi-grid scheme with pseudo time stepping was applied to solve the fully discrete equations, with no attempts at optimizing the smoother for the different schemes. The results where com-pared for a selection of higher order temporal schemes, including both classical and non-classical SBP-SAT schemes, as well as the popular and widely used diagonally implicit Runge-Kutta method ESDIRK4.
The numerical results indicate that some of the SBP-SAT time stepping schemes can compete with ESDIRK4 for efficiency. However, the classical SBP schemes using a larger number of stages did not perform well in this basic setting. Current work indicate that this disadvantage can be overcome with a more suitable choice of multi-grid scheme that works more efficiently for fully implicit methods.
Future work will aim at developing more efficient multi-grid schemes for fully discrete SBP-SAT schemes, making the convergence less dependent on the number of intermediate stages. This will include a modification of the dual time stepping technique used as smoother. Similar investigations for more realistic unsteady flow problems with boundary layers will also be considered, and possibly also combining the dual time stepping technique with Newton iteration for faster convergence.
Appendix A Boundary layer and coordinate stretching
The solid wall boundary at y= 0 in (3.1) causes the solution to change rapidly wihin
a thin layer close to this boundary. To study this, we consider a simplified steady version of (3.1):
ux=ε(uxx+ uyy) 0 < x, y <∞
u(x, 0) = 0 u(0, y) = 0.
Note that this problem has a discontinuity at x= y = 0. In order to study what happens
close to the solid boundary, we make the coordinate transformation y=√εY , which
gives
ux=εuxx+ uYY.
We now expand u with a regular perturbation in powers ofε as u= u0+√εu1+
εu2+ . . .. The equation then becomes u0x+
√
εu1x+εu2x+ . . . =εu0xx+ u0YY+ √
Collecting all terms of the same order and setting each to zero now yields the se-quence of equations u0x= u0YY u1x= u1YY u2x= u0xx+ u2YY . . .
In particular, the full equation for the leading zeroth order term in the expansion is
u0x= u0YY u0(0,Y ) = 1
u0(x, 0) = 0.
This problem has the following analytical solution:
u0=εr f( Y √ 4x) =εr f( y √ 4εx)
whereεr f is the so-called error function, defined by:
εr f(ξ) =√2π Z ξ
0
e−s2ds,
The error function becomes almost constant for values ofξ> 2, making u0 approxi-mately constant for y> 4qεx. For a fixed value of x, the boundary layer width is thus proportional to√ε, and the derivative of u0with respect to y is moreover proportional
to 1/√ε.
In order to resolve the boundary layer at y= 0 in, we can introduce a stretching
ηof the vertical coordinate such thatηy(0) is proportional 1/√ε. A possible choice with this property is
y= 1 +tanh(B(η− 1)) tanhB ,
where B= 9/4, leading toηy(0) = 1/√ε. The full stretching function is shown in Figure 3.1.
References
1. Berg, J., Nordstr¨om, J.: Stable Robin solid wall boundary conditions for the Navier-Stokes equations. Journal of Computational Physics 230(19), 7519–7532 (2011)
2. Bijl, H., Carpenter, M.: Iterative solution techniques for unsteady flow computations using higher order time integration schemes. International Journal for Numerical Methods in Fluids 47, 857–862 (2005)
3. Birken, P., Jameson, A.: On nonlinear preconditioners in newtonkrylov methods for unsteady flows. International Journal for Numerical Methods in Fluids 62, 565–573 (2010)
4. Boom, P., Zingg, D.: High-order implicit time-marching methods based on generalized summation-by-parts operators. Tech. rep., arXiv:1410.0201 (2014)
5. Boom, P., Zingg, D.: Runge-kutta characterization of the generalized summation-by-parts approach in time. Tech. rep., arXiv:1410.0202 (2014)
6. Carpenter, M., Gottlieb, D.: Spectral methods on arbitrary grids. Journal of Computational Physics
129, 74–86 (1996)
7. Carpenter, M., Nordstr¨om, J., Gottlieb, D.: A stable and conservative interface treatment of arbitrary spatial accuracy. Journal of Computational Physics 148, 341–365 (1999)
8. Del Rey Fernandez, D., Boom, P., Zingg, D.: A generalized framework for nodal first derivative summation-by-parts operators. Journal of Computational Physics 266, 214–239 (2014)
9. Kennedy, C.A., Carpenter, M.H.: Additive Runge-Kutta schemes for convection-diffusion-reaction equations. Applied Numerical Mathematics 44(1-2), 139–181 (2003)
10. Kleb, W., Wood, W., van Leer, B.: Efficient multi-stage time marching for viscous flows via local preconditioning. AIAA Paper 99-3267 (1999)
11. Knoll, D., Keyes, D.: Jacobian-free newton-krylov methods: A survey of approaches and applications. Journal of Computational Physics 193, 357–397 (2004)
12. Lundquist, T., Nordstr¨om, J.: The SBP-SAT technique for initial value problems. Journal of Compu-tational Physics 270, 86–104 (2014)
13. Nordstr¨om, J., Carpenter, M.H.: High-order finite difference methods, multidimensional linear prob-lems and curvilinear coordinates. Journal of Computational Physics 173, 149–174 (2001)
14. Nordstr¨om, J., Lundquist, T.: Summation-by-parts in time. Journal of Computational Physics 251, 487–499 (2013)
15. Sv¨ard, M., Carpenter, M., Nordstr¨om, J.: A stable high-order finite difference scheme for the com-pressible Navier-Stokes equations: far-field boundary conditions. Journal of Computational Physics
225(1), 1020–1038 (2007)
16. Sv¨ard, M., Nordstr¨om, J.: A stable high-order finite difference scheme for the compressible Navier-Stokes equations: No-slip wall boundary conditions. Journal of Computational Physics 227(10), 4805–4824 (2008)
17. Sv¨ard, M., Nordstr¨om, J.: Review of summation-by-parts schemes for initial-boundary-value prob-lems. Journal of Computational Physics 268, 17–38 (2014)
