10.2 Interface boundary conditions

Interface boundary conditions are needed when there are multiple grids in the computational domain, as discussed in Section 11 below. This applies to complex geometries, when multiple patches are used for computational efficiency, to adapt the domain decomposition to the wave zone, mesh refinement, or a combination of these.

A simple approach for exchanging the information between two or more grids is a combination of interpolation and extrapolation of the whole state vector at the intersection of the grids. This is the method of choice in mesh refinement, for example, and works well in practice. In the case of curvilinear grids and very–high-order FD or spectral methods, though, it is in general not only difficult to prove numerical stability, but even to find a scheme that exhibits stability from a practical point of view.

10.2.1 Penalty conditions

The penalty method discussed above for outer boundary conditions can also be used for multi-domain interface ones, including those present in complex geometries [115, 235, 118, 311, 312, 140]. It is simple to implement, robust, leads to stability for a very large class of problems and preserves the accuracy of arbitrary high-order FD and spectral methods.

Finite differences.
As an example, consider again an advection equation but now on the whole real line

where we assume that the initial data is smooth and has compact support. The domain of the real line is chosen just for simplicity, to focus on the interface procedure at . In the realistic case of compact domains, outer boundaries are also present and these can be treated by any of the methods discussed in Section 10.1.

Consider two grids, left and right, covering the intervals , and , respectively:

where the gridspacings need not agree, , and the difference operators and , which are not necessarily equal to each other either, and do not even need to be of the same order of accuracy, satisfy SBP with respect to scalar products given by the weights on their individual grids:

For the time, being assume that both SBP scalar products are diagonal, though. Then, the SAT semi-discrete approximation to the problem is

Notice that in this approach the numerical solution at any fixed resolution is bi-valued at the interface boundary , and in the same spirit as the penalty approach for outer boundary conditions; continuity of the fields at the interface is not enforced strongly but weakly.

Defining the energy

and using the approximation (10.38, 10.39) and the SBP property of the difference operators, its time derivative is
Then an estimate follows if two conditions are satisfied. One of them is . The other one imposes an additional constraint on the values of and :
• Positive :
The estimate is
• Negative : this is obtained from the previous case after the transformation
and
• Vanishing : this can be seen as the limiting case of any of the above two, with

The following remarks are in order:

• For the minimum values of allowed by the above inequalities, the energy estimate is the same as for the single grid case with outer boundary conditions, see Section 10.1.3, and the discretization is time-stable (see Section 7.4), while for larger values of there is damping in the energy, which is proportional to the mismatch at the interface.
• Except for the case of the most natural choice , the evolution equations for outgoing modes also need to be penalized in a consistent way in order to derive an energy estimate. However, as is always the case, the lack of an energy-type estimate does not mean that the scheme is unstable, since the energy method provides sufficient but not always necessary conditions for stability.
• The general case of symmetric hyperbolic systems follows along the same lines: a decomposition into characteristic variables is performed and the evolution equation for each of them is penalized as in the advection equation example. At least for diagonal norms, stability also follows for general linear symmetric hyperbolic systems in several dimensions. With the standard caveats for non-diagonal norms, the procedure follows in a similar way except that penalty terms are not only added to the evolution equations at the interface on each grid, but also near them. In practice, though, applying penalties just at the interfaces appears to work well in practice in many situations.

Spectral methods.
The standard procedure for interface spectral methods is to penalize each incoming characteristic variable, exactly as in the outer boundary condition case. Namely, as in Eq. (10.33) with lower bounds for the penalty strengths given by Eqs. (10.34) and (10.36) for Legendre and Chebyshev polynomials, respectively. We know from the FD analysis above, though, that in general this does not imply an energy estimate, and in order to achieve one, outgoing modes also need to be penalized, with strengths that are coupled to the penalty for incoming modes. However, the procedure of penalizing just incoming modes at interfaces appears to works well in practice, so we analyze this in some detail.

Figure 6 shows the spectrum of the Chebyshev penalty method as described, for an advection equation,

where there is an interface boundary at . In more detail: the figure shows the maximum real component in the spectrum for collocation points as a function of the penalty strength . There are no real positive values, and this remains true for different values of , in agreement with the fact that penalizing just the incoming mode works well in practice. The figure also supports the possibility that might actually be enough for stability, as in the outer boundary case.

The figure also shows the spectral radius as a function of the penalty strength. Beyond it grows very quickly, and even though as mentioned the timestep is usually determined by keeping the time integration error below that one due to spatial discretization, that might not be the case if the spectral radius is too large. Thus, it is probably good to keep .