10.1 Outer boundary conditions

10.1.1 Injection

Injecting boundary conditions is presumably the simplest way to numerically impose them. It implies simply overwriting, at every or some number of timesteps, the numerical solution for each incoming characteristic variable or its time derivative with the conditions that they should satisfy.

Stability of the injection approach can be analyzed through GKS theory [229], since energy estimates are, in general, not available for it (the reason for this should become more clear when discussing the projection and penalty methods). Stability analyses, in general, not only depend on the spatial approximation (and time integration in the fully-discrete case) but are in general also equation-dependent. Therefore, stability is, in general, difficult to establish, especially for high-order schemes and nontrivial systems. For this reason a semi-discrete eigenvalue analysis is many times performed. Even though this only provides necessary conditions (namely, the von Neumann condition (7.80View Equation)) for stability, it serves as a rule of thumb and to discard obviously unstable schemes.

As an example, we discuss the advection equation with “freezing” boundary condition,

ut = ux, x ∈ [− 1,1], t ≥ 0, (10.1 ) u(0,x ) = f (x), x ∈ [− 1,1], (10.2 ) u(t,1) = g(t), t ≥ 0, (10.3 )
where the boundary data is chosen to be constant, g(t) = f (1). As space approximation we use a Chebyshev collocation method at Gauss–Lobatto nodes. The approximation then takes the form
d --uN (t,xi) = (QuN )(t,xi), i = 0,1,...,N, (10.4 ) dt
where Q is the corresponding Chebyshev differentiation matrix; see Section 9.6. The boundary condition is then imposed by replacing the N th equation by
-d u (t,x ) = -dg(t) = 0. (10.5 ) dt N N dt
For this problem and approximation an energy estimate can be derived [209Jump To The Next Citation Point], but we present an eigenvalue analysis as a typical example of those done for more complicated systems and injection boundary conditions. Figure 4View Image shows the semi-discrete spectrum using a different number of collocation points N. As needed by the strong version of the von Neumann condition, no positive real component is present [cf. Eq. (7.82View Equation)]. We also note that, as discussed at the beginning of Section 9.8, the spectral radius scales as ∼ N 2.
View Image

Figure 4: Semi-discrete spectrum for the advection equation and freezing boundary conditions [Eqs. (10.1View Equation, 10.2View Equation, 10.3View Equation], a Chebyshev collocation method, and numerical injection of the boundary condition. The number of collocation points is (from left to right) N = 20,40,60. The scheme passes the von Neumann condition (no positive real component in the spectrum). In fact, as discussed in the body of the text, the method can be shown to actually be numerically stable.

There are other difficulties with the injection method, besides the fact that stability results are usually partial and incomplete for realistic systems and/or high-order or spectral methods. One of them is that it sometimes happens that a full GKS analysis can actually be carried out for a simple problem and scheme, and the result turns out to be stable but not time-stable (see Section 7.4 for a discussion on time-stability). Or the scheme is not time-stable when applied to a more complicated system (see, for example, [116, 1]).

Seeking stable numerical boundary conditions for realistic systems, which preserve the accuracy of high-order finite-difference and spectral methods has been a recurring theme in numerical methods for time-dependent partial differential equations for a long time, especially for nontrivial domains, with substantial progress over the last decade – especially with the penalty method discussed below. Before doing so we review another method, which improves on the injection one in that stability can be shown for rather general systems and arbitrary high-order FD schemes.

10.1.2 Projections

Assume that a given IBVP is well posed and admits an energy estimate, as discussed in Section 5.2. Furthermore, assume that, up to the control of boundary terms, a semi-discrete approximation to it also admits an energy estimate. The key idea of the projection method [317, 319, 318] is to impose the boundary conditions by projecting at each time the numerical solution to the space of gridfunctions satisfying those conditions, the central aspect being that the projection is chosen to be orthogonal with respect to the scalar product under which a semi-discrete energy estimate in the absence of boundaries can be shown. The orthogonality of the projection then guarantees that the estimate including the control of the boundary term holds.

In more detail, if the spatial approximation prior to the imposition of the boundary conditions is written as

ut = Qu (10.6 )
and [for example, as is many times the case when SBP holds] there is a semi-discrete energy,
E = ⟨u,u⟩, (10.7 ) d
with respect to some scalar product, for which an estimate holds up to boundary terms:
-d dtEd = ⟨Qu, u⟩ + ⟨u,Qu ⟩, (10.8 )
with
⟨Qu, u⟩ + ⟨uQu ⟩ ≤ 2αEd + boundary terms (10.9 )
for some constant α independent of the initial data and resolution. Due to the SBP property, the boundary terms are exactly those present in the continuum energy estimate, except that so far they cannot be bounded because the numerical solution is not yet required to satisfy the boundary conditions.

The latter are then imposed by changing the semi-discrete equation (10.6View Equation) to

ut = 𝒫Qu (10.10 )
where 𝒫 projects u to the space of gridfunctions satisfying the desired boundary conditions, is time-independent (which in particular requires the assumption that the boundary conditions are time-independent) and symmetric
⟨u, 𝒫v ⟩ = ⟨𝒫u, v ⟩ (10.11 )
(and, being a projection, 2 𝒫 = 𝒫). Since the projection is assumed to be time-independent, projecting both sides of Eq. (10.6View Equation) implies that
(𝒫u )t = 𝒫Qu (10.12 )
and for any solution of (10.6View Equation) satisfying the boundary conditions, 𝒫u = u,
ut = 𝒫Qu. (10.13 )
Then
d-Ed = ⟨𝒫Qu, u ⟩ + ⟨u,𝒫Qu ⟩ = ⟨Qu, 𝒫u ⟩ + ⟨𝒫u, Qu ⟩ = ⟨Qu, u⟩ + ⟨u, Qu ⟩ dt ≤ 2αEd + boundary terms.
Since the solution now satisfies the boundary conditions, the boundary terms can be bounded as in the continuum and stability of the semi-discrete IBVP follows. In principle, this method only guarantees stability for time-independent, homogeneous boundary conditions. Homogeneity can be assumed by a redefinition of the variables; however, the restriction of a time-independent boundary condition is more severe.

Details on how to explicitly construct the projection can be found in [227]. The orthogonal projection method guarantees stability for a large class of problems admitting a continuum energy estimate. However, its implementation is somewhat involved.

10.1.3 Penalty conditions

A simple and robust method for imposing numerical boundary conditions, either at outer or interpatch boundaries, such as those appearing in domain decomposition approaches, is through penalty terms. The boundary conditions are not imposed strongly but weakly, preserving the convergence order of the spatial approximation and leading to numerical stability for a large class of problems. It can be applied both to the case of FDs and spectral approximations. In fact, the spirit of the method can be traced back to Finite Elements-discontinuous Galerkin methods (see [147] and [28] for more recent results). Terms are added to the evolution equations at the boundaries to consistently penalize the mismatch between the numerical solution and the boundary conditions the exact solution is subject to.

Finite differences.
In the FD context the method is known as the Simultaneous Approximation Technique (SAT) [117]. For semi-discrete approximations of IBVPs of arbitrary high order with an energy estimate, both the order of accuracy and the presence of energy estimates are preserved when imposing the boundary conditions through it.

As an example, consider the half-space IBVP for the advection equation,

ut = λux, x ≤ 0, t ≥ 0, (10.14 ) u (0,x) = f(x), x ≤ 0, (10.15 ) u(t,0) = g(t), t ≥ 0, (10.16 )
where λ ≥ 0 is the characteristic speed in the − x direction.

We first consider a semi-discrete approximation using some FD operator D satisfying SBP with respect to a scalar product Σ, which we assume to be either diagonal or restricted full (see Section 8.3). As usual, Δx denotes the spacing between gridpoints,

xi = iΔx, i = 0,− 1,− 2,.... (10.17 )
In the SAT approach the boundary conditions are imposed through a penalty term, which is applied at the boundary point x 0,
d δi,0S --ui = λDui + -------(g − u0 ), (10.18 ) dt σ00Δx
where δi,0 is the Kronecker delta, S a real parameter to be restricted below, and σ00 is the 00-component of the SBP scalar product Σ. Introducing the semi-discrete SBP energy
Ed = ⟨u,u⟩Σ, (10.19 )
its time derivative is
-d 2 dtEd = (λ − 2S )u0 + 2gu0S. (10.20 )
For homogeneous boundary conditions (g = 0) both numerical and time-stability follow if S ≥ λ ∕2:
d-E = − (2S − λ )u2≤ 0. (10.21 ) dt d 0
Strong stability follows if S > λ∕2, including the case of inhomogeneous boundary conditions: for any 𝜀 such that
0 < 𝜀2 < 2S − λ, (10.22 )
Eq. (10.20View Equation) implies
( ) ( )2 ( )2 d-Ed = (λ − 2S )u2 + 2(𝜀u0) gS- ≤ (λ − 2S)u2 + (𝜀2u2) + gS- = − &tidle;𝜀u2 + S- g2,(10.23 ) dt 0 𝜀 0 0 𝜀 0 𝜀
where
&tidle;𝜀 := 2S − λ − 𝜀2 > 0. (10.24 )
Integrating Eq. (10.23View Equation) in time,
∫ ( ) ∫ 2 t 2 2 S 2 t 2 ∥u (t)∥Σ + &tidle;𝜀 u0(τ )d τ ≤ ∥u(0)∥Σ + -𝜀 g (τ )d τ, (10.25 ) 0 0
thereby proving strong stability.

In the case of diagonal SBP norms it is straightforward to derive similar energy estimates for general linear symmetric hyperbolic systems of equations in several dimensions, simply by working with each characteristic variable at a time, at each boundary. A penalty term is applied to the evolution term of each incoming characteristic variable at a time, as in Eq. (10.18View Equation), where λ is replaced by the corresponding characteristic speed. In particular, edges and corners are dealt with by simply imposing the boundary conditions with respect to the normal to each boundary, and an energy estimate follows.

The global semi-discrete convergence rate can be estimated as follows. Define the error gridfunction as the difference between the numerical solution u and the exact one u(e) evaluated at the gridpoints,

e(t) = e(t,x ) := u(t,x ) − u(e)(t,x ), i = 0, − 1,− 2,..., (10.26 ) i i i i
It satisfies the equation
d- -δi,0S--- dtei = λDei − σ00Δx e0 − Fi. (10.27 )
where Fi denotes the truncation error, here solely depending on the differentiation approximation:
( ) Fi(t) = F (t,xi) = λ u(xe)(t,xi) − Du (e)(t,xi) , (10.28 ) { r = 𝒪 (Δx )2p at and close to boundaries, (10.29 ) 𝒪 (Δx ) in the interior,
with r < 2p in general. For example, in the diagonal case r = p, and in the restricted full one r = 2p − 1 [see the discussion below Eq. (8.29View Equation)].

Using Eq. (10.26View Equation) and the SBP property, the norm of the error satisfies

d --∥e∥2Σ = (λ − 2S)e20 − 2⟨e,F ⟩Σ ≤ 2∥e∥ Σ∥F ∥Σ, (10.30 ) dt
where we have used S ≥ λ∕2 and the Cauchy–Schwarz inequality in the second step. Dividing both sides of the inequality by 2∥e∥ Σ and integrating we obtain
∫t ∥e(t)∥ ≤ ∥F (τ )∥ dτ, t ≥ 0. (10.31 ) Σ Σ 0
Since by Eq. (10.29View Equation) the norm of the truncation error is of the order (r + 1∕2), it follows that the Σ-norm of the error converges to zero with rate (r + 1 ∕2). This also implies that the error converges pointwise to zero with rate,
ei(t) = 𝒪 ((Δx )r) , t > 0, i = 0,− 1,− 2,... (10.32 )
In particular, this proves that the error at the boundary point x0 = 0 converges to zero at least as fast as r (Δx ), and therefore, the penalty term in Eq. (10.18View Equation) is bounded. Two remarks are in order:

Spectral methods.
The penalty method for imposing boundary conditions was actually introduced for spectral methods prior to FDs in [198, 199]. In fact, as we will see below, the FD and spectral cases follow very similar paths. Here we only discuss its application to the collocation method. Furthermore, as discussed in Section 9.4, when solving an IBVP, Gauss–Lobatto collocation points are natural among Gauss-type nodes, because they include the end points of the interval. We restrict our review to them, but the penalty method applies equally well to the other nodes. We refer to [236] for a thorough analysis of spectral penalty methods.

Like the FD case, we summarize the method through the example of the advection problem (10.14View Equation, 10.15View Equation, 10.16View Equation), except that now we consider the bounded domain x ∈ [− 1,1] and apply the boundary condition at x = 1. Furthermore, we first consider a truncated expansion in Legendre polynomials. A penalty term with strength τ is added to the evolution equation at the last collocation point:

d ′ --uN (t,xi) = λu N(t,xi) + τ λ(g(t) − uN (t,xN ))δiN i = 0,1,...,N, (10.33 ) dt
where {xi}Ni=0 are now the Gauss–Lobatto nodes; in particular, x0 = − 1 and xN = 1.

Using the discrete energy given by Eq. (9.95View Equation) and SBP property associated with Gauss quadratures discussed in Section 9.4, a discrete energy estimate follows exactly as in the FD case if

1 N (N + 1 ) τ ≥ ---N-≥ ----------, (10.34 ) 2A N 4
where AN N is the Legendre–Gauss–Lobatto quadrature weight at the right end point x = 1 N [cf. Eqs. (9.63View Equation, 9.65View Equation)]. Notice that this scaling of the penalty with the weight is the analog of the (σ00Δx ) SBP weight term for the FD case [Eq. (10.18View Equation)]. The similarity is not accidental: both weights arrive from the discrete integration formulae for the energy. Notice:

Devising a Chebyshev penalty scheme that guarantees stability is more convoluted. In particular, the advection equation is already not in the Chebyshev norm, for example (see [209]). Stability in the L2 norm can be established using the Chebyshev–Legendre method [144], where the Chebyshev–Gauss–Lobatto nodes are used for the approximation, but the Legendre ones for satisfying the equation. In this approach, the penalty method is global, because it adds terms to the right-hand side of the equations not only at the endpoint, but at all other collocation points as well.

A simpler approach, where the penalty term is only applied to the boundary collocation point, as in the Legendre and FDs case, is to show stability in a different norm. For example, in [137] it is shown that a penalty term as in Eq. (10.33View Equation) is stable for the Chebyshev–Gauss–Lobatto case in the norm defined by the weight, cf. Eq. (9.17View Equation),

∘ ------ ω&tidle;(x) = (1 + x)ωChebyshev(x) = 1 +-x-, − 1 < x < 1, (10.35 ) 1 − x
if
N 2 τ ≥ --- . (10.36 ) 2
In Figure 5View Image we show the spectrum of an approximation of the advection equation (with speed λ = 1 and g(t) = 0) using the Chebyshev penalty method, N = 20, 40,60 collocation points and S := τ N 2 = 0.3. There are no eigenvalues with positive real component. Even though this is only a necessary condition for stability, it suggests that it might be possible to lower the bound (10.36View Equation) on the minimum penalty needed for stability. On the other hand, for S = 0.2 (and lower values) real positive values do appear in the spectrum. Compared to the injection case (Figure 4View Image), the spectral radius for the cases shown in Figure 5View Image are about a factor of two larger. However, in most applications using the method of lines, the timestep is dictated by the condition of keeping the time integration error smaller than the spatial one, so as to maintain spectral convergence, and not the CFL limit.
View Image

Figure 5: Spectrum of the Chebyshev penalty method for the advection equation and (from left to right) N = 20,40,60 collocation points and S = 0.3.

  Go to previous page Go up Go to next page