High-fidelity Sound Propagation in a Varying 3D Atmosphere

A stable and high-order accurate upwind finite difference discretization of the 3D linearized Euler equations is presented. The discretization allows point sources, a varying atmosphere and curved topography. The advective terms are discretized using recently published upwind summation-by-parts (SBP) operators and the boundary conditions are imposed using a penalty technique. The resulting discretization leads to an explicit ODE system. The accuracy and stability properties are verified for a linear hyperbolic problems in 1D, and for the 3D linearized Euler equations. The usage of upwind SBP operators leads to robust and accurate numerical approximations in the presence of point sources and naturally avoids the onset of spurious oscillations.


Introduction
Numerical methods applicable to computational aeroacoustics is an active research field, requiring a deep knowledge of flow-physics, mathematical modeling, numerical analysis and computer science. A detailed review of the progress of computational aeroacoustics can be found in [10]. Commonly, different models are used for sound generation and sound propagation. Sound generation is a highly non-linear phenomenon while sound propagation over long distances is a mostly linear phenomenon. It is therefore natural to treat these two physical phenomena with different models and numerical techniques. In the present study the main focus is to capture sound propagation over long distances efficiently, under the assumption that the sound source can be accurately modeled as a point-source. Due to the linear nature of long-range sound propagation it can be modeled using the three dimensional (3D) linearized Euler equations. This model can give more accurate results than other commonly used models. Such models include for example, ray-tracing, and diffusion based models (see e.g. [37,52]). The Euler equations includes effects from variations in the atmosphere and wind that are usually neglected in these types of models [19].
It is well known that compared to first-and second-order accurate methods, higher order methods capture wave dominated phenomena such as sound propagation more efficiently since they, for a given error tolerance, allow a considerable reduction in the degrees of freedom. In particular, high-order finite difference methods (HOFDM) are ideally suited for problems of this type. See the pioneering paper by Kreiss and Oliger [22]. 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 [1,7,14,20,24,46]). Another difficulty concerns non-smooth data such as point sources [42].
The main focus in this study is to develop a stable HOFDM for the linearized Euler equations in a 3D varying atmosphere with curved topography in the presence of point sources. A well-proven HOFDM for well-posed initial boundary value problems (IBVP), is to combine summation-by-parts (SBP) operators for approximating derivatives [23,31,45], and either the simultaneous-approximation-term (SAT) method [9], or the projection method [32,[40][41][42], to impose boundary conditions (BC). Recent examples of the SBP-SAT approach can be found in [5,15,21,27,28,36,38]. The following review papers [13,49] are highly recommended for those interested in the SBP-SAT methodology. The SBP operators found in literature are essentially central finite difference stencils (see for example [4,12,23,25,26,31,45]). These SBP operators are closed at the boundaries with a careful choice of one-sided difference stencils, to mimic the underlying integration-by-parts formula in a discrete inner product. In this study, central-difference SBP operators defined on a regular equidistant grid will be referred to as standard SBP operators.
For linear IBVP with smooth data (here referring either to physical data or the underlying curvilinear grid) an SBP-SAT approximation constructed from standard SBP operators yields a stable and accurate approximation. For IBVP with non-smooth data or non-linear problems, the exclusive usage of standard SBP operators in combination with SAT or projection does not guarantee an accurate solution, even though the scheme is stable. To damp spurious oscillations, the addition of artificial dissipation is most often necessary. Adding robust and accurate artificial dissipation is far from trivial, and commonly involves tuning of parameters. It is imperative that the addition of artificial dissipation does not destroy the stability and accuracy properties or introduce stiffness, which in practice requires a careful boundary closure of the added artificial dissipation. A procedure for adding artificial dissipation to the standard SBP operators was introduced in [34]. In [27], SBP operators defined on a regular grid with non-central finite difference stencils in the interior are introduced, referred to as upwind SBP operators. An advantage of the upwind operators, compared to standard SBP operators, is that they naturally introduce artificial dissipation.
In acoustics, efficient techniques to truncate unbounded domains is an important topic that deserves special attention. For wave equations, two commonly used techniques to truncate unbounded domains are: local high-order Absorbing Boundary Conditions (ABC) and Perfectly Matched Layers (PML). The latest developments for the ABC and PML methods are nicely summarized in [18]. As mentioned in [6], the generalization of these methods to more general hyperbolic systems such as the linearized Euler equations has not been completely successful. In [3] the PML method was generalized towards the linearized Euler equations with constant background flow but it is unclear if that method can be extended to the case with variable coefficients. An alternative to ABC and PML is found in [11], where a buffer region technique for computing compressible flows on unbounded domains is proposed. An extension of this technique for the linearized Euler equations with variable coefficients is something we would like to address in a coming study. A well-proven, first order accurate, non-reflecting BC for Euler and Navier-Stokes equations is to impose homogeneous characteristic BC (see for example [47,50]). Since boundary treatment is not the main focus of the present study, we will employ this method.
The main result of this work is the formulation of a stable and high-order accurate upwind SBP-SAT discretization applicable to the linearized Euler equations in a 3D varying atmosphere and topography with a point source. We show that the usage of upwind SBP operators, as compared to standard SBP operators, leads to more accurate and robust approximations, in particular when modeling sound generation from point sources and in wave boundary interactions. The most significant difference between the schemes is observed on coarse grids, where spurious oscillations is the main error source in the standard discretization.
In Sect. 2 the SBP operators are introduced in one dimension (1D). To illustrate some of the difficulties related to the point source in a simple setting, the stability analysis for a 1D hyperbolic system is discussed in Sect. 3. The accuracy properties are verified by performing 1D numerical simulations, including point sources. The extension to the linearized Euler equations in a 3D varying atmosphere is introduced in Sect. 4 and the stability analysis is discussed in Sect. 5. Verification of accuracy and stability by numerical studies of the 3D Euler equations is performed in Sect. 6. Section 7 summarizes the work.

The SBP Operators
In this section we introduce the upwind and standard SBP operators and discuss their properties. Before we start, some definitions are needed. Let the inner product for real-valued functions u, v ∈ L 2 [x l , x r ] be defined by (u, v) = x r x l u v dx and let the corresponding norm be u 2 = (u, u). The domain (x l ≤ x ≤ x r ) is discretized using the following N equidistant grid points: The approximate solution at grid point x i is denoted v i . The discrete solution vector is given by We define an inner product for discrete real-valued vector functions In the definitions and discretizations the following vectors (and matrix) will be used frequently: (1)

Standard SBP Operators
The standard SBP operators can be found in earlier papers (see e.g. [4,25,[29][30][31]33,35]). For completeness we restate the definition for these first derivative SBP operators: approximating ∂/∂ x, using a pth-order accurate interior stencil, is said to be a pth-order diagonal-norm first-derivative SBP operator if the diagonal matrix H defines a discrete inner product, and Q + Q T = 0.

Upwind SBP Operators
In a recent study [27] diagonal-norm SBP operators with non-central interior stencils were introduced, referred to as upwind SBP operators, approximating ∂/∂ x, using pth-order accurate interior stencils, are said to be pth-order diagonal-norm upwind SBP operators if the diagonal matrix H defines a discrete norm, One important property of the upwind SBP operators, compared to the standard operators, is the symmetric part of Q ± that when combined with flux splitting techniques, naturally introduce artificial dissipation which damps spurious oscillations. The usage of diagonalupwind SBP operators in combination with the SAT technique of imposing BC allow for general stability proofs for linear hyperbolic and hyperbolic-parabolic systems of equations (see [27] for details).

Accuracy of SBP Operators
LetD denote either a standard SBP operator (i.e., D 1 , given by Definition 2.1) or an upwind SBP operator (i.e., D ± , given by Definition 2.2). The following definition is relevant for the convergence properties of the SBP operators.

Definition 2.
3 Let x q be the projection of the polynomial x q q ! onto the discrete grid-points, i.e., the vector x. We say thatD is pth-order ifDx q = qx q−1 for q = 0 . . . p, in the interior and for q = 0 . . . p/2 or q = 0 . . . (p − 1)/2 for odd p, at the boundaries.
It is important to keep in mind that the formal boundary accuracy alone does not dictate the expected convergence rate of the numerical approximation. The expected convergence rate can be shown to be higher. In [48] it is shown that a point-wise stable approximation of an IBVP involving derivatives up to order q yields a convergence rate of order q + r , where r is the order of accuracy at the boundaries. Let p denote the interior accuracy of a diagonal-norm SBP operator (both standard and upwind). The boundary accuracy for diagonal-norm SBP operators is restricted to ( p/2)th-order accuracy when p is even (see [31]), and ( p − 1)/2 when p is odd. The expected convergence rate for first order hyperbolic problems is (r +1)th, and (r + 2)th for parabolic problems and second order hyperbolic problems. As an example both the 5th order upwind SBP operator and the 4th order standard SBP operator has r = 2 which gives the expected convergence rate 3 for first order hyperbolic IBVP. In the present study point-wise stability is not proven, but the numerical convergence studies presented in Sect. 6 indicate that the expected convergence rates from the assumption of point-wise stability are obtained. tion for the proposed numerical method is more clearly seen by first analyzing a simplified 1D scalar hyperbolic problem with constant coefficients. This simpler problem shares some of the important numerical difficulties with the full 3D linearized Euler equations, such as the numerical approximation of point sources. After comparing the upwind and standard discretization for the simplified 1D problem we will do a similar comparison for the 3D linearized Euler equations in Sect. 6.
The simplified problem we will be looking at is the following well-posed scalar hyperbolic equation, including a point source placed at x c , where δ(x) is the delta distribution and the source term function g(t) is a smooth function with g(0) = 0. Here f (x) is the initial data and g l (t) the boundary data.

Point Source Discretization
A highly accurate discretization of point sources for hyperbolic problems, based on moment conditions, is introduced in [51]. For a moment condition of order M > 0, approximating where d is the discretization of δ. However, these moment conditions alone do not guarantee a high order approximation for hyperbolic problems discretized with standard SBP operators, since there is no damping mechanism for spurious oscillations related to the Nyquist mode, (sometimes referred to as the π-mode). Spurious oscillations are typically triggered by nonsmooth features, for example point sources. In [42] a set of additional smoothness conditions were introduced, to prevent spurious oscillations when combining point sources and standard SBP operators. To enforce the smoothness conditions of order P > 0 we introduce the following additional constraint In Fig. 1 we compare the numerical approximations of (2) using standard SBP operators with and without the additional smoothness conditions for the point source (this scheme (5) is introduced in Sect. 3.2). The necessity of the smoothness conditions is obvious from the figure.
Apart from added complexity, an issue is that the smoothness conditions are constructed for the periodic problem, which introduces some complications if the sources are placed near physical boundaries. In this special case a modified boundary approximation of the smoothness conditions obtained by investigating the singular values of the difference operator is required (See [42] for details). In the present study we introduce an approach which naturally works for sources placed at the boundary.
Our approach is to replace the standard SBP operators (5) with upwind SBP operators (6), recently introduced in [27]. These operators introduce artificial dissipation which removes spurious oscillations without the additional smoothness conditions. In this work, the discretizations of the delta distribution combined with the standard and upwind operators will be denoted d S and d U respectively. In the standard case the discrete delta distribution d S obeys both the smoothness and moment conditions up to the interior accuracy order of the

Two Different Schemes
In this section, comparison of the standard and the upwind SBP-SAT schemes discretizing (2) is presented. The semi-discrete SBP-SAT approximation of (2) using a standard first derivative SBP operator D 1 is given by, where τ is a penalty parameter chosen for stability and d S is the discrete approximation of the delta distribution obeying both the moment and the smoothness conditions introduced in 3.1. The order of the moment and smoothness conditions is chosen to be the same as the accuracy order of D 1 .
The semi-discrete SBP-SAT approximation of (2) using an upwind first derivative SBP operator D − is given by, In this case, the penalty parameter τ is the same as in the standard case but the discrete approximation of the delta distribution d U , only obeys the moment conditions. Hence, in this setting no special conditions are needed when the point source is put near the boundary.

Stability Analysis
For stability, it is enough to study the homogeneous versions of (5) and (6), i.e., setting g(t) to zero. Therefore, the only difference in the stability analysis of the schemes (5) and (6) arise from the differences between the operators. Multiplying the homogeneous version of (5) by v T H and adding the transpose yields, Stability follows by choosing τ ≤ − 1 By the same procedure the upwind discretization (6) yields the energy estimate Choosing τ = − 1 (as in (8)) yields The difference between (8) and (10) is the term 2v T Sv that introduces artificial dissipation.

Numerical Experiments in 1D
In this section, the accuracy of the standard (5) and upwind (6) discretization schemes are verified, both by a convergence study and by a comparison of the errors at the outflow boundary. In the numerical experiments we consider the source term function, . Both the initial data ( f (x)) and boundary data (g l (t)) are set to zero. The semi-discretizations are integrated in time using the classical 4th order Runge-Kutta method, with a time step dt = 0.1h.
The errors of the SBP-SAT approximations (5) and (6) are verified against an analytic solution. Considering the coordinate change τ = t − (x − x c ), the solution to (2) is Figure 2 shows this solution at 3 different times.
In the convergence study the solution is integrated to t = 1, which is well before the pulse hits the boundary. The convergence rate k is calculated as  where e (N ) = u re f − u (N ) l 2 , and dim is the dimension of the problem (in 1D dim = 1 and in 3D dim = 3). The convergence results from the standard and upwind SBP schemes are presented in Tables 1 and 2 respectively. Both SBP-SAT approximations converge as expected (since the pulse does not interact with the boundary the interior accuracy of the operators k = p is expected). This verifies that the smoothness conditions are not necessary when employing upwind SBP operators.
In Fig. 3 we compare the errors at the outflow of the upwind approximation (6) with a 5th order upwind operator and the standard approximation (5) with a 4th order standard operator both with N = 201. These operators have the same boundary accuracy and are therefore globally 3rd order accurate. Here, an error is created at the source in both cases but when the pulse leaves the domain, high frequent components of the error are reflected at the boundary in the standard case, while the upwind discretization damps the error completely. (At t = 2 the l 2 norm of the error in the upwind case is 4.23 × 10 −10 ).

Euler Equations
By using the ideal gas law as equation of state, the Euler equations can be written on quasilinear form as ∂u ∂t where u = (ρ, u, v, w, p) T are the primitive variables. Here ρ is the density u, v, w the velocity components in the x, y and z direction and p is the pressure. The coefficient matrices are given by,Â where γ = 1.4 is the heat capacity ratio for air [2].
To linearize the Euler equations (13) we consider a small perturbation εu 1 around a known steady state solution to the full non-linear equations, denoted u 0 . This is done by splitting the field into two parts as u = u 0 + εu 1 , where u 0 will be referred to as the background state. The vector εu 1 is a small perturbation around the background state. The forcing function is also split asF =F 0 + εF 1 . By linearizing the Euler equations as presented in "Appendix I" the following linear system for u 1 is obtained whereÊ The matricesÂ,B andĈ can be symmetrized simultaneously (see [2] for further details) by multiplication from the left with This gives the transformed system where q = S −1 sym u and A, B and C are symmetric matrices.

Background State
When linearizing the Euler equations (13) a known steady state solution to the equations is assumed. Different forcing functionsF are chosen depending on what kind of waves the model targets. For testing purposes, we have followed the sound propagation model in [19], and set the forcing term to zero. This choice does not change the stability analysis, so the model can easily be changed to target other types of waves.
To find simple analytical solutions to the Euler equations which can be used for atmospheric modeling, we restrict the background fields in this study to variations with altitude. In a setting where the bottom boundary (Bo) is parallel to the x y-plane and the altitude increases with z the following analytic solution is obtained. Here K > 0 is constant and ρ(z) > 0, u(z) and v(z) are real valued functions. With this solution, various cases of wind speeds and speed of sound profiles varying with altitude can be modeled. However, in the case of a curved bottom boundary, we need to add the constraint u(z) = 0 and v(z) = 0 when z ∈ Bo to obey the wall boundary condition. In this case a realistic analytic background state with wind is difficult to obtain. For the general case, we would need a background state consistent with the model. This can be achieved by an accurate steady state solution to the non-linear Euler equations. However, this is not within the scope of this project.

Curvilinear Transform
In this section a transform of the Euler equations from a general curvilinear geometry ∈ (x, y, z) to a cuboid domain ∈ (ξ, η, ζ ) is presented. Consider the mapping which has the Jacobian The derivatives in each spatial dimension can be written as J q x = q ξ (y η z ζ − y ζ z η ) + q η (y ζ z ξ − y ξ z ζ ) + q ζ (y ξ z η − y η z ξ ),

Fig. 4 Computational domain in 3D
This allows rewriting (15) as These matrices are symmetric since A, B and C are symmetric and the metric coefficients are scalar.

Well-Posed Boundary Conditions
To introduce BC for the Euler equations we denote our 6 boundaries We (west), Ea (east), So (south),N o (north), Bo (bottom) and To (top) as shown in Fig. 4 where the computational 3D domain is We < ξ < Ea, So < η < N o and Bo < ζ < To. To model the atmospheres interaction with the ground a wall BC is set at the bottom boundary. At the other boundaries, characteristic boundary conditions (CBC) are imposed. CBC is a first order accurate non-reflecting BC [47,50] that is well-proven for the compressible Euler and Navier-Stokes equations. In the framework of the linearized Euler equations, the CBC can be written using the following flux splitting, where A , B and C are diagonal matrices. The derivation of these matrices can be found in [43,44]. Using this splitting, the characteristic BC are written as At the bottom boundary the normal component of the velocity is set to zero as Here L = 0 ζ x ζ y ζ z 0 where the derivatives of ζ can be expressed in terms of derivatives of x, y and z as

The Continuous Energy Analysis
In this section an energy analysis is shown for the Euler equations (16) with the BC in (17) and (18). To simplify the notation in the coming energy analysis, we will assume constant coefficient matrices. We stress that the energy analysis can be performed also for the variable coefficient case (see for example [39]).
The following inner product is used: (1) , u (2) , . . . , u (m) ] T and v = [v (1) , v (2) , . . . , v (m) ] T are vector valued functions with m components. Further, let the inner product be defined as (u,v) = Multiplying (16) by q T , setting data to zero and assuming constant coefficients gives, By adding (19) and its transpose and integrating over the domain the growth rate To obtain an estimate for the wall BC first shown in [50], using the technique presented in [16] we arrange the flux splitting ofC so that where C+ > 0, C0 = 0 C− < 0 and T C is orthogonal. Note that in this case an orthogonal T C exists since the coefficient matrixC is symmetric. This allows rewriting (18) as which is negative if ( C− + R T C+ R) is negative definite. This has been checked for all simulations performed in this work.

Semi-Discrete Analysis in 3D
In this section a standard and an upwind energy stable 3D SBP-SAT discretzation of (16) is introduced. To simplify notation in the energy analysis, we will only consider the case with constant background flow. We stress that the energy analysis can be performed also for the variable coefficient case, by a skew-symmetric splitting of the equations (see for example [39]). The skew-symmetric splitting is straightforward and unrelated to the numerical difficulty of imposing well-posed BC and will therefore not be pursued here. Hence, to simplify the notation in the coming energy analysis, A, B and C are assumed constant coefficient matrices (here consistent with E = 0). Before we present the discretizations some necessary definitions and notation are introduced.

Definitions in 3D
The continuous domain is discretized by N ξ N η N ζ points where The numerical solution at each grid point 112 . . . , q For a more compact notation we introduce the Kronecker product: where C is a p ×q matrix and D an m ×n matrix. The Kronecker product fulfills the relations Let I n be the n × n identity matrix for any n. To extend the 1D difference operators to 3D the following notation is introduced. For any N ξ × N ξ matrix P let For any N η × N η matrix P let and for any N ζ × N ζ matrix P let Further let the 1D inner product be extended as We also extend the vectors e l and e r to 3D as Further let any m × m matrix G be extended by To P E To . Additionally, let any (m N 3 × 1) vector V be extended at the boundary by

Stability Analysis
The SBP-SAT discretization of (16) with the BC in (17) and (18) can be written with the SAT-term

Lemma 1 The scheme (24) is stable if the continuous problem is well-posed and the penalty parameters are chosen as
Proof Multiplying (24) by 2q T H ξηζ and adding the transpose leads to: The discretization is stable if BT ≤ 0. Since all characteristic boundary terms are treated in a similar manner it is enough to show this for the bottom and top boundary. Due to the SBP property of D 1ζ (Definition 2.1) we get and Since well-posedness requires the matrix For the upwind formulation the characteristic variables going left in the domain are discretized with D + and the characteristics going to the right with D − . Therefore, the coefficient matrices are once again split into a positive and a negative part, this time as where α A , α B and α C are the magnitudes of the globally largest eigenvalues of the matrices A,B andC. With this splitting, the following upwind discretization is obtained, where the SAT-term is the same as in the standard case (24). This gives the energy estimate where BT < 0 is already shown. The scheme is stable since the additional terms compared to the standard case (24) provides artificial damping.

Convergence
To verify the accuracy of the discretizations, two convergence tests against numerically obtained reference solutions are performed. One in a setting with a curved bottom boundary, and one with a wind speed that varies with altitude. These tests are not combined into one test since no simple background field with wind on the curved geometry is known to the authors (for this a steady state solution to the Euler equations (13) The initial condition in pressure is set as the Gaussian , with  Tables 3  and 4.    In the second test, performed on a cuboid geometry, wind is added to the background field. A convergence test on the domain x, y, z ∈ [0, 2000] is run against a reference solution computed with the standard scheme (24) with 6th order SBP operators on a grid with N ξ , N η , N ζ = 257. The simulation is integrated until t = 3s with the initial pressure (29) Figure 6 shows snapshots of the pressure component along xz-plane, at 4 different times in the reference simulation. The convergence results are presented in Tables 5 and 6. For both settings the scheme with the standard operators shows the expected convergence (k (2) = 2, k (4) = 3 and k (6) = 4), while the upwind scheme shows a convergence rate higher than the expected rates (k (3) = 2 and k (5) = 3). We can not explain this higher convergence, however it has been observed in earlier studies with the upwind operators [27].

Resolution on Coarse Grids
To allow large 3D-simulations it is important to get accurate results on coarse grids. In this case the main difficulties are to resolve non-smooth features, as for example point sources, and to keep the accuracy at the boundaries. To compare the standard (24) and upwind (25) discretizations on coarse grids we perform a test in a setting including a point sound source by adding forcing term F = δ 3 (x − x c )g(t) in the pressure component. The 3D source discretization δ 3 (x − x c ) is a product of the 1D point source discretization in each coordinate direction, see [42] for details. The source term function is chosen as on a domain bounded by x, y, z ∈ [0, 2000] with characteristic BC at all boundaries. The initial condition is set identically to zero in all components. In Fig. 7 we compare the errors in the pressure component when using the 5th order upwind discretization (25) and the 4th order standard discretization (24) with N ξ , N η , N ζ = 33. These operators have the same boundary accuracy, therefore they are both globally 3rd order accurate. The errors are calculated against a numerically obtained reference solution computed with a 6th order standard discretization with N ξ , N η , N ζ = 257 (see Fig. 8). Note that the errors in the upwind simulation are smaller than in the standard case, both before and after the pulse leaves the domain. This difference is related to spurious oscillations which are triggered at the point source and reflected at the boundaries. In the upwind case, these oscillations are damped by the built in artificial dissipation.
In order to illustrate the effect of the artificial dissipation at the boundaries, Fig. 9 displays the discrete energy ||q|| H , of the 5th order upwind discretization (25), the 4th order standard Fig. 9 Energy decay when the pressure pulse leaves the domain for (25) and (24) and the reference solution discretization (24) both with N ξ , N η , N ζ = 33 and the reference solution (6th order standard discretization with N ξ , N η , N ζ = 257) when the pressure pulse leaves the domain. To isolate the boundary behavior, the simulations are initiated at t = 3.5s with the reference solution at this time as initial condition (see Fig. 8). The parameter setting is the same as in the reference solution (30). The energy remaining in the simulation after the pulse leaves the domain at t ≈ 7s is an undesired artifact caused mainly by reflected errors at the boundaries and errors from using characteristic BC to model a non-reflecting BC. Figure 9 shows that the energy from the upwind discretization is damped at about the same rate as the reference solution which is dominated by the error from the characteristic BC. The error in the standard discretization is magnitudes larger, and dominated by reflected spurious oscillations.

Conclusions and Future Work
The main motivation in the present study has been to derive a provably stable high-order accurate SBP-SAT approximation of the linearized Euler equations in a 3D varying atmosphere and curved topography, including point sources. This is achieved by utilizing novel upwind SBP operators with built in artificial damping.
The upwind SBP-SAT discretization of the 3D Euler equations leads to highly robust and accurate approximations, verified through numerical computations in 1D and 3D. Numerical experiments show that in the presence of point sources, the usage of upwind SBP operators efficiently avoid the onset of spurious oscillations. At the boundaries, the upwind SBP-SAT approximation yield a more efficient artificial damping of reflected waves when imposing CBC, as compared to the exclusive usage of central-difference (standard) SBP operators. In a coming study we hope to extend the current model with more accurate non-reflecting boundary conditions. We also aim for an extension of the upwind SBP-SAT methodology towards more general and non-linear problems in 3D, such as the compressible Navier-Stokes equations.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.