Well-posed transparent boundary conditions for the shallow water equations

Applied Numerical Mathematics 38 (2001) 445–474 www.elsevier.com/locate/apnum

Well-posed transparent boundary conditions for the shallow water equations ✩ Ivar Lie Norwegian Meteorological Institute, P.O. Box 43 Blindern, N-0313 Oslo, Norway

Abstract We derive transparent boundary conditions for the shallow water equations often used as a simple model for both atmospheric and oceanic flows. Then well-posed results are presented when the time discretization is done by the characteristic-based semi-Lagrange method and the space discretization is done by two variants of the mixed finite element methods. Since the mixed finite element discretization is related to the staggered finite difference schemes extensively used in geophysical models, the results points to how the transparent boundary conditions can be used in the finite difference discretization.  2001 IMACS. Published by Elsevier Science B.V. All rights reserved. Keywords: Transparent boundary condition; Semi-Lagrangian time discretization; Mixed finite element

1. Introduction High-quality transparent boundary conditions are important in the use of limited-area models in numerical weather prediction in order to have accurate model results also near the boundary. More or less satisfactory boundary conditions are used in the numerical weather prediction community, see, for example, [14,15,17], and these conditions work reasonably well for large and medium-scale models. More accurate boundary conditions are required when running fine-scale models on small domains since the large-scale forcings will be completely determined by the boundary values. The ideal situation would be to find a high-quality transparent boundary condition which is easy to implement. Unfortunately (and fairly obviously) we are not in that situation. There are at least two problems: • The exact transparent boundary condition, in the cases where it can be found, is not a practical mathematical object to work with. Appropriate approximations must be used or one have to use other (approximate and heuristic) techniques to find reasonable transparent boundary conditions. • If one has derived a set of transparent boundary conditions, they may often be difficult to implement properly and efficiently into an existing discretization scheme. ✩

This work was done as a part of the project: Transparent boundary conditions for the HIRLAM model. E-mail address: [email protected] (I. Lie).

0168-9274/01/$ – see front matter  2001 IMACS. Published by Elsevier Science B.V. All rights reserved. PII: S 0 1 6 8 - 9 2 7 4 ( 0 1 ) 0 0 0 4 5 - 9

446

I. Lie / Applied Numerical Mathematics 38 (2001) 445–474

In this paper we will only discuss briefly how one can derive high-quality transparent boundary conditions for the equations used in numerical weather prediction. What we present are a few examples of approximate boundary conditions for a version of the shallow water equations. Work on transparent boundary conditions for the governing equations of numerical weather prediction will appear in forthcoming reports. We will concentrate on the second item above, namely how transparent boundary conditions can be incorporated into a discretization scheme. The background for this work is that many weather centers want to have well-posed transparent boundary conditions for their models. This means moving away from the current flow relaxation scheme, and other approximations. Therefore, new efficient transparent boundary conditions have to be derived and tested, first in an experimental setting, then in an operational setting. This work has started with the papers of Aidan McDonald, see [14–16], where he discusses the problems of incorporating a class of transparent boundary conditions into a semi-Lagrangian staggered finite-difference discretization, well known to the numerical weather prediction community. Using the classical difference schemes, he develops several variants of algorithms to compute boundary values and discusses the implementation aspects. In [23] we developed a multidomain technique for solving the shallow water equations where transparent boundary conditions were incorporated. The discretization we used was mixed finite elements. This technique works very well in practice, see also [25]. We will therefore develop this technique further in this paper, and provide details so one can compare the boundary treatment with the finite difference approach. It is to be emphasized that interface conditions between subdomains in a multidomain setting, as developed in [23], and transparent boundary conditions are two aspects of the same thing, in the latter case the neighboring subdomain is some exterior domain. Hence once we have appropriate transparent boundary conditions, we can formulate numerical weather prediction problems (and many other related problems) in a multidomain setting. This will provide greater flexibility in terms of models and model nesting. Imposing transparent boundary conditions based on characteristics is not straightforward using the classical finite difference schemes on a staggered grid. The different finite difference schemes tested in [14,15] show this clearly. We formulate the shallow water equations with transparent boundary conditions in the mixed finite element method. This means that the boundary conditions are imposed weakly, and the imposition comes directly out of the variational formulation, in contrast to the finite difference case where the boundary conditions have to imposed strongly by special finite difference formulas. The difference between strong and weak imposition of boundary conditions is usually small, and becomes negligible when the solution is smooth. Moreover, the mixed finite element method of course yields a linear system to solve, and that can be interpreted as a finite difference scheme, which may point to how proceed in the finite difference case. In this paper we will use the following simplified version of the shallow water equations: Du + ∇φ = 0, Dt Dφ
+ φ ∇ · u = 0. Dt

(1a) (1b)

I. Lie / Applied Numerical Mathematics 38 (2001) 445–474

447

Here u denotes the wind vector and φ the geopotential. For simplicity we have not included rotation here, but this is not conceptually different with respect to the numerical methods we will present. In the case of rotation, the momentum equation reads: 

Du + Mu + ∇φ = 0, Dt

M=



0 f

−f  . 0

(2)

The equations are linearized in the sense that φ
is a fixed value. We assume that we solve this PDE system in a domain Ω × [0, T ] with Lipschitz-continuous boundary ∂Ω. Initial conditions are given on Ω and our goal is to derive transparent boundary conditions on ∂Ω for t ∈ (0, T ]. Since (1) is a hyperbolic system the initial and boundary conditions have to satisfy certain compatibility conditions in order to produce a solution which has a certain smoothness, for details see [21]. This aspect, often ignored in applications, is actually very important if we want to prove some error estimates. The PDE system (1) is linear after semi-Lagrangian time discretization (see below), so that the discrete system can be solved by some linear solver. We emphasize that our approach can be used to solve more general problems:In [23] the nonlinear shallow water equations are solved in a multidomain setting. This means that Ω = ni=1 Ωi , Ωi ∩ Ωj = ∅ for i = j . Let Γij denote the interface between Ωi and Ωj , and let Ni denote the index set of neighbors to Ωi . The boundary for a subdomain is then partitioned into the boundary for Ω and interfaces: ∂Ωi = (∂Ω ∩ ∂Ωi ) ∪



Γij .

j ∈Ni

The nonlinear shallow water equations are then then solved on Ωi , i = 1, . . . , n, with transparent boundary conditions on ∂Ω ∩ ∂Ωi and interface conditions on Γij , j ∈ Ni . The transparent boundary conditions and the interface conditions between the subdomains are actually the same expression, the only thing that differs is that on Γij the values of the incoming characteristics for Ωi are taken from Ωj , while for transparent boundary the values are taken from the exterior domain. Because the subdomains are non-overlapping we have to use an iterative process to obtain a globally smooth solution, for details, see [20,23]. This paper is organized as follows: We will first present the derivation of the transparent boundary conditions for (1), using characteristics and also a higher order approximation to the exact transparent boundary condition. In, e.g., [23] we use a semi-Lagrangian algorithm which is a combination of the classical variant and a variant with time interpolation. This is discussed in Section 3. We will also present two variants of the mixed finite element method, the classical dual mixed method and a variant which gives better approximation of the geopotential. This will be discussed in some detail, since these methods are not commonly used for numerical weather prediction.

2. Transparent boundary conditions In this section we will derive characteristic boundary conditions for the PDE system (1) which will be used in the space discretizations below. To this end, we write (1) in the following way: Ut +

2  i=1

Ai (U)Uxi = 0,

(3)

448

I. Lie / Applied Numerical Mathematics 38 (2001) 445–474

where U = [u, φ]T is the vector of unknowns. Now from standard results for hyperbolic systems, see, e.g., [9], we find the characteristic variables of (3) by computing the left eigenvectors of the matrix 2 



Ai U
ni ,

(4)

i=1

where U
is a frozen value of U valid at the boundary and ni are the components of the unit normal vector n at ∂Ω. Evidently both U
and n will depend on where on ∂Ω we are. In our case we have 



u  A1 =  0  φ


0

1

u

 0 , 

0

u





v  A2 =  0 

0

0 v

0  

1 ,

φ
v

where u and v denote the frozen values of the wind vector u, for example, the results from the previous timestep. The eigenvalues of (4) are:

u·n+




u · n, u · n − φ,

  φ
,

and the corresponding characteristic variables: 

ψ1 = φ + φ
u · n, ψ2 = u × n,

(5a) (5b)

φ
u · n.

(5c)

ψ3 = φ −



We see that we  must distinguish between three cases:


In this case ψ1 is always outflow, ψ3 is always inflow and the direction of ψ2 depends (a) |u · n| < φ: on the flow.  So for well-posedness, ψ3 (and sometimes ψ2 ) must be specified.


In this case all characteristics are outflow and hence none of them have to be specified. φ:

(b) u · n





In this case all characteristics are inflow and all of them have to be specified. (c) −u · n > φ: The characteristic boundary conditions derived above is a special case of the boundary conditions that can be derived from approximations to the Neumann operator, also called the Dirichlet-to-Neumann map. We will describe the Neumann operator techniques briefly, following [4]. We consider the transparent boundary condition in the x1 -direction and write (3) in the form: ∂U ∂U ∂U − A2 =− . (6) A1 ∂x1 ∂t ∂x2 In [4], A1 is assumed to be diagonal and the PDE system is assumed to be strictly hyperbolic in the x1 -direction. Hence, 

A1 = 



A11

0

0

A21

,

Hence we can write ∂U ∂U ∂U +E =A , ∂x1 ∂t ∂x2

A11 = diag(λ1 , . . . , λk ) < 0,

A21 = diag(λk+1 , . . . , λn ) > 0.

(7)

I. Lie / Applied Numerical Mathematics 38 (2001) 445–474

449

−1 with A = −A−1 1 and E = −A1 A2 . Following [4,30], we now take the Fourier–Laplace transform in t and x2 and denote the transformed U variable by U
, we obtain

∂U
= (Aiξ + Eiω)U
= M(ξ, ω)U. (8) ∂x1 We now have to diagonalize the symbol M(ξ, ω). This is done by a smooth matrix V (ξ, ω) invertible under certain conditions, defining the symbol of a pseudodifferential operator so that if 



∂ ∂ , , x1 , x2 U, W =V ∂t ∂x2 the PDE system (8) can be written 



0 ∂W  Λ11 (ξ, ω) , = ∂x1 Λ21 (ξ, ω) Λ22 (ξ, ω)

(9)

where Λ11 (ξ, ω) is a k × k pseudodifferential operator where for the principal symbol we have −1 Λ11 (1, 0) = diag(λ−1 1 , . . . , λk ). In general the PDE system can have variable coefficients, and hence the pseudodifferential operators we use may also depend on x1 and x2 even if we don’t specify it explicitly. Note that Λ11 has all the negative eigenvalues corresponding to inflow. Thus the theoretical perfectly absorbing boundary condition annihilates these components at the boundary if we assume that the values for the variables outside Ω are zero. If these values are not zero we can perform a homogenization. With the notation πk (W1 , . . . , Wk , Wk+1 , . . . , Wn ) ≡ (W1 , . . . , Wk ) the perfectly absorbing boundary condition is: πk W = πk V (ξ, ω)U = 0,

(10)

at the boundary. This nonlocal and abstract condition can be made more practical by considering the asymptotic expansion of V : 1 1 V (ξ, ω)  V0 (ξ, ω) + V−1 (ξ, ω) + 2 V−2 (ξ, ω) + · · · , ξ ξ where V−j is homogeneous of degree zero. The pseudodifferential operators above are still non-local, so in order to have local conditions, these operators are expanded in Taylor series around (1, 0) (normal incidence). We can now write 

V−j (ξ, ω) = V−j

ω 1, ξ



= V−j (1, 0) +

 +1  ω V−j (1, 0) + O   . k k! ∂ω ξ

  k  ω 1 ∂k k=1

ξ

The first approximation to the boundary condition is πk V0 (1, 0)U = 0,

(11)

which is the characteristic boundary condition taking into account only the normal incidence flow. The second approximation taking into account also flow tangential to the boundary is: 



∂ ∂V0 ∂ V0 (1, 0) + (1, 0) + V−1 (1, 0) U = 0. πk ∂t ∂ω ∂x2 Other boundary conditions are discussed in [4].

(12)

450

I. Lie / Applied Numerical Mathematics 38 (2001) 445–474

The assumption that A1 should be diagonal can be relaxed, it is sufficient that the PDE system be strictly hyperbolic in the x1 -direction. So we can proceed directly from (8) to find the spectral decomposition of M(ξ, ω) for fixed ξ, ω sufficiently close to (1, 0) so that M is diagonalizable. Since −1 −1 M(ξ, ω) = −A−1 1 iξ − A1 A2 iω = −iξ A1 (I + A2 µ),

for µ = ω/ξ , we may use the scaled symbol M1 (1, µ) = −A−1 1 (I + A2 µ),

(13)

since we are primarily interested in the eigenvectors, and the true eigenvalues can be found easily. This is the approach used in [30] where in addition least-squares approximations of the eigenvectors are constructed. We will compute the spectral decomposition of M1 (1, µ) based on left eigenvectors: M1 (1, µ) = T −1 (1, µ)Λ(1, µ)T (1, µ),

(14)

by a perturbation method with perturbation parameter µ. The explicit expression for M1 (1, µ) is: 

     −1 M1 (1, µ) = −A1 (I + A2 µ) = −    

u(1 ¯ + µv) ¯

−µφ


u¯ 2 − φ



+ µv) φ(1 ¯

u¯ 2 − φ
1 + µv¯ u¯ uµ ¯ φ


u¯ 2 − φ


u¯ 2 − φ


0 −

−

1 + µv¯

 

u¯ 2 − φ
   µ  ,  u¯  u(1 ¯ + µv) ¯   u¯ 2 − φ


where u¯ and v¯ denote frozen values of the variables. Note that M1 (1, µ) = A(0) + A(1) µ, (1) (1) (0) = −A−1 with A(0) = −A−1 1 and A 1 A2 , so the term A µ can be considered as a perturbation of A . The eigenvalues of M1 (1, µ) are:

λ1 = −

1 + µv¯ u¯

and λ2,3 = −

u(1 ¯ + µv) ¯ ±




u¯ 2 + v¯ 2 ) − µ2 φ
2 + 2φµ
vv µ2 φ( ¯ + φ


. u¯ 2 − φ
We see that in order to determine the sign of the eigenvalues for µ = 0, we have to make certain ¯ we have to assume assumptions. For example, if the sign of λ1 is to be determined by the sign of u, that 1 + µv¯ > 0 which is a kind of glancing condition. Also, λ2,3 have to stay real. If we assume that φ
− u ¯ 2 > 0, the expression under the square root is positive if µ is within the interval 



φ
− u¯ 2 v¯ + φ
− u¯ 2 , , φ
− u ¯ 2 φ
− u ¯ 2

v¯ −



¯ then the flow which can be very small if v¯ and φ
− u¯ 2 are small. But if v¯ is small compared to uu, is almost normal to the boundary, and 1 + µv¯ > 0 will hold. Hence it seems that the restrictions on the

I. Lie / Applied Numerical Mathematics 38 (2001) 445–474

451

values of µ do not impose any problems in formulating a transparent boundary condition appropriate for the flow. If we expand λ2,3 in Taylor series around µ = 0, we get λ2,3 = −



φ


u¯ ±

− φ


u¯ 2





(1 + µv) ¯ + O µ2 .

(15)

Note that the eigenvalues in (15) have the structure: λ2,3 (µ) = λ(0) + µλ(1) ,

(16)

¯ (0) . with λ(1) = vλ We now compute the left eigenvectors for the eigenvalues of the form (15) and for λ1 . The eigenvectors are assumed to be of the form x(µ) = x (0) + µx (1), so we can write

x (0) + µx (1)

T







T

A(0) + µA(1) = λ(0) + µλ(1) x (0) + µx (1) .

(17)

This gives the relations x (0)TA(0) = λ(0) x (0)T ,



x (1)T λ(0) I − A(0) = −x (0)T λ(1)I − A(1) .

(18a) (18b)

The computational procedure is therefore: First find x (0) from (18a), and then find x (1) from (18b). This produces the left eigenvector x(µ). The spectral decomposition of A(0) is:

(0)

λ A







1 1 1  ,−  = − ,− , u¯ u¯ + φ
u¯ − φ


with eigenvectors: 



0  1 

0

  − φ
  . 0  

φ


0 1

1

For the eigenpair [−1/u, ¯ [0, 1, 0]T ] we find  





 

 

0

 u¯ 

0

1

   

x(µ) =  1  + µ  0  ,

but the second component of the vector cannot be determined uniquely by the perturbation method we use. The eigenvector corresponds to the ‘characteristic boundary condition’: v + µ(φ + uu) ¯ = 0,

452

I. Lie / Applied Numerical Mathematics 38 (2001) 445–474

or in operator form: ∂ ∂v + (φ + uu) ¯ = 0, ∂t ∂y

(19)

which may also be written 



∂ ∂v ∂ ∂ + v¯ ¯ − v¯ = 0. v + (φ + uu) ∂t ∂y ∂y ∂y

The exact eigenvector for the eigenvalue λ1 is 

µu¯



     1 + µv¯  ,  

µ which gives the boundary condition 



∂ ∂ ∂ + v¯ ¯ = 0. v + (φ + uu) ∂t ∂y ∂y

(20)

Hence we propose to use (20) instead of (19). Note that if we use the second component of the locally linearized momentum equation: ∂v ∂v ∂φ ∂v + u¯ + v¯ + = 0, ∂t ∂x ∂y ∂y we can write (20) as 



∂ u¯ ∂v ∂u − = 0. +u −u¯ ∂x ∂y ∂y For the eigenpair [−1/(u¯ +






[ φ,
0, 1]T ] we obtain φ),

    
φ
v¯ φ          x(µ) =  0  + µ  −u¯ φ
 ,    

v¯

1 which gives: 

(1 + µv) ¯ φ+ or










− µ φ
uv φu ¯ =0

      ∂  ∂
− ∂ u¯ φv
= 0. + v¯ φ + φu ∂t ∂y ∂y

Similarly for the eigenpair [−1/(u¯ − 




[− φ),

(21)




0, 1]T ], we obtain: φ,

     ∂  ∂  ∂

= 0. + v¯ φ − φu + u¯ φv ∂t ∂y ∂y

(22)

I. Lie / Applied Numerical Mathematics 38 (2001) 445–474

453

The boundary conditions (19), (21) and (22) are extensions of the characteristic-based method, and take into account that the flow may have a tangential component at the boundary. They can be implemented by using the semi-Lagrangian time discretization, but this time one-dimensional with direction along the 


has to be approximated boundary. The advection operator is (∂/∂t + v∂/∂y) ¯ and the term (∂/∂y)(u¯ φv) by some quadrature formula, see [13]. In the analysis leading up to the well-posedness results we will only use the characteristic boundary conditions, and we refer to [13] for more details on the higher order conditions.  We will now discuss what conditions to use in the different flow directions. With u − φ
< 0, λ3 (A(0) ) >0, so (22) is an inflow boundary condition corresponding to the ‘fast’ characteristic variable
When u < 0, λ1 (A(0) ) > 0 and (19) becomes an inflow boundary condition, an extension ψ3 = φ − φu. of the condition v = v0 , where v0 is the value for the exterior domain. 

If u − φ
> 0, then λi (A(0) ) < 0, i = 1, 2, 3, hence all characteristics are outflow and no boundary conditions have to be enforced.  The last case is u + φ
< 0, all characteristics are inflow and (19), (21) and (22) have to be enforced. Our reasoning of the direction of the flow has been based on the matrix λi (A(0) ), but will also hold for λi (M1 ) provided (1 + µ) > 0, i.e., the flow direction is the same. The boundary conditions obtained above are not directly comparable with those derived in [4], since the method in [4] is based on first diagonalizing the matrix A1 . Note also the structure of the boundary conditions (21) and (22): The first term represents the characteristic variable advected by the operator (∂/∂t + v∂/∂y), ¯ and this accounts for the non-normal incidence. The second term is a contribution involving the tangential variation of the tangential velocity, so that both velocity components are present. The condition (20) can be interpreted by using the momentum equation as described above. We see that if ∂ u/∂y ¯ = 0, the expression is the vorticity of the velocity field, a condition that has been proposed on a pure physical basis very early in the development of numerical models for ocean and atmosphere. The perturbation approach used in this section can evidently be generalized to take into account higher order expansions for λ2,3 : Assume that λi (µ) =

N 

µj λ(j ) ,

x(µ) =

j =0

N 

µj x (j )

j =0

and N  j =0





µj x (j )T A(0) + µA(1) =

N  j =0

µj λ(j )

N 

µk x (k) ,

(23)

k=0

which by comparing coefficients of equal powers of µ yields a recursive procedure to compute x (0), x (1) , . . . , x (N) , and then higher order transparent boundary conditions. However, the well-posedness of such higher order boundary conditions will have to be investigated in each case.

454

I. Lie / Applied Numerical Mathematics 38 (2001) 445–474

3. Semi-Lagrangian time discretization In this section we will present briefly the semi-Lagrangian time integration schemes we use in this paper. It is assumed that the reader is familiar with the terminology in this field, for an explanation of the basic principles, see, e.g., [27]. We will discuss two versions of the semi-Lagrangian techniques, namely the usual one, and the variant with time interpolation suggested for multidomain formulations, see [23,25]. The two variants of the semi-Lagrangian method, the usual one and the one with time interpolation are different in the case when the departure point of the trajectory is outside the computational domain, or subdomain in multidomain formulation. The notation is as follows: Superscript “0” means value at the departure point of the trajectory, while no superscript means value at the arrival point, i.e., the values we want to compute. With a 2-level semiLagrangian method our PDE system reads: h h (24a) u + ∇φ = u0 − ∇φ 0 = Ru , 2 2 h h (24b) φ + φ
∇ · u = φ 0 − φ
∇ · u0 = Rφ , 2 2 where the timestep is h. This PDE system will be the starting point for our space discretization described in the following two sections. Note that provided we have sufficiently regular u and φ, we can take the divergence operator of (24a) and insert into (24b) to get the following Helmholtz equation: −'φ +

4 h2 φ


φ=

4 h2 φ




 h
Rφ − φ ∇ · Ru .

2

(25)

In [23], it is also shown how (25) can be solved by a Galerkin finite element method with transparent boundary conditions in a multidomain setting. If we use piecewise linear basis functions, the resulting discrete value φh will be found in a subspace of H 1 (Ω), and hence the velocity obtained from (24a) will only have L2 -regularity. See [23] for a comparison of computational complexity of this method compared to the mixed finite element methods we will discuss below. In this paper we will not discuss the trajectory calculation, for details we refer to [27] and also to [1] for more recent algorithms. We will now describe the right hand sides of (24a) and (24b) in more detail since this will be important when we consider the space discretization below. Note first that the trajectories in the continuous case will be the same for all variables, since they are defined everywhere. In the discrete case when variables are staggered each variable will have its own trajectory, and hence the quantities appearing on the right hand sides of (24a) and (24b) will be different even if they are denoted by the same symbol. In order to simply the discussion below, we consider only the equation for passive advection ct + U (x, t)cx = 0.

(26)

Assume that we trace a particle backwards from an arrival point (xm , tn ) to a departure point (x ∗ , tn−1 ) using the conventional semi-Lagrangian algorithm. When determining the departure point, a spatial interpolation is carried out, using grid points at the departure time level tn−1 . In this way we form a continuous function c(x, tn−1 ) for x ∈ [xl , xl+k ] for some variable c.

I. Lie / Applied Numerical Mathematics 38 (2001) 445–474

455

Fig. 1. Semi-Lagrangian algorithm with interpolation in time.

Fig. 2. The complete algorithm.

We can of course also interpolate in time at a certain grid point xl and form a continuous function c(xl , t), provided we have points c(xl , tn−k ), . . . , c(xl , tn ) available to determine c(xl , t) with sufficient accuracy for t ∈ [tn−k , tn ]. Now the value of c(xm , tn ) can be found by tracing the particle backwards until it crosses the line x = xl . Assuming that this happens at time t = t ∗ we find c(xl , t ∗ ) by interpolation in time, using c(xl , tn−k ), . . . , c(xl , tn ). This is the idea behind the semi-Lagrangian method with interpolation in time, illustrated in Fig. 1 for the one-dimensional case. The complete algorithm is based on a combination of the algorithm with space interpolation and with time interpolation. It is illustrated in Fig. 2 for the one-dimensional case in a multidomain formulation (here: two domains) and works as follows: (1) Use space interpolation variant to compute the values at all grid points in a subdomain at the next time level, except when a grid point has its departure point at the previous time step located outside this subdomain.

456

I. Lie / Applied Numerical Mathematics 38 (2001) 445–474

(2) Use space interpolation variant to compute the value at the interface using only local information. In the figure this is done in Ω1 , for the interface points A, B and C at time level tn−2 , tn−1 and tn , respectively. Communicate the interface value to the neighboring domain, Ω2 in the figure. (3) Use time interpolation variant to compute those grid points that has their departure point at the previous time step located outside the subdomain. This is done by calculating the trajectory’s intersection with the interface and interpolating in time at the intersection point. Hence we find the value at the wanted grid point. This applies to grid point D and E in Ω2 in the figure. Remark. Note that this algorithm also gives a method for treating inflow boundaries accurately, since in the case where Ω1 is substituted by the exterior domain, the values at the inflow part of the interface between the subdomains Ω1 and Ω2 are given by the (time-dependent) boundary values. With this algorithm we do not have to use extra gridpoints outside the inflow boundaries, as done in implementations of the space interpolation variant. Remark. The stability of the numerical scheme depends both on the stability of the conventional semi-Lagrangian variant and on the stability of the variant with time interpolation. The stability of the conventional variant is well known, see, e.g., [27], and the time interpolation variant is analyzed in [24] and found stable for a family of interpolants. For the details on this, see [24]. The inflow characteristics are assumed to be known at the boundary at t ∗ from time interpolation, for example ψ3∗ is calculated as ψ3∗

=

n+1 

j

αj ψ3 ,

(27)

j =n−k

so the unknowns at the boundary will appear in the expression for ψ3∗ . This relation must be used to find the departure point values, one cannot interpolate every variable and use these values directly. More details will appear when we discuss the fully discrete case below. We outline the details of the algorithm for the case of passive advection in one space dimension. The algorithmic details for forced advection and for 2D and 3D are very similar, so the principles are adequately covered by the derivation below. Assume that we use a two-time-level [27] semi-Lagrangian scheme. We first consider space interpolation variant, and integrate (26) along the trajectory from a departure point (x ∗ , tn ) at time level n to the arrival point (xi+1 , tn+1 ). This gives c(xi+1 , tn+1 ) − c(x ∗ , tn ) = 0, (28) h where h is the time step and, assuming xi  x ∗  xi+1 , x ∗ = xi + αi . In order to obtain the displacement αi and thus x ∗ we take U (xi + 12 αi , tn+1 − 12 h) to be the approximate inverse slope at the midpoint of the trajectory between (x ∗ , tn ) and (xi+1 , tn+1 ). We then get 



αi h h. αi = U xi + , tn+1 − 2 2 We can obtain αi from (29) by for example fixed-point iteration [j +1] αi



[j ]

(29)



α h = U xi + i , tn+1 − h, 2 2

(30)

I. Lie / Applied Numerical Mathematics 38 (2001) 445–474

457

Fig. 3. Time displacement on the interface.

where j is the iteration index. When αi is known, we use spatial interpolation to obtain c(x ∗ , tn ) and thus (28) gives c(xi+1 , tn+1 ). We next consider the part of the algorithm containing the variant with interpolation in time. In this case we integrate (26) along the trajectory from the departure point (xl , t ∗ ) at the interface to the arrival point (xm , tn+1 ): c(xm , tn+1 ) − c(xl , t ∗ ) = 0, (31) β where β = tn+1 − t ∗ , see Fig. 3. So here we have to determine t ∗ , or equivalently, the “time displacement” β. Like in the SLT with space interpolation we can do this by using U (x, t) as the approximate inverse slope at the midpoint between t ∗ and tn+1 . We now get 



'xml β , tn+1 − · β, (32) hl = U xl + 2 2 where 'xml = xm − xl . This equation can be solved for β if U can be evaluated at an arbitrary point, if necessary by interpolation in time. In the special case where U is independent of time, we can solve directly for h. When U is time-dependent (32) becomes a nonlinear equation in β which has to be solved by an iterative method. If fixed-point iteration is used, we get hl . (33) β [j +1] = U (xl + 'xml /2, tn+1 − β [j ] /2) The convergence properties of this method, and other possible iterative methods are discussed in [12].

4. The dual mixed method Mixed finite element methods are extensively used in many application fields, for example electromagnetics, fluid flow and structural analysis. The theory is covered is many books, see, e.g., [3,6, 22]. The principle behind mixed FEM is to construct different approximations for the unknowns in the PDE system, in our case the vector u and the scalar φ. The starting point is always the weak form (variational form) of the PDE system. In order to use the weak form in a rigorous way, we use the space L2 (Ω) and the L2 -based Sobolev spaces H m (Ω). In addition to this, a special type of space is used, see, e.g., [6]: 



u ∈ H (div, Ω) = u: u ∈ L2 (Ω)2 , ∇ · u ∈ L2 (Ω) .

458

I. Lie / Applied Numerical Mathematics 38 (2001) 445–474

If we introduce the vectorial test function v and the scalar test function ψ (regularity will be discussed shortly), the weak form of (24a) and (24b) is: 

u · v dx +

Ω



Ω

h 2

h φψ dx + 2



∇φ · v dx =

Ω



 Ω

φ
∇ · uψ dx =

Ω

Ru · v dx,

(34a)



Rφ ψ dx.

(34b)

Ω

The (classical) dual mixed method is obtained by applying the Green’s theorem to the term  Ω



Ω

u · v dx −

h 2

h φψ dx + 2



∇ · vφ dx =

Ω



 Ω

φ
∇ · uψ dx =

Ω

Ru v dx −

h 2





φv · n ds,

 Ω

∇φ ·v dx: (35a)

∂Ω

Rφ ψ dx.

(35b)

Ω

For this formulation we must have u, v ∈ H (div, Ω), φ, ψ ∈ L2 (Ω), so the weak formulation is: Find u ∈ H (div, Ω), φ ∈ L2 (Ω) such that (35a), (35b) are satisfied for ∀v ∈ H (div, Ω) and ∀ψ ∈ L2 (Ω). We see immediately that the geopotential has only L2 -regularity, and therefore can be approximated for example by piecewise constant functions. A more regular φ can be found using the primal–dual formulation described in the next section. Consider now the boundary term in (35a). The value of φ is not known in the discrete case (see below)
and has to be expressed by other quantities. We will in this and the next section assume that |u · n| < φ, which means that the characteristic variable ψ3 is inflow, and we use it to find φ. Hence, 

∂Ω

φv · n ds =

 

ψ3 +





φ
u · n v · n ds =

∂Ω



ψ3 v · n ds +

∂Ω

 

φ
u · nv · n ds.

(36)

∂Ω

The first term is a known function since ψ3 has to be specified, while the second term contains the unknown u. In this way we have imposed the characteristic inflow boundary condition weakly, a method well suited to this mixed formulation. Note also that only the normal component of the wind vector enters into the weak formulation, the characteristic ψ2 cannot be imposed weakly.
the characteristic ψ3 is also inflow so we can use it as done above. In the case where −u · n > φ, However, in this case u · n is also known at the boundary since ψ1 is also inflow, hence the boundary term in (35a) is a known function.
we suggest that the boundary value of φ should be taken as the value in In the case where u · n > φ, the neighborhood of the boundary, in the discrete case that means the adjacent cell. This implies that the boundary term will contain the unknown φ. In order to prove the existence and uniqueness of the continuous problem (35) with the boundary integral evaluated as in (36), we use the Helmholtz equation (25) which is equivalent to the strong formulation of (35) under certain smoothness restrictions. The weak formulation of (25) is: Find φ ∈ H 1 (Ω) such that 

Ω

∇φ · ∇ψ dx +

4 h2 φ




Ω

φψ dx =



Ω

R1 ψ dx +



∂Ω

∇φ · nψ ds,

(37)

I. Lie / Applied Numerical Mathematics 38 (2001) 445–474

459

for all ψ ∈ H 1 (Ω). Here we have set R1 =

4



 h
Rφ − φ ∇ · Ru .

2 h2 φ
It remains to interpret the boundary integral. In order to find the equivalent to the inflow boundary condition given by the characteristic variable ψ3 , we use the momentum equation (24a) and take the inner product with n: h (38) u · n + ∇φ · n = Ru · n. 2 

Now using the characteristic ψ3 = φ − φ
u · n, we find u · n from (38) and insert it into the characteristic variable to obtain the boundary condition for (25):   h
(39) φ + φ ∇φ · n = φ
Ru · n + ψ3 , 2 which we recognize as a Robin boundary condition for φ. Hence if we find ∇φ · n from this condition and insert into the boundary integral of (37), we have the complete weak formulation. It is well known that the solution of an elliptic problem like (37) with Robin boundary conditions exists and is unique, see, e.g., [18]. We can now formulate the existence theorem: Theorem 1. The problem (35) has a unique solution (u, φ) ∈ H (div, Ω) × L2 (Ω). In addition, φ is the solution of (37), and we have that u = −(h/2)∇φ + Ru . Proof. The proof consists of obtaining an equivalence between (35) and (37), and follows exactly the proof of Theorem 1 in [26]. ✷ Now consider the discrete formulation. We assume that Ω is partitioned into square elements of size k. Let φh , uh , ψh and v h denote the discrete approximations to the respective continuous quantities. Also, let Vφ ⊂ L2 (Ω) be the space of piecewise constant functions on each element. Hence on such an element, φh is constant. Let Vu ⊂ H (div, Ω) denote the space of vectorial basis functions of the simplest Raviart–Thomas type over an element, for details see, e.g., [3]. The discrete formulation is then: • Find uh ∈ Vu , φh ∈ Vφ such that 

Ω



h uh · vh dx − 2 φh ψh dx +

Ω

h 2

∂Ω

∇ · uh φh dx =

Ω





Ω

φ
∇ · uh ψh dx =

Ω

are satisfied with 



φh v h · n ds =



h Ru v h dx − 2





φh v h · n ds,

(40a)

∂Ω

Rφ ψh dx

(40b)

Ω

ψ3 v h · n ds +

∂Ω

for all v h ∈ Vu and ψh ∈ Vφ .

 

∂Ω

φ
uh · n v h · n ds,

(41)

460

I. Lie / Applied Numerical Mathematics 38 (2001) 445–474

The requirement to obtain a nonsingular solution of the resulting linear system is: dim(Vu ) + dim(Vφ ) = dim(Vv ) + dim(Vψ ),

(42)

which in essence says that the resulting linear system is required to be square. In order to analyze the discrete problem we write it in abstract form using the following discrete bilinear forms: a(uh , v h ) =



uh · v h dx,

(43a)

∇ · uh φh dx,

(43b)

φh ψh dx,

(43c)

Ω

b(uh , φh ) =



Ω

c(φh , ψh ) =



Ω

so that we can write (40a) and (40b) as follows: h h (44a) a(uh , v h ) − b(uh , φh ) = (Ru , v h ) − φh v h · n, 2 2 h
φb(uh , ψh ) + c(φh , ψh ) = (Rφ , ψh ). (44b) 2 This form of the abstract problem is slightly different from the standard formulation with an elliptic problem, but we can use the theory for generalized saddle-point problems in [19]. We also refer to [3, Chapter IV], where non-homogeneous boundary conditions are discussed. We have the following existence theorem: Theorem 2. Assume that the condition (42) and the following inf-sup conditions (LBB conditions) hold: a(uh , v h )
β1 > 0, (45) (a) inf sup uh ∈Vu v h ∈Vu uh H (div) v h H (div) (b) (c)

b(uh , φh )
β2 > 0, φh ∈Vφ uh ∈Vu uh H (div) φh L2 inf sup

inf sup

φh ∈Vp ψh ∈Vφ

(46)

c(φh , ψh )
β3 > 0. φh L2 ψh L2

(47)

Then the problem (44) has a unique solution. Proof. The inf-sup conditions (the only non-trivial one is (b)) can be shown using the classical theory of mixed methods, see, e.g., [3,6,22]. ✷ We also have the following error estimate: Theorem 3. Let u and φ be the solution of (35), and let uh and φh be the solution of (44). With u ∈ H 1 (Ω)2 , φ ∈ H 1 (Ω) and ∇ · u ∈ H 1 (Ω) we have:



u − uh H (div) + φ − φh L2  C k |φ|H 1 + |u|H 1 + |∇ · u|H 1 , where C is a constant independent of the element size.

(48)

I. Lie / Applied Numerical Mathematics 38 (2001) 445–474

461

✷

Proof. See, e.g., [22].

If we use higher order elements, for example the Raviart–Thomas elements of order m + 1, a similar estimate is used to establish k-order convergence if u ∈ H m (Ω)2 , φ ∈ H m (Ω) and ∇ · u ∈ H m(Ω). As usual we first compute the linear system based on (40) for one element, and then perform the assembling of the linear system for all elements. The linear system for element K reads: h Ak UK − D T φK = R1 , 2 h
2 k φK + φD UK = R2 , 2

(49a) (49b)

where UK = [uK,1 , uK,2 , uK,3 , uK,4 ]T , D = [1, 1, 1, 1]T and 

−1  2  2 1  −1 AK =  6 0  0 

0

0



0 0 2 −1

0   0  

.

−1   

2

The integrals involving the right hand sides of (40a) and (40b) are computed using the explicit expressions for Rφ and Ru .  K

4 h  Rφ ψh dx = φK0 − φ
u0 , 2 j =1 K,j

(50)

where quantities with superscript 0 refers to the departure point. The integral 

Ru · wi dx

K



exists only in the sense of distributions, and we have to use Green’s theorem on the term K ∇φh0 · v dx in order to compute it in the ordinary sense. This means that a boundary term will appear and the computation have to proceed as described for the left-hand side of (40a). We introduce a global numbering of elements and edges, and relate the local numbers of each element to this global numbering. In the standard way, the velocity unknowns changes sign in order to have continuous fluxes over element edges. Also let the elements by organized in an nxe · nye array. Let ne denote the number of elements and nu the number of edges. The global linear system has the structure: h A U + BΦ = R1 , 2 h φ
C U + DΦ = R2 , 2 where A ∈ Rnu ×nu , B ∈ Rnu ×ne , C ∈ Rne ×nu and D ∈ Rne ×ne .

(51a) (51b)

462

I. Lie / Applied Numerical Mathematics 38 (2001) 445–474

The matrix A is block diagonal with tridiagonal blocks of size 4. The blocks have the structure: 



2 

1 1 Ak =  6 0 

0

1 0 0  4 1 0 

 + diag(α1 , 0, 0, α4 ),

1 4 1  

0 1 2

where the matrix diag(α1 , 0, 0, α4 ) represents the contributions from the boundary term  

φ
uu · n v h · n ds.

∂Ω

The matrix B consists of elements εik , i = 1, . . . , nu , k = 1, . . . , ne . For an element K and an edge i, ε is defined as: εKi =

     1,

−1,

    0,

i ∈ ∂K, nxy · nKi = 1, i ∈ ∂K, nxy · nKi = −1, i∈ / ∂K,

where nxy is the unit vector in either positive x- or y-direction, and nKi is the unit normal on edge i. The means that each row in B has exactly two elements, one for each element sharing the edge, except for ¯ a boundary edge where there is only one element. If we consider φ
a constant, we can write C = φ
C, T ¯ where C = −B . In the more general case the matrices B and C will be structurally transposes. The matrix D is a diagonal matrix with elements m(K), which in our case will all be equal so we can write D = k 2 I . The vector R1 contains contributions from the integrals 

Ru · v h dx

Ω



and

ψ3 v h · n ds,

∂Ω

while the vector R2 contains contributions from the integral 

Rφ ψ dx. Ω

Before we go into more detail about the boundary treatment, consider the solution of (51). Methods for the solution of linear systems for saddle-point problems has been addressed in many papers and we will not go into detail here, only mention a couple of methods that may be efficient in our case. We refer to, for example, [22, pp. 591–600] for a overview of methods and their performance. There are at least three ways to solve our linear system. Note first that the block matrix 



A  C

B D

is nonsymmetric so a nonsymmetric solver has to be used whether we use direct or iterative methods. The first method is to solve the system as it is by a sparse direct solver or a preconditioned iterative solver. This has been done in many applications of mixed methods, see, e.g., [22], and also in [23]. The second

I. Lie / Applied Numerical Mathematics 38 (2001) 445–474

463

method utilizes the simple structure of the matrix D, and solves Φ from (51b) and inserts into (51a). The result is a linear system in U : 

A + φ




h2 h BB T U = R1 − 2 BR2 , 2 4k 2k

(52)

so the (reduced) matrix has a Schur complement structure. The third method solves U from (51a) and inserted into (51b) to obtain an equation for Φ: 

 h2 T −1 h
k I + φ B A B Φ = R2 − φ
CA−1 R1 , 2

4

2

(53)

showing again the Schur complement structure. Since in the rectangular element case there we have nu = (nxe + 1)nye + (nye + 1)nxe , the linear system (53) is considerably smaller that (52), which again is smaller than the total system. The price we pay for solving (52) compared to the full system is only a matrix–matrix multiplication (which can destroy sparsity). However, the matrix in (52) is symmetric positive definite, and can therefore be solved more efficiently than the original nonsymmetric system solving (53) involves inverting A and 3 matrix–matrix products, which may be costly. However, inverting A can be done easily since A is block diagonal with tridiagonal blocks of relatively modest size. The matrix in (53) is also symmetric positive definite, and given the size of the linear system it is likely that solving (53) is faster than solving (52). Now returning to the detailed treatment of the boundary conditions in the discrete form, we first note that all edges and all elements occurs as unknowns in the linear system (51). Hence there are no values outside the domain required to evaluate the boundary values. This is in contrast to most finite-difference discretizations, where special formulas and often extrapolation are required at the boundary. To illustrate this, consider the the momentum equation for a boundary edge uK1 in an element K: h (a11 + α1 )uK1 + a12 uK2 + b11 φK = RuK1 + ψ3,1 , 2

(54)




The quantity ψ3,1 represents the value of the where a11 = 13 , a12 = 16 , b11 = (h/2)εK1 , α1 = φ/k. incoming characteristic ψ3 at edge K1. The contributions α1 and ψ3,1 comes from the boundary term as explained above, and we see how it enters into the coefficient matrix as well as the right hand side. The formula (54) can be read off directly as a finite difference formula, but the contributions from the boundary integrals cannot be directly interpreted as boundary conditions imposed in the strong form. The equation for element K contains the discrete divergence operator:

φK + φ




h −uK1 + uK2 −uK3 + uK4 + 2k k k



=

1 Rφ . k2 1

(55)

This formula can therefore be directly interpreted as a finite difference formula. We now proceed to the computation of the right hand sides, and we start with the continuity equation, which is computed in (50). We have to interpret the quantities φK0 and u0K,j . Consider the trajectory from φK at time tn+1 . We first assume that the departure point is located inside the domain, somewhere in an element K  at time tn . Since the value of φ is constant within an element, we take for φK0 the value φK  . Let uK  ,j represent the velocity components of the element K  at time tn . For u0K,j , take the values uK  ,j .

464

I. Lie / Applied Numerical Mathematics 38 (2001) 445–474

Consider the case when the departure point is located outside the domain. Assume that the trajectory crosses the boundary in an element K  and time t ∗ , tn  t ∗  tn+1 . Now φ is not known at the boundary, so we have to calculate the value using the inflow characteristic φ = ψ3 +



φ
u · n.

Let the normal component of u at the boundary edge at time t ∗ be given by time interpolation: n+1 

u∗ · n =

β j uj · n,

β j ∈ R.

(56)

j =n−k

Hence we can compute φ ∗ = ψ3∗ +



φ
u∗ · n,

(57)

for the value φ 0 . We see that the unknown un+1 appears in the formula and will therefore contribute to the coefficient matrix. So in the case of time interpolation, the departure point values contains terms that contribute to the coefficient matrix. This is very different from the finite difference schemes used in [14],
φ ∗ is known on the and makes it difficult to compare the two approaches. In the case where −u · n > φ, boundary, so no information from characteristics is needed, and un+1 does not appear in the expression. The velocity components u0K,j are the velocity components of the element where the trajectory crosses the boundary, which we can denote by u∗K  ,j . Hence the other velocity components will also have to be interpolated using the formula (56), and will contribute to the coefficient matrix as well. Consider now the right hand side in the momentum equation, for which we obtain 

Ru · v dx =

Ω

 



h u − ∇φ 0 · v dx = 2



0

Ω

h u · v dx + 2 0

Ω



0

h φ ∇ · v dx − 2



0

Ω

φ 0 v · n ds,

(58)

∂Ω

using Green’s formula. Recall that the departure point values now refer to the trajectories for the velocity components. Since the departure point will not in general be located on an edge, a departure point value, say u0i must be computed by interpolation: u0i =



βm um ,

(59)

m∈Ni

where the values um are all located at time level n and Ni denotes an index set of neighboring values from which the interpolation formula is built. Because of the special structure of the Raviart–Thomas basis, a departure point for horizontal edge values have to be interpolated from horizontal edge values only, and similarly for a vertical edge. For an element K at time tn+1 the four velocity components will have separate trajectories in general, therefore the coefficients in the interpolation formula (59) will be specific to each velocity component. It remains to interpret φ 0 . We here assume that φ 0 is the departure point value for the trajectory with arrival point in the center of element K. This means that φ 0 is the same quantity in the right hand sides for both the momentum equation and the continuity equation. The boundary integral in (58) again have to be computed using the inflow characteristic: φ = 0

ψ30

+



φ
u0 · n,

where now u0 refers to the velocity at the boundary edge at time tn , and this is known.

I. Lie / Applied Numerical Mathematics 38 (2001) 445–474

465

When the departure point is located outside the domain, let the trajectory cross the boundary in an element K  at time t ∗ , tn  t ∗  tn+1 . We can now interpolate the velocity component at the boundary edge using the formula (56) and find the departure point value. The value of φ 0 is taken to be φ interpolated in time (using (56)) in the element K  . The value for φ 0 needed in the boundary integral is computed from the incoming characteristic using (57). So for time interpolation the terms in (58) contains velocity values at time tn+1 which contribute to the coefficient matrix of the linear system. In
modifications similar to what was discussed above must be used. the case where −u · n > φ, The existence theorem implies the following stability result which expresses the well-posedness: Corollary 4. Let uh and φh be the solution of (44). Then there exist constants c11 , c12 , c21 , c22 independent of the discretization such that 







uh H (div)  c11 Ru L2 + c12 Rφ L2 ,

(60a)

φh L2

(60b)

   c21 Ru 

L2

  + c22 Rφ L2 .

Proof. See [19]. We have used the symbols Ru and Rφ as the right hand sides since we see from the results above for the time interpolation variant that the actual right hand sides contains u-values, more precisely boundary values u·n at time tn+1 . Since the trace operator γn : u → u·n is linear and continuous from H (div, Ω) to H −1/2 (∂Ω), see, e.g., [6, Theorem 2.5], we have u · n−1/2,∂Ω  uH (div,Ω) , which implies estimates of the form:   Ru   d11 Ru  + d12 u,   Rφ   d21 Rφ  + d22 u,

(61a) (61b)

where Ru and Rφ refer to quantities not containing values at time tn+1 . Hence the corollary is a true stability estimate. ✷ 5. The primal–dual mixed method The primal–dual mixed method is the form we get without applying Green’s theorem to the PDE system, i.e., 

Ω



Ω

u · v dx +

h 2

h φψ dx + 2



Ω



Ω

∇φ · v dx =



Ω

φ
∇ · uψ dx =

Ru · v dx,

(62a)



Rφ ψ dx.

(62b)

Ω

We see that for the integrals to have a meaning, the variables and the test functions must have the following regularity: u ∈ H (div, Ω), φ ∈ H 1 (Ω), v ∈ (L2 (Ω))2 and ψ ∈ L2 (Ω). Hence the primal– dual method gives a more accurate approximation of φ than the dual-mixed method, which has only φ ∈ L2 (Ω). The method is therefore attractive when the vector and scalar variables are to be approximated with the same accuracy.

466

I. Lie / Applied Numerical Mathematics 38 (2001) 445–474

The primal–dual mixed method was introduced by Thomas and Trujillo in the context of domain decomposition methods, see [28]. It is probably a bit puzzling that the original weak form has not been used so much, but we will see when we come to the discrete formulation that the primal–dual method uses some special vectorial test functions and special care has to be taken on the boundaries when constructing the approximation to φ. A detailed derivation of the primal–dual mixed method for a porous medium flow problem is presented in [7], and the formulation of the discrete system for (62) is using results from that paper. For the discrete formulation, define the following finite-dimensional spaces: Let Vu ⊂ H (div, Ω) be the space of vectorial basis functions of the simplest Raviart–Thomas type (the same as in the previous section). Also let Vφ ⊂ H 1 (Ω), the detailed description of this space will be given below. For the test spaces, let Vψ ⊂ L2 (Ω) be the space of piecewise constant functions on each element, and let Vv ⊂ (L2 (Ω))2 be a special vectorial test space, to be described below. The discrete formulation is then as follows: • Find uh ∈ Vu , φh ∈ Vφ such that 

Ω



Ω

h uh · vh dx + 2 h φh ψh dx + 2



∇φh · v h dx =

Ω





Ru · v h dx,

(63a)

Rφ ψh dx,

(63b)

Ω

φ
∇ · uh ψh dx =

Ω



Ω

for all v h ∈ Vv , ψh ∈ Vψ . We now discuss the spaces Vv and Vφ in more detail. Let v K ∈ Vv , and write vK =

4 

VK,j v j (x, y)

(64)

j =1

with v 1 = [−1 0]T , v 2 = [1 0]T , v 3 = [0 −1]T , and v 4 = [0 1]T . The supports of v j , j = 1, . . . , 4, are the triangles formed by the diagonals of an element such that v i points out of the element. We see that dim(Vu ) = dim(Vv ) and from (42) this means that Vφ and Vψ must have the same dimension. Consider now Vφ . The problem here is that the integral 

∇φh · v h dx

Ω

must be evaluated at a boundary edge, which means that φh must be differentiable there. We can solve this problem by the following construction: Let Pi be the barycenters of the elements, and let P
k be the midpoints of the boundary edges. We denote by χi the finite element “hat” functions with center at Pi . These functions are modified such that χi (P
k ) = 0, a restriction applying only to the “hat” functions for elements with one or more boundary edges. Moreover, let χ
k be “half-hat” functions centered at P
k such that χ
k (Pi ) = 0. Hence we can write: ∇φh =



φ(Pi )∇χi +

i



φ P
k ∇ χ
k ,

k

where we can use the inflow characteristic to determine φ(P
k ):



φ P
k = ψ3 P
k +





φ
uh · n P
k ,

I. Lie / Applied Numerical Mathematics 38 (2001) 445–474

467

which yields ∇φh =



φ(Pi )∇χi +

i

!



ψ3 P
k +



"



φ
uh · n P
k ∇ χ
k .

(65)

k

We have now apparently increased the dimension of Vφ by introducing the degrees of freedom P
k , but we will now show that these extra degrees of freedom are eliminated when we evaluate the integrals in the formulation. Let IP be the index set of χ -functions that intersect the element K with barycenter P . Similarly, let I
P be index set of χ
functions that intersect this element, and denote by v j the vectorial test function for an arbitrary edge. Then 

∇φh · v j dx =





φh (Pi )

∇χi · v j dx +





ψ3 P
k

k∈I
P K  

φ
uh · n P
k ∇ χ
k · vj dx. + k∈I
P K

i∈IP

K





∇ χ
k · v j dx

K

(66)

Consider now the different terms on the right hand side of (66): The first term contains the contributions from the neighboring elements of K, and hence only contributions from the ordinary degrees of freedom for φh . The second term does not contain any unknowns since ψ3 (P
k ) is assumed to be known and χ
k is given. In the third term we observe that due to the properties of the Raviart–Thomas elements the value of uh · n is constant along an edge, and hence uh · n(P
k ) is identical to the velocity unknown at the edge. This means that the third term will contribute to the coefficient matrix for velocities in much the same way as we saw for the boundary term in the dual mixed formulation. Thus we see that there are no explicit degrees of freedom for φ at P
k , and the requirement dim(Vφ ) = dim(Vψ ) is satisfied. The modifications of the above procedure in the case where u · n > φ
are as follows: We let φ(P
k ) = φ(Pj ), where Pj is the cell center of the cell for which P
k is a boundary point. This means that the extra degrees of freedom are eliminated without using information from the characteristics.
φ(P
k ) is known at the boundary since all characteristics are inflow. This means When −u · n > φ, again that the extra degrees of freedom are eliminated, and we do not have to use (66). The analysis of the primal–dual method is somewhat more involved, and we have use the following discrete bilinear forms: a(uh , v h ) =



Ω

b2 (uh , ψh ) =

uh · v h dx, 

b1 (φh , v h ) =



∇φh · v h dx,

(67a)

φh ψh dx,

(67b)

Ω

∇ · uh ψh dx,

c(φh , ψh ) =

Ω



Ω

so that (63a) and (63b) can be written: h a(uh , v h ) + b1 (φh , vh ) = (Ru , v), 2 h
c(φh , ψh ) + φb 2 (uh , ψh ) = (Rφ , ψh ). 2

(68a) (68b)

This formulation is within the framework of the generalized saddle-point problems discussed in [19,29].

468

I. Lie / Applied Numerical Mathematics 38 (2001) 445–474

We have the following existence theorem: Theorem 5. Assume that the condition (42) and the following inf-sup conditions (LBB conditions) hold: (a) (b) (c) (d)

a(uh , vh )
β1 > 0, uh ∈Vu v h ∈Vv uh H (div) v h L2 b1 (φh , vh )
β2 > 0, inf sup φh ∈Vφ v h ∈Vv φh H 1 v h L2 b2 (uh , ψh )
β3 > 0, inf sup u ψh ∈Vψ uh ∈Vu h H (div) ψh L2 c(φh , ψh )
β4 > 0. inf sup ψh ∈Vψ φh ∈Vφ φh H 1 ψh L2 inf sup

(69) (70) (71) (72)

Then the problem (68) has a unique solution. Proof. The only non-classical inf-sup condition here is (b), and this is shown in [7] using the so-called Fortin lemma and the properties of the Raviart–Thomas elements. The items (a), (c) and (d) can be shown using classical techniques in, e.g., [3,22]. ✷ We also have the error estimate: Theorem 6. Let u and φ be the solution of (62) and let uh and φh be the solution of (63). With u ∈ H 1 (Ω)2 , ∇ · u ∈ H 1 (Ω) and φ ∈ H 2 (Ω) we have:



u − uh H (div) + φ − φh H 1  C k |u|H 1 + |∇ · u|H 1 + |φ|H 2 ,

(73)

where C is a constant independent of the element size. Proof. See [29, Theorem 4].

✷

We see that for the primal–dual method the error estimate in φ is in the H 1 -norm compared to the dual mixed method where the error in φ was estimated only in the L2 -norm. The price we pay for this is the extra regularity requirement on φ: H 2 compared to H 1 . For the fully discrete system, we use the same global numbering of element and edge unknowns as for the dual mixed method, and we write the global linear system as follows: h A U + D1 Φ = R1 , 2 h φ
D U + C Φ = R2 , 2

(74a) (74b)

where A ∈ Rnu ×nu , D1 ∈ Rnu ×ne , D ∈ Rne ×nu , C ∈ Rne ×ne . As in the dual mixed case, the matrix A is block diagonal with tridiagonal blocks of size nxe + 1.

I. Lie / Applied Numerical Mathematics 38 (2001) 445–474

The blocks have the form: 

5

1

0

k   1 10 Ak =  24  0 1

1

 

0

0

469



0 

0 

 + diag(α1 , . . . , α4 ),

10

1 

1

5

(75)



where diag(α1 , . . . , α4 ) comes from the boundary terms in:  



φ
uh · n P
k ∇ χ
k · v j dx,

(76)

K

for details, see [7]. The matrix D is the discrete divergence matrix consisting of elements εki , k = 1, . . . , ne , i = 1, . . . , nu . Thus there will be four matrix elements in each row. The matrix D1 comes from the integral 

∇χi · v j dx,

(77)

K

and will be structurally symmetric to D, implying that each row will have two nonzero elements, except for a boundary edge which will have a single nonzero element. The values of the matrix elements in D1 will be different from those of D, as we see from the integral expression (77). The matrix C comes from the integral 

φh ψh dx,

(78)

K

and is diagonal since the support of the test function ψh is only a single element. The exact form of φh is given in [7, p. 9]. Hence the full linear system has a nonsymmetric, although structurally symmetric matrix. As for the dual mixed method, a nonsymmetric solver has to be used, and the different approaches to the solution of the linear system described briefly in the previous section can also be used here. Consider again an example of edge and element equations. The equation for a boundary edge 1 of element K is: (a11 + α1 )uK1 + a12 uK2 + d11 φK = RuK1 .

(79)

Note that in this case we do not have any boundary integrals, so the boundary contributions will only appear through the matrix element α1 . The equation for element K is: 

−uK1 + uK2 −uK3 + uK4 + 2k k k

h φ




+

c11 1 φ1 = 2 RφK , 2 k k

(80)

which is quite similar to (55). The strategy for the computation of the right hand sides in (63a) and (63b) is almost the same as for the dual-mixed method in the previous section, and will not be repeated here.

470

I. Lie / Applied Numerical Mathematics 38 (2001) 445–474

In the same way as for the dual mixed method, the existence theorem implies the following stability result which expresses the well-posedness: Corollary 7. Let uh and φh be the solution of (68). Then there exist constants c11 , c12 , c21 , c22 independent of the discretization such that 







uh H (div)  c11 Ru H 1 + c12 Rφ H 1 ,

(81a)

φh H 1

(81b)

   c21 Ru 

H1

  + c22 Rφ 

H1.

Proof. See [19]. We proceed as in the proof of Corollary 4, and hence   Ru   d11 Ru  + d12 u,   Rφ   d21 Rφ  + d22 u,

(82a) (82b)

where Ru and Rφ refer to quantities not containing values at time tn+1 , and where the norm in this case is the H 1 -norm provided ∇ · u ∈ H 1 (Ω) and φ ∈ H 2 (Ω). Hence we have the required stability result. ✷ Note that for the primal–dual method we can use the H 1 -norm throughout in the estimates, which is a definite improvement over the dual mixed method where L2 -norms were used in the estimates for φ.

6. Numerical illustration In this section we present a few numerical results in order to illustrate the well-posedness results presented above. The code used for the experiments is using a dual mixed method and with characteristic boundary conditions imposed weakly. The method and the computational details, along with many more experiments with transparent boundaries, is described in [8]. The results presented here is from an adjustment test where we consider a 10000 × 10000 km domain with gridsize either 100 km or 200 km. The reference geopotential is 5000 g, and a Gaussian peak of size 500 g is introduced at t = 0. The timestep is 3600 s. The system is initially completely out of geostrophic balance and will adjust itself by radiating gravity waves. What we should expect is that if the transparent boundary conditions are efficient, the geopotential field becomes smaller and the velocities in H (div)-norm decreases. In Fig. 4 we present the peak geopotential and the rms divergence for a 240 h run. We see clearly that the peak geopotential decreases almost monotonely (the value at 240 h is 11.04), and the divergence has a nice exponential decay with time. Both curves are nonsmooth because of the nonlinearities in the equations. In order to see if the performance of the transparent boundary is very sensitive to the gridsize, we compare runs with 50 × 50 and 100 × 100 resolutions. We see from Fig. 5 that the sensitivity of the results, in terms of maximal geopotential and velocity divergence, are very small with respect to resolution. Of course we cannot decrease the resolution too much, but within reasonable bounds the performance of the transparent boundary conditions seems to be robust against variations in the resolution.

I. Lie / Applied Numerical Mathematics 38 (2001) 445–474

Fig. 4. Maximum geopotential and rms divergence of wind field ·108 , vs. forecast time.

Fig. 5. Maximum geopotential and rms divergence of wind field ·108 , vs. forecast time for two resolutions.

471

472

I. Lie / Applied Numerical Mathematics 38 (2001) 445–474

7. Discussion and conclusion In this paper we have shown how two versions of the mixed finite element method for space discretization, combined with semi-Lagrangian time integration can incorporate characteristic boundary conditions to give a well-posed problem. The numerical implementation of this method, done in [23,25], gives good results for the tests run so far, despite the fact that the characteristic boundary conditions are only a fairly rough approximation to the exact transparent boundary condition. We have also seen another advantage of the mixed finite element formalism, namely that we have quite a large body of theory we can use in order to prove well-posedness results and even error estimates. Note that these error estimates express the error bounds for the discrete solution with given, approximate transparent boundary conditions. What we want to obtain is an error estimate which also contains the error we commit when we approximate the exact transparent boundary condition. This goal may be somewhat ambitious since there are no results in the literature that contains such complete error estimates. Although classical finite difference discretization on a staggered grid and the mixed finite element method use the same grid, the two discretizations are also of different character, for example in the way the boundary conditions are imposed. This implies that comparisons of the discrete equations for a boundary edge done in the difference scheme and the mixed method is not totally fair. It is tempting to suggest that mixed finite elements, or the related finite volume method should be used in shallow water models or models governed by related equations since a range of different boundary conditions can easily be incorporated. It is too early to draw this conclusion now, at least one will have to investigate how the mixed finite element method can be combined with more sophisticated approximate transparent boundary conditions. Moreover, changing the basic discretization in an existing model is not done overnight, so a more practical goal may be to see how the existing discretization may benefit from other approaches, implying changes that will affect only parts of a model. The methods for space and time integration described in this paper are of course not limited to PDE systems or models using shallow water-like equations. We have already mentioned models for flow in porous media, see [7], and full 3D models for geophysical flows like the one presented in [2] can also be treated.

Acknowledgements The author would like to thank Aidan McDonald for careful reading of the paper, resulting in several improvements. We have benefited from the works of Astrid Holstad, who has derived the details of the discrete primal– dual method and done a successful implementation of the method for a porous media flow problem. She also developed the dual mixed finite element code used for the numerical illustration.

References [1] T. Alsos, Effective ODE-solvers for trajectory calculations in semi-Lagrangian methods used in weather forecasting, M.Sc. thesis, Department of Mathematical Sciences, The Norwegian University of Science and Technology, Trondheim, Norway, 1998.

I. Lie / Applied Numerical Mathematics 38 (2001) 445–474

473

[2] L. Bonaventura, A semi-implicit semi-Lagrangian scheme for the compressible nonhydrostatic equations of the atmosphere, GKSS Forschungszentrum, External Report, Institut für Gewässerphysik, Geesthacht, Germany, 1998. [3] F. Brezzi, M. Fortin, Mixed and Hybrid Finite Element Methods, Springer, Berlin, 1991. [4] B. Engquist, A. Majda, Absorbing boundary conditions for the numerical simulation of waves, Math. Comp. 31 (1977) 629–651. [5] M. Falcone, R. Ferretti, Convergence analysis for a class of high-order semi-Lagrangian advection schemes, SIAM J. Numer. Anal. 35 (1998) 909–940. [6] V. Girault, P.-A. Raviart, Finite Element Methods for Navier–Stokes Equations, Springer, Heidelberg, 1986. [7] A. Holstad, Temperature-driven fluid flow in porous media using a mixed finite element method and a finite volume method, Adv. Water Resources, 2001, to appear. [8] A. Holstad, I. Lie, Transparent boundary conditions using a mixed finite element formulation of the shallow water equation, Research Report No. 120, Norwegian Meteorological Institute, Oslo, Norway, 2001. [9] A. Jeffrey, Quasilinear Hyperbolic Systems and Waves, Research Notes in Mathematics, Vol. 5, Pitman, London, 1976. [10] H.-O. Kreiss, J. Lorenz, Initial-Boundary Value Problems and the Navier–Stokes Equations, Academic Press, San Diego, 1989. [11] I. Lie, On the stability and accuracy of semi-Lagrangian schemes, Research Report No. 72, Norwegian Meteorological Institute, Oslo, Norway, 1998. [12] I. Lie, Variants of semi-Lagrangian schemes for parallel processing, FFI/Report-96/00884, FFI, Kjeller, Norway, 1996. [13] I. Lie, The use of higher order transparent boundary conditions in shallow water models, Manuscript, Norwegian Meteorological Institute, Oslo, 2001. [14] A. McDonald, Well-posed boundary conditions for semi-Lagrangian schemes, Part I: The one-dimensional case, HIRLAM Technical Report No. 43, Met Eirann, Dublin, Ireland, 1999. [15] A. McDonald, Well-posed boundary conditions for semi-Lagrangian schemes, Part II. HIRLAM Technical Report, Met Eirann, Dublin, Ireland, 2000. [16] A. McDonald, Well-posed boundary conditions for semi-Lagrangian schemes: The two-dimensional case, HIRLAM Technical Report No. 47, January 2001, HIRLAM-5 project, SMHI, Norrköking, Sweden; also: A. McDonald, A step toward transparent boundary conditions for meteorological models, Technical Note 57, Met Eirann, Dublin, Ireland. Submitted to Mon. Wea. Rev. [17] F. Mesinger, Forward-backward scheme, and its use in a limited-area model, Contrib. Atmos. Phys. 50 (1977) 200–210. [18] J. Necas, Les Méthodes Directes en Théorie des Équations Elliptiques, Masson, Paris, 1967. [19] R.A. Nicolaides, Existence, uniqueness and approximation for generalized saddle point problems, SIAM J. Numer. Anal. 19 (1982) 349–357. [20] A. Quarteroni, A. Valli, Domain Decomposition Methods for Partial Differential Equations, Oxford University Press, Oxford, 1999. [21] J.B. Rauch, F.J. Massey, Differentiability of solutions to hyperbolic initial-boundary value problems, Trans. Am. Math. Soc. 189 (1974) 303–318. [22] J.E. Roberts, J.M. Thomas, Mixed and hybrid methods, in: P.G. Ciarlet, J.-L. Lions (Eds.), Handbook of Numerical Analysis, Vol. II, North-Holland, Amsterdam, 1991. [23] I. Lie, R. Skålin, Characteristic-based methods and domain decomposition for atmospheric models, in: Proceedings ENUMATH’95 Conference, Paris, 1995. [24] R. Skålin, I. Lie, Stability of semi-Lagrangian schemes with time interpolation, manuscript, 1996; available from I. Lie. [25] R. Skålin, Selected parallel algorithms for numerical weather prediction, Ph.D. Thesis, Department of Mathematical Sciences, Norwegian University of Science and Technology, Trondheim, Norway, 1997.

474

I. Lie / Applied Numerical Mathematics 38 (2001) 445–474

[26] P.-A. Raviart, J.-M. Thomas, A mixed finite element method for 2nd order elliptic problems, in: I. Galligani, E. Magenes (Eds.), Mathematical Aspects of the Finite Element Method, Lecture Notes in Mathematics, Vol. 606, Springer, New York, 1977, pp. 292–315. [27] A. Staniforth, J. Coté, Semi-Lagrangian integration schemes for atmospheric models—A review, Mon. Wea. Rev. 119 (1991) 2206–2223. [28] D. Trujillo, J.M. Thomas, Finite volume variational formulation. Applications to domain decomposition methods, Contemp. Math. 157 (1994) 127–132. [29] D. Trujillo, J.M. Thomas, Analysis of finite volumes method, in: S. Luckhaus, A. Bourgeat, C. Carasso, A. Mikelic (Eds.), Mathematical Modeling of Flow through Porous Media, World Scientific, Singapore, 1995, pp. 318–336. [30] L. Wagatha, Approximation of pseudodifferential operators in absorbing boundary conditions for hyperbolic equations, Numer. Math. 42 (1983) 51–64.