A hybrid framework for coupling arbitrary
summation-by-parts schemes on general meshes
Tomas Lundquist, Arnaud Malan and Jan Nordström
The self-archived postprint version of this journal article is available at Linköping University Institutional Repository (DiVA):
http://urn.kb.se/resolve?urn=urn:nbn:se:liu:diva-145303
N.B.: When citing this work, cite the original publication.
Lundquist, T., Malan, A., Nordström, J., (2018), A hybrid framework for coupling arbitrary summation-by-parts schemes on general meshes, Journal of Computational Physics, 362, 49-68. https://doi.org/10.1016/j.jcp.2018.02.018
Original publication available at:
https://doi.org/10.1016/j.jcp.2018.02.018
Copyright: Elsevier
A hybrid framework for coupling arbitrary
summation-by-parts schemes on general meshes
Tomas Lundquista, Arnaud Malanb, Jan Nordstr¨omc
aIndustrial CFD Research Group, Department of Mechanical Engineering, University of
Cape Town, South Africa (tomas.lundquist@uct.ac.za).
bIndustrial CFD Research Group, Department of Mechanical Engineering, University of
Cape Town, South Africa (arnaud.malan@uct.ac.za).
cDepartment of Mathematics, Computational Mathematics, Link¨oping University,
Sweden (jan.nordstrom@liu.se).
Abstract
We develop a general interface procedure to couple both structured and un-structured parts of a hybrid mesh in a non-collocated, multi-block fashion. The target is to gain optimal computational efficiency in fluid dynamics sim-ulations involving complex geometries. While guaranteeing stability, the pro-posed procedure is optimized for accuracy and requires minimal algorithmic modifications to already existing schemes. Initial numerical investigations confirm considerable efficiency gains compared to non-hybrid calculations of up to an order of magnitude.
Keywords: Hybrid meshes, Numerical interfaces, High order methods, Stability, Summation-by-parts, Weak boundary/interface conditions, Finite volume, Finite difference
1. Introduction
Computational fluid dynamics (CFD) has become a key technology in en-gineering innovation and science investigation. Despite significant advances, CFD as a design tool is still expensive when complex geometries are preva-lent. In these circumstances, the optimal method of discretization as well as the optimal level of grid refinement often varies throughout the computa-tional domain. Hence, the ability to combine general computacomputa-tional blocks in a modular fasion can be a powerful way of improving the efficiency in many types of simulations. For example, finite volume (FV) methods are still the
most widely used in CFD due their applicability to unstructured meshes, even though at most second order spatial accuracy is routinely employed. Finite difference (FD) methods on the other hand represent the most efficient type of high order schemes, but their use is limited to structured meshes, and thus less complex geometries. The best attributes of each of these methods may however be harnessed in an innovative manner by developing a hybrid scheme. Such an approach involves utilizing both methods on a particular problem so as to ensure optimal computational efficiency. The key challenge to this involves interfacing the two schemes in a stable and efficient manner, with minimal modifications to either method.
A systematic approach to the formulation of energy stable schemes is provided by the concept of summation-by-parts (SBP) operators, originally introduced in [1], together with weak boundary [2] and interface conditions [3]. The summation-by-parts framework was for a long time restricted to high order accurate FD methods, gradually incorporating techniques to handle curvilinear and multi-block grids [4, 5, 6] in a stable manner. The initial hybrid development of interest to this paper began with the formulation a node-centered edge-based FV scheme in SBP form [7, 8, 9]. This was extended to hybrid FD-FV couplings in [10, 11, 12]. The stability of this particular hybrid technique relies on a collocated FD-FV interface where the dual grid cells on the FV side are modified to match the integration weights and metric variables of the FD scheme.
Apart from FD and FV schemes, significant steps have also been taken towards unifying a much broader class of nodal methods with similar stability properties into the SBP framework. High order methods in this category include spectral [13, 14], discontinuous Galerkin (dG) [15], flux reconstruction [16, 17], time integration [18, 19, 20] and space-time [21, 22] methods. A recent related development to the hybrid FD-FV method described above is the use of so-called summation-by-parts preserving interpolation operators to couple non-collocated numerical domains, i.e. interfaces where the node distributions are non-matching, see [23, 24, 25, 26]. So far this approach has been applied to non-collocated FD-FD as well as FD-dG [25] couplings, where different grid resolutions are used on each side of an interface.
In this paper we extend the concept of summation-by-parts preserving operators and formulate a hybrid framework for general SBP methods posed on arbitrary grids. We use the term general here in the sense that the formu-lation does not rely on the explicit application of either tensor products or curvilinear coordinate transformations. As such, the new framework is well
suited for both structured and unstructured methods. To exemplify the new framework, we develop interface operators between FD blocks of arbitrary order involving curvilinear transformation, as well as between FD and FV blocks. Both the FD-FD and FD-FV couplings are combined with interpola-tion to allow for non-collocated node distribuinterpola-tions, i.e. so that different grid resolutions can be used on both sides of a hybrid interface.
The rest of this article is organized as follows. In sections 2 and 3 we lay out a general theoretical framework for energy stable interface couplings between SBP schemes on arbitrary geometries. In section 4 we confirm that both the FV and high order FD methods conform to the general SBP frame-work. Once the theoretical foundation is in place, we proceed in section 5 to derive new FD-FD as well as FD-FV interface operators for curvilinear grids, including the case of non-collocated grid distributions. In section 6 we perform numerical calculations and study the convergence properties and efficiency of the new hybrid schemes. Finally, in section 7, we draw conclu-sions.
2. Summation-by-parts operators
The summation-by-parts methodology may be understood as a systematic approach to construct stable numerical schemes for linear problems. Stability for non-linear problems with smooth solutions then follows by application of the energy method to all frozen coefficient versions of the problem at hand, see e.g. [27] for a review of this technique. In the remainder, we will assume that the reader is familiar with the energy method for analyzing well-posedness and stability.
2.1. Continuous analysis
For simplicity of presentation, the analysis in this paper will be carried out for the advection-diffusion problem,
ut+ aT∇u = ∇T∇u, x ∈ Ω ⊂ Rr, t > 0 1 2(a Tn − |aTn|)u − nT∇u = g(t, x), x ∈ ∂Ω, t > 0 u = f (x), x ∈ Ω, t = 0 (1)
posed on an arbitrary domain x ∈ Ω ⊂ Rr in r space dimensions. In (1), a is a constant vector, u(x, t) is the solution and we use well-posed Robin
boundary conditions. Other forms of boundary conditions in (1) are also possible, including Dirichlet conditions [28]. However, since the main focus in this article is on interface conditions, it is sufficient for our needs to consider only one form of well-posed boundary conditions.
We introduce the continuous inner products for sufficiently smooth scalar functions φ and ψ, (φ, ψ)Ω = Z Ω φψ dV, (φ, ψ)∂Ω= I ∂Ω φψ dS, (∇φ, ∇ψ)Ω = Z Ω (∇φ)T(∇ψ) dV (2) where RΩ and H∂Ω are volume and surface integrals, respectively. The inte-gration rules required for analyzing (1) with the energy method are,
(φ, aT∇ψ)Ω = (φ, aTnψ)∂Ω− (aT∇φ, ψ)Ω (3)
(φ, ∇T∇ψ)Ω = (φ, nT∇ψ)∂Ω− (∇φ, ∇ψ)Ω (4)
where n is the outward pointing unit normal to the domain. With zero boundary data g = 0, the energy method applied to (1) yields the (bounded) energy rate
∂ ∂tkuk
2
Ω+ 2k∇uk2Ω = −(u, |aTn|u)∂Ω, (5)
which proves well-posedness. 2.2. Discrete integration operators
In order to enable the application of rules analogous to (3)-(4) on the dis-crete side, the following three components are central and should be present (either explicitly or implicitly) in any numerical scheme of SBP type:
• An integration operator P over the whole discrete domain. • An integration operator P over the domain boundary.
• An operator E that restricts solutions to the domain boundary. The integration operator P can be seen as analogous to the concept of a mass matrix from finite element methods. For example, in the case of node-centered FV methods it is defined by the control volumes, while high or-der FD methods in summation-by-parts form are associated with high oror-der Gregory-type quadrature rules [29, 30]. Both P and P define inner products (and hence norms) of discrete solutions, and we will often refer to these op-erators as the discrete ”norms”. In this paper we assume that the opop-erators introduced above satisfy the following two properties.
Property 1. The discrete norms P and P can be represented by diagonal matrices with positive integration weights positioned along the diagonal. Property 2. The boundary restriction operator E is exact in the sense that it contains a subset of rows from a unit matrix of the same dimension as P, corresponding to all nodes located on the domain boundary.
Remark 1. Both properties 1 and 2 are required for stability proofs involv-ing curvilinear coordinate transformations, variable coefficient and nonlinear problems [31, 32, 33].
It is often natural to define a splitting of the boundary operators into n disjoint parts, i.e.
E = E1 .. . En , P = P1 . .. Pn . (6)
Such splittings will e.g. occur automatically in the case of tensor product extensions of discrete operators in one space dimension.
We introduce the following inner products analogous to (2),
(Φ, Ψ)P = ΦTPΨ, (EΦ, EΨ)P = (EΦ)TP (EΨ), (DΦ, DΨ)Ir⊗P = (DΦ)
T (Ir⊗P)(DΨ). (7) Here, D = (DT x1, . . . , D T xr) T, where D
xi is an operator that approximates the
xi− derivative. Finally, Ir is a unit matrix of the same size as the number of
spatial dimensions, r.
2.3. Advection and diffusion operators in summation-by-parts form
We define an operator of the form P−1Q to approximate the advective differential operator aT∇ in (1), where P is the discrete integration operator
introduced in the previous section, and Q is essentially skew-symmetric (it roughly corresponds to the stiffness matrix in finite element methods). We assume the existence of a summation-by-parts rule corresponding to (3) in continuous space, and consider the following specific form,
Q = ETP ΛE − QT (8)
where Λ = Diag(aTn). An artificial dissipation term in the form of a postive
normal to Ω is not uniquely defined at sharp corners or edges, it is natural to split the boundary as in (6) at such points. In that case, the following rule applies,
ETP ΛE = E1TP1Λ1E1+ . . . + EnTPnΛnEn. (9)
Note that the same node may be involved in two or more of the partitions in (9), and hence the corresponding row may appear more than once in E. This is typically the case for nodes at corners and edges of the domain where a splitting like (6) and (9) occurs.
Similarly to (8), we define a diffusion operator P−1R approximating ∇T∇, where R satisfies the rule corresponding to (4),
R = ETP D − DT(Ir⊗ P)D. (10)
In (10), D is defined as in (7), and D = Pr
i=0Diag(ni)EDxi approximates
the normal derivative nT∇.
Note that (8) and (10) perfectly imitate the continuous integration-by-parts rules (3)-(4) Indeed, for any two vectors Φ and Ψ, the summation-by-parts rules (8) and (10) lead to
(Φ, P−1QΨ)P = (EΦ, ΛEΨ)P − (P−1QΦ, Ψ)P
(Φ, P−1RΨ)P = (EΦ, DΨ)P − (DΦ, DΨ)Ir⊗P.
We will hence refer to operators that satisfy (8) and (10) as SBP operators, and two specific examples will be given in section 4.
Consider a discretization of (1) using the SBP operators (8) and (10) Ut+ P−1QU = P−1RU + P−1Σ, (11)
where the boundary condition is imposed in weak form as
Σ = ETP 12(Λ − |Λ|)u − un− g(t, X) , (12)
and X denotes the nodal grid distribition at the boundary. With g = 0, the discrete energy method yields
∂ ∂tkUk
2
P + 2kDUk2Ir⊗P = − (u, |Λ|u)P. (13)
2.4. Summation-by-parts preserving interface operators
Let Ω be divided into two subdomains ΩLand ΩR, and let each domain be
discretized with a scheme in SBP form (8) and (10). To communicate nodal solution values between the two sides of the interface in a stable and accu-rate manner, an pair of interface operators will be required. The following definition is a direct extension of the concept of SBP preserving interpolation in [23] to general hybrid meshes.
Definition 1. Consider a pair of linear operators [IL2R, IR2L], where IL2R is
used to transfer nodal values from the ΩL side of the interface to the ΩR side,
and IR2L is used in the opposite direction. We say that [IL2R, IR2L] is a pair
of SBP preserving interface operators with respect to the discrete integration operators [PL, PR], if
PLIR2L= IL2RT PR. (14)
The significance of the SBP preserving property (14) will become ap-parent in section 3, where interface couplings involving such operators are analyzed with the energy method.
In many situations it is convenient to construct the interface operators in two steps with the help of an intermediate virtual grid ΩM, associated
with the norm PM. For example, by defining IL2R = IM 2RIL2M we indicate
that nodal values are first transfered from ΩL to ΩM using IL2M, and then
from ΩM to ΩR using IM 2R. This approach is similar to the so-called glue
grid method of [25] for FD-FD and FD-dG couplings, and uses the fact that the combination of two SBP preserving operations is in itself SBP preserving (see Lemma 3 in [25]). Formulated in general terms below, we will make extensive use of this result.
Lemma 1. Let [IL2M, IM 2L] and [IM 2R, IR2M] be two pairs of SBP
preserv-ing operators with respect to the discrete norms [PL, PM] and [PM, PR]
re-spectively, and define
IL2R = IM 2RIL2M, IR2L= IM 2LIR2M. (15)
Then the pair [IL2R, IR2L] is summation-by-parts preserving with respect to
[PL, PR].
Proof. We have PLIM 2L = IL2MT PM and PMIR2M = IM 2RT PR. Hence,
PLIR2L= PL(IM 2LIR2M) = IL2MT PMIR2M = IL2MT I T
M 2RPR= (IM 2RIL2M)TPR= IL2RT PR
3. Energy stable weak interface conditions
Consider a discretization of (1) on the two subdomains ΩL and ΩR with
boundary and interface conditions imposed in weak form, Ut+ PL−1QLU = PL−1RLU + PL−1(ΣBL + Σ I L+ Σ V L) Vt+ PR−1QLV = PR−1RRV + PR−1(Σ B R+ Σ I R+ Σ V R), (16) where ΣB
L,R, ΣIL,R and ΣVL,R are penalty terms that impose boundary
con-ditions as well as inviscid and viscous coupling concon-ditions. One attractive feature of the weak formulation (16) is that all three of these conditions can be analyzed independently from each other with the energy method. Hence we can ignore most of the details regarding the physical boundary conditions, and focus on the inviscid and viscous coupling conditions.
We apply (9) to separate out the contributions from the common interface to the SBP properties (8) and (10)
QL= ELTPLΛLEL+ (ELB) TPB LΛ B LE B L − Q T L QR= ERTPRΛRER+ (ERB) TPB RΛ B RE B R − Q T R RL= ELTPLDL+ (ELB) TPB LD B L − D T L(Ir⊗ PL)DL RR= ERTPRDR+ (ERB) TPB RD B R− D T R(Ir⊗ PR)DR, (17)
where the superscript B is used for terms associated with the outer bound-aries to Ω. Note that DL and DR are the interface parts of the normal
derivatives in (10). For future convenience we also define
u = ELU, v = ERV, un= DLU, vn = DRV.
3.1. The inviscid coupling
Similarly to (12), the role of the inviscid interface terms ΣILand ΣIRin (16) is to impose accuracy conditions involving nodal interface values weakly. It was shown in [15] that the weak interface conditions applied in classical SBP multi-block schemes are similar to the coupling conditions used in a standard nodal dG method. Following this observation, we consider a formulation based on numerical fluxes,
ΣIL = ELTPL(ΛLu − dΛLu), ΣIR = E T
RPR(ΛRv − dΛRv) (18)
where dΛLu and dΛLv are linear combinations of nodal interface values, to be
The formulation used in [23], where SBP preserving interpolation opera-tors for FD schemes were first introduced, did not include an option to use upwind terms in the coupling, which is generally desired. This issue was resolved in [25] by the explicit use of an intermediate so-called glue grid on which all upwind difference terms were calculated. By slightly modifying this technique, we consider the following ansatz for the numerical fluxes in (18),
d ΛLu =αLΛL u + IR2Lv 2 + (1 − αL) ΛLu + IR2L(−ΛRv) 2 + σ 2IM 2L|ΛM|(IL2Mu − IR2Mv) d ΛRv =αRΛR v + IL2Ru 2 + (1 − αR) ΛRv + IL2R(−ΛLu) 2 + σ 2IM 2R|ΛM|(IR2Mv − IL2Mu). (19)
Note that the central contributions above are expressed as weighted sums of two different forms of average terms. Moreover, the upwind differencing is carried out on an intermediate interface grid (or glue grid) of the same type as discussed in section 2.4. We note however that even if no such intermediate grid has been used to construct IL2Rand IR2L, one could simply define either
M = R or M = L in (19) and still retain the upwind terms (ILL, IRR would
then be identity matrices).
Remark 2. On a collocated interface where PR = PL and ΛR = −ΛL, all
interface operators in (19) can simply be replaced with unit matrices, since (14) then holds trivially. In this case the numerical fluxes reduce to the standard Riemann flux,
d ΛLu =ΛL u + v 2 + σ 2|ΛL|(u − v), ΛdRv = ΛR v + u 2 + σ 2|ΛR|(v − u). (20) With the fully upwind choice σ = 1 in (20), it is straightforward to show that (18) can be rewritten into
ΣIL= 1 2E T LPL(ΛL− |ΛL|)(u − v), ΣIR = 1 2E T RPR(ΛR− |ΛR|)(v − u) (21) i.e. the interface conditions are imposed only at the inflow parts of the re-spective subdomains.
Proposition 1. Let αR = 1 − αL and σ ≥ 0 in (19). Then the inviscid
interface treatment in (16) is energy stable, and the energy rate due to the interface becomes ∂ ∂t(kUk 2 PL+kVk 2
Note that the right hand side to (22) is a non-positive term which provides a small amount of damping to the system if σ > 0.
Proof. The energy method applied to (16) yields the following contribution from inviscid terms in (16) (ignoring all viscous terms as well as terms asso-ciated with the outer boundaries),
1 2 ∂ ∂t(kUk 2 PL+ kVk 2 PR) = U T(ΣI L− QLU) + VT(ΣIR− QRV) = UT(ΣIL−1 2(QL+ QL) TU) + VT(ΣI R− 1 2(QR+ Q T R)V) = UT(ΣIL−1 2E T LPLΛLELU) + VT(ΣIR− 1 2E T RPRΛRERV) = uT(PL(ΛLu − dΛLu) − 1 2PLΛLu) + v T (PR(ΛRv − dΛRv) − 1 2PRΛRv) = uTPL 1 2ΛLu − dΛLu + vTPR 1 2ΛRv − dΛRv . Inserting (19) and multiplying with 2 yields
∂ ∂t(kUk 2 PL+ kvk 2 PR) = − u v T
Acentral+ σAupwind
u v where Acentral = 0 αLPLΛLIR2L+ (αL− 1)PLIR2LΛR αRPRΛRIL2R+ (αR− 1)PRIL2RΛL 0 Aupwind = PLIM 2L|ΛM|IL2M −PLIM 2L|ΛM|IR2M −PRIM 2R|ΛM|IL2M PRIM 2R|ΛM|IR2M .
Now consider the following splitting of the symmetric part of Acentral,
Acentral+ ATcentral= 0 αLPLΛLIR2L+ (αR− 1)ΛLIL2RT PR αRPRΛRIL2R+ (αL− 1)ΛRIR2LT PL 0 + 0 αRIL2RT ΛRPR+ (αL− 1)PLIR2LΛR αLIR2LT ΛLPL+ (αR− 1)PRIL2RΛL 0 . Recall that PL and PR are diagonal (Property 1) and hence commute with
ΛL and ΛR, respectively. After applying the SBP preserving condition (14),
Next, by applying (14) to each of the blocks in Aupwind, we get Aupwind= IT L2MPM|ΛM|IL2M −IL2MT PM|ΛM|IR2M −IT R2MPM|ΛM|IL2M IR2MT PM|ΛM|IR2M =IL2M 0 0 IR2M T 1 −1 −1 1 ⊗ PM|ΛM| IL2M 0 0 IR2M . From this it immediately follows that the inviscid contribution from the in-terface is given by,
∂ ∂t(kUk 2 PL+ kVk 2 PR) = − σ IL2Mu IR2Mv T 1 −1 −1 1 ⊗ PM|ΛM| IL2Mu IR2Mv = − σ(IL2Mu − IR2Mv, |ΛM|(IL2Mu − IR2Mv))PM
which concludes the proof. 3.2. The viscous coupling
For the viscous interface condition in (16), we generalize the formulation used in [11, 23], as
ΣVL = γLELTPL(un− IR2L(−vn)) + δLDLTPL(u − IR2Lv)
ΣVR = γRERTPR(vn− IL2R(−un)) + δRDRTPR(v − IL2Ru).
(23) The following result is a direct generalization of the stability proof in [23]. Proposition 2. The weak viscous coupling defined by (23) is stable, with zero contribution to the system energy, if
γR = δL, δR = γL, γL+ δL = −1. (24)
Proof. The total contribution from the viscous interface terms in (17) and (23) to energy can be written as,
1 2 ∂ ∂t(kUk 2 PL+ kVk 2 PR) + (kDLUk 2 Ir⊗PL + kDRVk 2 Ir⊗PR) = UT(ELTPLDLU + ΣVL) + V T(ET RPRDRV + ΣVR) = uT(PLun+ γLPL(un− IR2L(−vn))) + δLuTnPL(u − IR2Lv) + vT(PRvn+ γRPR(vn− IL2R(−un))) + δRvTnPR(v − IL2Ru) = u v un vn T Aviscous u v un vn ,
where Aviscous = 0 0 (1 + γL)PL γLPLIR2L 0 0 γRPRIL2R (1 + γR)PR δLPL −δLPLIR2L 0 0 −δRPRIL2R δRPR 0 0 .
After applying (24) and the SBP preserving property (14), the symmetric part of Aviscous above vanishes.
We summarize this section by stating the total energy rate of the discrete system (16) resulting from boundary and interface terms. From propositions 1 and 2 and the analysis in section 2.3 we get, with g = 0,
∂ ∂t(kUk 2 PL+ kVk 2 PR) + 2(kDLUk 2 Ir⊗PL+ kDRVk 2 Ir⊗PR) = − (uB, |ΛBL|uB)PB L − (vB, |Λ B R|vB)PB R − σ(IL2Mu − IR2Mv, |ΛM|(IL2Mu − IR2Mv))PM,
where uB = ELBU and vB = ERBV. This energy rate is similar and
corre-sponds to the previous ones in (5) and (13).
4. Examples of schemes in generalized summation-by-parts form In this section we verify that both the FD and the FV method can be expressed using the generalized SBP operators defined in (8) and (10). Once this has been achieved, it only remains to derive SBP preserving interface operators between the two schemes in order to apply the general coupling technique described in the previous section. For simplicity of presentation we consider (1) in two space dimensions, but the extension to three dimensions is straightforward.
4.1. The finite difference method
Consider a curvilinear domain (x, y) = (x(ξ, η), y(ξ, η)) ∈ Ω, where the Cartesian reference grid is given by (ξ, η) ∈ ˆΩ = [0, 1]2, see Figure 1 for an illustration. Moreover, let
Dξ= Pξ−1Qξ⊗ Iη, Dη = Iξ⊗ Pη−1Qη
where Pξ−1Qξand Pη−1Qη are one-dimensional summation-by-parts operators
in the ξ- and η-directions, satisfying Property 2. In other words,
Qξ,η = (eξ,η)T1(eξ,η)1 − (eξ,η)T0(eξ,η)0− QTξ,η (25)
ˆ nS ˆ nW ˆ nN ˆ nE ˆ Ω nS nW nN nE Ω
Figure 1: Curvilinear grid transformation.
4.1.1. The metric terms
Recently, it has been shown that a splitting of the inviscid operator into conservative and non-conservative parts is required for stability proofs [31, 11, 21]. We transform the differential terms in (1) as (see e.g. [11]),
aT∇u = 1 2J
−1
(ˆaT∇u + ˆˆ ∇T(ˆau)), ∇T∇u = J−1∇ˆT(J C ˆ∇u) (26) where ˆ∇ = (∂ξ, ∂η)T, ˆa = (ˆa1, ˆa2)T, ˆa1 = a1yη− a2xη, ˆa2 = −a1yξ+ a2xξ, and
C = ξ2 x+ ξy2 ξxηx+ ξyηy ξxηx+ ξyηy η2x+ ηy2 =ξx ηx ξx ηx T +ξy ηy ξy ηy T . Note that we have split the inviscid operator into a conservative and a non-conservative part in (26), which is possible thanks to the geometric conserva-tion law (ˆa1)ξ+ (ˆa2)η = 0. For future reference we also state the well-known
metric relations,
J ξx = yη, J ηx = −yξ, J ξy = −xη, J ηy = xξ (27)
where J = xξyη− yξxη is the Jacobian determinant of the coordinate
trans-formation.
normal vectors to ˆΩ and Ω can be written as ˆ nN = 0 1 T , nN = JN−1 −yξ xξ T ˆ nE = 1 0 T , nE = JE−1 yη −xη T ˆ nS = 0 −1 T , nS = JS−1 yξ −xξ T ˆ nW = −1 0 T , nW = JW−1 −yη xη T (28) where JN = JS = q x2 ξ+ y 2 ξ and JE = JW =px2η+ yη2. 4.1.2. The discretization
Following [31, 11, 21], we apply one-dimensional first derivative operators in SBP form (25) to discretize the differential terms in (26) as
P = J (Pξ⊗ Pη) (29) P−1Q = 1 2J −1 [DξAˆ1+ ˆA1Dξ+ DηAˆ2+ ˆA2Dη] (30) P−1R = J−1 DξJ DηJ C Dξ Dη . (31)
Here we have defined J = Diag(J ), ˆA1 = Diag(ˆa1), ˆA2 = Diag(ˆa2), and
C = ξx ηx ξx ηx T + ξy ηy ξy ηy T (32) where ξx = Diag(ξx), ηx = Diag(ηx), ξy = Diag(ξy) and ηy = Diag(ηy).
In order to show that these operators satisfy SBP rules in the generalized sense (8) and (10), we also need to define a set of boundary operators. Using (28) as well as the definitions in (26), we introduce
EN = (Iξ⊗ e1η), PN = JNPξ, ΛN = Diag(aTn)N = J−1N (ENAˆ2ENT),
EE = (e1ξ⊗ Iη), PE = JEPη, ΛE = Diag(aTn)E = −J−1E (EEAˆ1EET),
ES = (Iξ⊗ e0η), PS = JSPξ, ΛS = Diag(aTn)S = −J−1S (ESAˆ2EST),
EW = (e0ξ⊗ Iη), PW = JWPη, ΛW = Diag(aTn)W = J−1W(EWAˆ1EWT )
(33) where JN,E,S,W = Diag(JN,E,S,W). In the following two propositions we show
that the operators defined by (29)-(33) conform to the general SBP frame-work of section 2.
Proposition 3. The discrete operator P−1Q defined by (29) and (30) satis-fies the summation-by-parts rule,
Q = ENTPNΛNEN + EETPEΛEEE + ESTPSΛSES+ EWT PWΛWEW − QT.
Proof. Using (25) we may write, Q + QT = 1 2[(Qξ⊗ Pη) ˆA1+ ˆA1(Qξ⊗ Pη) + (Pξ⊗ Qη) ˆA2+ ˆA2(Qη⊗ Pξ) + ˆA1(QTξ ⊗ Pη) + (QTξ ⊗ Pη) ˆA1+ ˆA2(QTη ⊗ Pξ) + (Pξ⊗ QTη) ˆA2] = 1 2[((e T 1ξe1ξ− eT0ξe0ξ) ⊗ Pη) ˆA1+ ˆA1((eT1ξe1ξ− eT0ξe0ξ) ⊗ Pη) +(Pξ⊗ (eT1ηe1η− eT0ηe0η)) ˆA2+ ˆA2(Pξ⊗ (eT1ηe1η− eT0ηe0η))] = ((eT1ξe1ξ − eT0ξe0ξ) ⊗ Pη) ˆA1+ (Pξ⊗ (eT1ηe1η− eT0ηe0η)) ˆA2
where in the last step we have used the fact that all matrices involved are diagonal due to Properties 1 and 2, and hence commute. Next we expand the term associated with the north boundary according to,
(Pξ⊗ eT1ηe1η) ˆA2 = (Iξ⊗ eT1η)(Pξ⊗ e1η) ˆA2 = (Iξ⊗ e1η)T(Pξ⊗ 1)(Iξ⊗ e1η) ˆA2 = ENTPξENAˆ2.
Now, since ENENT = I, and ENTEN is diagonal (Property 2), we can further
write,
ENTPξENAˆ2 = ENTPξ(ENENT)ENAˆ2 = ENTPξ(ENAˆ2ENT)EN
= ENT(J−1N PN)(JNΛN)EN = ENTPNΛNEN
where in the second to last step we have applied the definitions in (33). The remaining three faces can be treated in exactly the same way.
Proposition 4. The operator P−1R defined by (29) and (31) satisfies the summation-by-parts rule, R = ET NPNDN + EETPEDE + ESTPSDS+ EWT PWDW − (DxTPDx+ DyTPDy) (34) where, Dx = ξxDξ+ ηxDη, Dy = ξyDξ+ ηyDη and, DN = J−1N EN(−yξDx+ xξDy), DE = J−1E EE(yηDx− xηDy) DS = J−1S ES(yξDx− xξDy), DW = J−1WEW(−yηDx+ xηDy).
Proof. We begin by expanding R using the definitions in (29), (31) and (32), as well as the one-dimensional summation-by-parts rule (25),
R =PJ−1 DξJ DηJ C Dξ Dη = (Qξ⊗ Pη)J (Pξ⊗ Qη)J ξx ηx (ξxDξ+ ηxDη) + ξy ηy (ξyDξ+ ηyDη) = (Qξ⊗ Pη)J (Pξ⊗ Qη)J ξxDx+ ξyDy ηxDx+ ηyDy = ((−eT0ξe0ξ+ eT1ξe1ξ) ⊗ Pη)J (Pξ⊗ (−eT0ηe0η+ eT1ηe1η)J ) ξxDx+ ξyDy ηxDx+ ηyDy − (QT ξ ⊗ Pη)J (Pξ⊗ QTη)J ) ξxDx+ ξyDy ηxDx+ ηyDy . (35) Now consider the boundary terms in (35). We get
((−eT0ξe0ξ+ eT1ξe1ξ) ⊗ Pη)J (Pξ⊗ (−eT0ηe0η+ eT1ηe1η)J ) ξxDx+ ξyDy ηxDx+ ηyDy = (−eT 0ξe0ξ+ eT1ξe1ξ) ⊗ Pη Pξ⊗ (−eT0ηe0η+ eT1ηe1η) J ξxDx+ J ξyDy J ηxDx+ J ηyDy = (−eT0ξe0ξ+ eT1ξe1ξ) ⊗ Pη Pξ⊗ (−eT0ηe0η+ eT1ηe1η) yηDx− xηDy −yξDx+ xξDy , where in the last step we have used the metric relations (27). With exactly the same manipulations that were used for the north boundary term in the proof of Proposition 3, we can now write,
(Pξ⊗ eT1ηe1η)(−yξDx+ xξDy) = ENTPξEN(−yξDx+ xξDy)
= ENTPNJ−1N EN(−yξDx+ xξDy) = ENTPNDN.
The remaining three terms can be rewritten in the same way, and hence the boundary contributions in (35) are equal to those in (34).
contri-bution in (34), (QT ξ ⊗ Pη)J (Pξ⊗ QTη)J ) ξxDx+ ξyDy ηxDx+ ηyDy = (QTξPξ−1⊗ Iη)(Pξ⊗ Pη)J (Iξ⊗ QTηP −1 η )(Pξ⊗ Pη)J ξxDx+ ξyDy ηxDx+ ηyDy = DTξP DT ηP ξxDx+ ξyDy ηxDx+ ηyDy = DξTP(ξxDx+ ξyDy) + DηTP(ηxDx+ ηyDy) =(ξxDξ+ ηxDη)TPDx+ (ξyDξ+ ηyDη)TPDy = DxTPDx+ DTyPDy.
We can thus conclude that (35) leads to (34), which completes the proof. 4.2. The finite volume method
The FV method has previously been shown to satisfy rules similar to the summation-by-parts properties associated with FD schemes, and we refer to [7, 8, 10, 11, 12] for a more detailed derivation of these. In short, a node-centered, edge-based FV method can be shown to yield discrete first derivative operators in the following form,
Dx = P−1Qx, Dy = P−1Qy
where the weights in P are the control volume magnitudes of the dual mesh. Moreover, Qx+ QTx = E TDiag(∆y i)E, Qy + QTy = E TDiag(−∆x i)E,
where ∆xi and ∆yi are the local distances in x− and y−coordinates between
neighbouring dual grid points on the domain boundary. We further define, Q = a1Qx+ a2Qy, R = P(D2x+ D
2 y)
where (a1, a2)T = a. It is now straightforward to confirm that these operators
satisfy (8) and (10), with (n)i = −Diag( 1 ∆si (∆yi, −∆xi)T), P = Diag(∆si), ∆si = q ∆x2 i + ∆y2i. (36)
5. New interface operators for hybrid couplings
We start by discussing the case where a FD block is connected to an-other FD block through an interface. After this is done, we show that the technique can also be extended to hybrid FD-FV couplings. The extension of these techniques for couplings between three-dimensional domains follows immediately using a tensor product formulation.
5.1. FD to FD couplings
The use of SBP preserving interface operators is required whenever the discrete tangential norms on the two sides of the interface are not identical to each other. This may either be a result of using different grid resolutions (non-collocated grids), or from using different types of FD operators (with different orders of accuracy). We will refer to these two cases as h−refinement and p−refinement, respectively. Obviously, both of these could occur simul-taneously, which we refer to as hp-refinement. Note that an hp−refinement operator could either be obtained by combining an h− and a p−refinement operator through Lemma 1, or alternatively be constructed for each combi-nation of different orders and levels of grid refinement.
We will first construct the interface operators with respect to the Carte-sian reference domain ˆΩ, via the same basic methodology that was used in [23]. The transformation to curvilinear coordinates will be done later. The goal is to satisfy the SBP preserving property (14) while minimizing the dis-cretization errors from both ˆIL2R and ˆIR2L. Let the lowest order of the two
FD schemes be accurate to order 2p in the interior. Then the norms PL and
PR integrate grid polynomials up to order 2p − 1 exactly [29, 30], and as a
consequene, the SBP preserving interface operators can be exact for polyno-mials up to order p − 1. This bound follows from a general accuracy result which we prove in Appendix A. Thus, we solve (14) together with
ˆ IL2RxjL= x j L, j = 0, . . . , p − 1 ˆ IR2Lx j R= x j R, j = 0, . . . , p − 1 (37) where xL = (0, ∆xL, 2∆xL, . . .)T, xR = (0, ∆xR, 2∆xR, . . .)T. By Taylor’s
formula, this leads to an accuracy order of p when applying the interface operators to sufficiently smooth functions. Moreover, all free parameters that remain after satisfying these exact relations may be used to minimize
the leading order error terms in a least squares sense. In other words to minimize an expression like
J X j=p (k ˆIL2RxjL− x j Rk 2+ k ˆI R2LxjR− x j Lk 2), (38)
where the nomenclature is defined previously. 5.1.1. p-refinement operators
If the same grid resolution is used on both sides of an interface, then the two SBP FD norms will be identical to each other in the interior. For example, the norms used by a standard second and fourth order interior finite difference scheme in SBP form are respectively given by
ˆ PL= ∆x 1 2 1 1 1 1 . .. , PˆM = 2∆x 17 48 59 48 43 48 49 48 1 . ..
where ˆPL and ˆPR are of the same dimension. The lower right boundary
closures follow from symmetry.
A pair of SBP preserving interface operators with respect to these norms may have the following general form,
ˆ IL2R = × × × × × × × × × × × × × × × × 1 . .. (39)
where ˆIR2L is obtained from (14), and thus has the same structure and
di-mension as ˆIL2R. Note that the boundary closure in (39) must be at least
of dimension 4 by 4, corresponding to the positions in ˆPL and ˆPR where the
integration weights differ from each other. Since the left FD scheme is second order in the interior, we solve (37) for p = 1, and minimize (38) using the remaining free parameters. Note however that both ˆIL2R and ˆIR2L are exact
Remark 3. The same procedure as described above may be applied to finite difference schemes with any combination of interior orders. As in the above example, the required size of the boundary closure in (39) will depend on the scheme with the highest order, while the number of exact accuracy relations will be limited by the scheme with the lowest order.
5.1.2. h− and hp-refinement operators
To demonstrate the general methodology for constructing h− or hp−refinement operators, we consider again the example where a second order FD scheme is used on the left block, and a fourth order on the right. In addition we assume that the grid resolution is twice as fine on the left interface, i.e.
PL= ∆x 1 2 1 1 1 1 . .. , PR= 2∆x 17 48 59 48 43 48 49 48 1 . .. .
Again building on [23], we apply a pair of interpolation stencils in ˆIL2R and
ˆ
IL2R such that the SBP preserving property (14) is satisfied in the interior.
To match the accuracy on the coarse fourth order right block in the interior of the interface, we choose to apply fourth order interpolation stencils of minimal bandwidth, see Figure 2 for an illustration.
ˆ IL2R IˆR2L ˆ PL ˆ PR
Figure 2: Schematic illustraing the fourth order interior interpolation stencils used in the interior of ˆIL2R (40) and ˆIR2L (41). Note that every second row in ˆIR2L represents an
The structure of ˆIL2R, including a block of unknown boundary closure
coefficients, can be written as,
ˆ IL2R= × × × × × × × × × × × × × × × × × × × × −1 32 0 9 32 1 2 9 32 0 − 1 32 −1 32 0 9 32 1 2 9 32 0 − 1 32 . .. . .. . .. . .. . (40)
Note that the first 4 columns of ˆIR2Lwill be affected by the boundary closure
in PR due to (14). Hence the general structure of ˆIR2L becomes
ˆ IR2L= × × × × × × × × × × × × × × × × × × −1/16 × × 0 × × × × × 9/16 −1/16 × × 0 0 1 × × 0 × 9/16 9/16 −1/16 1 −1/16 9/16 9/16 −1/16 . .. . .. . .. . .. . (41)
As before, we solve (37) exactly in the boundary closures for p = 1, and optimize (38) for accuracy using the remaining free parameters.
Following the same basic approach as above we can construct SBP pre-serving interface operators between finite difference blocks of any interior order using non-collocated node distributions at the interface.
5.2. Curvilinear interfaces
The final step in the coupling of arbitrary FD blocks to each other in a SBP preserving way involves a suitable treatment of curvilinear coordi-nate transformations. Consider the case where two FD blocks are connected through an east-west interface, see Figure 3. Let JL = Diag(JE)L and
nL
nR
ΩL ΩR
Figure 3: Coupling of two curvilinear blocks.
two respective block boundaries. Recall that the norms in the transformed domain may be written, as in (33),
PL = JLPˆL PR= JRPˆR. (42)
In the next proposition we show that interface operators derived in a Carte-sian setting may be extended to curvilinear domains.
Proposition 5. Let ˆIL2R and ˆIR2L be SBP preserving interface operators
with respect to [ ˆPL, ˆPR], and let
IL2R= J −1 2 R IˆL2RJ 1 2 L, IR2L= J −1 2 L IˆR2LJ 1 2 R. (43)
Then IL2R and IR2L are SBP preserving with respect to [PL, PR] given in
(42).
Proof. Note that since all norms are diagonal (Property 1), they commute with other diagonal matrices. We can thus write,
PLIR2L= (JLPˆL)(J −1 2 L IˆR2LJ 1 2 R) = J 1 2 L( ˆPLIˆR2L)J 1 2 R = J 1 2 L( ˆI T L2RPˆR)J 1 2 R = J 1 2 L(J −1 2 L I T L2RJ 1 2 R)(J −1 R PR)J 1 2 R= I T L2RPR
where in the third step we have used (14).
After the transformation (43) to curvilinear coordinates, the interface op-erators IL2Rand IR2Lare no longer exact for the first p − 1 polynomials as in
(37), not even for constant functions. As a consequence of this, free-stream preservation to machine precision is not possible when applying these oper-ators in a multi-block scheme. However, the order of accuracy is unchanged
by application of (43). For smooth grid functions ΦL and ΦR it follows from
(37) by application of Taylor’s formula, IL2RΦL = J −1 2 R IˆL2R(J 1 2 LΦL) = J −1 2 R (J 1 2 RΦR) + O(∆xp) = ΦR+ O(∆xp) IR2LΦR = J −1 2 L IˆR2L(J 1 2 RΦR) = J −1 2 L (J 1 2 LΦL) + O(∆xp) = ΦL+ O(∆xp),
since the values contained in JL,R are smooth functions of the metric
coeffi-cients. Moreover, we have already shown in section 3 that the scheme (16) is stable. Hence, the lack of free-stream preservation does not negatively impact either the stability or accuracy of non-collocated interface couplings between curvilinear domains.
5.3. FD to FV couplings
Next, we consider the hybrid case where a FD block on the left side of an interface is facing a FV mesh on the right side. As before the node distributions may be non-collocated by using different grid resolutions on the two sides, and we separate out the interface contributions to the boundary operators on both sides using (9). We start out from the observation made in [10] that for the special case of a non-curved interface with a uniform node distribution, the norm on the FV side exactly coincides with that of a second order interior FD scheme in SBP form. Hence, the same type of FD to FD operators ˆIL2R and ˆIR2L that were discussed in the previous section can be
used in this case.
For a general curved interface however the situation is more complicated. Even if the node distribution on the FV side is given by a smooth trans-formation from a uniform Cartesian grid, the norm PR is no longer exactly
identical to that of a second order interior FD scheme. In particular, PR in
this case only satisfies the second relation in (42) approximately. This means that we can not apply Proposition 5 directly in order to extend operators like ˆIL2R and ˆIR2L to couplings in the hybrid curvilinear case.
To overcome this, we divide the coupling procedure into two separate steps by introducing an intermediate grid as in Lemma 1. This intermediate grid shares the same node distribution as the FV interface, but the norm PM
is associated with a second order interior curvilinear FD scheme, see Figure 4 for an illustration. In summary, we consider the three norms (see (42) and (36)),
ΩL Ω R PL PM PR IL2M IM 2R IM 2L IR2M
Figure 4: A general hybrid FD-FV coupling is divided into two stages using an intermediate grid with the same node distribution as the FV interface. The norm on the intermediate interface grid is associated with a second order accurate curvilinear FD scheme.
We can now immediately apply Proposition 5 to define a pair of summation-by-parts preserving interface operators with respect to PL and PM,
IL2M = J −12 M IˆL2MJ 1 2 L, IM 2L = J −12 L IˆM 2LJ 1 2 M (45)
where ˆIL2M and ˆIM 2L are constructed exactly as discussed in section 5.1.
For the second step, it remains to account and adjust for the fact that PR is
only approximately equal to PM. As we shall now prove, this can be achieved
through simple rescalings with the following pair of diagonal matrices, IM 2R = P −1 2 R P 1 2 M, IR2M = P −1 2 M P 1 2 R. (46)
Note that these operators are accurate simply from the fact that PM ≈ PR.
Lemma 2. The diagonal operators defined in (46) are SBP preserving with respect to [PM, PR]. Proof. We have PMIR2M = PMP −1 2 M P 1 2 R = P 1 2 MP −1 2 R PR= IM 2RT PR.
The following result now follows directly from Lemma 1.
Proposition 6. The hybrid interface operators IL2R = IM 2RIL2M and IR2L=
Hence we conclude that the hybrid interface operators derived in this section can be applied in the generalized interface treatment of section 3 in a way that leads to stability.
5.4. Further examples of interface couplings
In this paper we have derived new interface operators for the FD-FD and FD-FV coupling cases. However, we stress that the general procedure out-lined in section 3 can be used for any pair of methods that can be written in the SBP form (8) and (10). For example, the analysis in section 4.1 can be directly extended to methods based on any tensor product extensions of one-dimensional first derivative SBP operators. This includes spectral and dG operators as well as more general types of SBP operators in one dimen-sion [34]. Once it has been confirmed that (8) and (10) are satisfied, the only remaining step is to construct interface operators that satisfy the SBP preserving condition (14) for the specific norms at hand. In the dG-dG case the interface operators can be constructed directly from the Lagrange poly-nomial bases, see [35] for more details. The FD-dG coupling technique based on SBP preserving projections presented in [25] could also be incorporated into this general framework with minor modifications.
5.5. A final note on accuracy
As we have seen, the achievable accuracy in the boundary closures of SBP preserving operators is limited to p, if the lowest order scheme involved is a FD scheme of order 2p in the interior. After implementing these operators weakly in (16), the resulting truncation errors are reduced to order p − 1 for inviscid penalty terms due to the scaling of integration weights in PL and PR
with the grid size, and to order p − 2 for viscous terms due to the additional factors DLT and DTR. For comparison, the difference operators P−1Q and P−1R introduce truncation errors at block boundaries of order p and p − 1
respectively. In the absence of interfaces requiring SBP preserving interface operators, these latter truncation errors are known to yield convergence rates of order p + 1 for the discrete solution [36]. The larger truncation errors from the interface operators are however restricted to a finite number of points in a two-dimensional domain, and it remains to be investigated how these affect the converge rates of solutions. This is considered next.
6. Numerical calculations
We test the new interface operators for hybrid couplings by applying the general formulation (16) to solve the model problem (1).
-1 -0.5 0 0.5 1 x -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 y
Figure 5: A four-block structured grid where the lower left block is discretized with a lower order scheme than the remaining blocks.
Figure 6: Error distributions using a second order interior FD scheme in the lower left block, and sixth order interior schemes in the remaining blocks, using N = 96 cells in each coordinate direction on each block. Left: = 0, right: = 0.01.
6.1. FD-FD couplings
We begin by studying the convergence properties of a hybrid scheme involving four collocated FD blocks with different orders of accuracy on a Cartesian domain, see Figure 5. We employ the exact solution
u = e−π24 (a 2
1+a22)tsin π
2(a1x + a2y − (a1+ a2)t)
and measure the global error at t = 2. Examples of global error distributions can be seen in Figure 6 using second and sixth order interior schemes, for = 0 and = 0.01 respectively. In both cases, the low order truncation errors at block corners are associated with elevated error levels. In the purely inviscid case = 0, these errors are propagated in the wave direction essentially without damping, while in the viscous case they are more dispersed.
101 102 103 N 10-10 10-5 100 ||e|| ||e|| ||e|| L 2
reference order 2, 2.5, 3, 3.5 & 4
= 0 101 102 103 N 10-10 10-5 100 ||e|| ||e|| L 2
reference order 1, 2, 3 & 4
= 0.01
Figure 7: Mesh convergence using three FD blocks of interior order 8 coupled to a single block with interior order 2, 4 or 6.
Mesh refinement plots are shown in Figure 7 for both the inviscid and the viscous collocated case. For simplicity we use three eighth order blocks in all cases, coupled to a single block of order 2, 4 or 6. These results indicate that, in general, a FD block of interior order 2p yields convergence rates of order p in the maximum norm. Moreover, the asymptotic convergence in the L2 norm is p + 1/2 in the inviscid case and p + 1 in the viscous case. These
rates may be compared to the general rate of p + 1 obtained if the same order is used everywhere, which avoids the need for SBP preserving interface operators. Hence there is a small order reduction in the inviscid case due to
the limited accuracy of SBP preserving interface operators, but not in the viscous case. This order reduction is outweighed in terms of efficiency by the possibility of using non-collocated interfaces, as will be demonstrated next.
103 104 105 106 107 # grid points 10-10 10-5 ||e|| non-collocated collocated reference order 2,3,4 = 0 103 104 105 106 107 # grid points 10-10 10-5 non-collocated collocated reference order 2,3,4 = 0.01
Figure 8: L2error versus total number of grid points. A non-collocated FD block of order
2, 4 or 6 is coupled to three coarse 8:th order blocks. The fine block has 40 × 40 cells in the first measurement, 80 × 80 in the second and so on. These results are compared to the corresponding collocated fine meshes using the same order in all four blocks.
6.1.1. Efficiency
Next, we consider a non-collocated multi-block case obtained by increas-ing the grid spacincreas-ings in the high order blocks in Figure 5 by a factor of 2. In Figure 8 we plot the L2 errors as a function of the total number of grid points.
These results are compared to those obtained using the same lower order dis-cretization and the same grid resolution in all four blocks. Note that in the latter case there is no need for interface operators, leading to convergence rates of p + 1 instead of p + 1/2 in the inviscid non-collocated case. Hence, for = 0 there is always a treshold level at which the collocated multi-block discretization approach becomes more efficient. However, for wide ranges of error tolerance levels the non-collocated approach is the most efficient. In the viscous case = 0.01, there is no order reduction for the non-collocated mesh, while an order of magnitude improvement in efficiency is achieved. 6.2. FD-FV couplings on a structured mesh
In this section we compare the accuracy of the new hybrid FD-FV in-terface operators to that of the previous method of modifying the FV dual
101 102 103 N 10-8 10-6 10-4 10-2 100 ||e||
||e|| , modified cells ||e|| L
2
, modified cells ||e|| , interface operators ||e|| L
2
, interface operators Reference, order 1, 1.5 and 2
= 0, 6:th order FD 101 102 103 N 10-8 10-6 10-4 10-2 100
||e|| , modified cells ||e|| L
2
, modified cells ||e|| , interface operators ||e|| L
2
, interface operators Reference, order 1 and 2
= 0.01, 6:th order FD 101 102 103 10-8 10-6 10-4 10-2 100 ||e|| N
||e|| , modified cells ||e|| L
2
, modified cells ||e|| , interface operators ||e|| L
2
, interface operators Reference, order 1, 1.5 and 2
= 0, 8:th order FD 101 102 103 N 10-8 10-6 10-4 10-2 100
||e|| , modified cells ||e|| L
2
, modified cells ||e|| , interface operators ||e|| L
2
, interface operators Reference, order 1 and 2
= 0.01, 8:th order FD
Figure 9: Mesh convergence plots for a structured FV block coupled to three high order FD blocks using two different interfacing techniques.
grid cells, which requires a collocated mesh. Thus, we insert a structured FV mesh in the lower right block of Figure 5. In Figure 9 we compare the convergence rates using sixth and eigth order interior FD schemes in the re-maining blocks. Both methods converge with order 1 in the L∞ norm, and
1.5 (inviscid case) or 2 (viscous case) in the L2 norm. However, the new
procedure achieves much lower error levels. The asymptotic improvement in the L2 norm is about an order of magnitude in the inviscid case, and even
greater in the L∞ norm. This again points to the superior performance of
the new hybrid framework.
6.3. Curvilinear FD to unstructured FV
In order to demonstrate the flexibility and potential gains in efficiency from the new hybrid approach, we consider a slightly more complex case
Figure 10: Schematic of a hybrid mesh with curved interfaces and grid refinement across interfaces. The numbers indicate the order of the FD interior stencils
Figure 11: Left: initial solution. Right: exact solution at t = 5.
involving an unstructured FV grid surrounded by several curvilinear FD blocks of varying order and grid resolution. The geometry and the different levels of refinement in both order and grid resolution are illustrated in Figure
Figure 12: Left: global error distribution at t = 5 for the hybrid mesh shown in Fig-ure 10, using 1024 cells around the circumference of the inner disc-shaped FV domain. Right: corresponding error distribution for a pure FV mesh using the same grid resolution everywhere as in the FV part of the hybrid mesh. Note the different scalings of the errors.
105 106 107 # nodes 10-5 10-4 10-3 10-2 ||e|| hybrid mesh FV mesh
Figure 13: Comparison between the hybrid mesh and a pure FV mesh in terms of the total number of grid points required for different levels of accuracy in the L2norm.
10. At the interfaces between the FV and FD parts of the domain, the nodal grid resolution on the FV side is finer by a factor of 2.
We consider the more demanding purely inviscid case = 0, and impose an exact solution given by several gaussian pulses that propagate through the domain. The initial solution as well as the exact solution at t = 5 are shown in Figure 11. An example of global error distribution at t = 5 is given in Figure 12, comparing the hybrid mesh to a pure FV mesh with the same grid resolution across the domain. The enhanced error levels associated
with block corners at the hybrid FD to FV interfaces are visible also here. However these errors are too small to negate the superior efficiency gained by using coarse high order finite difference blocks throughout most of the domain. An estimate of the increased efficiency is given in Figure 13, where the mesh complexity (the total number of nodes) is measured against the L2
error for a few different mesh realizations, including both the hybrid and the pure FV mesh cases. For the same accuracy, the FV mesh requires about an order of magnitude more computational nodes.
7. Conclusions
In this paper we have developed a stable general procedure to couple hybrid multi-block meshes in computational fluid dynamics using SBP pre-serving interface operators. The procedure applies to a wide range of SBP based numerical schemes, and requires minimal modifications of the different schemes involved.
Accordingly, we have derived new optimized SBP preserving interface operators for the coupling of curvilinear finite difference blocks of varying order and grid resolution, as well as between finite difference blocks and an unstructured finite volume mesh. In particular, the new FD-FV couplings are more flexible, easier to implement and leads to smaller discretization errors compared to the previous hybrid method, while still maintaining linear stability.
Our numerical results confirm gains in efficiency from using a hybrid mesh, in particular from the ability of using a fine unstructured grid on a small part of the domain, and more coarse structured high order blocks on the rest. An order of magnitude increase in efficiency is demonstrated, which is deemed a significant achievement.
Acknowledgements
This work is based on the research supported by the South African Re-search Chairs Initiative of the Department of Science and Technology and Na-tional Research Foundation of South Africa. In addition, part of the research leading to this work was supported by the AEROGUST project (funded by the European Commission under grant agreement number 636053). The partners in AEROGUST are: University of Bristol, INRIA, NLR, DLR, Uni-versity of Cape Town, NUMECA, Optimad Engineering S.r.l., UniUni-versity of
Liverpool, Airbus Defence and Space, Dassault Aviation, Piaggio Aerospace and Valeol.
References
[1] H.-O. Kreiss, G. Scherer, Finite element and finite difference methods for hyperbolic partial differential equations, in: C. De Boor (Ed.), Math-ematical Aspects of Finite Elements in Partial Differential Equation, Academic Press, New York, 1974.
[2] M. H. Carpenter, D. Gottlieb, S. Abarbanel, Time-stable boundary con-ditions for finite-difference schemes solving hyperbolic systems: Method-ology and application to high-order compact schemes, Journal of Com-putational Physics 111 (1994) 220–236.
[3] M. Carpenter, J. Nordstr¨om, D. Gottlieb, A stable and conservative interface treatment of arbitrary spatial accuracy, Journal of Computa-tional Physics 148 (1999) 341–365.
[4] J. Nordstr¨om, M. H. Carpenter, High-order finite difference methods, multidimensional linear problems and curvilinear coordinates, Journal of Computational Physics 173 (2001) 149–174.
[5] M. H. Carpenter, J. Nordstr¨om, D. Gottlieb, Revisiting and extend-ing interface penalties for multi-domain summation-by-parts operators, Journal of Scientific Computing 45 (2010) 118–150.
[6] J. Nordstr¨om, J. Gong, E. van der Weide, M. Sv¨ard, A stable and conservative high order multi-block method for the compressible Navier-Stokes equations, Journal of Computational Physics 228 (2009) 9020– 9035.
[7] J. Nordstr¨om, K. Forsberg, C. Adamsson, P. Eliasson, Finite volume methods, unstructured meshes and strict stability for hyperbolic prob-lems, Applied Numerical Mathematics 45 (2003) 453–473.
[8] M. Sv¨ard, J. Nordstr¨om, Stability of finite volume approximations for the laplacian operator on quadrilateral and triangular grids, Applied Numerical Mathematics 51 (2004) 101–125.
[9] M. Sv¨ard, J. Gong, J. Nordstr¨om, Stable artificial dissipation operators for finite volume schemes on unstructured grids, Applied Numerical Mathematics 56 (2006) 1481 – 1490.
[10] J. Nordstr¨om, J. Gong, A stable hybrid method for hyperbolic problems, Journal of Computational Physics 212 (2006) 436–453.
[11] J. Gong, J. Nordstr¨om, A stable and efficient hybrid scheme for viscous problems in complex geometries, Journal of Computational Physics 226 (2007) 1291–1309.
[12] J. Nordstr¨om, F. Ham, M. Shoeybi, E. van der Weide, M. Sv¨ard, K. Mattsson, G. Iaccarino, J. Gong, A hybrid method for unsteady inviscid fluid flow, Computers & Fluids 38 (2009) 875 – 882.
[13] M. Carpenter, D. Gottlieb, Spectral methods on arbitrary grids, Journal of Computational Physics 129 (1996) 74–86.
[14] M. H. Carpenter, T. C. Fisher, E. J. Nielsen, S. H. Frankel, Entropy stable spectral collocation schemes for the navier–stokes equations: Dis-continuous interfaces, SIAM Journal on Scientific Computing 36 (2014) B835–B867.
[15] G. J. Gassner, A skew-symmetric discontinuous galerkin spectral ele-ment discretization and its relation to SBP-SAT finite difference meth-ods, SIAM Journal of Scientific Computing 35 (2013) 1233–1253. [16] P. Castonguay, D. Williams, P. Vincent, A. Jameson, Energy stable
flux reconstruction schemes for advection-diffusion problems, Computer Methods in Applied Mechanics and Engineering 267 (2013) 400 – 417. [17] H. Ranocha, P. ¨Offner, T. Sonari, Summation-by-parts operators for
cor-rection procedure via reconstruction, Journal of Computational Physics 311 (2016) 299 – 328.
[18] J. Nordstr¨om, T. Lundquist, Summation-by-parts in time, Journal of Computational Physics 251 (2013) 487–499.
[19] T. Lundquist, J. Nordstr¨om, The SBP-SAT technique for initial value problems, Journal of Computational Physics 270 (2014) 86–104.
[20] P. D. Boom, D. W. Zingg, High-order implicit time-marching methods based on generalized summation-by-parts operators, SIAM Journal on Scientific Computing 37 (2015) A2682–A2709.
[21] S. Nikkar, J. Nordstr¨om, Fully discrete energy stable high order finite difference methods for hyperbolic problems in deforming domains, Jour-nal of ComputatioJour-nal Physics 291 (2015) 82 – 98.
[22] S. Nikkar, J. Nordstr¨om, A fully discrete, stable and conservative summation-by-parts formulation for deforming interfaces, Journal of Computational Physics 339 (2017) 500 – 524.
[23] K. Mattsson, M. H. Carpenter, Stable and accurate interpolation oper-ators for high-order multiblock finite difference methods, SIAM Journal on Scientific Computing 32 (2010) 2298–2320.
[24] A. Nissen, K. Kormann, M. Grandin, K. Virta, Stable difference meth-ods for block-oriented adaptive grids, Journal of Scientific Computing 65 (2015) 486–511.
[25] J. E. Kozdon, L. C. Wilcox, Stable coupling of nonconforming, high-order finite difference methods, SIAM Journal on Scientific Computing (2016) A923–A952.
[26] S. Wang, K. Virta, G. Kreiss, High order finite difference methods for the wave equation with non-conforming grid interfaces, Journal of Scientific Computing (2016) 1–27.
[27] M. Sv¨ard, J. Nordstr¨om, Review of summation-by-parts schemes for initial-boundary-value problems, Journal of Computational Physics 268 (2014) 17–38.
[28] M. Sv¨ard, J. Nordstr¨om, A stable high-order finite difference scheme for the compressible Navier-Stokes equations: No-slip wall boundary conditions, Journal of Computational Physics 227 (2008) 4805–4824. [29] J. E. Hicken, D. W. Zingg, Summation-by-parts operators and high
order quadrature, Journal of Computational and Applied Mathematics 237 (2013) 111–125.
[30] V. Linders, T. Lundquist, J. Nordstr¨om, On the Order of Accuracy of Finite Difference Operators on Diagonal Norm Based Summation-by-Parts Form, LiTH-MAT-R 2017/11, 2017.
[31] J. Nordstr¨om, Conservative finite difference formulations, variable coef-ficients, energy estimates and artificial dissipation, Journal of Scientific Computing 29 (2006) 375–404.
[32] M. Sv¨ard, On coordinate transformations for summation-by-parts op-erators, Journal of Scientific Computing 20 (2004) 29–42.
[33] J. Nordstr¨om, A. A. Ruggiu, On conservation and stability properties for summation-by-parts schemes, Journal of Computational Physics 344 (2017) 451–464.
[34] D. C. Del Rey Fern´andez, P. D. Boom, D. W. Zingg, A general-ized framework for nodal first derivative summation-by-parts operators, Journal of Computational Physics 266 (2014) 214–239.
[35] T. Lundquist, J. Nordstr¨om, An analysis of non-conforming grid techniques for high order summation-by-parts methods, LiTH-MAT-R 2017/02, 2017.
[36] M. Sv¨ard, J. Nordstr¨om, On the order of accuracy for difference approx-imations of initial-boundary value problems, Journal of Computational Physics 218 (2006) 333–352.
Appendix A. Accuracy
We show that the formal accuracy of SBP preserving interface operators satisfying Definition 1 is limited by the following general result.
Proposition 7. Consider a pair of summation-by-parts preserving interface operators [IL2R, IR2L] (14) with respect to the discrete norms [PL, PR], and
denote the nodal grids at the interface with XL and XR. Let ΦL = φ(XL),
ΦR = φ(XR), ΨL = ψ(XL) and ΨR = ψ(XR) denote the nodal interface
values of two different polynomials φ and ψ. If both IL2R and IR2L are exact
for all grid polynomials of order q, then the norms PL and PR must also
satisfy
(ΦL, ΨL)PL = Φ
T
LPLΨL = ΦTRPRΨR= (ΦR, ΨR)PR, (A.1)
Proof. Since IL2R and IR2L are exact for all grid polynomials of order q, we
have
IL2RΦL = ΦR, IR2LΨR= ΨL. (A.2)
Next we multiply the second equation above with ΦT
LPL from the left, to
obtain ΦTLPLIR2LΨR = ΦTLPLΨL. The summation-by-parts preserving
prop-erty (14) now gives (IL2RΦL)TPRΨR = ΦTLPLΨL. Finally, inserting the first
equation in (A.2) yields,
ΦTRPRΨR= ΦTLPLΨL
which proves the proposition.
For example, if PL are PR are norms are associated with FD schemes of
at least order 2p in the interior, then (A.1) is satisfied for q = p − 1, since both PL are PR are exact for all polynomials up to order 2p − 1 [29, 30]. This
