Two-component domain decomposition scheme with overlapping subdomains for parabolic equations

TL;DR

This paper introduces a novel domain decomposition method for parabolic equations using overlapping subdomains and indicator functions, generalizing ADI schemes with proven unconditional stability and demonstrated through numerical experiments.

Contribution

It proposes a new overlapping subdomain scheme based on indicator functions and a three-component splitting, extending ADI methods with stability analysis.

Findings

Unconditional stability conditions established for the new schemes.

Numerical results confirm the effectiveness of the method.

Generalization of ADI schemes to three-component splitting.

Abstract

An iteration-free method of domain decomposition is considered for approximate solving a boundary value problem for a second-order parabolic equation. A standard approach to constructing domain decomposition schemes is based on a partition of unity for the domain under the consideration. Here a new general approach is proposed for constructing domain decomposition schemes with overlapping subdomains based on indicator functions of subdomains. The basic peculiarity of this method is connected with a representation of the problem operator as the sum of two operators, which are constructed for two separate subdomains with the subtraction of the operator that is associated with the intersection of the subdomains. There is developed a two-component factorized scheme, which can be treated as a generalization of the standard Alternating Direction Implicit (ADI) schemes to the case of a special three-component splitting. There are obtained conditions for the unconditional stability of regionally additive schemes constructed using indicator functions of subdomains. Numerical results are presented for a model two-dimensional problem.

Equations

$$\mathcal{A}u=-\nabla (k(\mathbf{x})\nabla u)+c(\mathbf{x})u,\mathbf{x}\in \Omega$$

$$u(\mathbf{x})=0,\mathbf{x}\in \partial \Omega .$$

$$k(\mathbf{x})\ge \kappa >0,c(\mathbf{x})\ge 0,\mathbf{x}\in \Omega .$$

$$\frac{du}{dt}+\mathcal{A}u=f(t),0<t\le T,$$

$$u(0)=u^0,$$

$$(u,v)={\int }_{\Omega }u(\mathbf{x})v(\mathbf{x})d\mathbf{x},\Vert u\Vert =(u,u{)}^{1/2}.$$

$$d(u,v)=d(v,u),d(u,u)\ge \delta \Vert u{\Vert }^2,\delta >0,$$

$$(u,v{)}_d=d(u,v),\Vert u{\Vert }_d=(d(u,u){)}^{1/2}.$$

$$H_0^1(\Omega )=\{v(\mathbf{x})\in H^1(\Omega ):v(\mathbf{x})=0,\mathbf{x}\in \partial \Omega \}.$$

$$\left(\frac{du}{dt},v\right)+a(u,v)=(f,v),\forall v\in H_0^1(\Omega ),0<t\le T.$$

$$a(u,v)={\int }_{\Omega }(k\nabla u\cdot \nabla v+cuv)d\mathbf{x},$$

$$a(u,v)=a(v,u),a(u,u)\ge \delta \Vert u{\Vert }^2,\delta >0.$$

$$(u(0),v)=(u^0,v),\forall v\in H_0^1(\Omega ).$$

$$\Vert u(t)\Vert \le \Vert u^0\Vert +{\int }_0^t\Vert f(\theta )\Vert d\theta .$$

$$\Vert u\Vert \frac{d}{dt}\Vert u\Vert +a(u,u)=(f,u).$$

$$(f,u)\le \Vert f\Vert \Vert u\Vert ,$$

$$\frac{d}{dt}\Vert u(t)\Vert \le \Vert f(t)\Vert .$$

$$(Ay,v)=a(y,v),\forall y,v\in V^h.$$

$$A=A^{\ast }\ge {\delta }_{h}I,{\delta }_h>0.$$

$$\frac{dw}{dt}+Aw=\phi (t),0<t\le T,$$

$$w(0)=w^0,$$

$$\left(\frac{dw}{dt},w\right)+(Aw,w)=(\phi ,w).$$

$$\left(\frac{dw}{dt},w\right)\le (\phi ,w).$$

$$(\phi ,w)\le \Vert \phi \Vert \Vert w\Vert .$$

$$\frac{d}{dt}\Vert w\Vert \le \Vert \phi \Vert .$$

$$\Vert w(t)\Vert \le \Vert w^0\Vert +{\int }_0^t\Vert \phi (\theta )\Vert d\theta .$$

$$y^{n+\sigma }=\sigma y^{n+1}+(1-\sigma )y^n.$$

$$\frac{y^{n+1}-y^n}{\tau }+A{y}^{n+\sigma }={\phi }^{n+\sigma },n=0,1,...,N-1,$$

$$y^0=w^0.$$

$$I+\tau \left(\sigma -\frac{1}{2}\right)A\ge 0$$

Full text

††journal: arXiv.org

1 Introduction

Schemes of domain decomposition are considered for numerical solving time-dependent problems to partial differential equations. Iteration-free algorithms of domain decomposition take into account specific features of time-dependent problems in the most efficient way. In some cases, it is possible (see, e.g., Kuznetsov (1988, 1991)) without loss of accuracy of the approximate solution make only one iteration of the Schwarz alternating method at a new time level in solving boundary value problems for a parabolic second-order equation. Iteration-free schemes of domain decomposition are associated with various variants of additive schemes (splitting schemes) – see regionally additive schemes in Samarskii et al. (2002).

Methods of domain decomposition for solving unsteady problems can be classified by (I) the method of decomposition for a calculation domain, (ii) the choice of decomposition operators (exchange boundary conditions), and (III) the splitting scheme (approximation in time) employed. For multidimensional boundary value problems, it is possible to use domain decomposition methods with or without overlapping subdomains (Quarteroni and Valli, 1999; Toselli and Widlund, 2005). Methods without overlapping subdomains are connected with an explicit formulation of exchange data conditions on subdomain interfaces.

To construct decomposition operators for solving unsteady initial boundary value (IBV) problems for partial differential equations, it is convenient to use a partition of unity for a computational domain (Dryja, 1991; Laevsky, 1987; Samarskii and Vabishchevich, 1995; Vabishchevich, 1989, 1994a, 1994b). In decomposition methods with overlapping subdomains, the functions are associated with individual subdomains and take a value between zero and one. Results of studies on domain decomposition methods for Cauchy problems for partial differential equations are summarized in the book Samarskii et al. (2002). Among more recent studies, we highlight the works Vabishchevich (2008, 2011), where schemes of domain decomposition more suitable for numerical implementation are presented.

The construction of regionally additive schemes and study of their convergence are carried out on the basis of the general theory of splitting schemes (Samarskii, 2001; Marchuk, 1990; Vabishchevich, 2014). We highlight the simplest case of two-component splitting. In this case, we obtain unconditionally stable factorized splitting schemes, such as classical methods of alternating directions, predictor-corrector schemes and so on. A more interesting for computational practice is the situation, when a problem operator is divided into a sum of three or more noncommutative non-self-adjoint operators. In the case of such a multicomponent representation, splitting schemes are constructed on the basis of the concept of summarized approximation. For parallel computers, additively averaged splitting schemes are of particular interest. In the class of splitting schemes with full approximation (Vabishchevich, 2014), we highlight vector additive schemes, when the original equation is transformed into a system of similar equations (Abrashin, 1990; Vabishchevich, 1996; Abrashin and Vabishchevich, 1998). The most convenient approach for constructing additive operator–difference schemes of multicomponent splitting is based on regularization of difference schemes (Samarskii, 1967), where stability is achieved via perturbations of operators of the difference scheme.

The regionally additive schemes constructed on the basis of a partition of unity have certain defects. The most important among them is connected with the use of the uniquely defined operators with the generating coefficients in the overlapping domains. This leads to the fact that, for example, in the problems with constant coefficients of the equations one needs to use an algorithm with varying coefficients. In the paper Vabishchevich (2008), we constructed unconditionally stable regionally additive schemes with overlapping subdomains, which are more comfortable for the practical use than traditional schemes designed on the basis of a partition of unity. The regularization principle for operator–difference schemes (Samarskii, 1967) makes possible to construct schemes of component-wise splitting. In the present work, domain decomposition methods with overlapping subdomains are designed using a two-component splitting, which generalize standard factorized schemes (Samarskii et al., 2002; Vabishchevich, 2014). Decomposition operators are constructed on the basis of the indicator functions of the subdomains.

The paper is organized as follows. A model problem for a parabolic equation with self-adjoint elliptic operator of second order is formulated in Section 2. Approximation in space is constructed using Lagrangian finite elements, whereas approximation in time is based on conventional two-level schemes with weights. In Section 3, a two-component scheme of domain decomposition with overlapping subdomains is constructed on the basis of a partition of unity for the computational domain. New schemes of domain decomposition developed using indicator functions of subdomains are proposed in Section 4. They are based on the generalization of classical factorized schemes (operator analogs of ADI schemes). In Section 5, numerical experiments on the accuracy of the domain decomposition schemes are discussed for the model IBV problem. The results of the work are summarized in Section 6.

2 Problem formulation

Let $\Omega$ be a bounded domain ($\Omega\subset\mathbb{R}^{m},m=2,3$) with a piecewise smooth boundary $\partial\Omega$. We define an elliptic operator $\mathcal{A}$ so that

[TABLE]

on the set of functions

[TABLE]

Let $k(\bm{x})$ and $c(\bm{x})$ be smooth functions in $\overline{\Omega}$ and

[TABLE]

The Cauchy problem

[TABLE]

[TABLE]

is considered with $f(\bm{x},t)\in L_{2}(\Omega)$, $u^{0}(\bm{x})\in L_{2}(\Omega)$ using the notation $u(t)=u(\bm{x},t)$.

Let $(\cdot,\cdot),\|\cdot\|$ be the scalar product and norm in $H=L_{2}(\Omega)$, respectively:

[TABLE]

A symmetric positive definite bilinear form $d(u,v)$ such that

[TABLE]

is associated with the Hilbert space $H_{d}$, where the scalar product and norm are, respectively:

[TABLE]

Define $H^{1}_{0}(\Omega)$ as a subspace of $H^{1}(\Omega)$ such that

[TABLE]

Multiplying (3) by $v(\bm{x})\in H^{1}_{0}(\Omega)$ and integrating over the domain $\Omega$, we arrive at the equality

[TABLE]

Here $a(\cdot,\cdot)$ is the following bilinear form:

[TABLE]

where

[TABLE]

In view of (4), we put

[TABLE]

The variational (weak) formulation of the problem (1)–(4) consists in finding $u(\bm{x},t)\in H^{1}_{0}(\Omega)$, $0<t\leq T$ that satisfies (5), (6) with the condition (2) on the boundary.

For the solution of the problem (5), (6), we have the a priori estimate

[TABLE]

To show this, we put $v=u$ in (5) and get

[TABLE]

Taking into account the positive definiteness of the form $a(\cdot,\cdot)$ and the inequality

[TABLE]

we obtain

[TABLE]

From this inequality, in view of (6), it follows that the estimate (7) holds.

It is easy to obtain less trivial and more interesting a priori estimates for the solution of the problem (5), (6). We confine ourselves to the elementary estimate (7) aiming at obtaining similar estimates using various approximations in time and in space for the problem under consideration.

For numerical solving the IBV problem (3), (4), approximation in space is constructed using the finite element method. The weak formulation (5), (6) is employed. Define the subspace of finite elements $V^{h}\subset H_{0}^{1}(\Omega)$ and the discrete elliptic operator $A$ as

[TABLE]

The operator $A$ acts on the finite dimensional space $V^{h}$ and

[TABLE]

where $I$ is the identity operator.

For the problem (3), (4), we put into the correspondence the operator equation for $w(t)\in V^{h}$:

[TABLE]

[TABLE]

where $\varphi(t)=Pf(t)$, $w^{0}=Pu^{0}$ with $P$ denoting $L_{2}$-projection onto $V^{h}$.

Let us multiply equation (9) by $w$ and integrate it over the domain $\Omega$:

[TABLE]

In view of the self-adjointness and positive definiteness of the operator $A$, we have

[TABLE]

The right-hand side can be evaluated by the inequality

[TABLE]

By virtue of this, we have

[TABLE]

The latter inequality leads us to the desired a priori estimate:

[TABLE]

The estimate (11) expresses the stability of the solution of the problem (9), (10) with respect to the initial data and the right-hand side and is quite similar to the estimate (7).

To solve numerically the problem (9), (10), we apply the implicit two-level scheme (Samarskii, 2001). Let $\tau$ be a step of the uniform grid in time so that $y^{n}=y(t^{n}),\ t^{n}=n\tau$, $n=0,1,...,N,\ N\tau=T$. For a constant weight parameter $\sigma\ (0<\sigma\leq 1)$, we define

[TABLE]

To approximate equation (9), let us consider the following two-level scheme:

[TABLE]

[TABLE]

For $\sigma=0$, the scheme (12), (13) becomes the explicit scheme, for $\sigma=1$, we get the fully implicit scheme, whereas $\sigma=0.5$ corresponds to the symmetric scheme (the Crank-Nicolson scheme). The condition

[TABLE]

is necessary and sufficient for the stability of the scheme (Samarskii, 2001; Samarskii et al., 2002). In particular, the following statement is true.

Theorem 1

The difference scheme (12), (13) for $\sigma\geq 0.5$ and (8) is unconditionally stable in $H$ and its solution satisfies the a priori estimate

[TABLE]

Proof 1

Let us estimate the transition operator. Rewrite the scheme with weights (12) in the form

[TABLE]

where $S$ is the operator of transition to a new time level:

[TABLE]

Let us formulate the restrictions on the weight $\sigma$ that guarantee the following two-sided inequality:

[TABLE]

where $\|S\|\leq 1$. Taking into account the commutativity of the operators $B$ and $S$, (16) is equivalent to

[TABLE]

The right inequality is fulfilled for all $\sigma\geq 0$. The left inequality gives

[TABLE]

For non-negative operators $A$, it always holds for $\sigma\geq 0.5$. For (15), in view of (16), we get

[TABLE]

which provides the estimate (14).

3 Two-component scheme based on a partition of unity

For the differential problem under the consideration, we select the domain decomposition

[TABLE]

with overlapping subdomains ($\Omega_{12}\equiv\Omega_{1}\cap\Omega_{2}\neq\varnothing$) (Quarteroni and Valli, 1999; Toselli and Widlund, 2005). To organize parallel computations, each subdomain consists of a set of disconnected subdomains.

To construct schemes of domain decomposition, we use a partition of unity for the computational domain $\Omega$ (Laevsky, 1987; Mathew et al., 1998). Each separate subdomain $\Omega_{\alpha},\ \alpha=1,2$ we associate with the function $\eta_{\alpha}({\bm{x}}),\ \alpha=1,2$ such that

[TABLE]

where

[TABLE]

The standard approach (Samarskii et al., 2002) is based on the additive representation of the operator of the problem (9), (10):

[TABLE]

where each individual operator term $A_{\alpha}$ is associated with the separate subdomain $\Omega_{\alpha},\ \alpha=1,2$. For instance, in view of (1), it is natural to put

[TABLE]

In this case, for the representation (18), we have

[TABLE]

The construction and investigation of domain decomposition schemes for the unsteady problems (9), (10), (18) involves the consideration of the appropriate splitting schemes (Vabishchevich, 2014). In the case of two-component splitting, we can apply the following additive operator–difference schemes of ADI type: the Douglas–Rachford or Peaceman–Rachford scheme and factorized schemes, which generalize them.

In the factorized scheme, an approximate solution at a new time level is evaluated from the equation

[TABLE]

where

[TABLE]

For $\sigma=0.5$, the factorized scheme (19), (20) corresponds to Peaceman–Rachford scheme, whereas $\sigma=1$ results in the Douglas–Rachford scheme.

For the numerical implementation of the factorized scheme, we can introduce the auxiliary value $\widetilde{y}^{n+\sigma}$, which is determined from the equation

[TABLE]

For the approximate solution at a new time level, we have

[TABLE]

Thus, explicit-implicit approximations are used for $Aw$ in (9) taking into account the splitting (18).

Now we provide an estimate for the stability of the factorized scheme, which is similar to the estimate (14) for the scheme with weights (12), (13).

Theorem 2

The factorized difference scheme (13), (18)–(20) is unconditionally stable for $\sigma\geq 0.5$ and its solution satisfies the a priori estimate

[TABLE]

Proof 2

The consideration is conducted by analogy with the proof of theorem 1. From (19), we have

[TABLE]

In the case of (19), (20), the transition operator can be written

[TABLE]

where

[TABLE]

Using Kellogg’s lemma (see, e.g., Kellogg (1964); Grossmann et al. (2007); Vabishchevich (2014)), we obtain

[TABLE]

For $\sigma\geq 0.5$, we have

[TABLE]

In view of this, from (24), we get

[TABLE]

This inequality leads to the estimate (23).

4 Domain decomposition scheme based on indicator functions of subdomains

A new variant of domain decomposition schemes is designed using indicator functions for the subdomains $\Omega_{\alpha},\ \alpha=1,2$. Let us introduce

[TABLE]

Define also the indicator function for the domain of overlap $\Omega_{12}$:

[TABLE]

Thus, we have

[TABLE]

Operators of decomposition are constructed using the indicator functions of the subdomains $\Omega_{\alpha},\ \alpha=1,2$ and $\Omega_{12}$:

[TABLE]

For equality ((25), we have a three-component representation of the problem operator:

[TABLE]

where

[TABLE]

It is necessary to construct splitting schemes for the problem (9), (10), (26), where a transition to a new time level is performed, as before, by solving problems for the operators $A_{\alpha},\ \alpha=1,2$.

Taking into account the splitting (26), similarly to (21), (22), an approximate solution is obtained from equations

[TABLE]

[TABLE]

To study the stability of the scheme (27), (28), we employ theorem 2.

Suppose

[TABLE]

In view of (26), we get

[TABLE]

With the notation introduced above, the scheme (27), (28) is written as

[TABLE]

[TABLE]

where

[TABLE]

Let

[TABLE]

[TABLE]

This allows to write (29), (30) in the form

[TABLE]

[TABLE]

where

[TABLE]

To study the stability of the scheme (31), (32), we use theorem 2. Under the restriction $\sigma\geq 0.5$, the following estimate holds:

[TABLE]

where

[TABLE]

This implies

[TABLE]

where

[TABLE]

The main result of our consideration is the following basic statement on the stability of schemes of domain decomposition constructed on the basis of the indicator functions of the subdomains.

Theorem 3

The splitting scheme (13), (26)–(28) is unconditionally stable for $\sigma\geq 0.5$ and its solution satisfies the a priori estimate (33).

5 Numerical example

Numerical experiments presented here are of comparative nature. We consider the two-level scheme with weights (12), (13) as the reference scheme that provides the benchmark numerical solution for a comparison with the results of two schemes of domain decomposition with overlapping subdomains. Namely, these are the scheme based on the partition of unity for the domain (13), (18)–(20) and the scheme that uses the indicator functions of subdomains (13), (26)–(28). In our consideration, we monitor the proximity of the above decomposition schemes to the reference scheme. We confine ourselves to the case of schemes with $\sigma=1$.

We consider the model problem (1)–(4) in the unit square ($\Omega=[0,1]\times[0,1]$), where

[TABLE]

and $T=0.1$.

Let $\bar{y}$ be the benchmark numerical solution obtained using the implicit scheme (12), (13), whereas $y_{1},y_{2}$ are the solutions from the decomposition scheme (13), (18)–(20) and (13), (26)–(28), respectively. We evaluate the deviation of the solution obtained using the domain decomposition methods from the benchmark solution derived without domain decomposition as the error:

[TABLE]

The computational domain is divided by the variable $x_{1}$ into two subdomains with the overlap width $2\delta$. For the basic variant $\delta=0.05$, the functions $\eta_{\alpha}(x_{1}),\ \chi_{\alpha}(x_{1}),\ \alpha=1,2$, which define the partition operators, are shown in Fig.1. The uniform spatial grid $51\times 51$ was used with piecewise-linear finite elements on triangles. The number of steps in time is $N=50$.

The error of the domain decomposition schemes is presented in Fig.2. In the example considered here, the accuracy of the scheme based on the indicator functions of the subdomains (13), (26)–(28) is higher. The benchmark solution $\bar{y}$ at the final time moment is shown in Fig.3. The deviation of the solutions obtained with two schemes of domain decomposition from the benchmark one is given in Fig.4, 5.

Reducing the overlap width by half (see Fig.6) results in decreasing the accuracy of the approximate solution. For $\delta=0$, both schemes under the consideration coincide and lead to the domain decomposition scheme with non-overlapping subdomains.

The use of a finer grid in space (compare Fig.2 and Fig.7) leads to a drop in accuracy. This fact demonstrates a conditional convergence of the domain decomposition schemes. Calculations with a smaller time step (see Fig.8) allow to increase the accuracy of the numerical solution.

6 Conclusions

The boundary value problem for a second-order parabolic equation is considered. The standard finite element approximation in space is employed. For approximation in time, the two-level scheme with weight ($\theta$-method) is used.

Iteration-free schemes of domain decomposition with overlapping subdomains are constructed on the basis of a partition of unity for a computational domain. In the case of a two-component domain decomposition, approximation in time is designed using factorized schemes, which generalize the standard ADI schemes.

A new class of domain decomposition schemes with overlapping subdomains based on indicator functions of subdomains is proposed for unsteady problems. In this case, for the problem operator, we have a three-component representation with decomposition operators connected with subdomains and their intersection. The unconditional stability of these schemes is established for two-component schemes of domain decomposition.

Numerical results for the model two-dimensional problem demonstrate the robustness of the new decomposition scheme with overlapping subdomains. A comparison of the conventional scheme based on a partition of unity with the scheme that used the indicator functions of the subdomains demonstrates its advantage in accuracy.

Acknowledgements

This work was supported by the grant of Russian Federation Government (agreement # 14.Y26.31.0013).