DRAFT
Link¨
oping Studies in Science and Technology.
Dissertations, No. 1886
Error analysis of summation-by-parts
formulations: Dispersion, transmission
and accuracy
Viktor Linders
Department of Mathematics, Division of Computational Mathematics Link¨oping University, SE-581 83 Link¨oping, Sweden
DRAFT
Link¨oping Studies in Science and Technology. Dissertations, No. 1886Error analysis of summation-by-parts formulations: Dispersion, trans-mission and accuracy
Copyright c Viktor Linders, 2017
Division of Computational Mathematics Department of Mathematics
Link¨oping University
SE-581 83, Link¨oping, Sweden
viktor.linders@liu.se www.liu.se/mai/ms
Typeset by the author in LATEX2e documentation system.
ISSN 0345-7524
ISBN 978-91-7685-427-3
DRAFT
DRAFT
Abstract
In this thesis we consider errors arising from finite difference operators on summation-by-parts (SBP) form, used in the discretisation of partial differential equations. The SBP operators are augmented with simultaneous-approximation-terms (SATs) to weakly impose boundary conditions. The SBP-SAT framework combines high order of accuracy with a systematic construction of provably sta-ble boundary procedures, which renders it suitasta-ble for a wide range of prosta-blems. The first part of the thesis treats wave propagation problems discretised using SBP operators on coarse grids. Unless special care is taken, inaccurate ap-proximations of the underlying dispersion relation materialises in the form of an incorrect propagation speed. We present a procedure for constructing SBP operators with minimal dispersion error. Experiments indicate that they out-perform higher order non-optimal SBP operators for flow problems involving high frequencies and long simulation times.
In the second part of the thesis, the formal order of accuracy of SBP operators near boundaries is analysed. We prove that the order in the interior of a di-agonal norm based SBP operator must be at least twice that of the boundary stencil, irrespective of the grid point distribution near the boundary. This gen-eralises the classical theory posed on uniform and conforming grids. We further show that for a common class of SBP operators, the diagonal norm defines a quadrature rule of the same order as the interior stencil. Again, this result is independent of the grid.
In the final contribution if the thesis, we introduce the notion of a transmis-sion problem to describe a general class of problems where different dynamics are coupled in time. Well-posedness and stability analyses are performed for continuous and discrete problems. A general condition is obtained that is nec-essary and sufficient for the transmission problem to satisfy an energy estimate. The theory provides insights into the coupling of fluid flow models, multi-block formulations, numerical filters, interpolation and multi-grid implementations.
DRAFT
Sammanfattning p˚
a svenska
Partiella differentialekvationer ligger till grund f¨or matematiska modeller av m˚anga ¨amnesspecifika problem inom t.ex. fysik, kemi, ekonomi och meteorologi. Vanligen angrips s˚adana problem numeriskt genom att ers¨atta de styrande ek-vationerna med en l¨amplig diskret approximation, vilken i sin tur kan l¨osas med hj¨alp av en dator. Det ¨ar emellertid oundvikligt att den numeriska proceduren introducerar fel i l¨osningen. I allm¨anhet ¨ar det utmanande att analysera dessa fel, samt att b˚ade begr¨ansa och minimera dem. I denna avhandling beaktas just s˚adana utmaningar.
En method som begr¨ansar tillv¨axten av numeriska fel s¨ags vara stabil. F¨or finita differensapproximationer p˚a partiell summationsform (eng. ”summation-by-parts”, SBP) kan stabilitet systematiskt visas om rand- och begynnelsevil-lkor implementeras p˚a svag form. Att minimera felen ¨ar emellertid ett separat problem och beror p˚a de styrande ekvationernas natur.
Avhandlingen inleds med en analys av l˚aguppl¨osta SBP-approximationer av ek-vationer vars l¨osningar antar en v˚agform. F¨or h¨ogfrekventa v˚agor uppst˚ar s˚a kallade dispersionsfel, vilka tar sig uttryck i felaktiga v˚aghastigheter. I avhan-dlingen utvecklas en metod f¨or att konstruera nya SBP-approximationer med minskat dispersionsfel s˚a att detta fenomen minimeras.
I avhandlingens andra del studeras SBP-approximationers noggrannhetsord-ning, dvs. den takt med vilken det numeriska felet minskar d˚a approximationens uppl¨osning f¨orfinas. H¨ar presenteras ett bevis f¨or att noggrannhetsordningen vid ber¨akningsdom¨anens gr¨anser som mest ¨ar h¨alften av ordningen i dom¨anens inre. Detta resultat g¨aller ¨aven om approximationen ¨ar finare uppl¨ost vid r¨anderna ¨
an i dom¨anens ¨ovriga delar.
Slutligen behandlar avhandligen s˚a kallade transmissionsproblem, vilka beskriver scenarier d˚a olika dynamiska system sammankopplas i tiden. B˚ade differen-tialekvationer och diskreta approximationer analyseras och ett allm¨ant sta-bilitetsvillkor h¨arleds. Transmissionsproblem ¨ar vanligt f¨orekommande, b˚ade som ekvationer och i samband med deras diskretiseringar. Den h¨arledda teorin f¨or dessa problem ger s˚aledes insikter i fysikaliska fenomen s˚av¨al som en bred arsenal numeriska tekniker, exempelvis filtrering och interpolation.
DRAFT
List of Papers
This thesis is based on the following papers, which henceforth will be referred to by their roman numerals:
I. Linders, V. & Nordstr¨om, J. (2015). Uniformly best wavenumber approxi-mations by spatial central difference operators. Journal of Computational
Physics, 300, 695-709.
II. Delorme, Y. T., Puri, K., Nordstr¨om, J., Linders, V., Dong, S. & Frankel, S. H. (2017). A simple and efficient incompressible Navier-Stokes solver for unsteady complex geometry flows on truncated domains. Computers
& Fluids, 150, 84-94.
III. Linders, V., Kupiainen, M. & Nordstr¨om, J. (2017). Summation-by-Parts operators with minimal dispersion error for coarse grid flow calculations.
Journal of Computational Physics, 340, 160-176.
IV. Linders, V., Lundquist, T. & Nordstr¨om, J. (2017). On the order of accuracy of finite difference operators on diagonal norm based
summation-by-parts form. Link¨oping University Press, LiTH-MAT-R–2017/11–SE
(Submitted).
V. Nordstr¨om, J. & Linders, V. (2017). Well-posed and stable transmission problems. Link¨oping University Press, LiTH-MAT-R–2017/15–SE
(Sub-mitted).
With the exception of paper II, I have contributed to the above listed papers by deriving most of the novel theoretical results as well as performing the numerical calculations. I have written the manuscripts myself with editorial support from
Jan Nordstr¨om, Marco Kupiainen and Tomas Lundquist. My contribution to
paper II includes deriving and providing the difference stencils used to obtain the numerical results as well as selecting appropriate numerical filters to match the stencils.
DRAFT
Acknowledgements
It goes without saying that the content of this thesis is the result of the efforts and engagements of many people. The list is longer than I have made space for, though I suspect that they know who they are.
First, I would like to express my gratitude to my advisor, Jan Nordtr¨om, without whom this venture would not have commenced, continued nor concluded. His determination, perseverance and patience has been key to any success I have enjoyed during the past few years. I must also thank my second advisor, Marco Kupiainen, for much the same reasons. He has taught me everything I know about fast computers, small planes, huge ships and cramped caves.
I am sincerely grateful to the many friends and colleagues that I have been fortunate enough to get to know; some through work, others through PhD student commitments, but most through the efforts you have made to turn Link¨oping into the wonderful place it is. This acknowledgement is incomplete until I have raised my glass to chant for the K¨oren!
Among those who have made especially important contributions to this thesis, I must mention Tomas Lundquist. His devotion to understanding is unprece-dented and has been a great inspiration.
I must also give my thanks to Cristina La Cognata. Simply put, she made me better; as a researcher and as a person.
My family deserves all the gratitude I can muster, but I suppose they already know that. This thesis is dedicated to them.
Last, but certainly not least, to my love and best friend, Anna. In some sense you are the antipole of this thesis; void of errors, yet of unbounded energy.
Viktor Linders Link¨oping, 2017
DRAFT
Contents
Abstract iii Sammanfattning p˚a svenska v List of Papers vii Acknowledgements ix 1 Introduction 1 2 Well-posedness and stability 3
2.1 Well-posed problems . . . 3
2.2 Stable discretisations . . . 4
3 The energy method 7 3.1 Strong boundary conditions . . . 7
3.2 Weak formulations . . . 8 3.3 Transmission . . . 8 4 Summation-by-parts 11 4.1 Semi-discretisations . . . 11 4.2 SBP in time . . . 12 4.3 Second derivatives . . . 13 4.4 Generalised SBP . . . 14
5 Finite difference operators 15 5.1 The structure of P and Q . . . . 15
5.2 Accuracy and convergence . . . 16
6 Dispersion 19 6.1 Analytic and numeric dispersion relations . . . 19
6.2 Filters . . . 20
DRAFT
xii CONTENTS
References 26 I. Uniformly best wavenumber approximations by spatial central
differ-ence operators . . . 29
II. A simple and efficient incompressible Navier-Stokes solver for
un-steady complex geometry flows on truncated domains . . . 47
III. Summation-by-parts operators with minimal dispersion error for
coarse grid flow calculations . . . 61
IV. On the order of accuracy of finite difference operators on diagonal
norm based summation-by-parts form . . . 81
DRAFT
1
Introduction
Partial differential equations (PDEs) form a fundamental tool in the mathe-matical modelling of a plethora of natural phenomena arising in a variety of scientific disciplines. Important examples include fluid flow, electromagnetic interaction and wave motion. For typical realistic models, analytic solutions are beyond reach. Instead, one resorts to numerical methods, by approximat-ing the governapproximat-ing equations at a discrete set of grid points. Such procedures are accompanied both by theoretical difficulties and computational challenges. Nonetheless, the development of computer technology during the past decades has spurred a spiralling demand for increasingly complex models to be attacked by numerical means. Consequently, the need of efficient and robust numerical methods has been steadily increasing.
A fundamental question pertaining to any PDE is whether or not it is well-posed, i.e. if a unique solution exist that depends continuously on available data. A well behaved model must answer this question in the affirmative. Similarly, corresponding numerical approximations should stand the test of stability; the discretisations must not support unphysical growth of the solution. The latter is generally difficult to achieve close to the boundaries of the computational domain, in particular for methods with high order of accuracy. Yet, high order is desirable for an efficient implementation.
Summation-by-parts (SBP) operators, augmented with simultaneous approxi-mation terms (SAT) for weak boundary treatments [4], offer a systematic ap-proach towards provably stable, yet high order accurate numerical approxima-tions. By mimicking certain properties of the underlying PDEs, an SBP-SAT discretisation yields a numerical solution that satisfies bounds analogous to the true solution of the well-posed governing equations. Although originally pro-posed for finite difference methods [9], the SBP framework has been extended to methods outside this paradigm, including finite volume methods [12, 13], spectral Galerkin and element methods [3, 6], correction procedures via recon-struction [15] and temporal discretisations [14].
In any numerical approximation, errors are inevitably present. In this thesis, we are primarily concerned with the errors of finite difference stencils on SBP form. The first part of the thesis develops a theory of dispersion errors for finite difference approximations of wave propagation problems. Such errors materialise in the form of incorrect propagation speeds, which has a pronounced
DRAFT
2 1 Introduction
impact on the quality of numerical solutions on coarse grids. New SBP operators are developed that minimise the dispersion errors in a certain sense. Numerical simulations suggest that these optimised operators perform well for problems posed on grids that are too coarse for standard techniques.
Of further interest is the order of accuracy of an SBP operator, which reveals the asymptotic behaviour of the error as the grid is progressively refined. A classic theorem states that on a uniform grid, the accuracy near domain boundaries is at most half that of the interior domain [9]. The second contribution of the thesis is the generalisation of this result to arbitrary grids, which may or may not coincide with the physical domain boundaries. Further, it is shown that the accuracy of the quadrature rule (or mass matrix), inherent in the definition of the SBP operator, is limited by the accuracy of the operator in the interior domain.
Practical implementations of numerical methods certainly encompass more than just discretising the PDE. Typically, various techniques are employed in order to reduce the errors or improve the efficiency of the numerical method at hand. Examples include interpolating multi-block techniques, multi-grid accelerators and numerical filters. For each technique, it must be guaranteed that the errors introduced are bounded such that the resulting implementation remains stable. The final part of the thesis treats a general class of problems, referred to as transmission problems, for which the above listed techniques arise as special cases. Conditions for well-posedness and stability are derived through energy estimates and bounds on the solutions are obtained via scaling arguments. The introductory chapters of this thesis are organised as follows: In Chapter 2, the notions of well-posedness and stability are formalised. We continue by introducing the energy method in Chapter 3 and apply it to a model problem. The SBP-SAT framework is introduced in Chapter 4 and its efficacy for ob-taining stable discretisations is demonstrated using the discrete energy method. Specific details about the structure and properties of finite difference stencils on SBP form are presented in Chapter 5. In Chapter 6, a brief theory of dispersion errors of finite difference stencils is introduced, together with a discussion of nu-merical filters. The papers constituting the major part of the thesis are briefly summarised in Chapter 7. We remark that for consistency of the introduc-tory chapters, some notational differences may occur relative to the subsequent papers.
DRAFT
2
Well-posedness and stability
As indicated in the previous chapter, the concepts of well-posedness and stability are fundamental in the theory of PDEs and their numerical approximations. In the following, we will formalise these notions. We further introduce the concept of a semi-bounded operator. In due course we will see that constructing semi-bounded derivative approximations provides a systematic way of ensuring stability of numerical schemes.
2.1
Well-posed problems
Consider the linear initial-boundary value problem (IBVP)
ut+ D(u) = F, t > 0, x ∈ Ω,
B(u) = g, t ≥ 0, x ∈ Γ,
u = f, t = 0, x ∈ Ω ∪ Γ,
(2.1)
where Ω is an open d-dimensional region and D is a differential operator with smooth coefficients. The operator B defines a set of boundary conditions on the boundary Γ of Ω. The functions F , g and f are given forcing, boundary and initial data, which we assume are smooth and compatible.
For two functions u and v defined on Ω, we introduce the inner product and norm (u, v)H= Z Ω u>Hvdx, kukH = (u, u) 1/2 H , (2.2)
where H is a positive definite matrix.
Definition 2.1.1. Let V be the space of differentiable functions satisfying the
boundary conditions B(v) = 0. The differential operator D is semi-bounded if for all v ∈V, D(v) satisfies the inequality
(v, D(v))H ≥ −αkvk2H,
for some real constant α independent of of v. It is maximally semi-bounded if it is semi-bounded in the space V but not in any space with fewer boundary conditions.
DRAFT
4 2 Well-posedness and stability
Note that if D is semi-bounded in (2.1), and F = g = 0, an estimate of the solution u is viable by considering
d dtkuk
2
H = 2(u, ut)H= −2(u, D(u))H≤ 2αkuk2H, (2.3)
which after integration gives ku(x, t)k2
H ≤ e
2αtkf (x)k2
H. (2.4)
Definition 2.1.2. The IBVP (2.1) is well-posed if, for sufficiently smooth and
compatible data, it has a unique smooth solution satisfying the estimate
kuk2
I ≤ K(t)kf k 2
II, (2.5)
where K(t) is independent of u and the data. It is strongly well-posed if it satisfies the estimate
kuk2
I ≤ K(t) kf k2II+ kF k2III+ kgk2IV . (2.6)
The norms involved in (2.5) and (2.6) are generally different. Comparing (2.4) and (2.5), it is clear that if a unique solution exists, semi-boundedness of D imples well-posedness of (2.1). For linear problems, maximal semi-boundedness imples existence and uniqueness, and is thus sufficient for well-posedness. We say that the solution u satisfies an energy estimate if it is bounded in terms of data such as in (2.5) or (2.6).
It follows from (2.6) that the solution of a well-posed problem depends continu-ously on the available data. Consider (2.1) with perturbed data F + δF , g + δg and f + δf , and solution v. Then, if the problem is well-posed, the difference
w = v − u satisfies the estimate
kwk2 I ≤ K(t) kδf k 2 II+ kδF k 2 III+ kδgk 2 IV .
Thus, small perturbations in the data correspond to small variations in the solution.
2.2
Stable discretisations
Stability may be viewed as the discrete analogue of well-posedness. Consider the semi-discretisation of (2.1) given by
ut+ D(u, g) = F, t > 0,
u = f , t = 0, (2.7)
where u> = (u>1, . . . , u>d) with u>j = (u(j)0 , . . . , u(j)nj), j = 1, . . . , d is a grid vector
approximating the function u(x, t) on the d-dimensional grid x>= (x>1, . . . , xd>). Here, x> j = (x (j) 0 , . . . , x (j) nj) and x (j) k = x (j) 0 + Pk l=1h (j) l , where h (j) l denotes local
DRAFT
2.2 Stable discretisations 5
The spatial discretisation D(u, g) approximates the differential operator D aug-mented with the boundary operator B. The grid vectors f , g and F are obtained by projecting the data f, g and F onto the spatial grid.
We introduce the discrete inner product and norm (u, v)h= u>P v, kukh= (u, u)
1/2
h , (2.8)
where the positive definite matrix P is such that (2.8) approximates the contin-uous inner product and norm in (2.2) on the grid x.
Definition 2.2.1. The semi-discretisation (2.7) is stable if, for sufficiently
smooth and compatible data and sufficiently fine grids, it has a unique and smooth solution satisfying the estimate
kuk2
Ih ≤ Kd(t)kf k
2
IIh, (2.9)
where Kd is independent of u, the grid and the data. It is strongly stable if it
satisfies the estimate
kuk2 Ih ≤ Kd(t) kf k 2 IIh+ kFk 2 IIh+ kgk 2 IIIh . (2.10)
Here, ’smooth data’ refers to the projection of a smooth function onto the grid. Note that if F = g = 0 and D(u, 0) in (2.7) is semi-bounded in the inner product (2.8), an energy estimate follows in the same way as for the semi-bounded continuous problem.
Much like for the continuous problem, the energy estimate (2.9) implies that small perturbations in the data correspond to small variations in the semi-discrete solution. This is of particular importance when the perturbations are representations of truncation errors in the numerical approximation. In this case, (2.9) guarantees that the errors inherent to a stable discretisation will not be uncontrollably amplified.
Of course we must also discretise time prior to implementation, whence we find ourselves with a fully discrete problem. We will use the same notation for fully discrete problems as we have done in the semi-discrete setting, and simply treat time as another coordinate. In this case, we say that the fully discrete scheme is stable if, under the conditions of Definition 2.2.1, the estimate (2.9) is satisfied at the final time.
DRAFT
3
The energy method
In many cases, the question of well-posedness may be addressed by applying the energy method. We introduce the method below by means of an example. Subsequently, we will define SBP operators and show how the discrete energy
method yields energy estimates that mimic those of the continuous PDE.
3.1
Strong boundary conditions
Consider the following special case of (2.1):
ut+ ux= 0, t > 0, x ∈ (0, 1),
u = g, t ≥ 0, x = 0, u = f, t = 0, x ∈ [0, 1].
(3.1)
Note that the boundary condition is defined only at the boundary point x = 0. We refer to formulations such as (3.1) as an IBVP with strong boundary conditions. As a special case of (2.2), we choose the inner product and norm
(u, v) = Z 1
0
u(x)v(x)dx, kuk = (u, u)1/2.
The energy method applied to (3.1) (multiplying the equation by the solution, integrating by parts and imposing the boundary condition) yields
d
dtku(x, t)k
2= g2(t) − u(1, t)2≤ g2(t). (3.2)
Comparing (3.2) with (2.3), we clearly have semi-boundedness with α = 0. Fur-ther, we cannot retain semi-boundedness if the boundary condition is removed, whence it is maximal. Temporal integration of (3.2) to time t = T gives
ku(x, T )k2≤ kf (x)k2+ Z T
0
g2(t)dt. (3.3)
Comparing (3.3) with (2.6) from Definition 2.1.2, we see that (3.1) is strongly well-posed.
Herein lies the strength of the energy method; the procedure results in an evo-lution equation for the norm of the soevo-lution (the ’energy’). Integrating by parts
DRAFT
8 3 The energy method
produces a set of boundary terms. If these have the proper sign, a semi-bounded operator is obtained and an energy estimate is readily available. However, if they are of the wrong sign, there is in general no chance of bounding them in terms of kuk, and semi-boundedness is thus beyond reach.
3.2
Weak formulations
As an alternative to the strong formulation (3.1), we may define the boundary condition on the entire domain and impose it weakly as follows:
ut+ ux= L(u − g), t > 0, x ∈ [0, 1],
u = f, t = 0, x ∈ [0, 1]. (3.4)
Here, L is a lifting operator [2], defined for smooth functions φ and ψ through the relation
Z 1
0
φL(ψ)dx = φψ|x=0.
The energy method applied to (3.4) yields d
dtku(x, t)k
2= g2(t) − u(1, t)2
− (u(0, t) − g)2≤ g2(t), (3.5)
and maximal semi-boundedness follows in the same way as for the strong for-mulation. Integrating (3.5) in time reproduces the energy estimate (3.3); hence (3.4) is strongly well-posed. Note that we may impose the initial condition weakly in a similar fashion.
3.3
Transmission
So far, we have considered a single PDE and seen that the energy method allows us to obtain an energy estimate, if appropriate boundary conditions are available. In other situations, we may have to solve two or more PDEs that are coupled across an interface. In this case, the energy method yields a set of interface terms that must be bounded. Coupling PDEs can be done both in space and in time. In the latter case, we refer to the resulting model as a
transmission problem.
As a simple example, consider
ut+ ux= 0, 0 < t < 1, x ∈ (0, 1), vt+ vx= 0, 1 < t < 2, x ∈ (0, 1), u = 0, 0 ≤ t ≤ 1, x = 0, v = 0, 1 ≤ t ≤ 2, x = 0, u = f, t = 0, x ∈ [0, 1], v = u, t = 1, x ∈ [0, 1]. (3.6)
DRAFT
3.3 Transmission 9
We seek an energy estimate for v(x, 2) in terms of the initial data f . The energy method followed by integration in time results in
ku(x, 1)k2≤ kf (x)k2,
kv(x, 2)k2≤ ku(x, 1)k2.
Adding the two equations immediately gives kv(x, 2)k2 ≤ kf (x)k2, and an
en-ergy estimate is obtained.
A more realistic problem arises when we wish to couple two problems of the general form (2.1) as follows:
ut+ D1(u) = F1, t1< t < t2, x ∈ Ω, vt+ D2(v) = F2, t2< t < t3, x ∈ Ω, B1(u) = g1, t1≤ t ≤ t2, x ∈ Γ, B2(v) = g2, t2≤ t ≤ t3, x ∈ Γ, u = f1, t = t1, x ∈ Ω ∪ Γ, v = Xu, t = t2, x ∈ Ω ∪ Γ. (3.7)
Here, we have introduced the coupling matrix X = X(x, t, u). We assume that in the decoupled situation, Xu = f2 where f2 is data, there exist matrices H1
and H2 inducing norms akin to (2.2) such that ku(x, t2)k2H1 and kv(x, t3)k
2 H2
can be estimated in terms of kf1k2H1 and kf2k
2
H2 respectively.
It is in general non-trivial to obtain an energy estimate for v(x, t3) in the
trans-mission problem (3.7). The topic of Paper V is the derivation of conditions on
H1, H2 and X such that this can be achieved. It turns out that the precise
condition for the availability of an energy estimate is the same for continuous problems with strong or weak boundary conditions, semi-discrete problems and fully discrete problems using summation-by-parts in time (see Chapter 4). An energy estimate can be obtained if and only if
H1− X>H2X ≥ 0 (3.8)
is satisfied at the transmission time t = t2, where the inequality denotes positive
semi-definiteness.
Evidently, (3.7) is a very general model, and it has a large number of impor-tant problems as special cases. Likewise, the number of possible approximation methods that fit the framework is very large. Therefore, (3.8) serves as a general condition with a wide range of applications.
DRAFT
4
Summation-by-parts
The SBP-SAT framework is designed to mimic the energy rates obtained for continuous IBVPs. Thus, if the governing equations are well-posed, the corre-sponding discretisation will be stable. In this chapter we introduce the SBP-SAT procedure for semi-discrete and fully discrete problems by discretising the model problem (3.1).
4.1
Semi-discretisations
An SBP operator may be defined as follows:
Definition 4.1.1. A matrix D = P−1Q is an SBP operator of order q if 1. Dxm= mxm−1, m = 0, . . . , q,
2. P = PT > 0,
3. Q + Q>= B = diag(−1, 0, . . . , 0, 1).
Thus, for a function u(x, t), Du approximates the spatial derivative ux on the
grid x. Here, xm should be interpreted as elementwise exponentiation. We re-mark that there are generalisations of SBP operators [5] which will be discussed in due course.
Recall from (2.8) that our discrete inner product and norm take the form (u, v)h= u>P v, kukh= (u, u)
1/2 h .
Let us now make a comparison with the continuous setting. Consider the con-tinuous inner product (u, vx) and integrate by parts:
(u, vx) = Z 1 0 uvxdx = uv|10− Z 1 0 uxvdx = uv|10− (ux, v).
Analogously, using SBP operators, consider the discrete inner product (u, Dv)h= u>Qv = u>(B − Q>)v = unvn− u0v0− (Du, v)h.
DRAFT
12 4 Summation-by-parts
Thus, the SBP property ensures that the discretisation reproduces integration by parts. This is the key property that will allow us to discretely mimic contin-uous energy estimates.
A semi-discretisation of (3.1) is given by
ut+ P−1Qu = σ0P−1(u0− g)e0, (4.1)
where the right-hand side consists of an SAT term that weakly imposes the boundary condition at the first grid point. Thus, the SAT term is an example of a discrete lifting operator. Here, e0is the first column of the identity matrix,
and σ0is a penalty parameter that will be chosen to obtain semi-boundedness.
The discrete energy method (multiplying by uTP , adding the transpose of the
result and using Definition 4.1.1) applied to (4.1) yields d dtkuk 2 h= u 2 0(2σ0+ 1) − u2n− 2σ0u0g. (4.2)
Choosing σ0= −1 gives, after adding and subtracting g2on the right-hand side
of (4.2), d dtkuk 2 h= g 2− u2 n− (u0− g)2≤ g2, (4.3)
which term for term mimics the continuous energy rate (3.5). As in the con-tinuous case, maximal semi-boundedness is thus obtained. Integrating (4.3) in time gives ku(T )k2 h≤ kf k2h+ Z T 0 g2(t)dt. (4.4)
Comparing (4.4) with (2.10) in Definition 2.2.1, it is clear that the discretisation (4.1) is strongly stable.
4.2
SBP in time
It is possible to utilise the SBP-SAT framework for temporal as well as spatial discretisations [14]. The definition of an SBP operator acting in time is com-pletely analogous to Definition 4.1.1. Before proceeding, we define the Kronecker product of two matrices A and B as
A ⊗ B = a11B . . . a1nB .. . . .. ... am1B . . . amnB .
The Kronecker product satisfies several interesting properties. Among them, we will need to use the facts that (A⊗B)>= (A>⊗B>), (A⊗B)−1= (A−1⊗B−1)
and (A ⊗ B)(C ⊗ D) = (AC ⊗ BD), assuming that the matrices C and D have appropriate dimensions.
DRAFT
4.3 Second derivatives 13
A full discretisation of (3.1) is given by
Pt−1Qt⊗ Ix u + It⊗ Px−1Qx u = Σt⊗ Px−1 ((ux0− g) ⊗ ex)
+ Pt−1⊗ Σx (et⊗ (ut0− f )) ,
(4.5)
where the subscripts t, x indicate that the matrix dimension corresponds to the temporal and spatial grid respectively. Here, ex,tdenote the first column of the
identity matrices It,x, and Σt,xare penalty matrices yet to be determined. The
elements of the vector ux0 contain the numerical solution at the first spatial
grid point projected onto the temporal grid, and similarly for ut0. Finally, f
and g contain the discrete initial and boundary data.
The energy method applied to (4.5) (multiplying by u>(Pt⊗ Px), adding the
transpose and using the properties of the SBP operators) yields, after some rearrangement, u>t nPxutn= u > t0Pxut0+ u > x0Ptux0− u > xnPtuxn + u>x 0PtΣt(ux0− g) + (ux0− g) >Σ> tPtux0 + u>t0PxΣx(ut0− f ) + (ut0− f ) >Σ> xPxut0, (4.6)
where utn and uxn denote the solution vector at the final temporal and spatial
grid point respectively (although the temporal and spatial grids of course may have different dimensions). Choosing Σt= −Itand Σx= −Ix, then adding and
subtracting f>Pxf and g>Ptg to (4.6) results in
u>tnPxutn= f >P xf − (ut0− f ) >P x(ut0− f ) + g>Ptg − (ux0− g) >P t(ux0− g) − u > xnPtuxn ≤ f>P xf + g>Ptg. (4.7)
The energy estimate (4.7) is term for term analogous to the continuous estimate (3.3). Note that the discretisation is slightly dissipative since x0− g and ut0− f
may be non-zero. This effect vanishes with grid refinement. Clearly the full discretisation (4.5) is strongly stable.
4.3
Second derivatives
The SBP framework may be extended to second (and higher) derivatives [10]. We have reason to do so in the numerical experiments in Paper III. To obtain the appropriate symmetry conditions for a second derivative SBP operator, we consider the diffusion equation, ut= uxx. The energy method results in
d dtkuk
2= 2uu
x|10− 2kuxk2, (4.8)
where we have used the same inner product and norm as before. We wish to mimic this energy rate using a second order SBP operator, D2. A suitable
DRAFT
14 4 Summation-by-parts
approximation, S is an approximation of the first derivative at the boundaries, and B is given in Definition 4.1.1. Semi-discretising gives ut = D2u, and the
energy method results in d dtkuk
2
h= 2un(Su)n− 2u0(Su)0− 2kD1uk2h. (4.9)
Clearly (4.9) mimics the continuous energy rate (4.8).
Letting D2 = D2, where D is a first derivative SBP operator in the sense of
Definition 4.1.1, yields a second derivative operator that satisfies the ansatz. To see this, note that
D2= P−1(QD) = P−1((B − Q>)D) = P−1(−D>P D + BD).
There are other ways of designing D2that enjoys certain advantages compared
to using the first derivative twice, however the above approach suffices for our needs.
4.4
Generalised SBP
Consider an IBVP posed on a spatial domain x ∈ [a, b]. Generalised SBP (gSBP) operators [5] are obtained by relaxing the property Q + Q>= B as follows:
Definition 4.4.1. A matrix D = P−1Q is a gSBP operator of order q if 1. Dxm= mxm−1, m = 0, . . . , q,
2. P = PT > 0,
3. Q + QT = eT
beb− eTaea,
4. eaxm= am, ebxm= bm m = 0, . . . , q,
that is, ea,b are boundary interpolation operators defined on the grid x.
Such operators enjoy the freedom of utilising grids that do not necessarily con-form with the physical domain boundaries. In case they do concon-form with the boundaries, Definition 4.4.1 and 4.1.1 are equivalent. We shall return to gSBP operators in Paper IV as well as in the next chapter, though under the less precise moniker ’SBP operators’.
DRAFT
5
Finite difference operators
The definition of an SBP operator specifies certain symmetry conditions on the matrices P and Q. This is sufficient to obtain energy estimates, as we have shown in the previous chapter. However, up to this point we have not given these matrices any further meaning. In the context of finite element methods,
P and Q may respectively be identified as the mass and stiffness matrices, and
their meanings are derived thence. For spectral methods, it is natural to think of P as a Gaussian quadrature rule. In Papers I-IV we assume that the SBP operator constitutes a finite difference approximation. This chapter outlines the properties of P and Q in this specific context, with particular focus on structure and accuracy.
5.1
The structure of P and Q
In this chapter, we consider a physical domain x ∈ [a, b] and let the spatial grid be given by
x = (x0, x1, . . . , xr, xr+ h, xr+ 2h, . . . ) T
, (5.1)
where x0 < x1 < · · · < xr and h is the grid spacing in the interior of the
domain. In papers I-III it is assumed that the grid is uniform, i.e. that xj =
x0+ jh for all j, and further that x0 = a, i.e. it coincides with the physical
domain boundary. In paper IV these conditions are relaxed, and in paper V, the structure of the grid is completely arbitrary. Here, we will let x0, . . . , xr be
arbitrarily distributed, and resort to Definition 4.4.1 for the generalised SBP operators.
We let an SBP operator be defined on the grid x in (5.1), and further let it be based on finite difference stencils. By this we mean that Q is a banded matrix which, away from boundaries, assumes the form of a repeated interior
stencil consisting of central finite differences. In order to close the operator near
DRAFT
16 5 Finite difference operators
The most general form of Q near the left boundary is
Q = q0,0 . . . q0,r . . . q0,r+n−1 .. . . .. ... ... qr,0 . . . qr,r . . . qr,r+n−1 an .. . ... . .. ... ... . .. qr+n−1,0 . . . qr+n−1,r . . . qr+n−1,r+n−1 a1 . . . an −an . . . −a1 0 a1 . . . an . .. . .. . .. . .. . .. .
Here, qij are the stencil coefficients of the boundary blocks. The parameters
a1, . . . , an are the coefficients of the repeated central difference stencil, defined
to be of order 2s, s > 0, through the conditions
∂u ∂x = 1 h n X k=1
ak(u(x + kh, t) − u(x − kh, t)) + O(h2s). (5.2)
Note that the interior stencil by definition operates on a uniform grid. However, recall from (5.1) that r determines the number of non-uniformly distributed grid points. Hence, the central difference stencil may not utilise any of the first r grid points, which is why n − 1 extra rows are required in the boundary block to close the SBP operator. Note that if r = 0 we recover a uniform grid. We let the matrix P take the form of a diagonal matrix,
P = h diag(p0, . . . , pr+n−1, 1, 1, . . . ).
From Definition 4.4.1 it is clear that in order to ensure positive definiteness, we must have pj > 0, j = 0, . . . , r + n − 1. Note that the diagonal elements of P are
equal to h everywhere except in the r + n rows corresponding to the boundary block of Q.
The definition of the SBP operator does not specify that P must be diagonal, and indeed it is possible to construct operators that use other structures [9, 16]. However, the stability of such operators may deteriorate for problems with variable coefficients [11]; hence we will not deliberate them further.
Throughout this thesis we assume that the right boundary is treated in the same way as the left one. In practice, this means that pN −j = pj and qN −i,N −j =
−qi,j. This assumption implies that the grid points are placed symmetrically
with respect to the physical domain, i.e. the distance from the left boundary a to xj is the same as the distance from xN −j to the right boundary b.
5.2
Accuracy and convergence
In general, an SBP operator D is less accurate near the boundary than in the interior. We let q be the order of accuracy of the boundary block, and 2s the
DRAFT
5.2 Accuracy and convergence 17
order of the interior stencil. Then, the operator differentiates polynomials up to order q exactly, as indicated by Definitions 4.1.1 and 4.4.1. For a smooth function f (x), this further implies that the error satisfies
Df − f0= O(hq),
as is verified by a Taylor expansion. Here, f and f0 are projections of f and its derivative onto the grid (5.1). A fundamental result states that the convergence rate of an SBP discretisation is q+1 for first order hyperbolic problems, and q+2 for first order parabolic or second order hyperbolic problems, if the discretisation is pointwise stable [17]. Naturally, it is desirable to have q ≥ 2s − 1 in order to preserve the convergence rate of the interior stencil. A classic result for uniform, conforming grids dictates that this is impossible when P is a diagonal matrix. In fact, it is in general impossible to choose q > s [9]. A natural question is whether it is possible to remedy this drawback by departing from uniformity and conformity, and defining the operator on the arbitrary grid (5.1). Such strategies are commonly used to increase the order of accuracy of spectral methods. A major result of Paper IV answers this question in the negative; it is not possible to choose q > s on any grid. In practice, the choice q = s is essentially always made.
Another important result derived in Paper IV is that if q = s, then P defines a quadrature rule of order 2q. This has previously been established on uniform grids [7]. A first indication to why this assertion is true is given by the so called
compatibility conditions [9]: Recall from Definition 4.4.1 that Dxj= jxj−1, j =
0, . . . q. We multiply from the left by (xi)TP to obtain
(xi)TQxj = j(xi)TP xj−1= j1TP xi+j−1,
where in the last equality we have used the facts that P is diagonal and that
x0 = 1, i.e. the grid vector of all ones. Now swap the indices i and j and add
the result to obtain the simplest form of the compatibility conditions; (xi)T(Q + QT)xj = (i + j)1TP xi+j−1, i, j = 0, . . . , q.
From the properties of Q specified in Definition 4.4.1, it follows that
(i + j)1TP xi+j−1= bi+j− ai+j, i, j = 0, . . . , q. (5.3)
In other words, P integrates polynomials up to order 2q − 1 exactly. It can be shown [1] that for a smooth function f (x), this implies that the error
Z b a f (x) − 1TP f = O(h2q),
whence P does indeed define a quadrature rule of order 2q as asserted. In Paper IV it is further shown that this order cannot be exceeded.
DRAFT
6
Dispersion
A fundamental property of PDEs describing wave motion is the dispersion
re-lation, which encodes properties such as the wave speed into the model.
Dis-cretisations have their own numeric dispersion relations, and in general they are different from the analytic one. This chapter briefly introduces the theory of dispersion errors of finite difference stencils and discusses their implications and remedies. We shall concern ourselves exclusively with spatial discretisa-tions, and thus leave the time derivatives intact. For simplicity of presentation, we restrict our attention to central difference stencils. Complications arising at domain boundaries and interfaces are considered in Paper III.
6.1
Analytic and numeric dispersion relations
Let us return to problem (3.1) but change the setting to a periodic domain. We ignore initial conditions such that we may concern ourselves exclusively with the PDE;
ut+ ux= 0, 0 ≤ t. (6.1)
Suppose that u(x, t) has a uniformly convergent Fourier expansion. We may then consider solutions of (6.1) taking the form of a general Fourier mode,
u(x, t) = exp (i(κx − ωt)), with wavenumber κ and angular frequency ω.
In-serting u into (6.1) and solving for ω gives
ω = κ. (6.2) The phase speed vp, and group speed vg, are respectively defined as
vp:=
ω
κ = 1, vg:= ∂ω ∂κ = 1.
Thus (6.1) describes a wave travelling to the right at unit speed. The relation (6.2) is known as the analytic dispersion relation. The problem (6.1) is non-dispersive in the sense that the phase speed is independent of ω and κ, i.e. all waves of any frequency propagate equally fast.
Recall from (5.2) that a central difference stencil of order 2s is given by
∂u ∂x x j ≈ D1uj:= 1 h n X k=1 ak(uj+k(t) − uj−k(t)) .
DRAFT
20 6 Dispersion
Discretising (6.1) in space thus gives duj
dt + D1uj = 0.
Again, inserting a general Fourier mode uj(t) = exp (i(κxj− ¯ωt)) and solving
for the frequency ¯ω gives the numeric dispersion relation
¯ ω = 2 h n X k=1 aksin (khκ). (6.3)
The corresponding phase and group speeds are ¯ vp:= ¯ ω κ, ¯vg:= ∂ ¯ω ∂κ. (6.4)
Note from (6.3) and (6.4) that ¯vpand ¯vgdepend on the wavenumber, κ. The
ap-proximation is thus inherently dispersive, which is what gives rise to dispersion errors and subsequent errors in phase and group velocity.
It is convenient to introduce the normalised wavenumber ξ = κh and similarly normalise (6.3) with h; ¯ ξ := ¯κh ≡ ¯ωh = 2 n X k=1 aksin (kξ). (6.5)
Here we have simply defined ¯κ = ¯ω in analogy with (6.2). Noting that the
smallest wavelength the scheme can resolve is λmin = 2h, we find that the
largest resolvable wavenumber is κmax= 2π/λmin= π/h. Thus we always have
|ξ| ≤ π.
The dispersion error is defined as the discrepancy between the analytic and numeric wavenumbers;
Ed(ξ) := ξ − ¯ξ. (6.6)
Generally, Ed scales as O(h2s+1) for a central difference stencil of order 2s,
and sufficiently small h. In Figure 6.1, ¯ξ is plotted as a function of ξ for some
common central difference stencils. Clearly, small wavenumbers are better ap-proximated than large ones, and high order methods are better than low order ones. High order methods can thus be used on coarser meshes than their low order counterparts, which is beneficial from an efficiency point of view. The topic of Papers I-II is the derivation and implementation of new, optimised fi-nite difference stencils with even smaller dispersion error than the highest order methods, for a given number of stencil coefficients. Extensions to SBP operators are presented in Paper III.
6.2
Filters
We generally say that a Fourier mode of frequency ω is resolved if the corre-sponding dispersion error Ed(ξ) is sufficiently small. Figure 6.1 clearly
DRAFT
6.2 Filters 21
Figure 6.1: Numerical wavenumbers associated with central difference stencils.
order ones. Nevertheless, no finite difference stencil has satisfactory resolution in the vicinity of ξ = π, which is the source of aliasing errors, causing insta-bilities for non-linear problems. A candidate strategy to handle this issue is to effectively remove high-frequency content from the numerical solution. One way of achieving this is by introducing a numerical filter. Filters are discussed both in Paper II and Paper V in different contexts. We briefly introduce their properties below.
The idea of a filter is to operate on the solution in such a way that its energy content is biased towards the resolved wavenumbers. The solution ujcan be
ex-pressed as a sum of Fourier modes of the general form uj =PνAνcos ν(xj− φν),
where Aν(t) is the amplitude of the mode and φν is the phase. We seek an
am-plitude response function, R(ξ), such that the filter renormalises the amam-plitudes
of the Fourier modes as RAν.
A useful response function must satisfy several criteria: It should damp or completely remove the unresolved high-frequency Fourier modes while leaving the well resolved low-frequency spectrum virtually untouched. It should have a slope that approaches zero as ξ → 0 with order 2s such that the accuracy of the stencil is not compromised. It should be everywhere positive and its magnitude should be less than unity in order to prevent the growth of any
DRAFT
22 6 Dispersion
Fourier component. A simple function that satisfies these requirements is
R(ξ) = 1 − sin2sξ
2, (6.7)
although there are other options.
A candidate approach to designing an explicit filter is to introduce an approxi-mation of the 2sth derivative operator;
∂2su ∂x2s x j ≈ D2suj := 1 h2s " b0uj+ n X k=1 bk(uj+k+ uj−k) # . (6.8)
Inserting the Fourier mode uj= exp (iκxj) in the right-hand side of (6.8) gives
the Fourier image
ψ(ξ) := b0+ 2 n
X
k=1
bkcos (kξ).
If the approximation is second order accurate and n = s in (6.8), it can be shown [8] that ψ(ξ) = (−1)s+1 2 sinξ 2 2s .
Thus, letting F = (1 + αD2s), where α = (−1)s2−2s, results in a filter F with
response function (6.7).
While we will consider filters derived as above in Paper V, we emphasise that this is one among many possible designs. In Paper II, we utilise an optimised filter, chosen to match the properties of an optimised central difference stencils derived in Paper I. The resulting scheme has a resolution capacity akin to very high order methods while retaining the efficiency of a relatively narrow difference stencil.
DRAFT
7
Summary of papers
Paper I: Uniformly best wavenumber
approxima-tions by spatial central difference operators
Consider the central finite difference stencil (5.2). From (6.5) we know that its associated dispersion error is given by
ξ − 2
n
X
k=1
aksin (kξ).
In Paper I we study the properties of the dispersion error associated with
clas-sical central differences, i.e. those with maximal order of accuracy for a given
bandwidth. It is shown that the error is monotonically increasing with respect to ξ, and consequently that the numerical wavespeed necessarily underestimates the true wavespeed associated with the governing equation.
A strategy for devising new stencils with improved dispersion properties is de-veloped. Rather than maximising the order of accuracy with respect to the bandwidth, a set of stencil parameters are used to minimise the dispersion er-ror in the L∞-sense. It is shown that for each order and bandwidth, a unique optimal stencil exists that may be characterised in terms of the extrema of its dispersion error. New stencils are derived and compared with classical stencils as well as other low-dispersion methods, showing favourable results for a linear and a non-linear test case.
Paper II: A simple and efficient incompressible
Navier-Stokes solver for unsteady complex
geom-etry flows on truncated domains
In Paper II, the stencils from Paper I are tested in a realistic setting. The entropically damped artificial compressibility formulation of the incompressible Navier-Stokes equations is solved using the optimised stencils in combination with an optimised filter of the form (6.8). This, in turn, is combined with an immersed boundary method to efficiently treat complex geometries, and a new robust outflow boundary condition to enable higher Reynolds number simula-tions on truncated domains. Validation studies for the Taylor-Green vortex
DRAFT
24 7 Summary of papers
problem and the lid driven cavity problem in both 2D and 3D are presented. An eddy viscosity subgrid-scale model is used to enable large eddy simulations for the 3D cases. The resulting scheme is shown to be able to capture strong velocity gradients on relatively coarse grids compared to previously published studies.
Paper III: Summation-by-parts operators with
min-imal dispersion error for coarse grid flow
calcula-tions
In Paper III, SBP operators with minimal dispersion errors are constructed in a way analogous to the central difference stencils from Paper I and II. Their performance is investigated for severely under-resolved problems in combina-tion with small amounts of artificial dissipacombina-tion. It is shown that performance benefits can be obtained with appropriately optimised operators for wave prop-agation and turbulent flows involving large wavenumbers, long simulation times and large ranges of resolution scales. Further, the optimised operators remain stable for non-linear problems with less artificial dissipation than high order, non-optimal operators.
Paper IV: On the order of accuracy of finite
differ-ence operators on diagonal norm based
summation-by-parts form
A finite difference SBP operator of order q near the boundaries and order 2s in the interior necessarily satisfies q ≤ s, if uniform and conforming grids are used. Paper IV investigates the possibility of increasing q by departing from uniformity and conformity, as is frequently done for spectral methods. It is proven that q ≤ s must hold for any grid, irrespective of the point distribution near the domain boundaries. Further, it is shown that if q = s, the SBP matrix
P defines a quadrature rule of order 2q irrespective of the choice of grid, and
that this order cannot be exceeded.
Paper V: Well-posed and stable transmission
prob-lems
In Paper V, we study the transmission problem (3.7), which describes the tem-poral coupling of two or more IBVPs by enforcing a transmission condition,
v(x, t2) = Xu(x, t2), at some time t = t2. Here, X is a completely general,
possibly solution-dependent matrix, describing the information flow from the first IBVP to the second one.
The problem of obtaining energy estimates for transmission problems is in-vestigated. Continuous problems with strong and weak boundary conditions,
DRAFT
25
semi-discretisations and fully discrete problems using SBP in time are analysed. It is shown that an energy estimate can be obtained if and only if the condition
H1− X>H2X ≥ 0 (7.1)
is satisfied at the transmission time t = t2, where the matrices H1,2 induce the
norms in which u and v are estimated (see (2.2)). Further, it is shown that the condition (7.1) can always be satisfied if the matrices H1,2 are appropriately
scaled; however this has a non-negligible impact on the sharpness of the result-ing energy estimate. The theory is applied to a range of realistic situations, including coupled fluid flows, multi-grid accelerators, interpolatory multi-block formulations and explicit numerical filters.
DRAFT
26 REFERENCES
References
[1] B. K. Alpert. Hybrid Gauss-trapezoidal quadrature rules. SIAM Journal
on Scientific Computing, 20(5):1551–1584, 1999.
[2] D. N. Arnold, F. Brezzi, B. Cockburn, and L. D. Marini. Unified analysis of discontinuous galerkin methods for elliptic problems. SIAM journal on
numerical analysis, 39(5):1749–1779, 2002.
[3] M. H. Carpenter and D. Gottlieb. Spectral methods on arbitrary grids.
Journal of Computational Physics, 129:74–86, 1996.
[4] M. H. Carpenter, D. Gottlieb, and S. Abarbanel. Time-stable boundary conditions for finite-difference schemes solving hyperbolic systems: method-ology and application to high-order compact schemes. Journal of
Compu-tational Physics, 111(2):220–236, 1994.
[5] D. C. D. R. Fern´andez, P. D Boom, and D. W. Zingg. A generalized
framework for nodal first derivative summation-by-parts operators. Journal
of Computational Physics, 266:214–239, 2014.
[6] G. J. Gassner. A skew-symmetric discontinuous Galerkin spectral element discretization and its relation to SBP-SAT finite difference methods. SIAM
Journal on Scientific Computing, 35(3):A1233–A1253, 2013.
[7] J. E. Hicken and D. W. Zingg. Summation-by-Parts operators and
high-order quadrature. Journal of Computational and Applied Mathematics,
237:111–125, 2013.
[8] C. A. Kennedy and M. H. Carpenter. Comparison of several numerical methods for simulation of compressible shear layers. In NASA technical
paper 3438, Langley research, 1997.
[9] H.-O. Kreiss and G. Scherer. Finite element and finite difference methods for hyperbolic partial differential equations. Aspects of Finite Elements in
Partial Differential Equations, Academic Press, Inc., 1974.
[10] K. Mattsson and J. Nordstr¨om. Summation-by-Parts operators for finite difference approximations of second derivatives. Journal of Computational
Physics, 199(2):503–540, 2004.
[11] J. Nordstr¨om. Conservative finite difference formulations, variable coeffi-cients, energy estimates and artificial dissipation. J. Sci. Comput., 29:375– 404, 2006.
[12] J. Nordstr¨om and M. Bj¨orck. Finite volume approximations and strict sta-bility for hyperbolic problems. Applied numerical mathematics, 38(3):237– 255, 2001.
DRAFT
REFERENCES 27
[13] J. Nordstr¨om, K. Forsberg, C. Adamsson, and P. Eliasson. Finite volume methods, unstructured meshes and strict stability for hyperbolic problems.
Applied Numerical Mathematics, 45(4):453–473, 2003.
[14] J. Nordstr¨om and T. Lundquist. Summation-by-Parts in time. Journal of
Computational Physics, 251:487–499, 2013.
[15] H. Ranocha, P. ¨Offner, and T. Sonar. Summation-by-Parts operators for correction procedure via reconstruction. Journal of Computational Physics, 311:299–328, 2016.
[16] B. Strand. Summation-by-Parts for finite difference approximations for d/dx. Journal of Computational Physics, 110:47–67, 1994.
[17] M. Sv¨ard and J. Nordstr¨om. On the order of accuracy for difference ap-proximations of initial-boundary value problems. Journal of Computational
