Energy stable and accurate coupling of finite element methods and finite difference methods

We introduce a hybrid method to couple continuous Galerkin finite element methods and high-order finite difference methods in a nonconforming multiblock fashion. The aim is to optimize computational efficiency when complex geometries are present. The proposed coupling technique requires minimal changes in the existing schemes while maintaining strict stability, accuracy, and energy conservation. Results are demonstrated on linear and nonlinear scalar conservation laws in two spatial dimensions.


Introduction
Finite element (FE) methods and finite difference (FD) methods are among the most studied and employed numerical methods for partial differential equations (PDEs), which serve numerous industrial applications and other research areas, e.g., economics, engineering, chemistry, physics. FD methods are often seen as simple and efficient schemes for high-order approximations, see e.g., Mattsson and Nordström [23], Shu [28], despite being restricted to structured grids. FE methods are favored in various computational areas due to their theoretical reliability, generality, and the ability to handle complex domain shapes. However, the obtained matrices can be storage demanding; and the method itself is often considered to be more computationally expensive. For instance, while solving time-dependent nonlinear diffusion dominated problems, one has to reconstruct FE matrices in each time step, which is time-consuming. The greatest strengths of both methods can be combined by developing a hybrid scheme that allows different methods and refinement levels to be involved simultaneously. Hence, computational efficiency can be enhanced by utilizing the most suitable treatment for each portion of the domain. The major complication is known to consist in designing a method-to-method interface that results in an accurate and stable approximation. A favorable method should also accommodate the demands for simplicity, freedom to choose domain discretization, and minimal modification in the existing schemes to preserve their advantageous native properties.
A well-designed numerical framework for energy stable approximations of time-dependent problems is the combination of the summation-by-parts (SBP) operators, see Kreiss and Scherer [15], and the simultaneous-approximationterm (SAT) technique, see Carpenter et al. [4]. The SBP operators are used to approximate the governing equations. By design, they mimic the continuous integrations-by-parts properties. The SAT technique then imposes the boundary conditions in a provably stable manner. We refer the reader for a detailed discussion on the theory and applications of SBP-SAT in the following papers and references therein, [5,9,14,17,18,23,24,31]. A crucial benefit of numerical schemes in this form is that a proper formulation strictly prevents nonphysical energy growth in the discrete approximation, a property often referred to as "strict stability". The vast majority of previous papers on the SBP-SAT framework have been committed to FD schemes due to high-order operators' apparent advantages. Nontrivial geometries are then usually handled with curvilinear grids, see Nordström and Carpenter [24], Shadpey and Zingg [27], curved elements, see Crean et al. [7], and multiblock coupling, see Carpenter et al. [5]. To carry out further practical geometries, successful attempts have been made to couple FD methods with unstructured finite volume (FV) methods, see Nordström and Gong [25], Nordström et al. [26] by modifying the FV scheme at interfaces. The coupling of mixed order schemes on nonconforming grids was, for the first time, proposed in an energy stable manner by Mattsson and Carpenter [22] introducing SBP-preserving interpolation operators. Fundamental properties of SBP schemes are extended to blockto-block coupling: strict stability, accuracy, and conservation. This technique paves the way for many possibilities to couple different methods which can be written in SBP form. So far, the coupling of FD-FD by Mattsson and Carpenter [22], FD-discontinuous Galerkin (dG) method by Kozdon and Wilcox [14], and FD-FV by Lundquist et al. [18] are available in the literature.
In practical computation, properties of the concerning problem, e.g., nature of the solution, desired precision, geometries, often vary over the domain. An impediment with discretizations by a single conventional method, especially when using uniform precision, is that sometimes the preferable scheme or the required resolution for a minor region may bring up extreme computational necessity over the whole domain, e.g., when complex geometries present. Despite being of great interest, attempts to combine FD methods and FE methods have rarely been successful. Independent techniques have been proposed for several applications, e.g., seismic wave propagation, see [11,20,12], electrical modeling, see [30]. For example, [11,20] require a transition region between the two methods; the interface solution is updated similarly to imposing Dirichlet boundary conditions with a suitably averaging solution. In [30], at regions of interest, each rectangular FD block is split into two FE triangular elements resulting in locally modified schemes; accuracy is observed in cases of sufficiently gradual topographic gradients. Most relevantly, Gao and Keyes [12] derive an energy stable coupling of SBP FD schemes and FE schemes with quadrilateral elements under shared interface vertices/nodal points. Inline with other hybrid SBP methods, although the continuous Galerkin FE also has SBP properties, see Zemui [31], one difficulty in deriving a hybrid method lies within its nondiagonal mass matrix (can be understood as an SBP norm matrix). Note that the most related references, such as Lundquist et al. [18], Gao and Keyes [12], necessitate diagonal norm matrices.
In this paper, our main focus is to develop a hybrid method to couple FD and FE methods in a nonconforming multiblock fashion. The proposed technique results in a provably stable, accurate, and energy-conserved unified SBP method. Interface continuity is weakly imposed by the SAT technique using nondiagonal norm SBP-preserving interpolation operators. The SAT treatment on the FD side is the same as in our FD-FD coupling [22]. The SAT treatment for FE that we propose in this current paper is formulated in a general expression in the sense that it does not depend on the coupled scheme or the coupled mesh/grid. Furthermore, we show that even when considering nonlinear conservation laws, nondiagonal norm matrices, namely the continuous FE mass matrices, could be applicable without further techniques involved, e.g., mass lumping [31], to the classical Galerkin formulation. We include two techniques to construct the required interpolation operators. In the latter one, we introduce a way of acquiring the target interpolation operators through existing diagonal norm techniques, for instance, [1,14,18,19,22]. This paper is organized as follows. In Section 2, the numerical schemes of both FE and FD are presented in the context of SBP framework. Section 3 describes the proposed coupling technique. We show numerical verification and computational results in Section 4 on an advection-diffusion equation and a nonlinear Burgers' equation. Section 5 is the conclusion. Two approaches to constructing the required interpolation operators are shown in Appendix A.

Summation-by-parts operators
This section formulates summation-by-parts schemes by FD and FE approximations. Boundary treatment is briefly mentioned.

The continuous problem to a viscous scalar conservation law
Let Ω be a bounded domain in ℝ 2 . Consider the following scalar conservation laws where the flux vector ( ) = 1 ( ), 2 ( ) consists of two linear or nonlinear mappings 1 , 2 ∈ 1 (Ω; ℝ); = ( , ) is a positive-valued function on Ω; and is the normal vector pointing outward at the boundary. The notation "∇" refers to the gradient operator ∇ = , . In analysis of FD schemes, it is more convenient to write (1) as where the subscript notations , , indicate the partial derivatives with respect to the corresponding variables. Let ′ = ( ′ 1 , ′ 2 ) be the first derivative of with respect to the conserved variable . We assume that there exists a skew-symmetric splitting with the weight functions ,̃ , such that by using the alternative flux form, the total energy over time of a solution to (1) only conserves or decays, on the boundary or due to viscous dissipation, yet to be quantified. It is worth noting that the above assumption is not too severe, in which the existence of such a splitting is shown for a broad class of fluxes by Tadmor [29] by symmetrizing variables. This splitting technique is necessary for obtaining nonlinear energy estimates to (1). In (3), the following notations are used ∇ = 1 ∕ , 2 ∕ ,̃ ∇ = ̃ 1 ∕ ,̃ 2 ∕ . Similar notations used below are ( ) = ( 1 1 , 2 2 ) , (̃ ) = (̃ 1 1 ,̃ 2 2 ) , = ( 1 , 2 ) . One option of specifying the boundary condition is where Ω denotes the domain boundary. Stability of the problem (1) subject to the boundary condition (4) is analyzed below.
Definition 2.1. Given two real-valued vector functions 1 , 2 , and a real-valued function ( , ) > 0, the following definitions of a continuous inner product and its corresponding norm are used When ( , ) ≡ 1, the subscript notation is skipped. Inner products and a norm of two scalar functions are defined analogously by considering 1 , 2 as scalar functions. One can investigate necessary conditions for well-posedness of (1) with the boundary condition (4) using the energy method, see e.g., Kreiss and Scherer [15], Fernández et al. [9]. For continuous problems that can be written in the form = ( , ), a basic idea of the energy method is to examine energy growth utilizing the relation ‖ ‖ 2 = ( , ) + ( , ). A well-posed problem implies that ‖ ‖ 2 is nonpositive given zero boundary data. For two smooth scalar functions 1 , 2 ∈ 1 (Ω) and a smooth vector function 3 , continuous integration rules read Applying (5) to (1) gives Adding up the two equalities above yields where all the interior terms are cancelled by the following skew-symmetry assumption If the splitting is defined pointwise, i.e., 1 , 2 are constant weights, a sufficient spatial-derivative-free argument to (7) reads Given a specific flux, one can either solve (7) or (8) to determine the value of . It becomes apparent that the boundary condition (4) imposes nonlinear stability to (1) since bounded energy is followed. Indeed, inserting (4) into (6) with zero boundary data ≡ 0 gives

Properties of the discrete operators
Considering a general numerical method, we denote N the total number of degrees of freedom of the discretization, and N Ω the number of degrees of freedom distributed on the boundary. Let  be an integration operator over the whole domain  ≈ ∫ Ω , ∈ ℝ N , ≈ ( , );  be an integration operator over the boundary  | Ω ≈ ∫ Ω ; and  be a boundary selection operator | Ω =  . Precise definitions of the three matrix operators , ,  are often required by conventional summation-by-parts schemes, see e.g., Kreiss and Scherer [15], Lundquist et al. [18]. In the current study, , ,  are assumed to have the following properties: ( ) ,  are symmetric, positive definite;  is of size N × N; and  is of size N Ω × N Ω ; ( )  comprises a nonsquare binary matrix of which the only "1" element in each row one-to-one corresponds to a boundary index. Definition 2.2. Given two vectors , ∈ ℝ N , the following definitions of a discrete -inner product and its corresponding discrete -norm are used, The matrix  is often referred to as the norm matrix. Advantages of an SBP scheme come from the fact that its discrete operators, by design, mimic the continuous integration-by-parts. The following definitions introduce discrete operators of which properties are analogous to the continuous integration rules (5), see Carpenter et al. [4], Nordström and Carpenter [24].
where Diag( ) denotes the square diagonal matrix of which diagonal is the discrete vector .

Definition 2.4.
A matrix operator  −1  ( ) approximating ∇ ⋅ ∇ is said to be an SBP diffusion operator if where  =  ≥ 0 and  is a consistent approximation of ⋅ ∇ at the boundary.
The discrete energy method can be analogously applied to SBP schemes utilizing the properties (10), (11). A necessary condition for strict stability is that ‖ ‖ 2  being nonpositive given zero boundary data. We introduce FD notations in the following subsection.

The FD scheme
Consider a one-dimensional domain [ , ] discretized by a grid of equally spaced points { ≡ 1 < 2 < ⋯ < ≡ }. We first look at one-dimensional SBP operators for FD. The following two definitions (first presented in [21,23]) are central.
A discrete solution vector approximating a continuous function ∶ Ω ⟶ ℝ is then presented by where , approximates ( , ). It's often convenient to interprete as a "vector of vectors". The total number of degrees of freedom is N = × .
Discrete operators for the two-dimensional approximation are listed below. Definition of the Kronecker product "⊗" can be found in e.g., Mattsson and Carpenter [22]. Denote the identity matrix by , and the one of size × by . The following FD operators will be frequently used in the analysis, where , are the one-dimensional norm matrices in − and − directions. For ease of reading, it is worth noting that the first components in all of the above Kronecker products are of size × , and the second components are of size × . Note that the matrix ( ) 2 in 2 , 2 varies for different rows and columns, namely, ( ) , being a vector of length , corresponds to the value of ( , ) at ( , ), = 1, 2, … , . The corresponding two-dimensional norm matrix is ∶= = ⊗ . A consistent FD semi-discretization of (2) using diagonal norm SBP operators is given by where 1 , ′ 1 , 2 , ′ 2 are, respectively, the vectors containing values of 1 , ′ 1 , 2 , ′ 2 evaluated at nodal points with respect to the discrete variable . The diagonal matrices ′ 1 and ′ 2 are defined as ′ 1 = Diag( ′ 1 ( )), ′ 2 = Diag( ′ 2 ( )). The following proposition can be straightforwardly derived from how the above operators are defined, but useful in connection with FE schemes in the next section.
where , , , are vectors containing values of respectively evaluated at North, East, South, West boundary grid points;  = Diag( ), where is the discrete analogue of evaluated at the grid points,  ,  ,  ,  are parts of  at the North, East, South, West boundaries; and The following theorem completes the formulation of the FD approximation to (1), (4). Proof. Let = Diag( ). Deriving a discrete energy estimate from (13) and (14), one has It can be seen that the portion corresponding to each boundary is nonpositive. For example, the one corresponding to the North boundary is −(Diag(sign( )) ) ( The vector Diag(sign( )) = | | is also nonnegative. Therefore, the right-hand side of the above estimate is less than or equal to zero. Strict stability is thus guaranteed. The above estimate is analogous to the continuous energy estimate (9).

The FE scheme
We formulate and prove the SBP properties of the FE methods in this section.

The Galerkin FE approximation
Let { ℎ } ℎ>0 be a triangulation mesh of Ω, and ℙ 1 ( ) is the space of all linear functions on each element of { ℎ }. Denote 0 (Ω) the space of continuous functions on Ω. We seek the solution of (1) in the space of weakly differentiable functions 1 (Ω). To this purpose, we define a scalar-valued Lagrange FE space Note that all the degrees of freedom of continuous Galerkin FE methods are located on the nodal points. More specifically, using ℙ 1 elements, they are the vertices of the triangulation  ℎ . A FE approximation of (1) subject to the boundary condition (4) can then be formulated as: find ℎ ∈ ℎ such that We can write ℎ = span{ } N =1 , where , = 1, … , N are piecewise linear Lagrange basis functions. Inserting the (16) gives where the over dot denotes the time derivative; , ( ), , , R ( ), are respectively called the mass matrix, the convective flux matrix, the stiffness matrix, the Robin boundary mass matrix, the nonlinear Robin boundary mass matrix, and the Robin boundary vector. A different notation style is used for R ( ) to indicate that each entry corresponds to a -weighted line integral. The matrices are of full size N × N. To be more precise, for , = 1, … , N, As in the standard continuous Galerkin method, given that , , present nodal weights, the above integrals can be evaluated exactly, see Remark 2.1. We write instead of ( ) when ≡ 1. Pseudo-code for the assembly of the above matrices can be found in Larson and Bengzon [16], and Dao [8].
Remark 2.1. To clarify, we use the following quadrature rule to calculate the -, -, -dependent matrices, for each element on  ℎ , and for each edge on the boundary Γ. The above integrals are exact if , , are considered as nodal weights. Later we will show that the stability proofs do not rely on these quadrature rules. Thus, any other consistent rules could have been used. The use of (18) merely simplifies the SATs for coupling FE.
Then, the advection operator −1 and the diffusion operator (17) satisfy the summation-by-parts rules (10) and (11), with the norm matrix .
Proof. We prove that the FE advection operator under the general notation  = satisfies (10). By integrationby-parts, we have + = R ( ⋅ ). We need to show that R ( ⋅ ) =   Diag( ⋅ )  with the FE formulation of  and . For a general boundary selection operator , each row of  must be a one-to-one mapping from its "one" element with a node on the boundary. The combination of a left multiplier  and a right multiplier  thus expands the local boundary operators into the full domain size (with zeros on the extended rows and columns).

Relation (10) becomes obvious with the local boundary integration operator inserted
For the diffusion operator  = − + , we need to show that − + Ω =   −  with the FE formulation of , , , and . Firstly, we identify( ℎ ) ∈ ℎ , ℎ ∈ ℎ approximating ⋅ ∇ ℎ . The operator satisfies into the above equation gives where = ∫ Ω , , = 1, 2, … , N and Ω is defined above. Since ∫ Ω = 0 if either or corresponds to a nonboundary node, The multiplication of  and  is because  is defined in the SBP framework of size N Ω × N Ω . The matrix  is obtained by projecting onto the boundary,  =  . Note that  is not square because of the matrix dimensions in (11). From (19) and (20), we obtain Ω =  .
What is left is to prove the positive semi-definiteness of .
for all ℎ ∈ ℎ . Therefore, the desired property of  is achieved.
Remark 2.2. By the energy method, conservation property of the discrete fluxes (13) and (17) approximating the continuous skew-symmetric flux (3) is fulfilled under the following condition which naturally holds if the condition (8) of the flux holds for every point in space.

Simultaneous-approximation-term technique for interface coupling
This section describes our proposed SAT technique for interface treatment. To make this paper self-contained, we derive the coupling of both FD-FD and FD-FE. The former one is not the main focus, and partially covered in the literature, see Lundquist et al. [18], Mattsson and Carpenter [22]. Unlike the mentioned references which investigate linear problems, the technique presented in this paper is extended to nonlinear problems, as can be nonlinear in (1). Description of compatible SBP-preserving interpolation operators for nondiagonal norm schemes is included.

Continuous analysis
For simplicity, the analysis is done on the coupling of two domains Ω , Ω illustrated in Figure 2, where and present the parts of solution on the left and right domains, respectively, Ω Ω L R Ω I Figure 2: Coupling of two adjoining blocks.
The North, South, West boundaries of Ω and the North, East, South boundaries of Ω will be referred to as the "outer boundary" and the shared edge of the two adjoining blocks is referred to the "interface", denoted Ω . With the setup in Figure 2, Ω can be seen both as the East boundary of Ω and the West boundary of Ω . Applying the energy method on (24), a total energy estimate of the whole solution on Ω ∪ Ω reads where contains all the outer boundary terms being bounded with proper boundary treatment as discussed in Section 2, and contains the interface terms given by Therefore, a sufficient condition for an accurate and energy-conserved coupling is

Necessary properties of interface interpolation for stability and conservation
Let , be some discrete approximate solutions of , . An SBP-SAT semi-discretization of (24) is given by where SAT , SAT present the weak boundary treatment, and SAT , SAT present the weak interface treatment.
Another advantage of this weakly forcing technique is that SAT , SAT , SAT , SAT each can be treated separately by the energy method. Therefore, in the following subsections, we assume stable outer boundary treatment and only analyze the estimate terms regarding the additional involvement of SAT , SAT . To prove stability and energy conservation, making use of the following class of interpolation operators to translate the shared-interface nodal values, first introduced in Mattsson and Carpenter [22], is necessary.
We call 2 , 2 a pair of SBP-preserving interpolation operators if where = min( , ), and , are the orders of accuracy in boundary closure of the two involving methods.
The matrix mappings 2 and 2 are illustrated in Figure 3. For FD,  is the -norm we defined earlier. For FV, this operator is a diagonal matrix containing the Euclid distances between two neighboring interface nodes, Lundquist et al. [18]. For dG, this operator is a diagonal matrix presenting integration weights for interface edges, Kozdon and Wilcox [14]. However, it is well-known that in continuous Galerkin FE schemes, the one-dimensional norm matrix, i.e., the mass matrix, is not diagonal. In Appendix A, we present two techniques of constructing the pair 2 , 2 which satisfy and (27) and (28). The first technique Appendix A.1 is similar to Mattsson and Carpenter [22]. Appendix A.2 describes a way of acquiring the target interpolation operators through existing diagonal norm techniques. This construction is beneficial due to the fact that various techniques to construct the necessary diagonal norm interpolation matrices can be found in the literature, e.g., Almquist et al. [1], Kozdon and Wilcox [14], Lundquist et al. [18,19], Mattsson and Carpenter [22].

The FD-FD coupling
We consider the case where FD schemes are used on both Ω and Ω . Let ∈ ℝ N be an arbitrary FD solution. We denote by ∈ ℝ the part of on Ω , where grid points are distributed on Ω . For example, in the single-domain notations (15), the orientation of Ω , Ω gives that coupling SAT treatment is given by where SAT , SAT weakly couple the flux ( ) = ( ) at Ω SAT = ( −1 ) ⊗ ( 1 1 ( )) − 2 ( 1 1 ( )) , The terms SAT ′ , SAT ′ weakly couple = at Ω with the weights of 1 ′ ( ), 1 ′ ( ), respectively contain the viscosity coefficients of the left and the right domain evaluated at the grid points on Ω . The terms SAT , SAT naturally allow interface jumps of the viscosity , i.e., when ≠ . This property can be useful in applications involving discontinuous media, or when a significant difference in the amount of artificial viscosity is needed to stabilize the coupled schemes.  Proof. Applying the energy method on each equation of (26) and adding them together, a total energy estimate can be written as where the outer boundary terms all contained in are assumed to be finely tuned as given in (14) and Theorem 2.7; and the interface terms consist of where , ′ estimate the energy growth contributed by the discrete operators approximating the flux ∇ ⋅ ( ); concerns the energy growth by the operators approximating ∇ ⋅ ( ∇ ); the terms , ′ , are the contributions of the added SATs (29). Concretely, the terms SAT , SAT , resulting in the estimate , are used to control , The estimate is derived from the SBP property in Definition 2.3. By substituting = 1∕2, = −1∕2, one can see that is canceled by . Accordingly, the terms SAT in which ′ is also canceled by ′ with = 1∕2, = −1∕2. The remaining mixed terms containing 2 , 2 are The SBP property (28) implies that 2 = 2 . Moreover, we have 1 1 ( ) = ̃ 1 ′ 1 ( ) and 1 1 ( ) = ̃ 1 ′ 1 ( ) by the discrete skew-symmetry property (23). The above choice of the penalty parameters makes = − and = − . Therefore, the sum of two terms in the first line of (30) is zero. The sum of the two remaining terms in the second line is also zero. and ′ together imply conservation to the interface interpolation due to the discrete skew-symmetry (23). The same can be said for the other pair and ′ . Similarly, the estimate of SAT , SAT are used to control One can see that the mixed terms, e.g., , are canceled as long as we choose = − , = − . To cancel , we also need + = −1. Therefore, it is sufficient to conclude that with the above choice of parameters, all the interface terms vanish ≡ 0. The interface treatment is thus stable and conservative.

The FD-FE coupling
This section contains the main contribution of this paper. The following definitions align FE schemes with the coupling to FD schemes. (i) Each row of contains exactly one "1" element uniquely corresponding to a node on the concerned part of the boundary, (ii) rearranges the boundary nodes to the same spatial order of the coupled solution.
The following properties of the above boundary selection operators are useful for the stability proofs in this section.

Property 3.4. Let be a FE boundary selection operator. The following statements hold
where , are matrices of appropriate sizes.
Proof. We prove (i). From = , right multiplying both sides by gives = .
Since each row of contains exactly one "1" element corresponding to the same position on each column of , it is clear that = . We have the first statement proven.
Similarly, we prove (ii) by right multiplying both sides by to obtain = . However, ≠ N . The additional constraint ( N − ) = 0 is necessary for getting = . Assuming that is a square matrix of size N × N, this additional constraint can be understood as, all the columns of of which indices corresponding the nonboundary nodes are zero.
We consider the setup in Figure 2 where FD is used on Ω and FE is used on Ω . The interface Ω is thus both the West boundary of Ω and the East boundary of Ω . Therefore, the FE operators activating on the West boundary, e.g., , and the FD operators activating on the East boundary, e.g., ( ⊗ ), are frequently used in this section. The notations , are used instead of , since there is only one domain involving FD. Definition 3.1 is restated below in the FD-FE coupling context.

Definition 3.5. { 2 , 2 } is called a pair of SBP-preserving FD-FE interpolation operators if
Notice that is analogous to the norm in the FD-FD coupling. Despite it not being diagonal, has nice properties including being symmetric and positive definite. For example, if the node distribution on Γ is equidistant with step length , and ℙ 1 is used, then is a banded matrix, Let be a discrete vector containing nodal values of a FE solution in ℎ . We denote = the vector contains values of at the nodes on Ω . We prove that the following FD-FE coupling SATs impose stability to (26), ( 1 1 ( )) − 2 ( 1 1 ( )) , and the terms SAT , SAT couple the diffusion term Other notations are defined in Section 2.4. The assembly of the above matrices is described in Larson and Bengzon [16], Dao [8]. Note that the minus signs coming with ( ) in ( ) and the third component of SAT are in order to make ( ) analogous to the FD operator ( 1 ⊗ ) since is a negative vector. ( 1 1 ( )) − 2 ( 1 1 ( )) , which is in the same form as the diagonal norm coupling (29). We show the relevance of the quadrature rule (18) with the SAT (34) in the Theorem 3.7. Strict stability for the FD-FE coupling is shown in the theorem below. Proof. A total energy estimate can be obtained by applying the energy method to both schemes The structure of this proof is similar to the proof of Theorem 3.2. We only focus on the key differences coming from FE. The interface energy estimate reads where , ′ come from the operators approximating the flux; comes from the operators approximating the diffusion term; , ′ , come from the corresponding SATs. The following estimate is derived from the SBP property of the discrete operators Property 3.4 is used to handle the transpose terms. Under Lemma 3.6, we have Thus, the choice of = = 1∕2, = = −1∕2 makes −( + ′ ) appear in the sum ( + ′ ). Notice that Lemma 3.6 is only necessary to show that the energy term R (̃ 1 ′ 1 ( )) coming from the discrete operator is canceled. Because R is symmetric, the energy terms are canceled with the same choice of , as long as the quadrature rule to compute R in the scheme and in the coupling SATs are the same. Similar to the proof of Theorem 3.2, the mixed terms containing the interpolation operators 2 , 2 are also canceled because of the interpolation property (31) and the discrete skew-symmetry (23). The coupling of the diffusion terms can be proved to be stable in a similar manner.
Due to the fact that contains only zeros on the rows which do not correspond to the West boundary, we have = . As long as = − , = − and + = −1 in , is canceled by . To show that the remaining terms in containing 2 , 2 also vanish, we need to show that Indeed, from the property (31) which is It it sufficient to conclude that the interface terms vanish ≡ 0 with the above choice of parameters. This result comes from the construction of the SATs. The interface treatment is thus stable and conservative.

Remark 3.2.
Similar to the FD-FD coupling, the use of SBP-preserving interpolation operators used for FD-FE can also be shown to reserve the full order of accuracy. It is known that when narrow stencil FD operators of th-order are used, the accuracy at boundary closures is of ∕2th-order resulting in ( ∕2 + 1)th-order of overall accuracy. However, since the solution space ℙ 1 is currently used for FE, we expect that with roughly the same number of degrees of freedom of each method, the overall accuracy is bounded by second-order even if higher-order FD operators are used.

Numerical examples
This section includes numerical experiments of the proposed coupling technique on the linear advection-diffusion equation, and a nonlinear viscous Burgers' equation.

The FD-FE coupling
A FD-FE semi-discretization of the coupling problem (24) under the governing equation (36) is given by where , , are obtained by limiting the integration of to Γ , Γ , Γ , respectively. The coupling penalty terms are given by where = − , = − , + = −1.

Eigenvalue analysis for validation of the interface treatment
In this subsection, we show a numerical validation of Theorem 3.7. Consider the model problem (24) = (1, 0) . Figure 4 shows another example using 41 2 grid points on the FD domain, and an unstructured mesh using Delaunay triangulation with maximum element diameter ℎ = 0.2 on the FE domain (nonmatching coupling with 41 FD grid points, and 51 FE nodal points on the interface). The reference solution, is initialized at ( 0 , 0 ) = (−1, 0) and transported from left to right by setting = (1, 0) . The diffusion parameter is set to zero = 0. The other parameters are 1 = 1, 2 = 0.2. Sixth-order FD operators and ℙ 1 finite elements are used. The solution is examined at = 0 and = 1.  A semi-discrete problem to (1) can be written in the form where is a square real-valued matrix. In Figure 5, the eigenvalues of are plotted for several choices of the penalty parameters imposing the interface conditions in Theorem 3.7. The horizontal axis indicates the real part of the eigenvalues. The imaginary axis is vertical. Figure 5(a) and 5(c) are with two different sets of penalty parameters taking into account a diffusion term = 0.01 both satisfying Theorem 3.7 while the parameters used in Figure 5(b) and 5(d) violate the theorem. It can be clearly seen that in the former cases, all the yielding eigenvalues have nonpositive real part. Therefore, with some explicit method to integrate in time, for instance Runge-Kutta 4, and a sufficiently small time step, the fully discrete scheme is stable.

Convergence study
We use the following two dimensional Gaussian as reference solution The solution satisfies the boundary condition (4) with zero boundary data = 0 if the initial solution is small enough at the inflow boundary. For convergence results in this section, we choose = 0.005. For the inviscid problem, we set = 0. The 2D Gaussian (42) is used as the reference solution. We set 1 = 1, 2 = 0.2. The 2 -error in each FD solution part and in each FE solution part are respectively calculated by where is the exact solutions at the time of evaluation. The convergence rate in each method block is calculated as for ℎ 1 ≠ ℎ 2 being two different mesh/grid sizes. Since multiple blocks with different methods are involved, the global convergence rate is calculated based on the sum of all 2 -errors, and the scale of mesh/grid sizes, e.g., ℎ 2 ∕ℎ 1 = 2 when refining all the meshes/grids of size ℎ 1 twice. We use the FD diagonal norm SBP operators formulated in Mattsson and Nordström [23]. To demonstrate a proposed technique of constructing nondiagonal norm interpolation operators, the blocks are coupled by the 'glue grid' FD-FD interpolation operators Kozdon and Wilcox [14], and our proposed procedure Appendix A.2.
We choose the multiblock setup in Figure 6(a) to investigate convergence with an intention of examining possible degradation at corner points where the normal vector is not defined. Tables 1, 2 In Tables 1 and 2, is equal to both grid size in the FD domain ( = = ), and the number of discretization points along each axis of the FE mesh (matching interface). In Table 3 and 4, we perform the same study but when there are twice as many FE mesh points as FD grid points being located on the interfaces. The obtained rates are expected between second-order FE approximation and fourth-order FD approximation. Order degradation by the edges/corner point is not observed. Simpler experiments involving equal distribution of the two methods can be found in Dao [8]. There one can see that the global convergence is bounded by the lowest second-order rate in ℙ 1 FE regions.

Superconvergence preservation
For the convergence results, we use regular meshes for the FE discretization, illustrated in Figure 7,   We test propagating the solution by both direction, from the FE domain to the FD domain ( 0 , 0 ) = (1, 0) , = (−1, 0) , and in the opposite direction ( 0 , 0 ) = (−1, 0) , = (1, 0) . The solution is evaluated at = 2. The convergence rates in fourth-order FD cases reported in Table 5 are noticeably close to a fourth-order-fourth-order coupling, cf., FD-FD fourth-order case in Mattsson and Carpenter [22]. We observe the so-called superconvergence phenomenon in FE approximations, see e.g., Andreev and Lazarov [2]. At some special points, the order of convergence is higher than the convergence rate of the total error if calculated analytically. Those points are called superconvergence points, and they happen to be the nodal points of the regular mesh Figure 7 with the advection operator ⋅ ∇ . In the viscous case, there also exist such points which give rise to this behavior. It turns out that (44) is not the most correct way to assess the discrete error in the sense of mimicking the continuous measurement ‖ ‖ 2 (Ω) = √ ∫ Ω 2 . A straightforward fix to avoid this problem is calculating the error with many more sample points added to the computation mesh. From a different angle, it cannot be denied that (44) is desirable once our FE approximation fits in the SBP framework. (Simply, the defined discrete norm was applied.) To that extent, interpreting the advection operator as fourth-order and the diffusion operator as second-order may explain all the rates obtained. The authors have done separate numerical experiments on the FE approximation to verify that the current advection operator is fourthorder on the regular mesh regarding the error assessment (44). This is strong evidence that the proposed coupling technique preserves the native properties of the existing schemes.

Nonlinear Burgers' Equation in 2D
We consider the following nonlinear problem in two spatial dimensions to (7) gives 1 = 2 = 2 3 . From (4), a set of stable boundary condition reads where = ( 1 , 2 ) > 0 is a constant vector and is a real constant.

The FD-FD coupling
A FD-FD semi-discretization of the coupling problem (24) under the governing equation (46) is given by where the interface SATs are

The FD-FE coupling
A FD-FE semi-discretization of the coupling problem (24) under the governing equation (46) is given by where R , R , R are obtained by limiting the integration of R to Γ , Γ , Γ , respectively. The coupling penalty terms are given by

Convergence study
An analytic solution of (46) subject to the boundary condition (47) is  Tables 6 and 7 are same to the ones used in Table 1. In both tables, we can observe the expected convergence rate of second-order. It can also be seen that the use of higher-order FD operators slightly lowers the error and improves convergence rates.

Conclusion
The multiblock coupling of FD schemes and FE schemes in this paper is presented within the SBP framework. The technique uses a SAT technique to impose the interface continuity weakly. SBP-preserving interpolation operators resolve the nonconforming interface distribution. Accuracy, stability, and conservation are investigated on general scalar conservation laws. Numerical results are shown for a linear hyperbolic/parabolic problem and a nonlinear problem.
Future works may concern extensions for curvilinear grids, optimal time-stepping, an option for upwind operators [17], high-order finite elements or a solution to the case when block corners are unmatching. The work can be structurally extended for nonlinear systems.

A. Appendix: Construction of the nondiagonal norm interpolation operators
We show two ways of constructing the necessary interpolation operators in the interface SATs in this section. The first approach is similar to the one introduced in Mattsson and Carpenter [22]. The second way requires constructing a pair of diagonal norm interpolation matrices.

FD
FE FE lumped Figure 8: Construction of nondiagonal norm interpolation operators using a "bridge" diagonal norm layer.
Example A.1. 2 order FD -ℙ 1 FE matching uniformly distributed interface Recall the notations , being the one-directional norms in the -direction of the FE scheme and the FD scheme, respectively, ℎ being the step size, "L", "R", "M" being the FD, FE layers, and the lumped FE layer. We have For the interpolation between FE lumped and FE, we use 2 being the identity matrix and Theorem A.1. Assume that the diagonal norm interpolation operators 2 , 2 fulfill the acuracy and SBP requirements (27), (28) Since the lumped ℙ 1 mass matrix corresponds to a second-order SBP FD diagonal norm, and the corresponding interpolation operators are of the form in Example A.1 reassembling the Simpson's rule naturally hold the properties (28), (27) to the desired accuracy, and (52). By the relation (53), Theorem A.1 always holds true.