A new Algorithm Based on Factorization for Heterogeneous Domain Decomposition

# A new Algorithm Based on Factorization for Heterogeneous Domain Decomposition

Martin J. Gander Section de mathématiques, Université de Genève, 2-4 rue du Lièvre, CP 64, CH-1211 Genève 4, Switzerland. Martin.Gander@unige.ch    Laurence Halpern LAGA, Université Paris XIII, 99 Avenue J.-B. Clément, 93430 Villetaneuse, France. halpern@math.univ-paris13.fr    Véronique Martin LAMFA UMR-CNRS 7352, Université de Picardie Jules Verne, 33 Rue St. Leu, 80039 Amiens, France. veronique.martin@u-picardie.fr
###### Abstract

Often computational models are too expensive to be solved in the entire domain of simulation, and a cheaper model would suffice away from the main zone of interest. We present for the concrete example of an evolution problem of advection reaction diffusion type a heterogeneous domain decomposition algorithm which allows us to recover a solution that is very close to the solution of the fully viscous problem, but solves only an inviscid problem in parts of the domain. Our new algorithm is based on the factorization of the underlying differential operator, and we therefore call it factorization algorithm. We give a detailed error analysis, and show that we can obtain approximations in the viscous region which are much closer to the viscous solution in the entire domain of simulation than approximations obtained by other heterogeneous domain decomposition algorithms from the literature.

Key words. Heterogeneous domain decomposition, viscous problems with inviscid approximations, transmission conditions, factorization algorithm

AMS subject classifications. 65M55, 65M15

## 1 Introduction

The coupling of different types of partial differential equations is an active field of research, since the need for such coupling arises in various applications. A first main area is the simulation of complex objects, composed of different materials, which are naturally modeled by different equations; fluid-structure interaction is a typical example, and many techniques have been developed for this type of coupling problems, see for example the book [33], or the review on the immersed boundary method [32], and [9] for domain decomposition coupling techniques. A very important area of application is the simulation of the cardiovascular system [16]. A second main area is when homogeneous objects are simulated, but the partial differential equation modeling the object is too expensive to solve over the entire object, and a simpler, less expensive model would suffice in most of the object to reach the desired accuracy; air flow around an airplane is a typical example, where viscous effects are important close to the airplane, but can be neglected further away, see the early publication [10], and also [7] and the references therein. An automatic approach for neglecting the diffusion in parts of the domain is the -formulation, see [27] [5], and there are also techniques based on virtual control, originating in [11], see [1] for the case with overlap, and [25] for the case without, and also [12] for virtual control with variational coupling conditions. A third emerging area is the coupling of equations across dimensions, for example the blood flow in the artery can be modeled by a one dimensional model, but in the heart, it needs to be three dimensional, see for example [15]. All these techniques have become known in the domain decomposition community under the name heterogeneous domain decomposition methods, a terminology sparked by the review [36], and the literature has become vast in this field.

We are interested in this paper in the second situation, where the motivation for using different equations comes from the fact that we would like to use simpler, less expensive equations in areas of the domain where the full model is not needed, and we use as our guiding example the advection reaction diffusion equation. We are in principle interested in the fully viscous solution, but we would like to solve only an advection reaction equation for computational savings in part of the domain. Coupling conditions for this type of problem have been developed in the seminal paper [23], but with the first situation described above in mind, i.e. there is indeed a viscous and an inviscid physical domain, and the coupling conditions are obtained by a limiting process as the viscosity goes to zero, see also [24], and [3, 8] for an innovative correction layer, and [6] for the steady case.

Dubach developed in his PhD thesis [13] coupling conditions based on absorbing boundary conditions, and such conditions have been used in order to define heterogeneous domain decomposition methods in [18]. A fundamental question however in the second situation described above is how far the solution obtained from the coupled problem is from the solution of the original, more expensive one on the entire domain. A first comparison of different transmission conditions focusing on this aspect appeared in [19]. In [20], coupling conditions were developed for stationary advection reaction diffusion equations in one spatial dimensions, which lead to solutions of the coupled problem that can be exponentially close to the fully viscous solution, and rigorous error estimates are provided. The coupling conditions are based on the factorization of the differential operator, see also [29], and the exact factorization can be used in this one dimensional steady case. We study in this paper time dependent advection reaction diffusion problems, where the exact factorization of the differential operator can not be used any more, due to the non-local nature of the factors, and new ideas are needed in order to obtain better coupling conditions than the classical ones developed for situation one, where the domains are really physically different.

We present in Section LABEL:Sec2 our new factorization algorithm. In Section LABEL:sec:wpalgo1 and LABEL:sec:wpalgo2 we give a detailed error analysis of our algorithm, and prove asymptotic error estimates when the viscosity is becoming small. In Section LABEL:SecNum we present numerical experiments which show that our theoretical error estimates are sharp, and that the new factorization algorithm gives approximate solutions which are one order of magnitude more accurate in the viscous region than the best heterogeneous domain decomposition methods known from the literature.

## 2 A new coupling algorithm based on factorization

We now explain how the factorization technique that led to coupling conditions of excellent quality for one dimensional problems in [20] can be used to obtain a new coupling algorithm for evolution problems which we will call factorization algorithm.

### 2.1 Model problem

We consider the time dependent advection reaction diffusion equation

 \@fontswitchLadu:=∂tu−ν∂2xu+a∂xu+cu=fin (−L1,L2)×(0,T),\@fontswitchB1u(−L1,⋅)=g1on (0,T),\@fontswitchB2u(L2,⋅)=g2on (0,T),u(x,0)=hin (−L1,L2), \hb@xt@.01(2.1)

where is the velocity field, is the viscosity, and is a reaction term. The , are suitable boundary operators, representing Dirichlet or absorbing boundary conditions. We consider the following choice of these operators, depending on the sign of the advection term ,

 \@fontswitchB1\@fontswitchB2a>0Id∂t+a∂x+ca<0IdId \hb@xt@.01(2.2)

In the case , the flow is given at the inflow boundary, and an absorbing boundary condition is prescribed at the outflow boundary. This can be compared to the situation of the tail of a wing, where the flow goes from the complicated model region into the simplified model region. For negative , the flow is prescribed at the inflow and outflow boundary, which can be compared to the situation of the front of the wing, where the flow goes from the simplified model region into the complicated model region, and a boundary layer forms.

We assume that the initial data is compactly supported in . The forcing term is compactly supported in , and the boundary values and are compactly supported in . Regularity and compatibility conditions on the data need to be enforced to have a sufficiently regular solution, see Section LABEL:sec:wpalgo1.

### 2.2 The new algorithm based on factorization

Using Nirenberg’s factorization, we can factor the advection-diffusion operator into a product of two evolution operators in opposite directions,

 −ν∂2x+a∂x+c+∂t=−ν(∂x+\@fontswitchP+(∂t))(∂x+\@fontswitchP−(∂t))(mod\@fontswitchC∞).

This factorization has been used to design absorbing boundary conditions and paraxial equations for hyperbolic problems, see [2]. For parabolic problems, Nataf and coauthors [29, 34] computed approximations of via a double sweep, and also obtained transmission conditions for Schwarz domain decomposition methods [35], which led to the new class of optimized Schwarz methods, see [17] for an overview. The same factorization can also be used to obtain incomplete LU preconditioners [21, 22], and is the underlying mathematical structure of the recently developed sweeping preconditioner [14]. We now use this factorization to define our new factorization algorithm: we define two subdomains, and , and want to couple the advection-diffusion equation in with an advection equation in , defined by the transport operator . Our goal is to obtain a coupled solution which is as close as possible to the fully viscous solution of the original problem.

We start with the case , where according to (LABEL:eq:defbo) the exterior boundary condition is . Suppose there exists a decomposition with a transport operator propagating to the right, and a transport operator propagating to the left. The original problem

 ⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩˜\@fontswitchLma\@fontswitchLau=f in Ω×(0,T),u(−L1,⋅)=g1 on (0,T),\@fontswitchLau(L2,⋅)=g2 on (0,T),u(⋅,0)=h in Ω

can then be solved by introducing , and solving the two problems

 ⎧⎪ ⎪⎨⎪ ⎪⎩˜\@fontswitchLmauma=f in Ω2×(0,T),uma(L2,⋅)=g2 on (0,T),uma(⋅,0)=\@fontswitchLau(⋅,0) in Ω2, ⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩˜\@fontswitchLma\@fontswitchLauad=f in Ω1×(0,T),uad(−L1,⋅)=g1 on (0,T),\@fontswitchLauad(0,⋅)=uma(0,⋅) on (0,T),uad(⋅,0)=h in Ω1, \hb@xt@.01(2.3)

which leads to . Unfortunately, the exact factorization is expensive, but we can use an approximation with a remainder,

 \@fontswitchLad=νa2(\@fontswitchLma\@fontswitchLa−\@fontswitchR) with \@fontswitchR=(∂t+c)2 and \@fontswitchLma=∂t−a∂x+c+a2ν. \hb@xt@.01(2.4)

The viscous solution satisfies , and the algorithm corresponding to (LABEL:algoparfait) is

Since is unknown to evaluate the remainder, we approximate it by solving an advection equation, and our new factorization algorithm is

 ⎧⎪ ⎪⎨⎪ ⎪⎩\@fontswitchLauka=f in Ω2×(0,T),uka(0,⋅)=uk−1ad(0,⋅) on (0,T),uka(⋅,0)=h in Ω2, ⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩\@fontswitchLmaukma=a2νf+\@fontswitchRuka in Ω2×(0,T),ukma(L2,⋅)=g2 on (0,T),ukma(⋅,0)=f(⋅,0)+νd2xh in Ω2, \hb@xt@.01(2.5) ⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩\@fontswitchLadukad=f in Ω1×(0,T),ukad(−L1,⋅)=g1 on (0,T),\@fontswitchLaukad(0,⋅)=ukma(0,⋅) on (0,T),ukad(⋅,0)=h in Ω1,

where we start with a given initial guess . We will prove well posedness of this algorithm in Section LABEL:sec:wpalgo1, and give precise error estimates when is small, which show that the new factorization algorithm gives one and a half orders of magnitude better solutions in the viscous subregion than the best other coupling algorithms from the literature.

When , we have the factorization with remainder in reverse order, , and now the operator propagates to the left, and to the right. The viscous solution satisfies , and introducing , the algorithm corresponding to (LABEL:algoparfait) is

Since is unknown to evaluate the remainder and the boundary conditions, we approximate it again by solving an advection equation, and our new factorization algorithm becomes

 ⎧⎪⎨⎪⎩\@fontswitchLau1a=f in Ω2×(0,T),u1a(L2,⋅)=g2 on (0,T),u1a(⋅,0)=h in Ω2, ⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩\@fontswitchLau2a=a2νf+\@fontswitchRu1a in Ω2×(0,T),u2a(L2,⋅)=\@fontswitchLmau1a(L2,⋅) on (0,T),u2a(⋅,0)=\@fontswitchLmau1a(⋅,0) in Ω2, \hb@xt@.01(2.6) ⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩\@fontswitchLaduad=f in Ω1×(0,T),uad(−L1,⋅)=g1 on (0,T),\@fontswitchLmauad(0,⋅)=u2a(0,⋅) on (0,T),uad(⋅,0)=h in Ω1,

where one could also directly compute and . There is no iteration for in the algorithm, because the boundary condition at in the first step can not be updated naturally from the viscous solution in . We will study this algorithm in detail in Section LABEL:sec:wpalgo2, and show that it gives an order of magnitude better solutions in the viscous subregion than the other coupling algorithms from the literature.

### 2.3 Well-posedness results for advection reaction diffusion problems

We work in the usual Sobolev spaces in time and space, and for , in the hyperbolic case, and the anisotropic spaces in the parabolic case. For clarity, we will add an index defining time or space in the Sobolev space, for instance . We introduce for any domain the anisotropic Sobolev spaces (see [28])

 Hr,s(Ω×(0,T))=L2(0,T;Hr(Ω))∩Hs(0,T;L2(Ω)). \hb@xt@.01(2.7)

If is in , then for any integer and , we have

 ∂j∂xj∂k∂tku∈Hμ,ν(Ω×(0,T)),whereμr=νs=1−(jr+ks). \hb@xt@.01(2.8)

We introduce the space of traces of functions in for the half-space (and similarly for ). Denoting by the trace of the k-th derivative in time on the initial line, , and by the trace of the j-th derivative in space on the boundary , , the trace space is defined by

 Vr,s:={(fk,gj)∈∏k
###### Theorem 2.1 ([28])

For positive real numbers and such that , the trace map

 u↦{{∂ku∂tk(x,0)}k

is defined and continuous from onto .

We start with well-posedness results for the advection equation, by stating a general result, applicable to in or , and in . To this end, we introduce and consider

 \@fontswitchLbv:=∂tv+b∂xv+ηv=p in \@fontswitchO×(0,T). \hb@xt@.01(2.11)

Let be the spatial part of the operator , i.e. . We denote the boundary point where the flux enters the domain by , the other boundary point by , and define the characteristic time . If , and , and if , and . Note that is a continuous function of . We therefore equip (LABEL:eq:genadv) with initial and boundary conditions

 v(⋅,0)=h,v(x−,⋅)=g. \hb@xt@.01(2.12)

The following well-posedness result can be found in [31].

###### Theorem 2.2

If , and , then the transport problem (LABEL:eq:genadv,LABEL:eq:genadvdata) has a unique weak solution , given by (the characteristic function of in is denoted by )

 v(x,t) =h(x−bt)e−ηt1t<τ(x)+g(t−τ(x))e−ητ(x)1t>τ(x) \hb@xt@.01(2.13) +∫t(t−τ(x))+p(x−b(t−s),s)e−η(t−s)ds.

If for some we have , and , with the compatibility conditions

 dktg(0)=(k−1∑j=0(−\@fontswitchMb)j∂k−1−jtp)(x−,0)+(−\@fontswitchMb)kh(x−)for0≤k≤γ−1, \hb@xt@.01(2.14)

then and . Furthermore, we have for the estimates

 η∥∂ktv∥2L2x,t+|b|∥∂ktv(x+,⋅)∥2L2t≤1η∥∂ktp∥2L2x,t+∥∂ktv(⋅,0)∥2L2x+|b|∥dktg∥2L2t. \hb@xt@.01(2.15)

Similarly, we also use well-posedness results for the advection reaction diffusion equation

 \@fontswitchLadu:=∂tu−ν∂2xu+a∂xu+cu=fin \@fontswitchO×(0,T),\@fontswitchB1u(x1,⋅)=g1on (0,T),\@fontswitchB2u(x2,⋅)=g2on (0,T),u(x,0)=hin \@fontswitchO, \hb@xt@.01(2.16)

with boundary operators according to (LABEL:eq:defbo). We define to be the spatial part of the operator , i.e. .

###### Theorem 2.3

For , let , , for negative , and for positive , , with the compatibility conditions for and for negative given by

then problem (LABEL:eq:genadvdiff) has a unique solution in .

Proof. Existence and regularity results are well-known for Dirichlet boundary conditions on both sides, see [28, 31], so we do not consider the case of negative advection further. In [31] more precise results with error bounds in for the hyperbolic equation (see Theorem LABEL:th:parabolic2) can be found. In the case where , due to the absorbing boundary, we need to modify the proof on the right boundary, and we use a Fourier transform in time. A weak solution is obtained by a variational formulation, like in [30, 4] for instance. The regularity is obtained as follows: we first modify the boundary condition in (LABEL:eq:genadvdiff) on the right at to Dirichlet data . Because of the compatibility conditions on the left, and imposing symmetric compatibility conditions on on the right, there is a unique solution , see [28]. The difference is solution of the homogeneous case of equation (LABEL:eqadvisqapos), but the boundary condition on the right becomes . By the regularity results above, is in . To estimate , we will make use of the Fourier transform. We extend all functions by in , and smoothly into , and define

 ^v(ω)=12π∫Re−iωtv(t)dt.

Since the initial value vanishes, the equation is Fourier transformed in time to

This is for each an ordinary differential equation, with characteristic roots

 λ+(ω)=12ν(a+√a2+4ν(c+iω)),λ−(ω)=12ν(a−√a2+4ν(c+iω)), \hb@xt@.01(2.19)

with and . The general solution is

 ^v(x,ω)=ℓ+(ω)eλ+x+ℓ−(ω)eλ−x.

Using the boundary conditions, we then get the solution

 ^v(x,ω)=^q2(ω)eλ+(x−x1)−eλ−(x−x1)νλ2+eλ+(x2−x1)−νλ2−eλ−(x2−x1), \hb@xt@.01(2.20)

where we have used the relation . The value at can be equivalently written as

 ^v(x2,ω)=^q2(ω)e−(λ+−λ−)(x2−x1)−1νλ2+((λ−λ+)2e−(λ+−λ−)(x2−x1)−1). \hb@xt@.01(2.21)

In order to estimate the regularity of , we need to estimate the multiplicative factor on the right for large . We can see from (LABEL:deflambda) that . Therefore Since , we conclude that . Then is solution of the advection-diffusion equation with Dirichlet boundary conditions, and the data has sufficient regularity to conclude.

## 3 Properties of the factorization algorithm for positive advection

We consider the advection-diffusion equation in with Dirichlet boundary condition on the left, and absorbing boundary condition given by the transport operator on the right (see [26]),

 \@fontswitchLadu:=∂tu−ν∂2xu+a∂xu+cu=fin Ω×(0,T),u(−L1,⋅)=g1on (0,T),\@fontswitchLau(L2,⋅)=g2on (0,T),u(⋅,0)=hin Ω. \hb@xt@.01(3.1)

We suppose in this section that and are compactly supported in , that is compactly supported in