Provably non-stiff implementation of weak coupling conditions for hyperbolic problems

In the context of coupling hyperbolic problems, the maximum stable time step of an explicit numerical scheme may depend on the design of the coupling procedure. If this is the case, the coupling procedure is sensitive to changes in model parameters independent of the Courant–Friedrichs–Levy condition. This sensitivity can cause artificial stiffness that degrades the performance of a numerical scheme. To overcome this problem, we present a systematic and general procedure for weakly imposing coupling conditions via penalty terms in a provably non-stiff manner. The procedure can be used to construct both energy conservative and dissipative couplings, and the user is given control over the amount of dissipation desired. The resulting formulation is simple to implement and dual consistent. The penalty coefficients take the form of projection matrices based on the coupling conditions. Numerical experiments demonstrate that this procedure results in both optimal spectral radii and superconvergent linear functionals.

electrodynamics, computational geophysics, and computational fluid dynamics. Typical examples include the Euler equations coupled to the elastic wave equation, and electromagnetic wave interactions with scattering across material interfaces. A difficulty that can arise when developing numerical methods for these types of coupled problems is that the resulting discretization may suffer from stiffness. In such a case, an explicit time stepping procedure may require a significantly smaller time step than what is expected from each separate subproblem. For example, within fluid-structure interaction, the impact of the coupling procedure is a well-known problem. The socalled added mass instability is an effect of coupling light fluids to thin, dense structures immersed in fluid. This problem has led to a plethora of different coupling strategies, e.g., see [2,6,32]. Potential stiffness issues related to large density contrasts also appear in linear settings when coupling acoustic and elastic wave equations. Stiffness can also appear in other settings, e.g., for frictional sliding interfaces [17], and narrow cracks [24].
The imposition of coupling conditions is closely related to the imposition of boundary conditions where stiffness can occur for similar reasons. In this paper, we will for ease of presentation switch between these two use cases. We focus on weakly imposed boundary and coupling conditions by including additional source terms in the governing equations, so-called penalty terms. The role of the penalty terms is to force the numerical solution on the interface towards satisfying the actual coupling conditions [4,30]. The penalty term can be interpreted as specific weightings of the coupling conditions in combination with the equations.
There is some flexibility in choosing penalty weights, and the choice can have a striking impact on a scheme's performance. Well-designed penalty terms can be the difference between having a scheme that works in practice or a suboptimal scheme that runs orders of magnitude slower. Unfortunately, it is quite easy to construct suboptimal penalty terms that cause these problems, even though they are semi-discretely stable. Coupling procedures that overcome these problems use either upwind numerical fluxes based on solving the Riemann problem [34], or a specific characteristic treatment [17,20,21]. In each of the above cases, the penalty terms introduce some amount of potentially unwanted additional dissipation.
Adding certain types of artificial dissipation may result in improved accuracy [16,21]. However, dissipation may be suboptimal for wave propagation problems since that may prevent the use of symplectic or staggered time-stepping methods that have improved accuracy and stability properties [33]. For instance, a benefit of the fourth order staggered Runge-Kutta time-stepping scheme introduced in [11] is that it has a factor 16 smaller truncation error and a factor of two larger stability region along the imaginary axis compared to the classical fourth order Runge-Kutta scheme. Unfortunately, the staggered Runge-Kutta scheme is unstable even for a small amount of artificial dissipation [11,33]. An energy conservative weak coupling procedure for coupling acoustic materials across nonconforming grids is presented in [8]. However, this coupling procedure causes stiffness for problems with large density contrasts. To the best of our knowledge, there are no known penalty procedures that are both energy conservative and non-stiff. In this paper, we present a systematic procedure for constructing both energy conservative and dissipative penalty terms that are provably non-stiff.
We argue that numerical schemes should be provably stable and non-stiff. A powerful method for proving semi-discrete stability of a scheme is the energy method [13]. This method give sufficient conditions for semi-discrete stability, but not necessary ones. However, when the energy method fails, it is usually an indication that there is an underlying instability present that may be triggered under certain circumstances. Without a proof of stability, there are no guarantees that the output produced by a numerical scheme is reliable. A proof of stability also helps to rule out bugs and can serve as a strict test of the implementation. For instance, one can randomly initialize all solution fields and model parameters. Afterwards, one numerically computes the discrete energy rate in the implementation and asserts that it is nonpositive.
Similar arguments can be made for having a provably non-stiff scheme. A proof of non-stiffness gives a sufficient condition in order to prevent stiffness from occurring due to sensitivity in model parameters. We will show that if the matrix norm of each penalty matrix can be bounded by the CFL condition, it is provably non-stiff. Since this is only a sufficient condition, a lack of proof does not imply that the scheme is stiff. However, failure to bound the matrix norm could be an indication that the numerical scheme is sensitive to certain parameter values. A proof of non-stiffness also helps to rule out bugs in the implementation. For instance, if the simulation of a provably stable and non-stiff scheme becomes unstable due to a change in a model parameter independent of the CFL-condition, it must be caused by an incorrect implementation. In addition, since the proof of non-stiffness constrains the parameters to some degree, this analysis can limit the potential parameter space that one must probe to select optimal parameter values.
The key developments in this work are applicable to any numerical scheme in Summation-By-Parts (SBP) form that weakly imposes boundary or coupling conditions via Simultaneous-Approximation-Terms (SAT). Numerical schemes that are in SBP-SAT form lead to a proof of stability via the energy method and apply to finite difference [5,19,28], finite volume [22], spectral collocation [3], discontinuous Galerkin [9], and flux reconstruction schemes [26]. See [7,30] for reviews. When developing SBP-SAT schemes, one can focus on formulating and analyzing the continuous problem. Once this analysis is complete one can proceed by discretizing the continuous formulation and obtain a provably stable scheme without much difficulty, as shown in [20]. This aspect allows us to construct penalty terms in the continuous setting, leaving out specific implementation details for the variety of numerical schemes that belong to the SBP-SAT family. Hence, the focus of this paper is on constructing and analyzing penalty terms in a continuous setting without placing emphasis on particular discretization or implementation details. This paper is organized as follows. We begin in Sect. 2 by discussing a motivating example: the coupling of wave equations. This example demonstrates how naively selected penalty weights cause stiffness. Section 3 presents the general theory for constructing non-stiff penalty terms. In particular, it is shown that the penalty terms are connected to a projection matrix. Section 4 shows that penalty terms constructed via a projection matrix formula automatically result in a dual consistent formulation. Section 5 introduces a diagnostic test to flag penalty formulations suffering from stiffness, and shows that certain penalty terms obtained by the projection matrix formula are provably non-stiff. Section 6 revisits the motivating example and applies the general theory to construct both non-stiff, energy conservative, and energy dissipative penalty terms. Section 7 introduces numerical experiments that exemplify our theoretical developments. In particular, we present a challenging 2D air-water interface problem for the wave equation. Finally, we discuss our findings in Sect. 8.

A motivating example
Artificial stiffness can easily emerge in coupled problems [10,17]. In the decoupled case, each problem has some maximum stable time step Δt (i) max , for i = 1, 2 when using a specific explicit time integrator. On the other hand, if the two subdomains are coupled together, one might obtain a different Δt max that can be much larger than either Δt (1) max or Δt (2) max . This phenomenon stems from a suboptimal selection of penalty weights in the coupling procedure. We demonstrate this problem by coupling two wave equations at an interface Γ (see Fig. 1).
The wave equation for each side of the interface is written as a first order hyperbolic system of equations in Cartesian coordinates, where The vector q (i) collects the pressure p (i) and the components of the velocity field v (i) = (v (i) x , v (i) y ), decomposed with respect to the Cartesian basis. Moreover, The material parameters ρ 1 > 0 and ρ 2 > 0 are the respective densities on each side of the interface. Similarly, c 1 and c 2 are the respective wave speeds. The forcing/penalty terms f (i) are responsible for weakly imposing coupling conditions in each subdomain.

Penalty terms
A general formulation of the penalty terms is where u = (q (1) , q (2) ). The role of the lifting operator L is to make sure that the penalty term acts on the boundary only [1,27]. In one dimension, it is analogous to the Dirac delta function δ(x − x 0 ), where x 0 is the boundary point. In higher dimensions, it is defined by for smooth and vector-valued functions u and v, and where ∂Ω is an appropriate boundary part of the domain Ω. This relationship is similar to the divergence theorem in that it relates a volume/surface integral (left-hand side) to a surface/line integral (right-hand side). The matrix L T is the boundary operator, determined by the coupling conditions, and Σ (i) are penalty matrices. The penalty matrices describe how to form linear combinations of the coupling conditions for each equation.

Well-posedness
For the problem to be well-posed, there are three general requirements. See [20] for further details. For existence, there cannot be too many boundary conditions, for uniqueness there cannot be too few boundary conditions, and for a bound, the boundary condition must have appropriate form. In this work, we weakly impose the minimum number of appropriate boundary/coupling conditions that bound the solution.

Coupling conditions
One set of well-posed and physically relevant coupling conditions for the wave equation are In (4), n = (n x , n y ) is the unit normal on the interface Γ , directed towards Ω (2) . The coupling conditions (4) can be expressed in matrix vector form L T u, as done in (2), where and v (1) (2) · n. As previously mentioned, the penalty terms are formed by taking a linear combination of the coupling conditions (4) and distributing them among the equations, yielding

The energy method
The weights in the penalty matrices Σ (1) and Σ (2) are determined by bounding the solution via the energy method [13].
Let u = Ω u T udΩ denote the L 2 -norm. By differentiating q (i) 2 with respect to t and substituting ∂ t q (i) using (1) it follows that By applying Gauss's theorem and the definition of the lifting operator (3), the volume integral in (7) can be converted into the following surface integral, x n x + A (1) y n y ,Ã (2) = −A (2) x n x − A (2) y n y . Equation (8) shows that the energy q (i) is conserved in the interior, and changes in the energy rate are purely determined by fluxes across the interface. Note that in (8), the energy rate contribution from exterior boundaries have been neglected.

Proposition 1
The solution of the problem (1) is bounded and the energy (9) is conserved if the penalty matrices in (6) are chosen as where k i and l i are real parameters that satisfy the constraints Proof A direct calculation yields, u TÂ u = α(q (1) ) TÃ(1) q (1) + β(q (2) ) TÃ (2) (2) n .

Stiffness
Proposition 1 states sufficient conditions for obtaining an energy conservative coupling. However, the constraints (12) imply that there is one free parameter. The problem we deal with in this paper is: how should this parameter be chosen? It might be tempting to try to make all parameters share the same magnitude. For instance, as was done in [8], one can choose However, as we will see next, this choice causes stiffness (as do many other choices).
Before we proceed, we need to discuss the spatial discretization of the coupled wave equation problem. The discretization of (1) in space leads to linear system of ordinary differential equations: where p and v discretize the pressure and velocity fields, and M h is a sparse matrix that holds the discretization of the spatial operators and all penalty terms; applied to each side of the interface (see Sect. 5.2 for more details). Since we are only interested in understanding how the penalty parameters influence Δt max , the specific spatial discretization technique is not important (e.g., finite difference, finite volume, discontinuous Galerkin, etc.). The maximum stable time step of an explicit time integrator can be estimated as In (17), h is the grid spacing (taken to be constant and the same in each subdomain), and c max = max(c 1 , c 2 ). To prevent stiffness from occurring, Δt max must be no smaller than either Δt (1) max and Δt (2) max , associated with the respective uncoupled cases. The naive parameter choice (15) violates this condition. For this choice, the spectral radius of M h , denoted (M h ), depends on the density ratio and that causes stiffness (see Fig. 2). For example, if ρ 1 /ρ 2 ≈ 100 then Δt max must be reduced by a factor of ≈ 10 compared to min Δt (1) max , Δt (2) max . Our ambition in this paper is to find a procedure that removes the kink in Fig. 2 and minimizes the spectral radius.

General treatment
As made apparent by the previous example, the design of penalty terms can strongly influence the efficacy of the numerical scheme. One option is to perform numerical experiments to find optimal penalty weights, but that can be cumbersome and imprac- Fig. 2 The maximum stable time step Δt max as a function of the density ratio ρ 1 /ρ 2 (on a log-log scale). This plot was generated by discretizing the wave equation in (18) using c = c 1 = c 2 , fourth order SBP finite difference operators [28], and 100 grid points in each subdomain tical. In this section, we present a general approach that can automatically leads to non-stiff penalty weights.
We again consider the problem (1) but for simplicity focus on one subdomain, and now with general and symmetric A x and A y matrices: Since we are only considering one subdomain, all superscripts have been dropped. Our procedure for deriving penalty terms is inspired by the procedure developed in [17,18]. We start with a modified ansatz of the penalty term, In this formulation, δq = δq(q) is an unknown penalty vector that depends on q such that a minimum number of boundary conditions are weakly imposed.

The energy method
Next, we introduce the bilinear form Then, the energy method applied to (18) can now be written as We can reformulate the above by completing the square, yielding The solution is bounded if the right-hand side of (21) is non-positive. However, here we enforce the stronger requirement that each term on the right-hand side is individually non-positive (not just their sum): We will see that the two conditions in (22) uniquely determine δq.

Diagonalization
Since A is symmetric, it is diagonalizable by a complete set of orthogonal eigenvectors, Here, X is a matrix that contains the eigenvectors arranged as column vectors, and Λ is diagonal matrix containing the corresponding eigenvalues. The eigenvalue matrix Λ is split into positive Λ + > 0 and negative parts Λ − < 0, and X = [X + , X − ]. Without loss of generality, we assume that A contains no zero eigenvalues (if it does, they are ignored). We also define the transformation w = √ |Λ|X T q, leading to This way of defining w + and w − simplifies the upcoming presentation. For the other variable, δq, we define similar transforms, i.e., When the eigendecomposition is applied to the quadratic form Φ(q, q), we find

Boundary conditions
Recall that the goal is to determine the penalty vector δq such that (22) holds subject to the boundary conditions.

Strong boundary conditions
For strong boundary conditions, it is natural to impose (see [20]) where R is a rectangular matrix that imposes a minimum number of boundary conditions. From (25), we get For strong boundary conditions, δq = 0, the energy rate (21) becomes It follows that the right-hand side in (27) is non-positive if R T R ≤ I + , where I + is the identity matrix related to the size of w + . Clearly, this condition is satisfied when R is sufficiently small.

The weak primal condition
To satisfy Φ(q − δq, q − δq) ≥ 0 in (22), we proceed in the same way and specify If R T R ≤ I + , we obtain

The weak dual condition
As seen in (22), we also need Φ(δq, δq) ≤ 0. To obtain that, we specify where δ R is an another rectangular matrix.
Interestingly, the condition (δ R) T (δ R) ≤ I − corresponds to the well-posedness condition for the dual problem in (18). This is further discussed in Appendix 4.

Well-posedness
The conditions (28) and (30) determine δq uniquely under the conditions given in the following Proposition.
Proof The starting point is the eigenvector transformation (24) . Since X T = X −1 , the inverse transformation is Inserting (30) into the above yields Also, inserting (30) into (28) yields Inserting this final expression into δq results in the formula (32).
Proposition 2 shows that δq is formulated in terms of the boundary condition w − − Rw + . Hence, this formulation imposes a minimum number of boundary conditions, as required for well-posedness, and we have proven the following result.

Projection matrix formula
There is an alternative procedure to determine δq, obtained via δq = Pq, where P is a projection matrix. By definition, a projection matrix satisfies the idempotency property P = P P.

Proposition 4 Let
Then δq = Pq, where P is the projection matrix Proof We start by deriving the formula for P and then show that this matrix is a projection matrix. Note that δq in (32) can be written as Hence, it remains to show that L T δL = I − − Rδ R. The result immediately follows from eigenvector orthogonality: For P to be a projection matrix it must satisfy P = P P. We have By applying Proposition 4, the penalty term (19) can be formulated as

The dual problem
We show how the general formulation presented in the previous section is connected to the dual problem.

The strong dual problem
Define the functional The dual of the primal problem (18) where f * is the PDE describing the dual problem (to be determined). By inserting (18) into J (q * , f ) and applying integration by parts in time and space, we get where Φ(q, q * ) is defined in (20), Ψ (q, q * ) = q T q * T 0 , and The terms Φ and Ψ arising from integration by parts must vanish for the dual problem to be well-defined. The term Ψ vanishes due to the initial condition u(x, 0) = 0 and end condition u * (x, T ) = 0. Next, we derive the dual boundary conditions. Following Sect. 3.2, by diagonalizing and splitting A into positive and negative parts, we get where w + * = (X + ) T q * , and w − * = (X − ) T q * . The strong primal boundary condition (26) leads to Hence, to make Φ vanish, the dual boundary conditions must be By introducing the time reversal transformation τ = T − t, the dual problem becomes subject to the boundary conditions (38). Proof The energy rate for the dual problem (39) with f * = 0 is By imposing the strong dual boundary conditions (38), we get Hence, the dual problem is well-posed if R R T ≤ I − . On the other hand, the primal problem is well-posed if R T R ≤ I + (see Proposition 3). We can show that the two conditions are the same using the singular value decomposition (SVD) R = U Σ V T , where and U and V are real orthogonal matrices, and Σ is a diagonal matrix. By Thus, both the primal and dual problems are well-posed if Σ ≤ I .

The weak dual problem
Since the analysis of the weak dual problem is very similar to the weak primal problem, some details are omitted. To weakly impose the boundary conditions (38), we set the right-hand side of (39) to In (40), δq * is an unknown penalty vector. This vector is determined in the same manner as the penalty vector δq of the primal problem. When the dual problem (39) uses the penalty term (40), the energy method leads to To bound the solution, we require First, consider Φ(q * − δq * , q * − δq * ) and specify By inserting (42) into (37) and requiring R R T ≤ I − , we get Next, consider Φ(δq * , δq * ) and specify By inserting (43) into (37) and requiring (δ The results established thus far are summarized in the following proposition.

Dual consistency
To achieve superconvergent linear functionals, it is important that the weak primal and dual problems are related to each other via an appropriate choice of penalty terms, f and f * [21]. By substituting f and f * in (36) using the penalty terms (19) and (40), we get To arrive at (45), we have neglected the integrals over t and set the initial and end conditions to zero. As shown in the following proposition, (45) establishes a relationship between the penalty parameters δ R and δ R * that is analogous to the relationship between the primal (26) and dual (38) boundary conditions.

Proposition 7
The weak primal and dual problems (see Proposition 3 and 6) satisfy Proof See Appendix A for proof.

Bounded matrix norms
For spatial discretizations with all eigenvalues located in the correct half plane, the maximum stable time step is determined by the stability region of the particular time stepping scheme. For the spatial discretization to be robust, it must handle any choice of model parameters without causing unwanted growth of the spectral radius, which can be costly. To avoid computing the spectral radius, our goal in this section is to a priori analyze how the matrix norm will scale (up to a constant factor) for the spatial discretization of a particular model problem. Since this analysis relies on symbolic computation, it can identify model parameters that cause suboptimal scaling of the matrix norm. If the scaling is the same as the expected CFL-condition obtained from the uncoupled problem, we say that discretization is non-stiff. This definition will be made precise soon. In more detail, we focus on the following objectives: 1. To develop an easy-to-use diagnostic test that indicates how the penalty parameters influence the matrix norm of the spatial discretization. 2. To prove that the projection matrix formula (35) results in penalty terms that are non-stiff, in the sense discussed below.

Matrix analysis
Before proceeding, we briefly summarize key concepts and results from matrix analysis. For a matrix A ∈ R m×m , the 2-norm is defined by the induced vector norm where

The Frobenius norm is
The 2-norm and Frobenius norm satisfy the following properties: where ⊗ is the Kronecker product, and (A) is the spectral radius of A. See e.g., Golub and Van Loan [12] for further properties and proofs of (49)-(51).

Model problem
As the model problem, we consider the 1D hyperbolic PDE and discretize it in space with SBP operators, Here , where x j = jh and h is the grid spacing for an equidistant grid. The matrix D x is defined by D The matrix Q x satisfies the summation-by-parts property: This property is key for proving energy stability of a SBP-based scheme. For further details about SBP-based discretizations, see [7,30], and the references therein.
To simplify the forthcoming analysis, we assume that we only need to apply a penalty term on the left boundary, such that L h in (53) becomes Moreover, we assume that D x is a consistent approximation of the first derivative such that it differentiates a constant exactly: Before we can prove that the discretization (53) is non-stiff for appropriately chosen penalty parameters, we give the following definition of non-stiffness.

Definition 2
The semi-discrete approximation (53) of (52) is non-stiff if there exists a positive constant γ , independent of the model parameters, such that

Diagnostic test
Our objective is to decide how h M h F (or the 2-norm) scales for a particular choice of penalty terms without performing any numerical computations. This capability is useful because it enables one to perform a diagnostic test and assess the robustness of a particular weak coupling procedure without having to implement it and compute the spectral radius for a range of model parameters. If the test passes, then the weak coupling procedure is provably non-stiff. If the test fails, it indicates that the coupling procedure is sensitive to certain choices of model parameters. Our test consists of the following four steps.
1. Determine Σ such that the problem is energy stable and the penalty parameters σ i j are parameterized by the model parameters (e.g., a i j in A x , and l i j in L).

Compute the spectral radius (A x )
3. Compute the norm of the penalty matrix Σ L T F . 4. Check if there exists a constant γ , independent of the model parameters, such that: Σ L T F ≤ γρ(A x ). 5. If the test fails, repeat it for the 2-norm. 6. If the test still fails, redesign the coupling procedure to be provably non-stiff.
We recommend checking the Frobenius norm first, because it is much easier to compute than the 2-norm. To compute the 2-norm one must symbolically compute the singular values of the penalty matrix.
The diagnostic test is motivated by the following assumption and Proposition.
Assumption 1 For finite sized (fixed N ) matrices, D x and L h , there exists constants s d,l , r d,l > 0 such that

Proposition 8
The matrix M of the spatial discretization (53) satisfies Proof See Appendix B for proof.

Remark 1
If there are multiple penalty terms present in M h , we use the reverse triangle inequality on each term to obtain a lower bound, i.e., Likewise, we use the triangle inequality to obtain an upper bound.

Suboptimal scaling
As an illustration, consider the motivating example presented in Sect. 2.
1. From (6), the penalty matrix on each side of the interface is 2. The spectral radius is (A

The Frobenius norm of each penalty matrix is
4. For the naive choice of parameters (15), we get The presence of ρ 2 /ρ 1 causes linear growth in h M F for ρ 2 /ρ 1 1 and ρ 2 /ρ 1 1. 5. Since we have failed to bound the penalty matrices in the Frobenius norm, we proceed to check the 2-norm: The presence of ρ 1 /ρ 2 also causes linear growth in h M 2 . Hence, we cannot rule out that the spectral radius scales suboptimally as a function of the model parameters. In fact, the numerical experiment presented in Fig. 2 shows unwanted growth of the spectral radius in ρ 1 /ρ 2 . 6. Since we have failed to prove that this coupling is non-stiff, we will redesign it using the projection matrix formula.

Provably non-stiff penalty treatment
As is shown in the following proposition, the projection matrix formula yields a useful estimate for well-posed linear hyperbolic problems. For certain choices of the penalty parameters δ R, the method is provably non-stiff.

Proposition 9
Let Σ L T = A P, where A is given in (19) and P is computed using the projection matrix formula (35). Then Proof See Appendix C for proof.
If the denominator I − − Rδ R 2 approaches zero, the norm of h M 2 scales suboptimally. There are at least two choices of δ R that prevent the denominator from vanishing for all R. For the energy dissipative choice δ R = 0, (56) yields For the energy conservative choice δ R = −R T , (56) yields In each case, we used R 2 ≤ 1, which is required by well-posedness.

Revisiting the motivating example
To demonstrate how to derive dual consistent and non-stiff penalty terms by example, we apply it to the motivating example presented in Sect. 2. We summarize our procedure and carry out each step at a time.
1. (The penalty formulation) Formulate the penalty terms such that the energy rate can be written similarly to (21). Since we are studying the coupled problem (1), we need to scale the energy rate as done in (9), where u = (q (1) , q (2) ), A = diag(αÃ (1) , βÃ (2) ), and δu is to be determined.

(Projection matrix formula)
Determine δu = Pu using the projection matrix formula defined in Proposition 4. 6. (Matrix norms) Following Sect. 5, compute the matrix norms AP/α F and AP/β F to verify that the formulation is provably non-stiff.
5. We use the projection matrix formula in Proposition 4 for δ R = −R T , to obtain the penalty vector After inserting (60) into (57), the final penalty terms become where 6. Following the steps laid out in the diagnostic test in Sect. 5, we arrive at since ρ 1 > 0, ρ 2 > 0. Hence, this formulation is provably non-stiff in the sense of Definition 2. Alternately, if we repeat steps 1-5 for the energy dissipative choice δ R = 0, the penalty terms become

General penalty term formulation
By inspecting (61)-(62) we can formulate the following general penalty terms.
In (63), ζ controls the amount of energy dissipation, and k 1 , k 2 are listed in Proposition 1, and they are here defined as where m is a free parameter. This formulation recovers all other formulations.
To illustrate the results, we revisit the spectral radius investigation presented in Sect. 2.4. The energy conservative and dissipative choices (66)-(67) each cause the spectral radius to be independent of the density ratio, as shown in Fig. 3.

Numerical experiments
In this section, we conduct several numerical experiments. These experiments demonstate the importance of having a scheme with bounded matrix norm for robustness, and show the predictive ability of Propositions 8-9. In addition, we show that linear functionals are superconverging. Finally, we investigate advantages and disadvantages of developing a near-energy conservative scheme (in both space and time) of the wave equation in two spatial dimensions.

Matrix norm behavior
Consider the wave equation in a single domain on the unit interval 0 ≤ x ≤ 1, subject to the boundary conditions The left boundary condition (68) is implemented by taking R = 0 and δ R = 0 in the projection matrix formula (Proposition 4). By Proposition 3, the problem is well-posed and energy dissipative. We consider two different well-posed implementations of the right boundary condition, described next.

Naive choice
The right boundary condition can be implemented naively by applying the penalty term (1, t). The penalty matrix is chosen as Σ = [0, 1/ρ] T , which results in the energy dissipation rate By Proposition 8, this penalty term causes supoptimal scaling, since: Hence, for α 1 the growth rate of h M h is linear in α. We confirm this prediction by numerically computing M h 2 as well as lower and upper bounds on it by applying Proposition 8. We discretize in space by applying the fourth order SBP operators [28] using 100 grid points. The constants in Proposition 8 are numerically computed to be: Figure 4 show that the norm h M h 2 grows linearly as predicted when α increases.

Projection matrix formula
The right boundary condition (69) can also be implemented by applying the projection matrix formula (Proposition 4). In this case, R = (α − 1)/(1 + α), whereas δ R is a free, but bounded parameter. We investigate how δ R influences h M h 2 for a fixed α = 10 and vary δ R within its stability limit −1 ≤ δ R ≤ 1. Figure 5 shows that any

Superconvergence and dual consistency
SBP-SAT schemes that are dual consistent exhibit superconvergence for linear functionals. As explained in Sect. 4, a necessary condition for dual consistency is that the weak primal boundary conditions are related to the weak dual boundary conditions in a certain way. Interestingly, Proposition 7, shows that all penalty terms constructed using the projection matrix formula satisfy this condition. We numerically verify this result by constructing both energy conservative and energy dissipative penalty treatments that result in superconvergent linear functionals.
Consider again the wave equation in a single domain 0 ≤ x ≤ 1. We use the method of manufactured solutions to construct the solution p(x, t) = cos(kt) sin(kx), v(x, t) = − sin(kt) cos(kx).
The manufactured solution (70) satisfies the wave equation subject to the initial and boundary conditions: p(0, t) = 0, v ( 1, t) = − sin(kt) cos(k). These boundary conditions correspond to R = −1 at both x = 0 and x = 1. The manufactured solution (70) results for example in the linear functionals As before, the boundary conditions are implemented using the projection matrix formula. We compare the energy dissipative choice δ R = 0 against the energy conservative choice δ R = −1. Besides this parameter choice, all other settings are the same. In (70), we use k = 8π . We advance in time using a fourth order Runge-Kutta method until the final time T = 1.2 using the timestep Δt = h/4 (the wave speed c is set to c = 1 in this case). This experiment uses fourth order SBP operators, and the expected convergence rate is third order in the variables p, and v according to [29,31]. The linear functionals (72) should superconverge at fourth order rate [14]. Tables 1  and 2 lists the errors and convergence rates for the energy dissipative and energy conservative penalty parameter choices, respectively. As expected, both choices result in superconverging functionals.

Energy conservation versus energy dissipation
For energy conservative problems, it may be beneficial to design a numerical scheme that conserves energy in a semi-discrete sense. Weak boundary conditions that have dissipation (γ > 0 in (63)) cause the spectrum of the spatial discretization to have eigenvalues with a real part. The amount of dissipation controls the distance of the real part of the eigenvalues to the imaginary axis. Since some Runge-Kutta methods have a stability region that is wider along the imaginary axis compared to the real axis, the presence of dissipation can have a negative impact on the maximum stable time step. Staggered time-stepping methods, such as leap-frog and the fourth order staggered Runge-Kutta method [11] have stability regions that are primarily confined to the imaginary axis. On the other hand, they have smaller truncation error and larger stability region along the imaginary axis compared to their non-staggered counterparts. For instance, the fourth order staggered Runge-Kutta method has about a factor of two larger stability region along the imaginary axis, and a factor of 16 smaller truncation error, compared to the classical fourth order Runge-Kutta scheme. However, to allow for dissipation may improve the accuracy of a numerical scheme. To understand how these benefits and drawbacks influence performance, we compare the computational efficiency of energy conserving discretizations advanced in time with the fourth order staggered Runge-Kutta (SRK4) scheme against energy dissipating discretizations advanced in time using the fourth order, 5-4 solution 3, 2N low-storage Runge-Kutta method (RK4) [15].
To investigate the performance and demonstrate the importance of how the coupling parameters are chosen, we present a challenging application problem featuring a lightdense media interface. The interface is located at y = 0 and the light medium rests on top of the dense medium. The material properties in the light medium are ρ 1 = 1, c 1 = 1, and they are ρ 2 = 800, c 2 = 4 in the dense medium. These parameter values result in an impedance contrast of Z 2 /Z 1 = 3200, which is representative of a water-air interface. This large impedance contrast makes the problem computationally challenging. Figure 6 shows a coarse grid simulation of the pressure field for a fixed time. An explosive source initiated in the air medium sends out waves that both reflect against and transmit across the air-water interface (black horizontal line).
The air and water media are each discretized by the fourth order SBP staggered finite difference operators presented in [23]. Initially, all fields are set to zero. To initiate the simulation, we use a singular source term with a prescribed source time function. This source term is written as δ(x − x s )δ(y − y s )g(t) and acts in the right-hand side of the pressure equation. The source is positioned in close vicinity of the interface, at (x s , y s ) = (0.0, 0.05). The delta functions δ(x − x s ) and δ(y − y s ) are discretized in a line by line manner to fourth order accuracy by imposing moment conditions, see [25] for details. As the source time function, g(t), we use the Ricker wavelet where f p = 6.4 and t 0 = 0.25. These settings result in λ min /h ≈ 6 grid points per minimum wavelength in the air medium. This estimate defines λ min = c min / f max , where f max is the frequency at 5% of the peak amplitude of the Ricker wavelet in the Fourier domain.
Since we are mostly interested in understanding how the coupling treatment influences accuracy at the interface, we measure the error along it. Due to the discontinuous nature of the weak coupling conditions, the error is defined by the jump in a quantity that is continuous across the interface. We choose the jump in pressure. To approximate the solution anywhere along the interface, we use cubic Lagrange interpolation. In particular, we focus on measuring the error in the middle of the interface, at (x r , y r ) = (0.0, 0.0).

Accuracy
Before investigating the computational efficiency of each method, we compare the accuracy of the energy conservative and energy dissipative couplings obtained by the projection matrix formula (35). When solving the air-water interface problem on the coarse grid with the energy conservative coupling (66), numerical artifacts emerge. These artifacts are large-amplitude spurious oscillations that propagate along the interface. Figure 7a shows the presence of these oscillations at the same time and on the same grid as shown in Fig. 6. While these spurious oscillations vanish with grid refinement and appear to be confined to the interface, they prevent one from performing coarse grid computations that demand a reasonably (less than 10% error) accurate solution at the interface. Fortunately, we can suppress these spurious oscillations by selecting a different δ R in the projection matrix formula.
To appreciate what causes the spurious oscillations to develop, recall that the interface problem couples a light medium to a dense medium. This problem has a large impedance contrast Z 2 /Z 1 1 across the interface. In this limit, we study the reflection and transmission of plane wave solutions that are incident normal to the interface. The following solution solves (1), wherek 1 ,k 2 are wavenumbers that satisfy the dispersion relation ω =k 1 c 1 =k 2 c 2 , and a R , a T are reflection and transmission coefficients, determined by the coupling conditions (4). The reflection and transmission coefficients are: In the limit when Z 2 /Z 1 → ∞, we get a R = −1, a T = 0, and there is no transmission into the dense medium. An interpretation of this limit is that a wave propagating in the light medium senses the interface as a rigid wall. On the other hand, a wave propagating in the dense medium senses the interface as a free surface, and transmits into the light medium with double the amplitude. If we neglect any transmission, the problem decouples into two problems with the boundary conditions v (1) The only way to implement these boundary conditions in an energy conserving manner is to set m → ∞ and ζ = 0, in (63)-(64). The choice (73) results in k 1 = 1, k 2 = 0, and (63) yields where v (2) n = p (1) = 0 in the decoupled case. In the limit when Z 2 /Z 1 → ∞, the energy conservative coupling (66) approaches the penalty terms of (74) at a linear rate. We speculate that this rate is too slow when there is no artificial dissipation present to suppress the spurious oscillations. Instead, we choose δ R so that it corresponds to (74) by directly taking the limit (we do not modify Z 2 , Z 1 in the problem itself) Obviously, if Z 1 > Z 2 we would have taken the opposite limit Z 2 /Z 1 → 0 instead. Note that the choice (75) preserves the energy conservation property because δ R is still an orthogonal matrix. Since we used the projection matrix formula to derive this coupling, it is also both dual consistent and provably non-stiff. This simple modification of the penalty parameter results in a dramatic improvement in accuracy of the coarse grid simulation (Fig. 7b). However, the dissipative choice δ R = 0 in Fig. 6 is the most accurate option. Figure 8 shows the error in pressure at the center grid point of the interface as a function of time for each type of coupling. The modified energy conservative coupling (73), SRK4-EC2, reduces the maximum amplitude of the error by a factor of 18 more than the coupling (66), SRK4-EC1. While this simple modification resulted in a dramatic improvement for this particular problem, we have not investigated other problems with moderate impedance contrasts.

Time step selection
We determine the maximum stable time step Δt max = CFLh/c max for each type of coupling and time stepping method. In each case, we maximize CFL by performing several coarse grid computations using the bisection method. Table 3 lists CFL numbers for each method up to three decimal places.
As expected from the 1D investigations, the naive coupling, RK4-NA, results in more than an order magnitude reduction in ν compared to the other couplings. Interestingly, SRK4-EC1 and SRK4-EC2 result in a factor of four larger CFL compared  Table 3 for label descriptions to RK4-ED. We believe this improvement is due to two reasons: 1. SRK4 has a factor of ≈ 2 larger stability region along the imaginary axis compared to RK4. 2. The dissipative part of (67) leads to eigenvalues with real parts that are ≈ 2 smaller than the RK4 stability envelope along the imaginary axis. By decreasing ζ in (63) it should be possible to increase CFL, but most likely at the expense of decreasing the accuracy in the solution.

Computational efficiency
With each method tuned to run at its time step stability limit, we compare their computational efficiency (error as a function of computational time). The norm of the error is the maximum amplitude of the error in time in the middle of the interface (see Fig. 8). Figure 9 shows the computational efficiency of each method. The naive method, RK4-NA, performs the worst, running more than an order of magnitude slower than the other methods. The computational efficiency of SRK4-EC is limited by the large errors caused by the inaccurate coupling. On coarser grids, SRK4-EC2 performs significantly faster than RK4, but the trend reverses on sufficiently fine grids. RK4-ED is the most accurate method on all grids. While not shown, we also tested RK4-EC2 and that resulted in an almost indistinguishable difference in error compared to SRK4-EC2. To model the speedup of SRK4-EC2 compared to RK4-ED, we apply a cubic polynomial least square fit to each data series. Using this model, we observe that the maximum performance increase is about 150 % in favor of SRK4-EC2 for 10% error. Due to the improved accuracy of the dissipative coupling, this performance increase rapidly diminishes. At an error level of 0.01% the performance increase is completely saturated. In conclusion, the energy conservative method outperforms the energy dissipative method on coarse grids due to better stability properties. On the other hand, the energy dissipative method outperforms the energy conservative method on fine grids due to better accuracy properties (see Fig. 10).
To determine computational timings, we implemented each method in a similar manner using C++ and CUDA. The computational timing experiments were measured on a Nvidia GTX 2080 Ti card. No efforts were made to optimize the performance of any of the implementations.

Conclusions
When weakly coupling hyperbolic partial differential equations, one must select certain penalty parameters that describe how the coupling conditions are weighted. Within the SBP community, it is well-established how to constrain these parameters to give a proof of semi-discrete stability of a numerical scheme via the energy method. However, Fig. 9 Computational efficiency of the SRK4 and RK4 schemes Fig. 10 Modeled speedup of SRK4-EC2 compared to RK4-ED (speedup > 1 SRK4-EC2 is faster than RK4-ED) the energy method alone cannot fully constrain all of the parameters. The remaining parameters must be carefully selected because they can have a striking impact on stiffness and accuracy of the numerical scheme.
We showed the importance of this claim by simulating the interaction of waves across an air-water interface in an energy conserving manner. For one set of parameter values, the coupling treatment caused the numerical scheme to be an order of magnitude stiffer than expected from the CFL condition from the uncoupled case. Another set of parameter values prevented stiffness, but instead caused a degradation in accuracy. This accuracy degradation showed the development of spurious waves in the vicinity of the interface. In this study, we explained what causes stiffness and developed a general coupling procedure that is provably non-stiff and accurate.
To overcome stiffness in an automated manner, we presented a general formulation for designing penalty terms. This general formulation results in stable and dual consistent schemes that are provably non-stiff. In this formulation, the penalty terms are related to projection matrices that are guaranteed to exist as long as a determinant condition is satisfied. In the limit when the determinant approaches zero, the matrix norm grows without bounds. We gave two examples for which the determinant never vanishes; these examples are associated with either an energy conservative or an energy dissipative coupling.
A potential advantage of the energy conservative coupling compared to the energy dissipative coupling is that it is compatible with energy conservative time stepping schemes. For this reason, we compared the computational efficiency of a fourth order staggered Runge-Kutta with the energy conservative coupling against a fourth order Runge-Kutta with the energy dissipative coupling. We showed that energy conserving methods can outperform energy dissipative methods on coarse grids because they allow larger time steps. However, energy dissipative methods outperform energy conservative methods on fine grids because they have better accuracy properties. For the air-water interface problem, we observed a factor of 2.5 speedup of the energy conservative method compared to the energy dissipative method on a coarse grid.
There is one important aspect that we did not address: what mechanism of the penalty parameters controls accuracy? For the energy-conservative penalty parameters for the air-water interface problem, we saw that a simple change in a parameter value resulted in an order of magnitude reduction in the error level on the coarsest grid. We linked this accuracy improvement to what the coupling treatment must become in the limit that the impedance contrast approaches infinity. Ideally, we would like to have a simple diagnostic test, analogous to the stiffness test, that can establish if a set of parameters will cause an accurate interface treatment or not. It would also be useful to have an automated and general procedure that can select these parameters, thereby avoiding the cumbersome work of manually deriving and testing parameters on a case by case basis.
Funding Open access funding provided by Linköping University. Ossian O'Reilly was supported by the Southern California Earthquake Center (Contribution No. 10148), with contributions from National Science Foundation Cooperative Agreements EAR-1600087, United States Geological Survey Cooperative Agreement G17AC00047, NSF-ACI Award 1450451, and Pacific Gas and Electric. Jan Nordström was supported by Vetenskapsrådet (Award Number 2018-05084 VR), Sweden.

Declarations
Competing interests The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Code availability Codes for reproducing the results in this work are available at github.com/ooreilly/sbp.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
We diagonalize each side of (77) individually. By the primal (28) In general, the only possibility for (78) to be equal to (79) is if (46) holds.

B Proof of Proposition 8
Proof (i) Upper bound: By the triangle inequality and Kronecker product property (51), we obtain Since A x is symmetric, A x 2 = (A x ), and from (54), the upper bound in (55) follows.
(ii) Lower bound: By definition, any vector z = 0 yields a lower bound on M h 2 : In particular, if D x is a consistent approximation of the derivative, the choice z = 1⊗x implies that D x e = 0, and therefore By (51), (54), and since x = 0 is arbitrary,

C Proof of Proposition 9
Proof By the projection matrix formula (35) and the submultiplicative norm property (50), we have We will bound each of three factors in the right-hand side of (80) individually. By the definition of P in (35), ζ 1 satisfies By eigenvector orthogonality and orthonormality, ζ 2 satisfies Similarly, ζ 3 satisfies By inserting (81)-(83) into (80), the bound (56) follows.