(PDF) Entropy Stable Staggered Grid Discontinuous Spectral Collocation Methods of any Order for the Compressible Navier--Stokes Equations

In the spirit of making high-order discontinuous Galerkin (DG) methods more competitive, researchers have developed the hybridized DG methods, a class of discontinuous Galerkin methods that generalizes the Hybridizable DG (HDG), the Embedded DG (EDG) and the Interior Embedded DG (IEDG) methods. These methods are amenable to hybridization (static condensation) and thus to more computationally efficient implementations. Like other high-order DG methods, however, they may suffer from numerical stability issues in under-resolved fluid flow simulations. In this spirit, we introduce the hybridized DG methods for the compressible Euler and Navier-Stokes equations in entropy variables. Under a suitable choice of the numerical flux, the scheme can be shown to be entropy stable and satisfy the Second Law of Thermodynamics in an integral sense. The performance and robustness of the proposed family of schemes are illustrated through a series of steady and unsteady flow problems in subsonic, transonic, and supersonic regimes. The hybridized DG methods in entropy variables show the optimal accuracy order given by the polynomial approximation space, and are significantly superior to their counterparts in conservation variables in terms of stability and robustness, particularly for under-resolved and shock flows.

The purposeof the research is to study and introduce into scientific circulation a complex of finds of decorative elements of horse harness from the outskirts of the village of Rykhtychi near Drohobych in the context of cultural and historical processes that took place in the territory of Eastern Prykarpattia and Western Volyn during the Roman period. Their participants were a multi-ethnic barbarian population that came to these lands from the Middle Danube during or after the end of the Marcomannic Wars (last quarter of the 2nd century AD) and disappeared in the phase of migration of the Gothic tribes to the Northern Black Sea region (the 3rd century AD). The methodology of the research is based on various principles of scientific study: historicism, objectivity, systematicity, which determined the application of both general scientific (analytical, statistical, historical retrospective and perspective, synthesis and analysis) and special (comparative typological, stratigraphic, cartographic, material science analysis and etc.) methods. This enabled a versatile and critical approach to the study of the problems of emergence, nature and significance of the complex of finds from the vicinity of Rykhtychi and ensured the reliability f the obtained results and conclusions. The scientific novelty of the research work consists in the study and introduction into scientific circulation of a set of decorative elements of horse harness for riding from the outskirts of the village of Rykhtychi, unique in the lands of Ukraine, and determining its place among similar complexes in the territory of Central Europe. The Conclusions. It was determined that the find of decorative elements of horse harness for riding from the vicinity of Rykhtychi reflects the complex events of the last quarter of the 2d century AD, related to the Marcomannic Wars in the Middle Danube. Based on a number of analogies from the territory of Central and Northern Europe, it was determined that the complex can be a votive offering to the gods for success in solving certain cases, a positive result in achieving the set goal, etc. In barbarian societies, they were made in water environments (rivers, lakes, swamps), on the tops or gentle slopes of hills, in crevices of rocks. The find from the vicinity of Rykhtychi may mark one of the migration routes of a certain group of people from the Middle Danube region through the Carpathians to Western Pobuzhzhia.

Entropy Stable Staggered Grid Discontinuous Spectral Collocation Methods of any Order for the Compressible Navier--Stokes Equations Item Type Article Authors Parsani, Matteo; Carpenter, Mark H.; Fisher, Travis C.; Nielsen, Eric J. Citation Parsani M, Carpenter MH, Fisher TC, Nielsen EJ (2016) Entropy Stable Staggered Grid Discontinuous Spectral Collocation Methods of any Order for the Compressible Navier--Stokes Equations. SIAM Journal on Scientific Computing 38: A3129– A3162. Available: http://dx.doi.org/10.1137/15M1043510. Eprint version Publisher's Version/PDF DOI 10.1137/15M1043510 Publisher Society for Industrial & Applied Mathematics (SIAM) Journal SIAM Journal on Scientific Computing Rights Archived with thanks to SIAM Journal on Scientific Computing Download date 25/04/2020 05:36:51 Link to Item http://hdl.handle.net/10754/621846 Downloaded 11/20/16 to 109.171.137.210. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php SIAM J. SCI. COMPUT. Vol. 38, No. 5, pp. A3129–A3162 c 2016 Society for Industrial and Applied Mathematics ENTROPY STABLE STAGGERED GRID DISCONTINUOUS SPECTRAL COLLOCATION METHODS OF ANY ORDER FOR THE COMPRESSIBLE NAVIER–STOKES EQUATIONS∗ MATTEO PARSANI† , MARK H. CARPENTER‡ , TRAVIS C. FISHER§ , AND ERIC J. NIELSEN¶ Abstract. Staggered grid, entropy stable discontinuous spectral collocation operators of any order are developed for the compressible Euler and Navier–Stokes equations on unstructured hexahedral elements. This generalization of previous entropy stable spectral collocation work [M. H. Carpenter, T. C. Fisher, E. J. Nielsen, and S. H. Frankel, SIAM J. Sci. Comput., 36 (2014), pp. B835–B867, M. Parsani, M. H. Carpenter, and E. J. Nielsen, J. Comput. Phys., 292 (2015), pp. 88–113], extends the applicable set of points from tensor product, Legendre–Gauss–Lobatto (LGL), to a combination of tensor product Legendre–Gauss (LG) and LGL points. The new semidiscrete operators discretely conserve mass, momentum, energy, and satisfy a mathematical entropy inequality for the compressible Navier–Stokes equations in three spatial dimensions. They are valid for smooth as well as discontinuous flows. The staggered LG and conventional LGL point formulations are compared on several challenging test problems. The staggered LG operators are significantly more accurate, although more costly from a theoretical point of view. The LG and LGL operators exhibit similar robustness, as is demonstrated using test problems known to be problematic for operators that lack a nonlinear stability proof for the compressible Navier–Stokes equations (e.g., discontinuous Galerkin, spectral difference, or flux reconstruction operators). Key words. high-order accurate discontinuous methods, entropy stability, SBP-SAT, compressible Navier–Stokes, staggered grid, conservation AMS subject classifications. 65N12, 65N35, 65P40 DOI. 10.1137/15M1043510 1. Introduction. Numerous strategies exist for constructing high-order accurate discontinuous Galerkin (DG) spectral element methods. Popular variants adopt either the weak (integral) or strong (differential) form of the governing equations derived by integrating the equations once or twice against a test function. Various interior and interface flux approximations are used (e.g., quadrature free fluxes [1], skew-symmetric [21]), as are various quadrature rules (e.g., Legendre–Gauss (LG),1 Legendre–Gauss–Lobatto (LGL), or Legendre–Gauss–Radau points). Each design choice is motivated by desirable goals, such as efficiency, accuracy, flexibility, etc. ∗ Submitted to the journal’s Methods and Algorithms for Scientific Computing section October 27, 2015; accepted for publication (in revised form) July 27, 2016; published electronically October 4, 2016. Some content of this paper appeared in 54th AIAA Aerospace Sciences Meeting, AIAA 2016-1058, Curran, Red Hook, NY, 2016 [13]. http://www.siam.org/journals/sisc/38-5/M104351.html Funding: This work was partially supported by King Abdullah University of Science & Technology (KAUST) in Thuwal, Saudi Arabia. † King Abdullah University of Science and Technology (KAUST), Extreme Computing Research Center (ECRC), Computer, Electrical and Mathematical Sciences & Engineering (CEMSE), Applied Mathematics and Computational Science (AMCS), Thuwal, 23955–6900, Saudi Arabia (matteo.parsani@kaust.edu.sa). ‡ Computational AeroSciences Branch (CASB), NASA Langley Research Center (LaRC), Hampton, VA 23681 (mark.h.carpenter@nasa.gov). § Computational Thermal and Fluid Mechanics, Sandia National Labs, Albuquerque, NM 87185 (tcfishe@sandia.gov). ¶ CASB, NASA LaRC, Hampton, VA 23681 (eric.j.nielsen@nasa.gov). 1 These points are also referred to as Gauss points in the literature. A3129 Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. Downloaded 11/20/16 to 109.171.137.210. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php A3130 PARSANI, CARPENTER, FISHER, AND NIELSEN Kopriva and Gassner [27] presented a survey of design decisions made when constructing nodal DG algorithms, as well as their advantages and disadvantages. A popular design philosophy for the incompressible Navier–Stokes equations is to stagger the velocity and pressure fields on independent point sets. For example, Bernardi and Maday [2] and Malay, Patera, and Rønquist [32] develop a staggered approach for the Stokes problem, while Fischer and Rønquist [16] uses a similar formulation in large-scale computational studies. Kopriva and Kolias [28] extend the staggered ideas to the compressible Navier–Stokes equations. They collocate the solution variables at the Legendre–Gauss–Chebyshev quadrature points (0, . . . , N − 1), and evaluate the fluxes at the Legendre–Gauss–Chebyshev points (0, . . . , N ). This conservative combination of over-collocated fluxes proves to be more robust in practical problems than conventional LGL techniques [28]. Staggered collocation operators have evolved over the past two decades to include a rich set of approaches. For example, the spectral difference (SD) [31, 33] and flux reconstruction (FR) [25, 45] approaches, discretize the compressible Navier–Stokes equations in strong form on a staggered set of solution (order p) and flux points (order p + 1). The observed convergence rate is reported to be (p+1) for a sequence of nested uniform and quasi-uniform grids. An alternate design strategy based on a summation-by-parts (SBP), simultaneousapproximation-term (SAT) framework (i.e., SBP-SAT operators), is used in references [8, 39] to construct discontinuous collocation (DC) spectral element methods of any order, and are referred to as entropy stable (SS) DC algorithms (i.e., SSDC algorithms). Therein, the primary design motivation is a semidiscrete spatial operator that supports a nonlinear (entropy) stability proof for the compressible threedimensional (3D) Navier–Stokes equations, on curvilinear unstructured hexahedral elements. The SSDC operators discretize the governing equations in strong form at the 3D tensor product LGL points, and adjoining elements are coupled using a provably nonlinearly stable SAT penalty approach technique [38]. The resulting algorithm is similar to the strong form nodal DG method reported in reference [23], although it substantially differs in the treatment of the nonlinear Euler fluxes and the interface couplings. A novel choice of nonlinear fluxes ensures conservation of mass, momentum, and energy as well as the entropy inequality within each element; hence elementwise entropy stability.2 The SSDC operators are distinguished from conventional finite element methods (FEMs) (e.g., DG, SD, FR) in several important ways. The nonlinear L2 stability proof is sharp; indeed, elementwise entropy conservation and entropy conservative interface fluxes guarantee global entropy conservation (neutral nonlinear stability). This sharpness is achievable because hyperviscosity dissipation, dealiasing, over-integration, or filtering of the fluxes/solution are not needed.3 (Strong form SBP operators do not require assumptions of integral exactness, a common obstacle confronting weak formulations such as conventional FEM.) A second distinguishing feature is in the treatment of solid wall boundary conditions. Nonlinearly stable solid wall BCs that maintain the nonlinear stability properties of the interior operator exist for SBP-SAT operators of any order; for details see references [39, 41]. The sharp L2 stability properties of the SSDC operators provide enhanced robustness for under-resolved smooth flows and even enable capturing shocks of moderate strength (e.g., normal shock 2 Dissipation 3 Aliasing is required for shocked flows to enforce a physical entropy inequality. is still present but does not directly lead to nonlinear instability. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. Downloaded 11/20/16 to 109.171.137.210. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php ENTROPY STABLE STAGGERED GRID DISCONTINUOUS METHOD A3131 strengths M ≤ 1.8). The captured shocks are fully consistent with the Lax–Wendroff theorem [29] for weak solutions.4 Although the formulations presented in references [8, 39] are a significant step towards operational SSDC spatial discretizations of any order, noteworthy challenges remain: (1) arbitrary collocation points, (2) spatial- (h−) and order- (p−) adaptive refinement of hexahedral elements, and (3) triangular, prismatic and tetrahedral elements. Herein, an SBP-SAT framework is used to develop a staggered grid, SS spectral element formulation that includes a broader selection of collocation points. The SS mechanics developed in [8, 39] are extended to include solutions collocated at the LG points with fluxes computed at the LGL points. The new staggered SSDC operators based on the LG points have several advantages relative to the original SSDC LGL operators. First, the staggered SSDC algorithms store the solution at the LG points, thus breaking the link between the solution storage points and the associated (suboptimal: 2p−1) LGL quadrature rules. A solution polynomial of order p can (in principle) be used with flux polynomials of arbitrary order constructed on the LGL points. Second, operators with desirable properties exist for prolonging and restricting data between the LG and LGL points, thus allowing many possibilities to move data around an element while maintaining entropy stability. Developing additional element types and connectivities can benefit from this important capability. For example, the closely related problems of h-refinement at a 2 : 1 element interface compression or a p-refinement interface, require data movement along nonconforming interfaces onto a common intermediate mortar [3]. The quadrature points do not, in general, coincide on either side of the interface, and thus, the nonlinear stability proofs presented in [8, 39] do not immediately extend to this extremely important capability. Recently, however, formulations based on the LG points were developed that naturally extend to nonconforming interfaces [11]. Third, the LG points are an interior point distribution and as such, the solution variables are not collocated at the boundaries of the element. Thus, geometric boundary discontinuities at corners are handled in an integral sense without explicit knowledge of the boundary singularity. Neither is it necessary to evaluate formulas that are known to be singular on the boundaries (e.g., turbulence source terms). Extensive numerical tests, presented herein, reveal that the new staggered SSDC operators are significantly more accurate than the LGL SSDC operators [8, 39] of equivalent polynomial order. They are, however, more costly to implement. Simple counting arguments based on inviscid and viscous flux evaluations indicate that the cost of the staggered algorithm for a solution polynomial order p is comparable to that of an LGL operator [8, 39], with a solution polynomial order of (p + 1). Preliminary studies using an explicit temporal integrator indicate that the increased accuracy of the staggered approach partially offsets the additional cost, particularly at low polynomial order. Further investigation is needed to convincingly establish whether overcollocating the fluxes (including the SD [28, 31] or FR [25, 45] operators) can be justified from a cost perspective. An ongoing investigation continues that includes the effects of implicit temporal integrators as well as the impact of data movement; computationally intensive yet extremely accurate algorithms will be competitive in the future [24]. The paper is organized as follows. Section 2 includes the general theory of SBP operators, and differentiation and interpolation spectral collocation operators. Section 3 4 The authors are not aware of a proof that density is bounded away from zero at the continuous level! Thus, departure from the real manifold is still problematic for strong shocks. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. Downloaded 11/20/16 to 109.171.137.210. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php A3132 PARSANI, CARPENTER, FISHER, AND NIELSEN includes a brief survey of entropy stability at the continuous and semidiscrete level. In section 4, the data mechanics of the staggered approach, an entropy stability proof for the one-dimensional (1D) Burgers’ equation, and its extension to the 3D compressible Navier–Stokes equations, including the interface SAT coupling approach, are presented. The comparison approach used to capture shocks and maintain nonlinear stability at the semidiscrete level is briefly presented in section 5. In section 6, a theoretical cost comparison is made between the conventional collocated [8, 39] and the staggered SSDC operators. Section 7 presents numerous numerical studies. The Euler vortex and viscous shock propagation problems are used to demonstrate the superior accuracy of the staggered algorithm. The Taylor–Green vortex problem and supersonic flow past a 3D square cylinder (Mach number = 1.5 and Reynolds number = 104 ), are used to demonstrate the robustness in the limit of order one discretization errors. Sod’s problem and the sine-shock interaction problem are used to test the staggered operators on exact discontinuous solutions. The square cylinder in supersonic crossflow is used to demonstrate algorithmic robustness in three dimensions. Section 8 summarizes the results of the paper. 2. Methodology. 2.1. Governing equations. Consider the (nonlinear) 3D compressible Navier– Stokes equations for a calorically perfect gas expressed in the form (V ) (2.1) ∂f ∂fi ∂q = i , x ∈ Ω, t ∈ [0, ∞), + ∂t ∂xi ∂xi Bq = g (B) (x, t), x ∈ ∂Ω, t ∈ [0, ∞), q(x, 0) = g (0) (x), x ∈ Ω, ⊤ where the Cartesian coordinates, x = (x1 , x2 , x3 ) , and time, t, are independent (V ) variables, and index sums are implied. The vectors q, fi , and fi are the conserved variables, and the conserved inviscid and viscous fluxes, respectively. A 3D box i i h i h h H H H × xL × xL Ω = xL 3 , x3 2 , x2 1 , x1 is chosen as our computational domain with ∂Ω representing the boundary of the domain. The boundary vector g (B) is assumed to contain linearly well-posed Dirichlet and/or Neumann data. Herein, we have omitted a detailed description of the 3D compressible Navier–Stokes equations because it can easily be found in the literature. 2.2. SBP Operators: General properties. 2.2.1. First derivative. First derivative operators that satisfy the SBP convention, discretely mimic the integration-by-parts property (2.2) Z xH φ xL ∂q xH dx = φ q|xL − ∂x Z xH xL ∂φ q dx ∂x with φ an arbitrary scalar test function. At the discrete level, this mimetic property is achieved by constructing the first derivative approximation, Dφ, with an operator in the form (2.3) D = P −1 Q, ⊤ P = P ⊤, Q = B − Q, ζ ⊤ Pζ > 0, ζ 6= 0, B = Diag (−1, 0, . . . , 0, 1) , Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. Downloaded 11/20/16 to 109.171.137.210. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php ENTROPY STABLE STAGGERED GRID DISCONTINUOUS METHOD A3133 where ζ is an arbitrary vector. The matrix P can be thought of as a mass matrix (or integrator) much like in the finite element framework, or a volume that contains local grid information in the context of finite volume or finite difference numerical methods. The nearly skew-symmetric matrix Q, is an undivided differencing operator; all rows sum to zero, as do all columns save the first and last, which sum to −1 and 1, respectively. While the matrix P need not be diagonal, the class of diagonal norm SBP operators plays a crucial role in the development of SS SBP SAT operators (see [18, 8, 39, 40, 41]). Integration in the approximation space is conducted using an inner product with the integration weights contained in the norm P, (2.4) Z xH φ xL ∂q dx ≈ φ⊤ PDq, ∂x where (2.5) ⊤ q(x) = (q(x1 ), q(x2 ), . . . , q(xN )) , with x = (x1 , . . . , xN ) , x1 = xL , xN = x H , is the projection of continuous variables q onto the grid x. Substituting (2.3) into (2.4), the mimetic SBP property is demonstrated, (2.6) φ⊤ PP −1 Qq = φ⊤ B − Q⊤ q = φN qN − φ1 q1 − φ⊤ D⊤ Pq. 2.2.2. Complementary grids. Entropy analysis is frequently performed in indicial notation on a staggered set of solution and flux points. Conventional SBP operators are not directly applicable to a staggered form of analysis. Reusing existing entropy analysis mechanics necessitates recasting (not modifying) conventional SBP operators into a staggered form. Define on the interval xL ≤ x ≤ xH , the vectors of discrete solution points (2.7) ⊤ x = (x1 , x2 , . . . , xN −1 , xN ) ; xL ≤ x1 , x2 , . . . , xN −1 , xN ≤ xH . Since the approximate solution is constructed at these points, they are denoted as the solution points. Furthermore, create a set of (N + 1) intermediate points x defining control volume bounds around each solution point. These points are denoted flux points as they are similar in nature to the control volume edges employed in the finite volume method and are used to prove the nonlinear stability (entropy stability) as briefly shown in section 3.2 (see [8, 39] for a more detailed discussion). An example of the complete discretization operator for the fourth-order accurate polynomial interpolation (p = 4) using the LGL points in the standard 1D element (i.e., xL = −1, xH = 1) is illustrated in Figure 1. In this figure, the solution points are identified with • and the flux points are identified with |. The distribution of the flux points depends on the discretization operator D = P −1 Q. The spacing between the flux points is implicitly defined by the norm P = (pij ); the diagonal elements of P (i.e., pii , i =, 1, 2, . . . , N ) are equal to the spacing between flux points and their values depend on the LGL quadrature weights, ω,5 (2.8) ⊤ x = (x̄0 , x̄1 , . . . x̄N ) , x̄0 = x1 , x̄N = xN , x̄i − x̄i−1 = pii , i = 1, 2, . . . , N. 5 See section 2.3.1 for the functional form of the elements of the matrix P and the their dependence on the quadrature weights, ω. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. Downloaded 11/20/16 to 109.171.137.210. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php A3134 PARSANI, CARPENTER, FISHER, AND NIELSEN f¯2 f¯0 f¯1 f1 u1 −9 10 f4 u4 f3 u3 x3 x2 x1 x̄0 x̄1 −1 f¯3 f2 u2 x4 x̄3 x̄2 − p3 −16 45 7 +16 45 0 + p3 7 f¯4 f¯5 f5 u5 x5 x̄4 x̄5 +9 10 +1 Fig. 1. The 1D discretization for the fourth-order accurate polynomial interpolation (p = 4) LGL collocation in the standard element is illustrated. Solution points are identified with • and the flux points are identified with |. In operator notation, this is equivalent to (2.9) ⊤ ∆x = P1; 1 = (1, 1, . . . , 1) , where the N × (N + 1) matrix ∆ is defined as (2.10)  −1  0   ∆= 0   0 0 1 −1 0 0 0 0 1 .. . 0 0  0 0    , 0 0   −1 1 0  0 −1 1 0 0 .. . 0 0 and calculates the undivided difference of the two adjacent flux points. Note that in (2.8), the first and last flux points are coincident with the first and last solution points. This fact enables the endpoint fluxes to be consistent, (2.11) f¯0 = f (q1 ), f¯N = f (qN ). This duality is needed to define unique operators and is important in proving entropy stability. Remark 2.1. Herein, all vectors are represented in bold font while matrices (when possible) are represented using calligraphy font. The “overbar” symbol is reserved for flux point vectors; the (N + 1)-dimensional flux point vectors are numbered from 0 ≤ i ≤ N . In general, matrix ranks are inferred by the dimension of the vectors on which they operate. If the rank of a matrix is ambiguous, then the rank is provided in the text. 2.2.3. Telescopic flux form. [19] proves that any SBP matrix Q can be expressed as Q = ∆ (N +1) IN : the product of the N × (N + 1) difference operator ∆ and a unique (N + 1) × N interpolation operator (N +1) IN . Thus, all SBP derivative operators D = P −1 Q can be manipulated into the telescopic flux form, (2.12) ∂ f (q) = P −1 Q f + T(p+1) = P −1 ∆f + T(p+1) , ∂x where f = (N +1) IN f , and T(p+1) is the truncation error of the approximation of ∂f (q)/∂x. Conservation follows immediately for telescopic derivative operators in the form of (2.12). Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. Downloaded 11/20/16 to 109.171.137.210. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php ENTROPY STABLE STAGGERED GRID DISCONTINUOUS METHOD A3135 2.2.4. The semidiscrete operator. Based on the previous discussion of SBP operators and their equivalent telescopic form, the semidiscrete form of (2.1) becomes ∂q (B) (Int) g i + gi = Dxi −fi + cij Dxj q + Px−1 i ∂t (V ) (B) (Int) −1 g + P = Px−1 ∆ + f −f + g (2.13) , xi i xi i i i i q(x, 0) = g(0) (x), x ∈ Ω, where the subscript xi indicates the coordinate direction to which the operators ap(B) (Int) ply. The source terms gi and gi with i = 1, 2, 3, are constructed using a penalty approach and enforce boundary and interface conditions between elements,6 respectively; and g(0) represents the initial condition. Full implementation details of the interior operators, including the definition of the viscous coefficient tensors, cij , are included in previous works [19, 20, 8, 39, 38]. 2.3. Spectral collocation operators. 2.3.1. Differentiation. Consider the SBP operators constructed at the LGL points [9], which include the end points of the interval xL = −1 and xH = 1 (see Figure 1). Define the Lagrange basis polynomials relative to the N discrete LGL points, x, as N Y x − xk , Lj (x) = xj − xk (2.14) 1 ≤ j ≤ N. k=1 k6=j A representation of the differentiation operator D, which satisfies all the requirements for being an SBP operator, is given in the following theorem. Theorem 2.1. The derivative operator that exactly differentiates an arbitrary pth-order polynomial (p = N − 1) at the collocation points, x, can be expressed as D = P −1 Q (2.15) with (2.16) P = X ℓ ⊤ L(ηl ; x)L(ηl ; x) ωℓ , Q = X ℓ L(ηl ; x) ⊤ dL (ηl ; x) ωℓ , dx where ηℓ and ωℓ , 1 ≤ l ≤ N , are the abscissaes of the LGL points and their quadrature weights, respectively. L(x; x) is a column vector whose components are the Lagrange basis polynomials relative to the discrete nodes x (i.e., Lj (x) in (2.14)) Proof. The proof to this theorem can be found in [12, Appendix B]. The matrix P in (2.16) is symmetric and positive definite for any vector x [9]. Furthermore, because in the SSDC algorithms the discrete points x are the LGL points, the matrix P is a diagonal approximation (i.e., the so-called “mass lumped” approxR1 ⊤ b b = (b imation) of the full P-norm, which is defined as P pij ) = −1 L(x; x)L(x; x) (see Appendix B in [12]). To the best of our knowledge, diagonal norm SBP operators are necessary to prove strict entropy conservation or entropy stability [17, 20, 12]. 6 Herein, the solution between adjoining elements is allowed to be discontinuous. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. Downloaded 11/20/16 to 109.171.137.210. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php A3136 PARSANI, CARPENTER, FISHER, AND NIELSEN 2.3.2. Interpolation. Define on the interval −1 = xL ≤ x ≤ xH = 1 the vectors of discrete points, ⊤ e = (e x x1 , x e2 , . . . , x eM −1 , x eM ) , (2.17) ⊤ x = (x1 , x2 , . . . , xN −1 , xN ) , −1 ≤ x e1 , x e2 , . . . , x eM −1 , x eM ≤ 1, −1 ≤ x1 , x2 , . . . , xN −1 , xN ≤ 1. e and x are the LG points (i.e., the so-called Gauss points) Herein, the discrete points x and the LGL points, respectively. All the scalars, vectors, and matrices associated with the LG points are denoted with a “tilde” symbol. Next, define the interpolation e and x: operators that move data between x e−1 RLG−LGL , ILGL→LG = P ILG→LGL = P −1 R⊤ LG−LGL , (2.18) ⊤ e ILGL→LG = ILG→LGL P P, where (2.19) RLG−LGL = Z 1 −1 ⊤ e) L(x; x) dx. L(x; x In [12], it is shown that these polynomial interpolation operators exist and satisfy the relations (2.18), provided that the LGL points are of higher polynomial orders than the LG points (i.e., N > M ). e e are the LG points, the full P-norm Note that because the x is already a positive definite diagonal matrix. 3. Entropy consistent and SS SBP operators for the compressible Navier–Stokes equations. 3.1. Overview of the continuous entropy analysis. The nonlinear compressible Navier–Stokes equations given in (2.1) have a convex extension of their original form, that when integrated over the physical domain, Ω, depends only on boundary data and dissipative terms. This convex extension is motivated by an entropy function and provides a mechanism for proving the stability of the nonlinear system in the L2 norm. Indeed, Dafermos [14] shows that if a system of conservation laws is endowed with a convex entropy function, S = S(q), a bound on the global estimate of S can be converted into an a priori estimate of the solution vector q (e.g., the solution of system (2.1)). Interested readers should also refer to the work of Svärd [42] for the proof of the stability of system (2.1) in the L2 norm. Definition 3.1. A scalar function S = S(q) is an entropy function of system (2.1) if it satisfies the following conditions: • Differentiation of the convex function S(q) simultaneously contracts all the inviscid spatial fluxes as follows: (3.1) ∂S ∂fi ∂q ∂Fi ∂q ∂Fi ∂S ∂fi = = = , ∂q ∂xi ∂q ∂q ∂xi ∂q ∂xi ∂xi i = 1, 2, 3. The components of the contracting vector, ∂S/∂q, are the entropy variables denoted as w⊤ = ∂S/∂q. The Fi (q) are the entropy fluxes in the i-direction. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. A3137 Downloaded 11/20/16 to 109.171.137.210. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php ENTROPY STABLE STAGGERED GRID DISCONTINUOUS METHOD • The entropy variables, w, symmetrize system (2.1) if w assumes the role of a new independent variable (i.e., q = q(w)). Expressing (2.1) in terms of w yields (3.2) (V ) ∂w ∂fi ∂f ∂q ∂w ∂fi ∂w ∂ ∂q b cij = 0, i = 1, 2, 3, + − i = + − ∂t ∂xi ∂xi ∂w ∂t ∂w ∂xi ∂xi ∂xj ⊤ ⊤ with the symmetry conditions: ∂q/∂w = (∂q/∂w) , ∂fi /∂w = (∂fi /∂w) , and b cij = b cij⊤ .7 An extensive entropy analysis of the Euler and compressible Navier–Stokes equations can be found, for instance, in [43, 44, 8, 20, 39, 40, 34], and the references therein. Contracting system (2.1) with the entropy variables, w, results in the differential form of the (scalar) entropy equation, (3.3) ⊤ (V ) ∂S ∂Fi ∂S ∂fi ∂ ⊤ (V ) ∂w ∂S ∂q ∂S ∂fi (V ) fi + = + = = w fi − ∂q ∂t ∂q ∂xi ∂t ∂xi ∂q ∂xi ∂xi ∂xi ⊤ ∂ ⊤ (V ) ∂w ∂w b cij . w fi − = ∂xi ∂xi ∂xj Integrating (3.3) over the domain Ω and accounting for local sources of entropy dissipation yields a global conservation statement for the entropy in the domain Z Z h i d ⊤ (V ) wx⊤i b cij wxj dx. − S dx ≤ w fi − Fi (3.4) dt Ω ∂Ω Ω Equation (3.4) is sharp (i.e., =) for smooth flows. Physical arguments require that the five-by-five matrices b cij be positive semidefinite [17, 20]. Thus, the sign of the entropy change from viscous dissipation is always negative. An increase in entropy within the domain can only result from data that convect and/or diffuse through the boundaries, ∂Ω. Remark 3.1. The density and temperature are assumed to be strictly positive, i.e., ρ, T > 0. This ensures the convexity of the entropy function S = S(q) and a one-to-one mapping between w and q. (The proof can be found in Appendix B.1 of [18]). An additional numerical mechanism must be employed to bound ρ and T away from zero to guarantee positivity; positivity preservation is not considered herein. Remark 3.2. The scalar equation (3.4) is not a sharp integral statement of entropy conservation in the presence of shocks, as it does not precisely account for the dissipation of entropy at a shock. All that is known a priori is the sign of the jump in entropy. 3.2. Overview of the semidiscrete entropy analysis for the LGL points. The semidiscrete, nonlinear entropy estimate is achieved by mimicking term by term the continuous estimate given in (3.3) and (3.4). The analysis begins by contracting the discrete entropy variables, w⊤ with the semidiscrete equation (2.13). The resulting global equation that governs the semidiscrete decay of entropy is given by [8, 39] (3.5) 7 Using w⊤ P ∂q (V ) (B) (Int) , + w⊤ ∆f i = w⊤ ∆f i + w⊤ gi + w⊤ gi ∂t i = 1, 2, 3, (V ) the entropy variables, the viscous fluxes in the i-direction are defined as fi ∂w =b cij ∂x . The explicit form of the b cij matrices can be found in [17, 37]. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. j A3138 PARSANI, CARPENTER, FISHER, AND NIELSEN Downloaded 11/20/16 to 109.171.137.210. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php where w = w(q(1),(1),(1) )⊤ , w(q(2),(1),(1) )⊤ , . . . , w(q(N ),(N ),(N ) )⊤ (B) (Int) ⊤ , and the source terms gi and gi account for the enforcement of boundary and interface conditions, respectively. The entropy variables, w, are defined at the solution (V ) points whereas the quantities with an overbar, i.e., f and f i , are defined at the flux points (see Figure 1). Simplifying (3.5) into a semidiscrete convection-diffusion equation for the entropy, proceeds as follows. The semidiscrete time term is manipulated for diagonalnorm SBP operators as w⊤ P ∂q/∂t = 1⊤ P ∂S/∂t, (or the pointwise statement: wi⊤ Pi ∂qi /∂t = Pi ∂Si /∂t ∀i ). The viscous terms are first symmetrized as in (3.2) and then manipulated as (3.6) (V ) w⊤ ∆f i cij (Dxj w). = w⊤ B b cij Dxj w − (Dxi w)⊤ P b cij Dxj w account for diffusive fluxes on boundaries while the terms The terms w⊤ B b cij (Dxj w) account for viscous dissipation. See [8] for a complete deriva−(Dxi w)⊤ P b tion of the temporal and viscous manipulations. The inviscid terms of (3.5) are entropy conservative if in each coordinate direction i, they satisfy (3.7) w⊤ ∆f i = F i (qN ) − F i (q1 ) = Fi (qN ) − Fi (q1 ) = 1⊤ ∆Fi . A design-order nonentropy conservative discrete flux f i will in general only satisfy (3.7) to design order, and could lead to numerical instability if it is not SS. A general (SC) procedure for constructing design-order, entropy conservative fluxes, f , appears (SC) in [17, 18], where it is shown that a pointwise, entropy conservative flux f m can be constructed in each coordinate direction i by using a linear combination of qℓk weighted, two-point entropy conservative fluxes f S = f S (qℓ , qk ). The implementation formula is (3.8) (SC) fm = N m X X k=m+1 ℓ=1 2 qℓk f S (qℓ , qk ) , 1 ≤ m ≤ N − 1, where the coefficient qℓk corresponds to the l row and k column in the SBP matrix Q. Herein, the two-point entropy conservative flux f S of Ismail and Roe [26] is used in (3.8). Discontinuous formulations do not enforce solution continuity across adjoining el(Int) ement interfaces. Thus, the interface penalty terms gi in (3.5) account for interface connections at adjoining elements and are formulated using an SAT approach [38]. Inviscid coupling terms are first constructed that ensure entropy conservation across the interface. Then a characteristic based dissipation term is added, that results in an upwinding of the interface fluxes. The viscous coupling terms are patterned after a local discontinuous Galerkin-type (LDG-type) approach with additional interior penalty (IP) terms to damp neutral eigenvalues. An extensive analysis of all semidiscrete terms can be found in [18, 8, 7, 38, 37, 39]. Remark 3.3. The semidiscrete entropy stability does not necessarily lead to fully discrete stability as is usually the case for linear partial differential equations. However, as noted by Tadmor [43], entropy stability is enhanced by fully implicit time discretization. In contrast, explicit time discretization, leads to entropy production. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. Downloaded 11/20/16 to 109.171.137.210. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php ENTROPY STABLE STAGGERED GRID DISCONTINUOUS METHOD f¯2 f¯0 f¯1 f1 u1 f2 u2 ũ2 ũ1 −9 10 x̃3 p3 7 x̃4 −16 45 x4 x̄3 x̄2 0 +16 45 f¯4 f¯5 f5 u5 ũ4 x3 x2 − f4 u4 ũ3 x̃2 x̃1 x1 x̄0 x̄1 −1 f¯3 f3 u3 A3139 + p3 7 x5 x̄4 x̄5 +9 10 +1 Fig. 2. The 1D discretization for the fourth-order accurate polynomial interpolation (p = 4) e are identified with × and auxiliary points with the staggered approach is illustrated. Solution points x x are identified with •. Flux points x (used to prove the entropy stability) are identified with |. Thus, the entropy stability of explicit schemes hinges on a delicate balance between temporal entropy production and spatial entropy dissipation. Fully discrete entropy conservation is offered, for instance, by Crank–Nicolson time differencing [43]. 4. Staggered SBP operators. 4.1. Data mechanics of the staggered grid approach. Define a staggered grid algorithm (Sta-Grd-Alg) for building discrete differentiation operators using two e and x of dimension M and N , respectively (see Figure 2). sets of collocation points: x The proposed algorithm is similar to that proposed by Kopriva and Kolias [28]. Ase. Furthermore, assume sume that the time-dependent solution is stored at the points x that the extrema of x coincide with the endpoints of the domain, x1 = xL , xN = xH , to facilitate the imposition of interface or boundary data. Discrete differentiation of first- or second-order spatial terms (e.g., ∂f /∂x or ∂f (V ) /∂x), by using the Sta-Grd-Alg is accomplished as follows: e to x. • Interpolate the discrete entropy variables from x • Build the nonlinear fluxes f and f (V ) on the set of points x. • Build the interface and/or boundary penalties at the extrema of x. • Differentiate the fluxes on x, and impose the penalties by using the SAT approach. e. • Interpolate the discrete flux derivatives and penalties back to x • Advance the solution with a time integration scheme by using the interpolated e. flux derivative on x Tensor product arithmetic extends the approach directly to three spatial dimensions.8 An SBP-SAT stability proof is now presented for the Sta-Grd-Alg. It is valid for all diagonal-norm SBP operators. e and e (e.g., u e , P, Define tilde variables and operators that act on the set of points x e ), and the pair of interpolation operators ILG→LGL and ILGL→LG introduced in D e to x, while the ILGL→LG section 2.3.2. The ILG→LGL operator transfers data from x e. transfers data from x to x An entropy analysis of the 1D Burgers’ equation is presented before that of the compressible Navier–Stokes equations in multiple dimensions. The matrix manipulations used to prove entropy stability for Burgers’ equation, are structurally equivalent to those used in the compressible Navier–Stokes equations, although they are much easier to follow. 8 The Sta-Grd-Alg is valid for other grid distributions that do not support tensor product arith- metic. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. Downloaded 11/20/16 to 109.171.137.210. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php A3140 PARSANI, CARPENTER, FISHER, AND NIELSEN 4.2. Entropy stability of Burgers’ equation. Consider the 1D viscous Burgers’ equation ∂ εf (V ) ∂u ∂u ∂f (u) u2 ∂x + = , f (u) = , ∂x 2 ∂t ∂x ∂u ∂u = , x ∈ [xL = −1, xH = 1], t ∈ [0, ∞), f (V ) ∂x ∂x (4.1) u xL , t + u xL , t ∂u L u xL , t − ε x , t − g xL , t = 0, ∂x 3 u xH , t − u xH , t ∂u H u xH , t − ε x , t − g xH , t = 0, 3 ∂x u (x, 0) = g (0) (x) , where u = u(x, t) is the continuous solution vector, f = f (u) is a nonlinear flux function. The boundary conditions in (4.1) are constructed such that the semidiscrete energy only increases with respect to the imposed data and maintains the same form in the inviscid limit ε → 0. Define the interpolated solution vector u, and the diagonal velocity and viscosity matrices [⊓] and [ǫ] as (4.2) e; u = ILG→LGL u e ]; [⊓] = Diag[ILG→LGL u and the boundary operator nomenclature ⊤ e xL = (1, 0, . . . , 0)M , (4.3) (Du) xL = (Du) |x=−1 = (Du)⊤ e xL , [ǫ] = Diag[ILG→LGL e ǫ]; u xL = u|x=−1 = u⊤ e xL , g xL = g|x=−1 = g⊤ e xL xH , and g xH . The subscript M with similar definitions for e xH , u xH , (Du) in (4.3) denotes the size of the vector e xL . An entropy-entropy flux pair, (S, F ), and the entropy variable, w, for Burgers’ equation are [43] 2 3 u u (4.4) (S, F ) = ; u = w. , 2 3 Note that the entropy is guaranteed to be convex for all u (i.e., ∂ 2 S/∂u2 = 1).9 Apply the Sta-Grd-Alg to construct the quadratic inviscid and linear viscous fluxes as well as the boundary penalties. The entropy analyses for the St-Grd-Alg proceeds as expected for the semidiscrete temporal, and discrete viscous, and SAT terms (e.g., see the collocated analysis given in reference [8]). Differences, however, appear in the analysis of the quadratic flux term. The quadratic flux is discretized using a diagonal norm SBP operator, D = P −1 Q, with the Q matrix rearranged into telescoping form [12, 30]. The resulting discretized expression for the quadratic term in (4.1) is 1 ∂u2 ∂ f (u) = ≈ P −1 ∆f (u). ∂x 2 ∂x The resulting semidiscrete staggered grid operator is then [12] 9 Note that the entropy for the Burgers’ equation is not unique. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. ENTROPY STABLE STAGGERED GRID DISCONTINUOUS METHOD A3141 Downloaded 11/20/16 to 109.171.137.210. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php (4.5) d e + ILGL→LG P −1 ∆ f (u) = ILGL→LG D [ǫ] D u u dt ! u xL + u xL L L L − ILGL→LG P −1 e xL u x − ǫ (Du) x − g x 3 ! u xH − u xH H H H + ILGL→LG P −1 e xH u x − ǫ (Du) x − g x 3 with the initial data u(x, 0) = g(0) (x). To simplify the notation in (4.5), the dependence on the time, t, has been omitted. What remains is to construct an entropy (SC) . conserving flux f (u) = f i In [20] it is shown that an entropy conserving flux for the delta form operator P −1 ∆f (u) can be constructed, provided there exists a two-point entropy flux relation and a diagonal-norm SBP operator is used for D. Using the definition for the entropy variable u given in (4.4) and Tadmor’s integral relation [43], Z 1 f S (uk , uℓ ) = g (w(uk ) + ξ (w(uℓ ) − w(uk ))) dξ, g(w(u)) = f (u), (4.6) 0 yields a two-point entropy conservative flux (4.7) f S (uℓ , uk ) = The entropy conserving flux is given by (4.8) (SC) fi = i N X X 1 2 uℓ + uℓ uk + u2k . 6 2 qℓk f S (uℓ , uk ) k=i+1 ℓ=1 =2 N i X X k=i+1 ℓ=1 qℓk 1 2 uℓ + uℓ uk + u2k , 6 1 ≤ i ≤ N − 1, and is of design-order accuracy [12]. Contracting the quadratic term ILGL→LG P −1 ∆f ⊤ eu e yields the telescoping condition vector P u⊤ ∆f (4.9) (SC) = 1⊤ [⊓]∆f (SC) (SC) in (4.5) with the discrete = 1⊤ ∆F = F N − F 1 with Fi = (4.10) N i X X k=i+1 ℓ=1 qℓk 1 3 3 u + uk , (uℓ + uk ) f S (uℓ , uk ) − 6 ℓ 1 ≤ i ≤ N − 1, 1 3 u . 3 N d eu e⊤P e ≤ −ε(Du)⊤ PDu Collecting all the terms in this entropy analysis yields dt u F0 = 1 3 u , 3 0 FN = which establishes boundedness of the L2 norm of u provided that the boundary data 2 2 g xL and g xH are well behaved [12]. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. Downloaded 11/20/16 to 109.171.137.210. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php A3142 PARSANI, CARPENTER, FISHER, AND NIELSEN (-1,+1) (+1,+1) (-1,-1) (+1,-1) Fig. 3. Fully staggered two-dimensional (2D) tensor product elements. Solution LG points are identified with ×, flux LGL points are identified with • and . Note that black square symbols are used only for corner LGL flux points. 4.3. Staggered grid in two dimensions. Figure 3 shows a popular staggered data structure in two spatial dimensions, with the solution stored at the tensor product LG points (i.e., the so-called Gauss points) of order p (blue crosses). Herein, this data structure is referred to as the fully staggered approach. The fully staggered approach moves the entropy variables via 3D tensor product interpolations from LG (supporting a polynomial of order p) to LGL (supporting a polynomial of order p+1) points (black circles). The discrete operators reported in [8, 7, 38, 39] are used on the LGL points to construct the spatial residual. The temporal updates needed on the LG points are obtained by restricting the LGL residuals back to the LG points. Extension to general curvilinear coordinates follows immediately on the LGL points [8, 38]. 4.4. Tensor operators in three dimensions. Consider a single tensor product element and an SSDC discretization with M = p + 1 LG solution points in each coordinate direction10 [8, 7, 39, 37]; the following elementwise matrices are needed: (4.11) e = P eM ⊗ P eM ⊗ P eM ⊗ I5 , P I LG→LGL = (ILG→LGL )x1 x2 x3 = (ILG→LGL ⊗ ILG→LGL ⊗ ILG→LGL ⊗ I5 ) , I LGL→LG = (ILGL→LG )x1 x2 x3 = (ILGL→LG ⊗ ILGL→LG ⊗ ILGL→LG ⊗ I5 ) , Dx1 = (DN ⊗ IN ⊗ IN ⊗ I5 ) , ... Dx3 = (IN ⊗ IN ⊗ DN ⊗ I5 ) , Px1 = (PN ⊗ IN ⊗ IN ⊗ I5 ) , ... Px3 = (IN ⊗ IN ⊗ PN ⊗ I5 ) , Px1 x2 = (PN ⊗ PN ⊗ IN ⊗ I5 ) , ... Px2 x3 = (IN ⊗ PN ⊗ PN ⊗ I5 ) , 10 Recall from section 2.3.2 that with the staggered algorithm the number of LG and LGL points in one dimension are denoted by M and N , respectively. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. Downloaded 11/20/16 to 109.171.137.210. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php ENTROPY STABLE STAGGERED GRID DISCONTINUOUS METHOD A3143 P = Px1 x2 x3 = (PN ⊗ PN ⊗ PN ⊗ I5 ) , Bx1 = (BN ⊗ IN ⊗ IN ⊗ I5 ) , ... Bx3 = (IN ⊗ IN ⊗ BN ⊗ I5 ) , ∆x1 = (∆N ⊗ IN ⊗ IN ⊗ I5 ) , ... ∆x3 = (IN ⊗ IN ⊗ ∆N ⊗ I5 ) , eM is the norm of the LG points, while DN , PN , ∆N , and BN are the 1D where P SBP operators [37] defined on the LGL points, and IN is the identity matrix of dimension N . I5 denotes the identity matrix of dimension five.11 The subscripts in (4.11) indicate the coordinate directions to which the operators apply (e.g., Dx1 is the differentiation matrix in the x1 direction). The symbol ⊗ represents the Kronecker product. When applying these operators to the scalar entropy equation in space at the LG points, a hat is used to differentiate the scalar operator from the full vector operator. For example, b e = P b = (PM ⊗ PM ⊗ PM ) . eM ⊗ P eM ⊗ P eM ; P (4.12) P The vector of conservative variables of each element is ordered as (4.13) ⊤ ⊤ ⊤ ⊤ ⊤ ⊤ e = qe x e(M )(M )(M ) e(1)(1)(2) , . . . , qe x q e(1)(1)(1) , qe x = qe(1) , qe(2) , . . . , qe(M , 3) where the subscripts denote the ordering of the solution points in the coordinate e directions. Assume an equivalent definition and order for the entropy variables w, and an analogous definition for the variables at the LGL points.12 e to x, while the ILGL→LG transfers The ILG→LGL operator transfers data from x e data from x to x. Define (4.14) e w = ILG→LGL w; b cij = ILG→LGL e cij ; [b cij ] = Diag[ILG→LGL e cij ]. 4.5. SS semidiscrete operator on one domain. Using these definitions and the SBP operators (4.11), system (2.1) is discretized locally on an isolated element as [37, 38] i h de q (V ) (Int) (4.15) = ILGL→LG Px−1 gi , + ILGL→LG Px−1 ∆x i f i − D x i f i i i dt where Einstein notation is used to express the coordinate directions. The penalty (Int) interface terms gi with i = 1, 2, 3 are used to connect neighboring elements (see section 4.6). The inviscid fluxes f i of any order in (4.15) are computed using linear combinations of qℓk -weighted, two-point entropy conservative fluxes (f S (uℓ , uk ) as in (3.8). Thus, they are entropy conservative. The interpolated entropy variables e are used to build the fluxes on the LGL points. The viscous fluxes are w = ILG→LGL w also computed using interpolated entropy variables and the operators Dxi , i = 1, 2, 3, defined in (4.11). The viscous coefficient matrices b cij are again formed using interpolated data on the LGL points. 11 The 3D compressible Navier–Stokes equations form a system of five nonlinear partial differential equations. 12 With the staggered algorithm, the number of LGL points in one dimension is denoted by N (see section 2.3.2). Thus, N 3 is the number of LGL points in a 3D tensor-product element. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. Downloaded 11/20/16 to 109.171.137.210. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php A3144 PARSANI, CARPENTER, FISHER, AND NIELSEN Remark 4.1. The interpolations from and to the LG points are carried out in computational space by using an efficient tensor-product algorithm that require only the knowledge of the 1D ILG→LGL and ILGL→LG operators. The extension to general curvilinear coordinates follows immediately on the LGL points [8, 39]. 4.6. SS interface coupling. Consider two cubic tensor product elements by extending (4.15) to two adjoining elements denoted with the subscripts l (left) and r (right), respectively. Without loss of generality assume that all their faces are orthogonal to the three coordinate directions and are not boundary faces, i.e., they are not part of the boundary surface ∂Ω. The resulting expressions become [37, 38] (4.16a) i h de qℓ −1 (Int),q f − D [b c ] Θ , ∆ + ILGL→LGℓ Px−1 xi ,ℓ ij,ℓ j,ℓ = ILGL→LGℓ Pxi ,ℓ gi,ℓ xi ,ℓ i,ℓ i ,ℓ dt (Int),Θ Θi,ℓ − Dxi wl = Px−1 gi,ℓ i ,ℓ (4.16b) , (4.16c) de qr (Int),q gi,r , cij,r ] Θj,r = ILGL→LGr Px−1 + ILGL→LGr Px−1 ∆xi ,r f i,r − Dxi ,r [b i ,r i ,r dt (Int),Θ Θi,r − Dxi wl = Px−1 gi,r i ,r (4.16d) . In (4.16), Θi,ℓ and Θi,r are the vectors of the gradient of the entropy variables on the (Int),q (Int),Θ left and right elements in the i direction, whereas gi,(·) and gi,(·) are the penalty interface terms on the conservative variable and the gradient of the entropy variable, respectively [37, 38]. As indicated in (4.14), the matrices [b cij ] are block diagonal matrices with N 3 five-by-five blocks corresponding to the viscous coefficients of each (V ) cij Θj . LGL point.13 Note that (4.16) is obtained by using fi = b cij wxj = b To obtain an equation for the entropy of the system, we follow the entropy stability analysis presented in [8, 37, 38]. Therefore, multiplying the two discrete equations ⊤ e l and ([b e l⊤ P in the left element by w cij,l ]Θj,l ) P l , respectively, and the two discrete ⊤ e r and ([b e r⊤ P equations in the right element by w cij,r ]Θj,r ) P r , respectively, the expression for the time derivative of the entropy function S in each element is (4.17a) 2 q d e⊤ e b e [b cij,l ] Θj,l 1 Pl S l+ 2 dt Pl bx x ,l Bbx ,l F2,l + P bx x ,l Bbx ,l F3,l bx x ,l Bbx ,l F1,l + P + 1⊤ P 1 3 2 1 2 3 2 3 1 = wl⊤ Px2 x3 ,l Bx1 ,l [b c1j,l ]Θj,l + Px1 x3 ,l Bx2 ,l [b c2j,l ]Θj,l + Px1 x2 ,l Bx3 ,l [b c3j,l ]Θj,l (Int),q (Int),q (Int),q + Px1 x3 ,l g2,l + Px1 x2 ,l g3,l + wl⊤ Px2 x3 ,l g1,l ⊤ (Int),Θ ⊤ (Int),Θ ⊤ (Int),Θ + ([b c1j,l ]Θj,l ) Px2 x3 ,l g1,l + ([b c2j,l ]Θj,l ) Px1 x3 ,l g2,l + ([b c3j,l ]Θj,l ) Px1 x2 ,l g3,l , 13 In the staggered algorithm framework, N 3 is the number of LGL points in a 3D tensor-product element. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. Downloaded 11/20/16 to 109.171.137.210. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php ENTROPY STABLE STAGGERED GRID DISCONTINUOUS METHOD A3145 (4.17b) 2 q d e⊤ e b e [b cij,r ] Θj,r 1 Pr S r + 2 dt Pr ⊤ b b b bx x ,r Bbx ,r F3,r + 1 Px2 x3 ,r Bx1 ,r F1,r + Px1 x3 ,r Bbx2 ,r F2,r + P 1 2 3 c1j,r ]Θj,r + Px1 x3 ,r Bx2 ,r [b c2j,r ]Θj,r + Px1 x2 ,r Bx3 ,r [b c3j,r ]Θj,r = wr⊤ Px2 x3 ,r Bx1 ,r [b (Int),q (Int),q (Int),q + Px1 x3 ,r g2,r + Px1 x2 ,r g3,r + wr⊤ Px2 x3 ,r g1,r ⊤ (Int),Θ ⊤ (Int),Θ (Int),Θ ⊤ + ([b c1j,r ]Θj,r ) Px2 x3 ,r g1,r + ([b c2j,r ]Θj,r ) Px1 x3 ,r g2,r + ([b c3j,r ]Θj,r ) Px1 x2 ,r g3,r , e and 1 represent vectors with M 3 and N 3 elements, respectively where the vectors 1 e = (1, 1, . . . , 1)⊤ 3 ). (i.e., 1 M To simplify the notation, assume that the interface between the two tensorproduct cells lies at x1 = 0. We also assume that all the points that lie on the other faces of the two cubes are treated in an SS fashion; their contributions can then be neglected without loss of generality. Then, for our analysis, we can just focus on a pair of LGL interface nodes at x1 = 0. We then introduce the operators e(−) and e(+) , which “extract” from the cellwise solution and flux vectors only the variables associated with these LGL points (i.e., the node at the left, (−), and at the right, (+), of the interface).14 Therefore, (4.17) reduces to (4.18a) d e⊤ b elS e l + 1⊤ P bx x ,l F1,l + 2 1 P 2 3 dt = wl⊤ Px2 x3 ,l q [b cij,l ]Θj,l [b c1j,l ]Θj,l e (−) + (4.18b) = −wr⊤ Px2 x3 ,r (Int),Θ ⊤ q (Int),q g1,l , 2 [b cij,r ]Θj,r [b c1j,r ]Θj,r e ⊤ Pl wl⊤ Px2 x3 ,l + ([b c1j,l ]Θj,l ) Px2 x3 ,l g1,l d e⊤ b e r − 1⊤ P er S bx x ,r F1,r + 2 1 P 2 3 dt 2 Pr (+) (Int),Θ + ([b c1j,r ]Θj,r ) Px2 x3 ,r g1,r (Int),q + wr⊤ Px2 x3 ,r g1,r . The interface penalty terms are constructed as a combination of a LDG-type approach and an IP technique [38]: (Int),q g1,l (4.19a) 14 These i h (−) = +f 1 − f (SS) q (−) , q (+) e(−) h i h i 1 (−) (−) (+) (+) c1,j Θj − b c1,j Θj + − (1 + α) b e(−) 2 1 (−) (+) + [L] w −w e(−) , 2 two vectors are zero at all points except at the “(−), (+)” interface points. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. A3146 Downloaded 11/20/16 to 109.171.137.210. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php (4.19b) PARSANI, CARPENTER, FISHER, AND NIELSEN (Int),Θ g1,l (Int),q g1,r (4.19c) (4.19d) (Int),Θ g1,r 1 = − (1 − α) w(−) − w(+) e(−) , 2 h (+) i = −f 1 + f (SS) q (−) , q (+) e(+) h i h i 1 (+) (+) (−) (−) c1,j Θj − b + + (1 − α) b c1,j Θj e(+) 2 1 + [L] w(+) − w(−) e(+) , 2 1 (+) (−) = + (1 + α) w e(+) . −w 2 The LDG penalty terms involve the coefficients 12 (1 ± α) and act only in the normal direction to the face. The IP terms involve the block diagonal parameter matrix, [L] = Diag [L], with N 3 five-by-five blocks, L, which are left unspecified for the moment.15 Herein, the solution between adjoining elements is allowed to be discontinuous. An inviscid interface flux that preserves the entropy consistency of the interior high order accurate spatial operators [8] on either side of the interface f (SS) q (−) , q (+) is constructed as (4.20) f (SS) q (−) , q (+) = f (SC) q (−) , q (+) + Λ w(+) − w(−) , where f (SC) q (−) , q (+) is the entropy conservative inviscid interface flux of any order [18, 8, 7, 37, 39] constructed as in (3.8). Λ is a negative semidefinite interface matrix with zero or negative eigenvalues. The superscripts (−) and (+) denote the collocated values onthe left and right sides of the interface, respectively. The SS flux f (SS) q (−) , q (+) is more dissipative than the entropy conservative inviscid flux f (SC) q (−) , q (+) , as can be easily verified by contracting f (SS) q (−) , q (+) against the entropy variables [8, 39]. Note that in [41], grid interfaces for SS finite difference schemes are studied and interface fluxes similar to (4.20) are proposed. Substituting expressions (4.19) and (4.20) into (4.18) and summing all the contributions of the two elements result in (4.21) # " 2 2 q q d e⊤ b d e⊤ b e e e e = Υ(I) + Υ(V ) , + [b cij,l ]Θj,l [b cij,r ]Θj,r 1 P l Sl + 1 P r Sr + 2 dt dt Pl Pr where Υ(I) and Υ(V ) are the inviscid and the viscous interface terms. At the two interface nodes, these terms are ⊤ Υ(I) = w(+) − w(−) f (SS) q (−) , q (+) − ψ (+) − ψ (−) (4.22) ⊤ = w(+) − w(−) Λ w(+) − w(−) , (4.23) 15 In ⊤ Υ(V ) = w(+) − w(−) L w(+) − w(−) . the staggered algorithm framework, N 3 is the number of LGL points in a 3D tensor-product element. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. Downloaded 11/20/16 to 109.171.137.210. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php ENTROPY STABLE STAGGERED GRID DISCONTINUOUS METHOD A3147 Clearly, the interface contributions are dissipative if both the five-by-five matrices Λ and L are negative semidefinite. The matrix Λ can be constructed using different approaches. Herein, an upwind operator that dissipates each characteristic wave based on the magnitude of its eigenvalue is used [34, 8, 37, 39]. We are then left with the viscous interface term, Υ(V ) . The parameter values α = 0 and α = ±1 yield a symmetric LDG and a “flip-flop” narrow stencil (nearest neighbor) LDG penalty, respectively [37, 39]. An LDG value of α = 0 produces a global discrete operator that has a neutrally damped spurious eigenmode. Herein, when α = 0, this mode is damped using the IP dissipation. For a Reynolds number that approaches ∞, we would like the five-by-five matrix L to go to zero so that only the inviscid SS penalty contributions in (4.19) are recovered. To achieve that, the matrix L is constructed as (−) L = −β (Int) (4.24) (−) (+) (+) b c11 + b c11 , 2 px1 ,11 β (Int) > 0, where b c11 and b c11 are the positive semidefinite viscous coefficient matrices at the left and right side of the interface in the normal direction. The coefficient β (Int) can be used to modify the strength of the IP penalty term, although excessively large values of β (Int) reduce the maximum stable time step. In references [37, 39, 38], β (Int) is selected to be the maximum value for which the explicit stability constraint remains unaffected. The factor px1 ,11 in the denominator represents the normal local grid spacing which is incorporated in the diagonal SBP operator P.16 This term is introduced to get the correct dimension, and, as for the standard IP finite element approach, it increases the strength of L with increased resolution. 5. The comparison approach using entropy conservative schemes. Entropy conservative formulations suffer breakdown when used without dissipation to capture shocks. The entropy discontinuity at the shock is inconsistent with the underlying premise of an isentropic algorithm, and large amplitude oscillations around the shock are the result. More problematic, however, is that entropy conservative formulations cannot converge to the weak solution; there is no mechanism to admit the dissipation physically required at the shock. Herein the comparison approach presented in [8] is used to facilitate the addition of dissipation to the entropy conservative baseline formulation. A comparison approach uses a companion algorithm in conjunction with an entropy conservative formulation (or other entropy datum). At every point in the domain, the entropy generated by the companion scheme is compared with that of the entropy datum. If the entropy condition is violated (i.e., the inequality derived in (3.4)), then more dissipation is added locally. (A popular example of schemes developed using a comparison approach are the E-schemes of Osher [36] that use a Godunov scheme as the entropy datum.) The sufficient local conditions for entropy stability are [8] (SC) ⊤ ≤ 0, i = 1, 2, . . . , N − 1, (5.1) (wi+1 − wi ) f¯i − f¯i where f¯i is the inviscid Euler flux constructed using a high-order accurate operator (SC) which might not (always) be SS (e.g., a nodal DG method) and f¯i is the entropy 16 Herein, pxi ,11 with i = 1 appears in the definition of L because the interface between the two elements is orthogonal to the x1 direction [37, 38]. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. Downloaded 11/20/16 to 109.171.137.210. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php A3148 PARSANI, CARPENTER, FISHER, AND NIELSEN conservative flux constructed as given in (3.8). These conditions can be enforced by using a limiter function of the form (5.2) (SS) (SC) f¯i = f¯i + δ(f¯i − f¯i ), δ = √ b2 + c 2 − b (SC) ⊤ √ , b = (wi+1 − wi ) (f¯i − f¯i ), c = 10−12 2 b2 + c 2 (SS) with f¯i the entropy stabilized flux. The pointwise entropy stability conditions (SC) given in (5.1) and (5.2) are valid for any pair of comparison fluxes f¯i and f¯i , provided that both can be expressed in telescoping flux form [8]. Note, however, that the local conditions in (5.1) do not provide insight into the magnitude of dissipation required to achieve a nonoscillatory shock. 6. Discussion: A theoretical cost analysis. A pth-order 3D tensor-product 3 LGL element has (p + 1) collocation points. The total cost of the fully discrete, explicit-in-time, operator predominantly results from the inviscid and viscous differentiation operators. Thus, the operation count for a conventional nonstaggered 4 element [8, 39] scales as (p + 1) . By a similar accounting, the operation count for the 4 staggered operator scales as (p + 2) . Assuming that the explicit time stepping stability constraint is equivalent for both schemes, the work ratio between the conventional and staggered algorithm is given by (6.1) p+2 p+1 4 = Wst . Note that this theoretical cost estimate does not account for hardware effects such as cache efficiency and data movement. At present, the current implementation of the staggered SSDC algorithm is approximately 1.7–2.0 times more costly than that of conventional SSDC algorithms. Optimization of both SSDC algorithms for current and future computing hardware architectures is the goal of future work. After optimization, a thorough comparison between the two approaches will be performed. The staggered formulation is indeed more costly than the conventional SSDC algorithm. Alternatively, the results shown in section 7 demonstrate that the fully staggered SSDC approach is significantly more accurate than the conventional LGL SSDC operator for the same solution polynomial and grid density. Thus, a significantly finer mesh is needed by the conventional LGL SSDC operator to achieve a comparable error level. Simple accounting arguments that compare the respective number of elements needed to achieve comparable error tolerances, suggest that the staggered formulation will be competitive, particularly for low polynomial orders (e.g., 0 ≤ p ≤ 3). We reiterate that more testing is needed before establishing the viability of the staggered approach. 7. Numerical results: Accuracy and robustness. 7.1. Isentropic Euler vortex propagation. The test case considered in this section is the (inviscid) propagation of a vortex for which an exact solution is known. This is an excellent test problem for verifying the accuracy and functionality of the inviscid components of a compressible Navier–Stokes solver. It is fully described by Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. ENTROPY STABLE STAGGERED GRID DISCONTINUOUS METHOD A3149 Downloaded 11/20/16 to 109.171.137.210. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php (7.1) h i 2 2 f (x1 , x2 , x3 , t) = 1 − (x1 − x1,0 − U∞ cos(α) t) + (x2 − x2,0 − U∞ sin(α) t) , 2 γ−1 exp (f (x , x , x , t)) , T (x1 , x2 , x3 , t) = 1 − ǫ2v M∞ 1 2 3 8π 2 1 ρ(x1 , x2 , x3 , t) = T γ−1 , f (x1 , x2 , x3 , t) x2 − x2,0 − U∞ sin(α)t , exp u1 (x1 , x2 , x3 , t) = U∞ cos(α) − ǫv 2π 2 x1 − x1,0 − U∞ cos(α)t f (x1 , x2 , x3 , t) u2 (x1 , x2 , x3 , t) = U∞ sin(α) − ǫv , exp 2π 2 u3 (x1 , x2 , x3 , t) = 0, where U∞ , M∞ , and (x1,0 , x2,0 , x3,0 ) are the module of the freestream velocity, the Mach number, and the coordinates of the vortex center, respectively. In this study, the values U∞ = M∞ c∞ , ǫv = 5.0, M∞ = 0.5, γ = 1.4, and α = 45o are used; the domain is described by x1 ∈ (−5, 5), x2 ∈ (−5, 5), x3 ∈ (−1, 1), (x1,0 , x2,0 , x3,0 ) = (0, 0, 0), t ≥ 0. The boundary conditions are prescribed by penalizing the numerical solution against the exact solution by using an SAT approach [10], whereas the interface coupling between two adjoining elements is imposed using the SS treatment proposed in references [8, 38]. The accuracy of the following SS schemes is investigated using uniform Cartesian and unstructured nonuniform grids: • Fully staggered SSDC algorithm with pLG = 1, 2, 3, 4, 5, 10 and pLGL = 2, 3, 4, 5, 6, 11. • Conventional LGL SSDC algorithm [8, 39] with pLGL = 1, 2, 3, 4, 5, 10. 7.1.1. Uniform Cartesian grid. Different grid resolutions are examined, and the vortex is halfway out of the domain when the error measure is evaluated. This measures the effect of the penalty boundary conditions and the interior scheme. Figure 4 shows the convergence study for a sequence of maximum nine nested grids (2 × 2 × 2, 4 × 4 × 2, . . . 512 × 512 × 2)17 . The L2 norms of the error decay asymptotically towards the designed rate in each case (i.e., second order, third order, fourth order, fifth order, sixth order, and eleventh order, respectively), for both fully staggered and the conventional SSDC approaches. We highlight that the fully staggered approach is significantly more accurate than the conventional algorithm for the same grid resolution. However, simple theoretical counting arguments based on inviscid flux evaluations (see section 6), seem to indicate that the cost of the staggered algorithm for a solution polynomial order p is comparable to that of an LGL operator [8, 39] with a solution polynomial order of (p + 1). 7.1.2. Unstructured nonuniform grids. The goal of this study is to investigate whether the fully staggered algorithm is effectively (p + 1)-order accurate for more realistic meshes (p is the order of the polynomial supported by the LG points). As for the case of uniform Cartesian grids, different grid resolutions are examined, and the vortex is halfway out of the domain when the error measure is evaluated. A nested 17 The number of grid cells is doubled every time in each coordinate direction. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. 10-1 10-2 10-3 10-4 10-5 10-6 10-7 10-8 10-9 10-10 10-11 10-12 PARSANI, CARPENTER, FISHER, AND NIELSEN Grid convergence study: Euler vortex on Cartesian uniform grids Fully-staggered, Fully-staggered, Fully-staggered, Fully-staggered, Fully-staggered, Fully-staggered, Conventional, Conventional, Conventional, Conventional, Conventional, Conventional, L2 Error ( -norm) Downloaded 11/20/16 to 109.171.137.210. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php A3150 10-2 1/(DOFs) 1/3 10-1 pLG =1, pLGL =2 pLG =2, pLGL =3 pLG =3, pLGL =4 pLG =4, pLGL =5 pLG =5, pLGL =6 pLG =10, pLGL =1 pLGL =2 pLGL =3 pLGL =4 pLGL =5 pLGL =10 pLGL =11 100 Fig. 4. Uniform grid convergence study for the inviscid vortex test case. Fully staggered SSDC algorithm (solid lines); conventional LGL SSDC algorithm [8, 39]: Dashed lines; reference curves (p + 1): Red solid lines. Fig. 5. “Unstructured grid kernel” used to construct a sequence of nested grids for the inviscid vortex test case. family of grids is constructed by replicating N times the “unstructured grid kernel” shown in Figure 5. The grids used are: 1 × 1 × 2, 2 × 2, . . . , 16 × 16 × 2, where the first integers represent the number of times that the kernel is repeated in the x1 and x2 directions; and the last integer indicates the number elements in the x3 direction. Figure 6, show the convergence study for the sequence of maximum sixteen nested grids. For the fully staggered SSDC algorithm, the L2 norms of the error decay asymptotes towards the designed rate (i.e., (p + 1)) and, in each case, is more accurate than the conventional LGL SSDC algorithm for the same grid resolution. The conventional path converges instead to p. 7.2. Viscous shock. The second test case considered in this section is the propagation of a viscous shock for which an exact time-dependent solution is known. The compressible Navier–Stokes equations support an exact solution for the viscous shock profile, under the assumption that the Prandtl number P r = 34 . Mass and total enthalpy are constant across a shock. Furthermore, if P r = 34 then the momentum and Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. A3151 10-1 10-2 10-3 10-4 10-5 10-6 10-7 10-8 10-9 10-10 10-11 10-12 10-13 Grid convergence study: Euler vortex on unstructured nonuniform grids Fully-staggered, Fully-staggered, Fully-staggered, Fully-staggered, Fully-staggered, Fully-staggered, Conventional, Conventional, Conventional, Conventional, Conventional, Conventional, L2 Error ( -norm) Downloaded 11/20/16 to 109.171.137.210. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php ENTROPY STABLE STAGGERED GRID DISCONTINUOUS METHOD 10-1 1/(DOFs) pLG =1, pLGL =2 pLG =2, pLGL =3 pLG =3, pLGL =4 pLG =4, pLGL =5 pLG =5, pLGL =6 pLG =10, pLGL =1 pLGL =2 pLGL =3 pLGL =4 pLGL =5 pLGL =10 pLGL =11 100 1/3 Fig. 6. Nonuniform grid convergence study for the inviscid vortex test case. Fully-staggered SSDC algorithm (solid lines); conventional LGL SSDC algorithm [8, 39]: Dashed lines; reference curves (p + 1): Red solid lines. energy equations are redundant. If we assume that the shock is propagating along the x1 coordinate direction,18 the single momentum equation across the shock is given by (7.2) ∂v − (v − 1)(v − vf ) = 0; −∞ ≤ x1 ≤ ∞ , t ≥ 0; αv ∂x 1 u1,right γ−1 µ u1 ; vf = ;α = v = , u1,lef t u1,lef t 2γ P r ṁ where ṁ is the constant mass flow across the shock. An exact solution is obtained by solving the momentum equation (7.2) for the velocity profile, v: 1+v v−1 . (7.3) x1 = 21 α Log |(v − 1)(v − vf )| + 1−vff Log v−v f A moving shock is recovered by applying a uniform translation to the solution. A full derivation of this solution appears in the thesis of Fisher [17]. In this study, the values U∞ = M∞ c∞ , M∞ = 2.5, Re∞ = 10, and γ = 1.4 are used. The viscous shock, which at t = 0 is located at the center of the computational domain, is propagated in the direction parallel to the horizontal axis. The domain is described by x1 ∈ (−10, 10), x2 ∈ (−10, 10), x3 ∈ (−1, 1), t ≥ 0. The boundary conditions on the faces perpendicular to the horizontal axis are prescribed by penalizing the numerical solution against the exact solution; periodic boundary conditions are used on the remaining four boundary faces of the computational domain. This is an excellent test problem for verifying the accuracy and functionality of the inviscid and viscous components of a compressible Navier–Stokes solver. The accuracy of the following SS schemes is investigated using uniform Cartesian and unstructured nonuniform grids: 18 This assumption is just used to simplify the notation. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. PARSANI, CARPENTER, FISHER, AND NIELSEN 100 10-1 10-2 10-3 10-4 10-5 10-6 10-7 10-8 10-9 10-10 10-11 10-12 10-13 Grid convergence study: Viscous shock on Cartesian uniform grids Fully-staggered, Fully-staggered, Fully-staggered, Fully-staggered, Fully-staggered, Fully-staggered, Conventional, Conventional, Conventional, Conventional, Conventional, Conventional, L2 Error ( -norm) Downloaded 11/20/16 to 109.171.137.210. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php A3152 10-2 1/(DOFs) 1/3 pLG =1, pLGL =2 pLG =2, pLGL =3 pLG =3, pLGL =4 pLG =4, pLGL =5 pLG =5, pLGL =6 pLG =10, pLGL =1 pLGL =2 pLGL =3 pLGL =4 pLGL =5 pLGL =10 10-1 pLGL =11 100 Fig. 7. Uniform grid convergence study for the viscous shock test case. Fully staggered SSDC algorithm (solid lines) versus a conventional LGL SSDC algorithm [8, 39] (dashed lines). The red lines show the (p + 1) slopes. • Fully staggered SSDC algorithm with pLG = 1, 2, 3, 4, 5, 10 and pLGL = 2, 3, 4, 5, 6, 11. • Conventional LGL SSDC algorithm [8, 39] with pLGL = 1, 2, 3, 4, 5, 10. 7.2.1. Uniform Cartesian grids. Different grid resolutions are examined, and the viscous shock has been propagated halfway out of the domain when the error measure is evaluated. This measures the effect of the penalty boundary condition and the interior scheme. Figure 7 shows the convergence study for a sequence of maximum nine nested grids (2 × 2 × 2, 4 × 4 × 2, . . . 512 × 512 × 2).19 The L2 norms of the error decay asymptotes towards the designed rate in each case, for both fully staggered and conventional LGL SSDC schemes. As for the Euler vortex test case, we note that the fully staggered approach is more accurate than the conventional algorithm for the same grid resolution. However, simple theoretical counting arguments based on inviscid flux and viscous flux evaluations (see section 6), seem to indicate that the cost of the staggered algorithm for a solution polynomial order p is comparable to that of an LGL operator [8, 39] with a solution polynomial order of (p + 1). 7.2.2. Unstructured nonuniform grids. A family of nested grids is constructed for the grid convergence study, by replicating N times the “unstructured grid kernel” shown in Figure 8. Figure 9 shows the convergence study for the family of maximum sixteen nested grids. As was the case in the previous study performed with the inviscid propagation of the Euler vortex, the fully staggered SSDC approach is more accurate than the conventional SSDC algorithm for the same grid resolution. Furthermore, for the fully staggered algorithm, the L2 norms of the error decay asymptotically towards the designed rate (i.e., (p + 1)) whereas the conventional path converges to order p. 19 The number of grid cells is doubled every time in each coordinate direction. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. A3153 Fig. 8. Unstructured grid kernel used to construct a sequence of nested grids for the viscous shock test case. 10-1 10-2 10-3 10-4 10-5 10-6 10-7 10-8 10-9 10-10 10-11 10-12 10-13 Grid convergence study: Viscous shock on unstructured nonuniform grids Fully-staggered, Fully-staggered, Fully-staggered, Fully-staggered, Fully-staggered, Fully-staggered, Conventional, Conventional, Conventional, Conventional, Conventional, Conventional, L2 Error ( -norm) Downloaded 11/20/16 to 109.171.137.210. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php ENTROPY STABLE STAGGERED GRID DISCONTINUOUS METHOD 10-2 10-1 1/(DOFs) 1/3 pLG =1, pLGL =2 pLG =2, pLGL =3 pLG =3, pLGL =4 pLG =4, pLGL =5 pLG =5, pLGL =6 pLG =10, pLGL =1 pLGL =2 pLGL =3 pLGL =4 pLGL =5 pLGL =10 pLGL =11 100 Fig. 9. Nonuniform grid convergence study for the viscous shock test case. Fully staggered SSDC algorithm (solid lines) vs. conventional LGL SSDC algorithm [8, 39] (dashed lines). The red lines show the (p + 1) slopes. 7.3. Taylor–Green vortex. The goal of this section is to demonstrate that nonlinearly stable schemes do not require stabilization techniques (e.g., dealiasing, filtering, limiting, overintegration procedures) for successful computations of underresolved turbulent flows. The Taylor–Green vortex test case is used as a benchmark model problem, and is solved on the periodic cube [−πL ≤ x, y, z ≤ +πL]. The initial condition is given by the following analytical expressions, x x x 2 3 1 (7.4) cos cos , u1 = V0 sin L L L x x x 1 2 3 u2 = −V0 cos sin cos , L L L u3 = 0, ρ0 V02 2x1 2x2 2x3 p = p0 + cos + cos cos +2 , 16 L L L Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. A3154 PARSANI, CARPENTER, FISHER, AND NIELSEN Re , =800 M =0.08 0.012 0.01 0.008 dke/dt Downloaded 11/20/16 to 109.171.137.210. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php Taylor-Green vortex, 0.006 Fully-staggered, Fully-staggered, DNS, Brachet et al. 0.004 pLG pLG 0.002 0.0 0 2 4 6 Time 8 , =16, =15 pLGL pLGL , grid: 43 =17, grid: 43 =16 10 Fig. 10. Evolution of the time derivative of the kinetic energy for the Taylor-Green vortex at Re = 800, M = 0.08; fully-staggered SSDC algorithm. where u1 , u2 , u3 are the components of the velocity in the x1 -, x2 -, and x3 -directions, p is the pressure, and ρ is the density. The flow is initialized to be isothermal, i.e., p/ρ = p0 /ρ0 = RT0 . The Reynolds number for this flow is defined as Re = (ρ0 V0 L)/µ, where µ is the dynamic viscosity. Starting from the initial condition, the nonlinear interactions of different flow scales yield vortex breakdowns. This nonlinear process is initially anisotropic and laminar, but subsequently, it develops into fully anisotropic turbulence that decays with the typical spectral energy distribution. Because our framework solves the compressible Navier–Stokes equations, we choose a Mach number of M = 0.08, which leads to a flow field that is essentially incompressible. This allows for a fair comparison with some of the incompressible simulations reported in literature. Herein, two Reynolds numbers are considered: Re = 800 and Re = 1, 600. The Prandtl number is set to P r = 0.71. 7.3.1. Re = 800. In this section we report the results for Re = 800 on a uniform Cartesian grid with four hexahedrons in each coordinate direction. Therefore, the total number of elements in the grid is also N e = 43 = 64. We compute the numerical solution with the sixteenth- (pLG = 15, pLGL = 16), and seventeenth-order (pLG = 16, pLGL = 17) accurate fully staggered algorithms. The numerical solutions presented herein are compared with the direct numerical simulation (DNS) of Brachet et al. [5]. Figure 10 shows the kinetic energy dissipation rate of our computations and the reference data. The computation with a formally seventeenth-order accurate scheme compares very well with the DNS results, even on this very relatively coarse grid. 7.3.2. Re = 1, 600. The simulation is run using a fully unstructured grid which contains 42 hexahedrons. The distribution of the element in the plane is shown in Figure 11. Figure 12 shows the kinetic energy dissipation rate of our computations with pLG = 15 and pLGL = 16, and pLG = 25 and pLGL = 26, and the reference data of Carton de Wiart et al. [15]. The DOFs for the computation with pLG = 15 and pLGL = 16 are too few to accurately resolve the flow field, however, the computation with a formally seventeenth-order accurate scheme is stable throughout all the Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. A3155 Fig. 11. Plane distribution of the elements used for the Taylor-Green vortex at Re = 1, 600, M = 0.08. Taylor-Green vortex, 0.014 Re 0.012 0.01 dke/dt Downloaded 11/20/16 to 109.171.137.210. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php ENTROPY STABLE STAGGERED GRID DISCONTINUOUS METHOD , =0 08 Fully-staggered, =15, Fully-staggered, =25, DNS, Carton de Wiart et al. =1600 M . , , pLG pLGL =16 42 pLG pLGL =26 42 cells cells 0.008 0.006 0.004 0.002 0.0 0 5 10 Time 15 20 Fig. 12. Evolution of the time derivative of the kinetic energy for the Taylor–Green vortex at Re = 1, 600, M = 0.08; fully staggered SSDC algorithm. simulation. This is a feat unattainable with alternative approaches based on highorder accurate linear stable schemes without filtering, dealiasing, limiters, etc. [6, 22]. The robustness and accuracy of the SS fully staggered algorithm is attained also for higher-order discretization as indicated in Figure 12 for the solution computed with pLG = 25 and pLGL = 26. In this case, the kinetic energy dissipation rate computed with a formally twenty-seventh-order accurate scheme compares very well with the DNS results. However, the dissipation rate typically converges relatively fast. Thus, to assess if the numerical solution is effectively converged it is good practice to compare enstrophy [15] as reported in Figure 13. At this resolution the enstrophy is not converged yet (spectral methods require a resolution of about 2563 to converge at this Reynolds number). Thus, a finer grid or a larger polynomial order is required. 7.4. Sod’s shock tube. Sod’s shock tube problem is a classical Riemannian problem that evaluates the behavior of a numerical method when a shock, expansion, and contact discontinuity are present. Of particular interest is smearing in the shock and contact, or oscillations at any of the discontinuities. The governing equations are the time-dependent 1D Euler equations. Sod’s problem is initialized with Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. A3156 PARSANI, CARPENTER, FISHER, AND NIELSEN Re , =0 08 Fully-staggered, =25, DNS, Carton de Wiart et al. =1600 10 M . pLG pLGL , =26 42 cells 8 Enstrophy Downloaded 11/20/16 to 109.171.137.210. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php Taylor-Green vortex, 6 4 2 0 0 5 10 Time 15 20 Fig. 13. Evolution of the enstrophy for the Taylor–Green vortex at Re = 1, 600, M = 0.08; fully staggered SSDC algorithm. (7.5) ρ(x) = 1 1/8 x < 0.5, x ≥ 0.5, u(x) = 0, p(x) = x ∈ (0, 1), 1, 1/10, t≥0 . x < 0.5, x ≥ 0.5, All simulations use the value γ = 7/5 for the ratio of specific heats. The governing equations are integrated in time up to t = 0.2. Neither the corrected nor the uncorrected schemes experienced instability on this test problem for elements of polynomial order 1 ≤ pLG ≤ 9. Figures 14 and 15 show the density distribution computed with the fully staggered SSDC algorithm with pLG = 3, pLGL = 4 and pLG = 9, pLGL = 10, respectively. Resolution of the expansion and contact discontinuity is adequate with either approach. Both the provably SS and baseline algorithms exhibit oscillations at the shock, although they capture the location correctly. Note, however, that when the entropy flux comparison is activated, the amplitudes of the oscillations are significantly reduced. The only mechanisms for dissipation in the fully staggered SSDC operators scheme are (1) boundary conditions, (2) element interface upwinding, (3) pointwise entropy correction of the fluxes; this is a remarkably low level of dissipation. 7.5. Sine-shock interaction. Numerical results are presented for the shock entropy-wave interaction problem. The solution of this benchmark problem contains both strong discontinuities and smooth structures and is well suited for testing highorder shock-capturing schemes. The governing equations are the time-dependent 1D Euler equations subject to the following initial conditions: ( (3.857143, 2.629369, 10.33333) if 0 ≤ x < 4.5, (7.6) (ρ, u, p) = (1 + 0.2 sin 5x, 0, 1) if 4.5 ≤ x ≤ 9. The governing equations are integrated in time up to t = 1.8. The exact solution to this problem is not available. Therefore, a numerical solution that is obtained using the conventional SSDC LGL algorithm [8, 39] with pLGL = 3 on a uniform grid with 1024 elements is used as the reference solution. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. ENTROPY STABLE STAGGERED GRID DISCONTINUOUS METHOD A3157 1.0 Exact Fully-staggered, pLG =3, pLGL =4, without comparison Fully-staggered, pLG =3, pLGL =4, with comparison Density 0.8 0.6 0.4 0.2 0.0 0.2 0.4 x 0.6 0.8 1.0 Fig. 14. Density profiles for Sod’s shock tube problem at t = 0.2. Fully-staggered SSDC algorithm with pLG = 3, pLGL = 4 on 128 uniform elements. Line markers are plotted every 0 solution LG points Sod problem 1.0 Exact Fully-staggered, pLG =9, pLGL =10, without comparison Fully-staggered, pLG =9, pLGL =10, with comparison 0.8 Density Downloaded 11/20/16 to 109.171.137.210. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php Sod problem 0.6 0.4 0.2 0.0 0.2 0.4 x 0.6 0.8 1.0 Fig. 15. Density profiles for Sod’s shock tube problem at t = 0.2. Fully staggered SSDC approach with pLG = 9, pLGL = 10 on 32 uniform elements. Line markers are plotted every fifteen solution LG points. Figures 16 and 17 show the density distribution computed with the fully staggered SSDC algorithm with pLG = 3, pLGL = 4 and pLG = 9, pLGL = 10, respectively. The algorithmic observations made in Sod’s problem are equally valid for the Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. Downloaded 11/20/16 to 109.171.137.210. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php A3158 PARSANI, CARPENTER, FISHER, AND NIELSEN Fig. 16. Density profiles for sine-shock interaction problem at t = 1.8. Fully staggered SSDC algorithm with pLG = 3, pLGL = 4 on 128 uniform elements. Line markers are plotted every fifteen solution LG points. Fig. 17. Density profiles for sine-shock interaction problem at t = 1.8. Fully staggered SSDC approach with pLG = 9, pLGL = 10 on 32 uniform elements. Line markers are plotted every fifteen solution LG points. sine-shock interaction problem. The coarse grid pLG = 9, pLGL = 10 simulation shows inadequate resolution as well as oscillations in the vicinity of the shock. None of the simulations exhibited numerical instability triggered by negative pressures. 7.6. Computation of a square cylinder in supersonic freestream. The development of a high-order accurate SS discretization aims to provide the next generation of robust high-fidelity numerical solvers for complex fluid flow simulations, for Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. Downloaded 11/20/16 to 109.171.137.210. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php ENTROPY STABLE STAGGERED GRID DISCONTINUOUS METHOD A3159 which standard suboptimal algorithms suffer greatly or fail completely. By computing the flow past a 3D square cylinder at Re∞ = 104 and M∞ = 1.5, we provide numerical evidence of such robustness for the complete SS high-order spatial discretization. This supersonic flow is characterized by a very large range of length scales, strong shocks, and expansion regions that interact with each other, leading to complex flow patterns. During the past three decades, this fluid flow problem has been thoroughly investigated by several researchers for aerodynamic applications (see for instance, [35, 4]). The domain of interest spans one square cylinder edge in the x3 direction, and at the two planes perpendicular to this coordinate direction, periodic boundary conditions are used. The flow is computed using an unstructured grid with 43, 936 hexahedrons. A fourth-order accurate (pLG = 3, pLGL = 4) fully staggered SS discretization without any stabilization technique is used. The body surface is considered adiabatic, and the wall BCs are imposed using the SS approach presented in [39]. The solution is initialized using a uniform flow at M∞ = 1.5 with zero angle of attack. This test case was also used in reference [39] to assess the robustness of the conventional SSDC algorithm. A strong shock is formed in front of the bluff body at the beginning of the simulation. Subsequently, the discontinuity moves upstream until it reaches a “stationary” position that is about 2.15 square cylinder edges from the frontal surface of the body. During this phase, additional weaker shocks, which originate from the four sharp corners of the body, interact with the subsonic regions formed near the walls. This complicated flow pattern yields the formation of shocklets in the wake of the square cylinder. A global view of the “high-order grid,” the Mach number, density, temperature, and entropy contours at t = 100 are shown in Figure 18. The shock has already reached a stationary position at t = 100, and the flow past the square cylinder is completely unsteady, characterized by subsonic and supersonic regions. The formation of shocklets in the near wake region are clearly visible. An extensive parametric study of Mach numbers (1.1 < M < 1.8) and grids (10K − 40K elements) was performed using the supersonic square cylinder. All high Mach number simulations performed to date suggest that the staggered and conventional approaches have comparable robustness. The failure mode for both approaches is a negative density that develops in the vicinity of a strong shock. 8. Conclusions. An SBP-SAT framework is used to develop a generalized SS spectral element formulation that includes a broader selection of collocation points. A necessary condition for an SS, staggered operator is a restriction/prolongation interpolation pair that satisfies a precise Galerkin constraint (see (2.18)). This constraint is automatically satisfied when interpolating between the LG and LGL points, provided the polynomial order of the LGL points exceeds that of the LG points (as well as other less desirable combinations). Thus, the primary goals herein are to (1) extend the SS mechanics to a staggered configuration that collocates the solution variables at the LG points and the fluxes at the LGL points, and (2) begin to investigate the competitive advantages (if any) of the new SS staggered algorithm relative to the algorithms presented in [8, 39]. The entropy analysis techniques presented in [8, 39] are used to develop an entropy conservative, staggered grid, spectral collocation operator for the 3D compressible Navier–Stokes equation. Entropy conservative/stable inviscid interface operators are then developed for the staggered formulation. The viscous interface coupling is based Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. Downloaded 11/20/16 to 109.171.137.210. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php A3160 PARSANI, CARPENTER, FISHER, AND NIELSEN Mach 1 0 2 2.42 Fig. 18. Unsteady flow past a 3D square cylinder at Re∞ = 104 and M∞ = 1.5; fourth-order accurate fully staggered SSDC method (pLG = 3, pLGL = 4) without stabilization technique; t = 100. on an LDG/IP approach analogous to the coupling operators presented in [8, 39]. The resulting spectral collocation operators are design-order accurate for arbitrary order, conservative on the LGL points, and satisfy the additional secondary constraint of entropy stability. Extension to curvilinear meshes [17, 20] and SS solid wall BCs [37, 39] follow immediately. The new operators are tested on both smooth and discontinuous test problems. The 1D Sod’s and sine-shock interaction discontinuous test problems are simulated with and without entropy corrections. It is shown that the base entropy conservative algorithm with upwinded interface fluxes is remarkably robust despite large oscillations at flow discontinuities. Remarkably little dissipation is needed to achieve satisfactory solutions for strong shocks. Extensive numerical tests for the 3D Euler and compressible Navier–Stokes equations reveal that the new staggered grid SS algorithms are significantly more accurate than the collocated LGL operators of equivalent polynomial order, that are presented in [8, 39]. They are however, more costly to implement from a theoretical point of view. The cost of a flux evaluation for the staggered algorithm of solution polynomial order p is comparable to that of an LGL points operator [8, 39] of solution polynomial order (p + 1). Preliminary studies using an explicit temporal integrator indicate that the increased accuracy of the staggered approach nearly offsets the additional cost. Further, thorough investigation is required to fully establish the cost efficacy of overcollocated fluxes (including the SD and the FR operators), relative to conventional collocated fluxes. Such an analysis will undoubtedly require a deep code optimization for current and emerging data-centric computing architectures. An ongoing investigation continues and includes (1) the effects of implicit temporal integrators and (2) data movement on current and next generation hardware; computationally intensive yet extremely accurate, low memory footprint algorithms could be competitive in the future. This work provides an essential step towards an operational SS framework of arbitrary order. An obvious extension of this work is the development of SS, h- and/or p-refinement operators. Both scenarios require data interpolation from adjoining interfaces onto a common intermediate mortar, with quadrature points that do not in general coincide. A long-term goal includes SS operators for other element types including triangles, prisms, and tetrahedra. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. Downloaded 11/20/16 to 109.171.137.210. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php ENTROPY STABLE STAGGERED GRID DISCONTINUOUS METHOD A3161 Acknowledgments. For computer time, this research partially used the resources of the Extreme Computing Research Center at King Abdullah University of Science & Technology. Special thanks are extended to Dr. Mujeeb R. Malik for partially funding this work as part of NASA’s Transformational Tools and Technologies (T 3 ) project. REFERENCES [1] H. L. Atkins and C.-W. Shu, Quadrature-free implementation of discontinuous Galerkin method for hyperbolic equations, AIAA J., 36 (1998), pp. 775–782. [2] C. Bernardi and Y. Maday, A collocation method over staggered grids for the Stokes problem, Internat. J. Numer. Methods Fluids, 8 (1988), pp. 537–557. [3] C. Bernardi, Y. Maday, and A. T. Patera, A new nonconforming approach to domain decomposition: The mortar element method, vol., in Nonlinear Partial Differential Equations and their Applications, Collège de France Seminar, Vol. XI Pitman Res. Notes Math. Ser., Longman Scientific and Technical, Harlow, England 299 (1994), pp. 13–51. [4] T. Birch, S. A. Prince, and G. M. Simpson, An Experimental and Computational Study of the Aerodynamics of a Square Cross-Section Body at Supersonic Speeds, Technical report RTO-MP-069(I), DERA, United Kingdom, 2003. [5] M. E. Brachet, D. Meiron, S. Orszag, B. Nickel, R. Morf, and U. Frisch, The Taylor– Green vortex and fully developed turbulence, J. Stat. Phys., 34 (1984), pp. 1049–1063. [6] J. R. Bull and A. Jameson, Simulation of the compressible Taylor-Green vortex using highorder flux reconstruction schemes, in 7th AIAA Theoretical Fluid Mechanics Conference, AIAA 2014-3210, Curran, Red Hook, NY, 2014. [7] M. H. Carpenter and T. C. Fisher, Entropy Stable Spectral Collocation Schemes for the Navier-Stokes Equations: Discontinuous Interfaces, Technical report NASA TM 218039, 2013. [8] M. H. Carpenter, T. C. Fisher, E. J. Nielsen, and S. H. Frankel, Entropy stable spectral collocation schemes for the Navier–Stokes equations: Discontinuous interfaces, SIAM Journal on Scientific Computing, 36 (2014), pp. B835–B867. [9] M. H. Carpenter and D. Gottlieb, Spectral methods on arbitrary grids, J. Comput. Phys., 129 (1996), pp. 74–86. [10] M. H. Carpenter, D. Gottlieb, and S. Abarbanel, Time-stable boundary conditions for finite-difference schemes solving hyperbolic systems: Methodology and application to highorder compact schemes, J. Comput. Phys., 111 (1994), pp. 220–236. [11] M. H. Carpenter, M. Parsani, T. C. Fischer, and E. J. Nielsen, Towards an entropy stable spectral element framework for computational fluid dynamics, in 54th AIAA Aerospace Sciences Meeting, AIAA 2016-1058, Curran, Red Hook, NY, 2016. [12] M. H. Carpenter, M. Parsani, T. C. Fisher, and E. J. Nielsen, Entropy Stable Staggered Grid Spectral Collocation for the Burgers’ and Compressible Navier–Stokes Equations, Technical report NASA TM 218990, 2015. [13] M. H. Carpenter, M. Parsani, T. C. Fischer, and E. J. Nielsen, Towards an entropy stable spectral element framework for computational fluid dynamics, in 54th AIAA Aerospace Sciences Meeting, AIAA 2016-1058, Curran, Red Hook, NY, 2016. [14] C. M. Dafermos, Hyperbolic Conservation Laws in Continuum Physics, Springer-Verlag, Berlin, 2010. [15] C. de Wiart, K. Hillewaert, M. Duponcheel, and G. Winckelmans, Assessment of a discontinuous Galerkin method for the simulation of vortical flows at high Reynolds number, Internat. J. Numer. Methods Fluids, 74 (2014), pp. 469–493. [16] P. F. Fischer and E. M. Rønquist, Spectral element methods for large scale parallel NavierStokes calculations, Comput. Methods Appl. Mech. Engrg, 116 (1994), pp. 69–76. [17] T. C. Fisher, High-order L2 Stable Multi-Domain Finite Difference Method for Compressible Flows, PhD thesis, Purdue University, West Lafagette, IN, 2012. [18] T. C. Fisher and M. H. Carpenter, High-order entropy stable finite difference schemes for nonlinear conservation laws: Finite domains, J. Comput. Phys., 252 (2013), pp. 518–557. [19] T. C. Fisher, M. H. Carpenter, J. Nordström, N. K. Yamaleev, and R. C. Swanson, Discretely Conservative Finite-Difference Formulations for Nonlinear Conservation Laws in Split Form: Theory and Boundary Conditions, Technical report NASA TM 2011-217307, 2011. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. Downloaded 11/20/16 to 109.171.137.210. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php A3162 PARSANI, CARPENTER, FISHER, AND NIELSEN [20] T. C. Fisher, M. H. Carpenter, J. Nordström, N. K. Yamaleev, and R. C. Swanson, Discretely conservative finite-difference formulations for nonlinear conservation laws in split form: Theory and boundary conditions, J. Comput. Phys., 234 (2013), pp. 353–375. [21] G. J. Gassner, A kinetic energy preserving nodal discontinuous Galerkin spectral element method, Internat. J. Numer. Methods Fluids, 76 (2014), pp. 28–50. [22] G. J. Gassner and A. D. Beck, On the accuracy of high-order discretizations for underresolved turbulence simulations, Theoretical and Computational Fluid Dynamics, 27 (2012), pp. 221–237. [23] J. S. Hesthaven and T. Warburton, Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications, Texts Appl. Math. 54 New York, Springer, 2008. [24] M. Hutchinson, A. Heinecke, H. Pabst, G. Henry, M. Parsani, and D. Keyes, Efficiency of high order spectral element methods on petascale architectures, in Proceedings of the International Supercomputing Conference, ISC’16, 2016. [25] H. T. Huynh, A flux reconstruction approach to high-order schemes including discontinuous Galerkin methods, in 18th AIAA Computational Fluid Dynamics Conference, AIAA 20074079, AIAA, Recton, VA, 2007. [26] F. Ismail and P. L. Roe, Affordable, entropy-consistent Euler flux functions II: Entropy production at shocks, J. Comput. Phys., 228 (2009), pp. 5410–5436. [27] D. A. Kopriva and G. J. Gassner, On the quadrature and weak form choices in collocation type discontinuous Galerkin spectral element methods, J. Sci. Comput., 44 (2010), pp. 136–155. [28] D. A. Kopriva and J. H. Kolias, A conservative staggered-grid Chebyshev multidomain method for compressible flows, J. Comput. Phys., 125 (1996), pp. 244–261. [29] P. Lax and B. Wendroff, Systems of conservation laws, Comm. Pure Appl. Math., 13 (1960), pp. 217–237. [30] P. D. Lax, Hyperbolic systems of conservation laws and the mathematical theory of shock waves, CBMS-NSF Regional Conf, Ser. in Appl. Math 11, SIAM, Philadelphia, 1973. [31] Y. Liu, M. Vinokur, and Z. J. Wang, Spectral difference method for unstructured grids I: Basic formulation, J. Comput. Phys., 216 (2006), pp. 780–801. [32] Y. Maday, A. T. Patera, and E. M. Ronquist, The pn pn2 Method for the Approximation of the Stokes Problem, Technical report 4, Laboratoire d’Analyse Numerique, Paris, 1992. [33] G. May and A. Jameson, A spectral difference method for the Euler and Navier–Stokes equations on unstructured meshes, in 44th AIAA Aerospace Sciences Meeting and Exhibit, AIAA 2006-304, AIAA, Reston, VA, 2006. [34] M. L. Merriam, An Entropy-Based Approach to Nonlinear Stability, Technical report NASA TM 101086, 1989. [35] T. Nakagawa, Vortex shedding behind a square cylinder in transonic flows, J. Fluid Mech., 178 (1987), pp. 303–323. [36] S. Osher, Riemann solvers, the entropy condition, and difference, SIAM J. Numer. Anal., 21 (1984), pp. 217–235. [37] M. Parsani, M. H. Carpenter, and E. J. Nielsen, Entropy Stable Wall Boundary Conditions for the Compressible Navier–Stokes Equations, Technical report NASA TM 218282, 2014. [38] M. Parsani, M. H. Carpenter, and E. J. Nielsen, Entropy stable discontinuous interfaces coupling for the three-dimensional compressible Navier–Stokes equations, J. Comput. Phys., 290 (2015), pp. 132–138. [39] M. Parsani, M. H. Carpenter, and E. J. Nielsen, Entropy stable wall boundary conditions for the three-dimensional compressible Navier–Stokes equations, J. Comput. Phys., 292 (2015), pp. 88–113. [40] M. Svärd and S. Mishra, Entropy stable schemes for initial-boundary-value conservation laws, Z. Angew. Math. Phys., 63 (2012), pp. 985–1003. [41] M. Svärd and H. Özcan, Entropy-stable schemes for the Euler equations with far-field and wall boundary conditions, J. Sci. Comput., 58 (2014), pp. 61–89. [42] M. Svärd, Weak solutions and convergent numerical schemes of modified compressible Navier– Stokes equations, J. Comput. Phys., 288 (2015), pp. 19–51. [43] E. Tadmor, Entropy stability theory for difference approximations of nonlinear conservation laws and related time-dependent problems, Acta Numer., 12 (2003), pp. 451–512. [44] E. Tadmor and W. Zhong, Entropy stable approximations of the Navier–Stokes equations with no artificial numerical viscosity, J. Hyperbolic Differ. Equ., 3 (2006), pp. 529–559. [45] P. E. Vincent, P. Castonguay, and A. Jameson, A new class of high-order energy stable flux reconstruction schemes, J. Sci. Comput., 47 (2011), pp. 50–72. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

RELATED PAPERS

RELATED TOPICS

Log In

Entropy Stable Staggered Grid Discontinuous Spectral Collocation Methods of any Order for the Compressible Navier--Stokes Equations

Entropy Stable Staggered Grid Discontinuous Spectral Collocation Methods of any Order for the Compressible Navier--Stokes Equations

Related Papers

RELATED PAPERS

RELATED TOPICS