Diagonal-norm upwind SBP operators

doi:10.1016/j.jcp.2017.01.042

Journal of Computational Physics

Volume 335, 15 April 2017, Pages 283-310

https://doi.org/10.1016/j.jcp.2017.01.042 Get rights and content

Abstract

High-order accurate first derivative finite difference operators are derived that naturally introduce artificial dissipation. The boundary closures are based on the diagonal-norm summation-by-parts (SBP) framework and the boundary conditions are imposed using a penalty (SAT) technique, to guarantee linear stability for a large class of initial boundary value problems. These novel first derivative SBP operators have a non-central difference stencil in the interior, and come in pairs (for each order of accuracy). The resulting SBP-SAT approximations lead to fully explicit ODE systems. The accuracy and stability properties are demonstrated for linear first- and second-order hyperbolic problems in 1D, and for the compressible Euler equations in 2D. The newly derived first derivative SBP operators lead to significantly more robust and accurate numerical approximations, compared with the exclusive usage of (previously derived central) non-dissipative first derivative SBP operators.

Introduction

It is well known that higher order methods (as compared to first- and second-order accurate methods) capture wave dominated phenomena more efficiently since they allow a considerable reduction in the degrees of freedom, for a given error tolerance. In particular, high-order finite difference methods (HOFDMs) are ideally suited for problems of this type. (See the pioneering paper by Kreiss and Oliger [13] concerning hyperbolic problems.) The major difficulty with HOFDM is to obtain a stable boundary treatment, something that has received considerable past attention concerning hyperbolic and parabolic problems. (For examples, see [14], [33], [31], [1], [3], [10].) Roughly speaking, the numerical difficulties increase with the order of the spatial (and temporal) derivatives. For wave-dominated problems, in particular, it is imperative to use finite difference approximations that do not allow growth in time—a property termed “strict stability” [9].

A robust and well-proven high-order finite difference methodology, for well-posed initial boundary value problems (IBVP), is to combine summation-by-parts (SBP) operators [12], [32], [18] and either the simultaneous approximation term (SAT) method [4], or the projection method [28], [29], [19], [30] to impose boundary conditions (BC). Recent examples of the SBP-SAT approach can be found in [11], [16], [2], [26], [8], [25].

The SBP operators found in literature (see for example [12], [32], [18], [23], [24], [5]) are essentially central finite difference stencils closed at the boundaries with a careful choice of one-sided difference stencils, to mimic the underlying integration-by-parts formula in a discrete norm. SBP operators for various derivative orders can be constructed by repeated application of a central-difference first derivative SBP operator, here denoted $D_{1}$ . For example, $D_{1} D_{1}$ is a second derivative SBP operator, and is denoted a wide stencil second derivative SBP operator in the present study. For linear IBVP with smooth data (here referring either to physical data or the underlying curvilinear grid) an SBP-SAT approximation constructed from $D_{1}$ operators yields a strictly stable and accurate approximation. For IBVP with non-smooth data (or nonlinear problems), the exclusive usage of $D_{1}$ operators in combination with SAT (or projection) does not guarantee strict stability. For nonlinear first order hyperbolic or hyperbolic–parabolic IBVP (such as the Navier–Stokes equation) the addition of artificial dissipation (AD) is most often necessary, to damp spurious oscillations. Adding robust and accurate AD is however far from trivial, and most often involves tuning of “free” parameters. It is imperative that the addition of AD does not destroy the linear stability and accuracy properties or introduce stiffness, which in practice requires a very careful boundary closure of the added AD. How to add AD to traditional (central difference) SBP-SAT approximations based on $D_{1}$ was first analysed in [21] and later applied in [37], [20], where the notion of odd upwind SBP operators was first introduced.

For second order hyperbolic IBVP a well-known cure to suppress spurious oscillations (when non-smooth features are present) is to employ compatible narrow-stencil second derivative SBP operators (see for example [22], [23]), instead of using a wide-stencil approximation based on $D_{1}$ . The narrow-stencil second derivative SBP operators are very accurate, and have been successfully implemented in large-scale 2D and 3D problems involving mixed variable coefficient second derivative terms (see for example [39], [7]). However, compatible narrow-stencil second derivative SBP operators do not yet exist beyond 6th order accuracy.

The main motivation in [6] for introducing dual-pair SBP operators was to introduce damping of spurious oscillations when discretising second order hyperbolic equations. The dual-pair SBP operators are first derivative operators with non-central finite difference stencils in the interior. In [6] up to 8th order accurate dual-pair SBP operators are mentioned, but only the 4th order accurate case is presented explicitly. To obtain stability the dual-pair SBP operators are combined with the projection method in [6], to impose BC. In the present study, dual-pair SBP operators with an additional stability constraint are derived, and further combined with the SAT method for imposing BC. These novel SBP operators are here referred to as upwind SBP operators since they naturally introduce AD when combined with flux-splitting techniques for hyperbolic systems. (The dual-pair SBP operators in [6] do not have this property.) The novel upwind SBP operators lead to highly robust and accurate discretisations of hyperbolic problems, corroborated through numerical computations of both linear and nonlinear problems. The usefulness for hyperbolic–parabolic systems such as the Navier–Stokes equations, is indicated through stability proofs.

In Section 2 the novel upwind SBP definition is presented, including details of the necessary steps in the derivation. The relationship (and differences) to the previously derived (here refereed to as traditional) upwind SBP operators in [21] is discussed in Section 3, where the SBP-SAT method is introduced for linear first- and second-order hyperbolic problems in 1D. In Section 4 the accuracy and stability properties of the newly developed upwind SBP operators are verified and compared to previously derived (non-dissipative) central difference SBP operators, by performing numerical simulations. The stability analysis is extended to 2D hyperbolic–parabolic systems in Section 5. Verification of accuracy and stability by numerical long-time simulations of an analytic 2D Euler vortex on a multiblock grid is presented in Section 6. Section 7 summarizes the work. The upwind SBP operators are presented in the Appendix.

Section snippets

The finite difference method

Previously derived SBP operators (see for example [12], [32], [18], [23], [16], [15], [24], [5]) are essentially central finite difference stencils closed at the boundaries with a careful choice of one-sided difference stencils, to mimic the underlying integration-by-parts formula in a discrete norm. In the present paper the SBP operators are addressed by the accuracy of the interior finite difference stencil. Finite difference SBP operators may be further categorised by the structure of their

Flux-splitting

Consider $u_{t} = A u_{x}, x_{l} \leq x \leq x_{r}, t \geq 0,$ where $A = A^{T}$ is a constant coefficient matrix of size $k \times k$ . (Here also assuming initial data $u (x, 0) = f (x)$ .) Multiplying Eq. (15) by $u^{T}$ , integrating by parts and adding the transpose leads to, $\frac{d}{d t} {‖ u ‖}^{2} = u_{r}^{T} A u_{r} - u_{l}^{T} A u_{l} .$ Here $u_{l, r}^{T} = [u_{l, r}^{(1)}, \dots, u_{l, r}^{(k)}]$ are the unknowns at the left and right boundary, respectively.

Before proceeding to the semi-discrete approximation the flux vector $A u$ is split into two parts with positive and negative running characteristics, respectively (often

First order hyperbolic system

Consider the hyperbolic system defined by (24) and (27). The accuracy of the SBP-SAT approximation (29) will be verified against an analytic solution, where a narrow Gaussian pulse is propagated and reflected at the boundaries.

The following Gaussian profiles, $θ^{(1)} (x, t) = \exp {(- \frac{x - t}{r_{⁎}})}^{2}, θ^{(2)} (x, t) = - \exp {(- \frac{x - t}{r_{⁎}})}^{2},$ are introduced where $r_{⁎}$ defines the width of the Gaussian. Let $L = x_{r} - x_{l}$ denote the width of the domain. The initial data is set to: $u_{t}^{(1)} (x, 0) = 0$ , $u_{t}^{(2)} (x, 0) = 0$ , and $u^{(1)} (x, 0) = θ^{(2)} (x, 0) - θ^{(1)} (x, 0)$

Definitions

The following definition will be used in subsequent sections:

Definition 5.1

Let $\bar{x} = (x, y)$ denote grid coordinates in two dimensions. Denote the 2D bounding box $x_{l} \leq x \leq x_{r}$ , $y_{b} \leq y \leq y_{t}$ by $\bar{x} \in Ω$ , and the line $x = x_{r}$ , $y_{b} \leq y \leq y_{t}$ by $\bar{x} \in δ Ω_{E a s t}$ . A grid-function u restricted to $\bar{x} \in δ Ω_{E a s t}$ will be denoted $u_{E}$ .

The domain Ω is discretized with an $(m_{x} + 1) \times (m_{y} + 1)$ -point equidistant grid defined as $\begin{matrix} x_{i} = x_{l} + (i - 1) h_{x}, & i = 0, 1, \dots, m_{x}, & h_{x} = \frac{x_{r} - x_{l}}{m_{x} - 1}, \\ y_{j} = y_{b} + (j - 1) h_{y}, & j = 0, 1, \dots, m_{y}, & h_{y} = \frac{y_{t} - y_{b}}{m_{y} - 1} . \end{matrix}$

In the following, k denotes the number of unknowns in the

The Euler vortex problem

To test the accuracy of the novel upwind SBP operators, an Euler-vortex that satisfies the 2D Euler equations, under the assumption of isentropy is run across conforming multiblock interfaces. The problem setup consists of two blocks ( $10 \times 10$ unit area) having matching gridlines (see Fig. 5). The two blocks (left and right) are patched together at both ends using the SBP-SAT method, i.e., the east boundary of the right block is coupled to the west boundary of the left block, and the east boundary

Conclusions

The main focus has been to construct novel diagonal-norm upwind SBP operators with built in artificial dissipation when combined with flux-splitting techniques for hyperbolic systems. This is achieved by allowing non-central difference stencils and introducing a modified SBP definition. The boundary and interface conditions are treated with the SAT technique. The novel upwind SBP operators lead to highly robust and accurate discretizations of hyperbolic problems, corroborated through numerical

References (39)

Saul Abarbanel et al.
Asymptotically stable fourth-order accurate schemes for the diffusion equation on complex shapes
J. Comput. Phys.
(1997)
M. Almquist et al.
High-fidelity numerical solution of the time-dependent Dirac equation
J. Comput. Phys.
(2014)
M.H. Carpenter et al.
Time-stable boundary conditions for finite-difference schemes solving hyperbolic systems: methodology and application to high-order compact schemes
J. Comput. Phys.
(1994)
David C. Del Rey Fernández et al.
A generalized framework for nodal first derivative summation-by-parts operators
J. Comput. Phys.
(June 2014)
Leonid Dovgilovich et al.
High-accuracy finite-difference schemes for solving elastodynamic problems in curvilinear coordinates within multiblock approach
Appl. Numer. Math.
(July 2015)
Kenneth Duru et al.
Stable and high order accurate difference methods for the elastic wave equation in discontinuous media
J. Comput. Phys.
(2014)
B.A. Erickson et al.
High order accurate adaptive schemes for long time, highly intermittent geophysics problems
J. Comput. Appl. Math.
(2014)
J.E. Hicken
Output error estimation for summation-by-parts finite-difference schemes
J. Comput. Phys.
(2012)
S.K. Lele
Compact finite difference schemes with spectral-like resolution
J. Comput. Phys.
(1992)
K. Mattsson et al.
A solution to the stability issues with block norm summation by parts operators
J. Comput. Phys.
(2013)

K. Mattsson et al.

Optimal diagonal-norm SBP operators

J. Comput. Phys.

(2014)

K. Mattsson et al.

Stable and accurate wave propagation in discontinuous media

J. Comput. Phys.

(2008)

K. Mattsson et al.

Summation by parts operators for finite difference approximations of second derivatives

J. Comput. Phys.

(2004)

K. Mattsson et al.

High order finite difference methods for wave propagation in discontinuous media

J. Comput. Phys.

(2006)

K. Mattsson et al.

Stable and accurate schemes for the compressible Navier–Stokes equations

J. Comput. Phys.

(2008)

Ken Mattsson

Diagonal-norm summation by parts operators for finite difference approximations of third and fourth derivatives

J. Comput. Phys.

(2014)

Ken Mattsson et al.

High-fidelity numerical simulation of solitons in the nerve axon

J. Comput. Phys.

(2016)

N. Anders Petersson et al.

Discretizing singular point sources in hyperbolic wave propagation problems

J. Comput. Phys.

(2016)

B. Strand

Summation by parts for finite difference approximations for $d / d x$

J. Comput. Phys.

(1994)

Cited by (34)

Adaptive energy stable artificial dissipation for preserving scalar boundedness in turbulent flows
2023, Computers and Fluids
An adaptive dissipation scheme is developed that preserves scalar boundedness in a high-order finite difference framework satisfying the summation-by-parts (SBP) property. A sensor is introduced that switches from an SBP dissipation operator based on a high derivative order that is efficient at absorbing the highest energy modes to a second order derivative that targets energy across scales. Conditionally- and unconditionally-stable formulations are presented. Numerical experiments of a one-dimensional advection equation and a three-dimensional turbulent round jet demonstrate that scalar boundedness can be achieved within an acceptable threshold while retaining overall high-order accuracy. Although not demonstrated here, the same approach could be employed for shock capturing as well.
Implicit summation by parts operators for finite difference approximations of first and second derivatives
2023, Journal of Computational Physics
Implicit finite difference approximations are derived for both the first and second derivatives. The boundary closures are based on the banded-norm summation-by-parts framework and the boundary conditions are imposed using a weak (penalty) enforcement. Up to 8th order global convergence is achieved. The finite difference approximations lead to implicit ODE systems. Spectral resolution characteristics are achieved by proper tuning of the internal difference stencils. The accuracy and stability properties are demonstrated for linear hyperbolic problems in 1D and the 2D compressible Euler equations.
Dual-pairing summation by parts finite difference methods for large scale elastic wave simulations in 3D complex geometries
2022, Journal of Computational Physics
Citation Excerpt :
We call these operators Dual-Pairing SBP (DP SBP) operators. The main benefit for these operators [13,14] is that they have the potential to suppress poisonous spurious oscillations from unresolved wave-modes, which can destroy the accuracy of numerical simulations. However, these operators are asymmetric and dissipative, can potentially destroy symmetries that exist in the continuum problem.
High-order accurate summation-by-parts (SBP) finite difference (FD) methods constitute efficient numerical methods for simulating large-scale hyperbolic wave propagation problems. Traditional SBP FD operators that approximate first-order spatial derivatives with central-difference stencils often have spurious unresolved numerical wave-modes in their computed solutions. Recently derived high order accurate dual pair SBP operators based non-central (upwind) FD stencils have the potential to suppress these poisonous spurious wave-modes on marginally resolved computational grids. In this paper, we demonstrate that not all high order dual pair SBP FD operators are applicable. Numerical dispersion relation analysis shows that odd-order dual pair SBP FD operators also support spurious unresolved high-frequency wave-modes on marginally resolved meshes. Meanwhile, even-order dual pair SBP FD operators (of order $2, 4, 6$ ) do not support spurious unresolved high frequency wave-modes and also have better numerical dispersion properties.
For all the dual pair SBP FD operators we discretise the three space dimensional (3D) elastic wave equation on boundary-conforming curvilinear meshes. Using the energy method we prove that the semi-discrete approximation is stable and energy-conserving. We derive a priori error estimate and prove the convergence of the numerical error for smooth solutions. Numerical experiments for the 3D elastic wave equation in complex geometries corroborate the theoretical analysis. Numerical simulations of the 3D elastic wave equation in heterogeneous media with complex non-planar free surface topography are given, including numerical simulations of community developed seismological benchmark problems. Computational results show that even-order dual pair SBP FD operators are more efficient, robust and less prone to numerical dispersion errors on marginally resolved meshes when compared to the odd-order dual pair and traditional SBP FD operators. Finally, scaling tests demonstrate nearly perfect strong scaling and verify the efficiency of our parallel implementation of the high order upwind SBP methods.
A residual-based artificial viscosity finite difference method for scalar conservation laws
2021, Journal of Computational Physics
Citation Excerpt :
To overcome this issue, we combine the upwind SBP finite difference operators presented in Mattsson [28] with the proposed RV method. The upwind SBP operators naturally provide high-order accurate dissipation, suppressing high-frequency oscillations when paired with standard flux-splitting methods, see e.g., [28,39,25,22]. Our numerical tests show that the proposed method captures shocks accurately, is high-order accurate, does not invoke any small scale oscillations in smooth regions, and can be applied to both convex and non-convex fluxes.
In this paper, we present an accurate, stable and robust shock-capturing finite difference method for solving scalar non-linear conservation laws. The spatial discretization uses high-order accurate upwind summation-by-parts finite difference operators combined with weakly imposed boundary conditions via simultaneous-approximation-terms. The method is an extension of the residual-based artificial viscosity methods developed in the finite- and spectral element communities to the finite difference setting. The three main ingredients of the proposed method are: (i) shock detection provided by a residual-based error estimator; $(i i)$ first-order viscosity applied in regions with strong discontinuities; $(i i i)$ additional dampening of spurious oscillations provided by high-order dissipation from the upwind finite difference operators. The method is shown to be stable for skew-symmetric discretizations of the advective flux. Accuracy and robustness are shown by solving several benchmark problems in 2D for convex and non-convex fluxes.
A high-order finite-difference scheme to model the fluid-structure interaction in pneumatic seismic sources
2021, Journal of Computational Physics
A high-order accurate finite-difference scheme modeling the fluid-structure interaction inside a pneumatic seismic source is presented. The model consists of two deforming fluid compartments separated by a moving shuttle. The fluid is governed by the 1D Euler equations. Well-posedness of the continuous problem is analyzed and proven in the frozen coefficient case. A stable discretization is derived using summation-by-parts operators with the boundary conditions imposed weakly using the simultaneous-approximation-term method. The theoretical convergence rate of the numerical scheme is verified by numerical experiments. Simulation results are compared to pressure measurements from inside a pneumatic seismic source and capture many of the features observed in the data.
An efficient finite difference method for the shallow water equations
2020, Journal of Computational Physics
Citation Excerpt :
A robust and well-proven high-order finite difference methodology, for well-posed initial boundary value problems (IBVPs), is to combine summation-by-parts (SBP) operators [12,13] and the simultaneous approximation term (SAT) method [14] to impose boundary conditions (BCs). Recent examples of the SBP-SAT approach can be found in [15–20]. An added benefit of the SBP-SAT method is that it naturally extends to multi-block geometries while retaining the essential single-block properties: strict stability, accuracy, and conservation [21].
A high-order explicit finite difference scheme is derived solving the shallow water equations. The boundary closures are based on the diagonal-norm summation-by-parts (SBP) framework and the boundary conditions are imposed using a penalty (SAT) technique. Flux-splitting combined with upwind SBP operators is used to naturally introduce artificial dissipation. The scheme is tested against various benchmark problems where high-order convergence is verified for smooth solutions. A particular discretization of the source term is used leading to a well-balanced scheme. We also present an application: A simplified incident wave simulation with wave-channel interaction using a multi-block setup. Experiments suggest that a bathymetry consisting of many spikes could provide a dispersing effect on an incoming wave.

View all citing articles on Scopus

¹: Fax: +46 18 523049.

View full text

Diagonal-norm upwind SBP operators

Abstract

Introduction

Section snippets

The finite difference method

Flux-splitting

First order hyperbolic system

Definitions

The Euler vortex problem

Conclusions

J. Comput. Phys.

J. Comput. Phys.

J. Comput. Phys.

J. Comput. Phys.

Appl. Numer. Math.

J. Comput. Phys.

J. Comput. Appl. Math.

J. Comput. Phys.

J. Comput. Phys.

J. Comput. Phys.

J. Comput. Phys.

J. Comput. Phys.

J. Comput. Phys.

J. Comput. Phys.

J. Comput. Phys.

J. Comput. Phys.

J. Comput. Phys.

J. Comput. Phys.

J. Comput. Phys.