Lobachevskii Journal of Mathematics Vol. 8, 2001, 167–184 \begin{matrix}c \\ ◯ \end{matrix}E. Laitinen, A.V. Lapin, J.Pieskä

E. Laitinen, A.V. Lapin, J.Pieskä

LARGE SPLITTING ITERATIVE METHODS AND PARALLEL SOLUTION OF VARIATIONAL INEQUALITIES

Abstract. Splitting iterative methods for the sum of maximal monotone and single-valued monotone operators in a finite-dimensional space are studied: convergence, rate of convergence and optimal iterative parameters are derived. A two-stage iterative method with inner iterations is analysed in the case when both operators are linear, self-adjoint and positive definite. The results are applied for the mesh variational inequalities which are solved using a non-overlapping domain decomposition method and the splitting iterative procedure. Parallel solution of a mesh scheme for continuous casting problem is presented and the dependence of the calculation time on the number of processors is discussed.

1. Introduction

We solve the inclusion

Au + Cu ∋ F  (1)

in a finite-dimensional space V equipped by the inner product (., .) and the norm ||.|| = (., .)^{1/2}. The operators A, C satisfy the following assumptions:

A is bounded, hemicontinuous and uniformly monotone in V,  (2)

C : V → 2^{V}  is maximal monotone.  (3)

We recall, that A is hemicontinuous if the function t → (A(u + tv), w) is continuous for any u, v, w ∈ V and by the uniform monotonicity we mean that there exists a positive constant m such that (Au_{1} − Au_{2}, u_{1} − u_{2}) ≥ m||u_{1} − u_{2}||^{2}, ∀u_{1}, u_{2} ∈ V. Further, domain of C : V → 2^{V} is the subset domC ⊂ V of those u that C(u) ≠ ∅. The operator C is maximal monotone if 1) it is monotone, i.e. for each pair u_{1}, u_{2} ∈ domC we have

(γ_{1} − γ_{2}, u_{1} − u_{2}) ≥ 0 ∀γ_{1} ∈ C(u_{1}) and ∀γ_{2} ∈ C(u_{2}),

2) the condition

(γ_{1} − γ_{2}, u_{1} − u_{2}) ≥ 0 ∀u_{1} ∈ domC and ∀γ_{1} ∈ C(u_{1})

implies u_{2} ∈ domC and γ_{2} ∈ C(u_{2}). Properties (2),(3) ensure the existence of a unique solution u to equation (1) for any F ∈ V [1, 8].

Let further B be one more bounded, hemicontinuous and uniformly monotone operator and D be a linear self-adjoint and positive definite operator. We study the following iterative method for solving (1):

D^{−1}(Bu^{k+1} − Bu^{k}) + Au^{k} + C(Bu^{k+1} − Bu^{k} + u^{k}) ∋ F, k = 0, 1, 2, ... (4)

starting from any u^{0} ∈ V.

To implement this iterative method we can split (4) into the system of noncoupled inclusion and equation:

D^{−1}(u^{k+1/2} − u^{k}) + Au^{k} + Cu^{k+1/2} ∋ F,
Bu^{k+1} − Bu^{k} = u^{k+1/2} − u^{k}.

Due to the properties of D and B the iterations u^{k+1/2} and u^{k+1} are correctly defined for any k.

The reason for such kind of splitting is that we choose the operator D^{−1} in such a way that the solution procedure for the inclusion in this system is much easier than for problem (1). We can mention several particular cases of iterative method (4). Namely, if B = I, where by I we denote the identity mapping, then we derive preconditioned implicit iterative scheme for (1):

u^{0} ∈ V; D^{−1}(u^{k+1} − u^{k}) + Au^{k} + Cu^{k+1} ∋ F, k = 0, 1, 2, ...  (5)

If D = λI, B = I + λA, then method (4) can be written in the splitting form:

(u^{k+1/2} − u^{k})/λ + Au^{k} + Cu^{k+1/2} ∋ F,  (6)
(u^{k+1} − u^{k+1/2})/λ + Au^{k+1} − Au^{k} = 0,  (7)

which is the Douglas - Rachford scheme for problem (1). In the case D = λI and B = 1/2(I + λA) method (4) becomes the Peaceman- Rachford scheme. For other examples of splitting methods for variational inequalities with linear operator A see [7].

Solution of (5) for fixed k may be of the same complexity as for problem (1), at least in the case of linear operator A. On the other hand, to solve the Douglas-Rachford scheme (6),(7) we need to solve equation with singlevalued operator I +λA (7) and the equation with multivalued operator I +λC (6). As V is finite dimensional space then we can use well-known iterative procedures for the solution of system (7) or a direct method if A is linear operator. Operator I + λC is easy invertible in many practical cases (see examples below).

The splitting algorithms for finding a zero of the sum of two maximal monotone operators have been studied in [5, 6]. Some estimates for the rate of convergence for these methods have been obtained in the case when A is single-valued monotone and Lipschitz-continuous operator. More sharp estimates have been derived in [7] for Douglas - Rachford and Peaceman-Rachford schemes as well as for some other splitting algorithms (alternating direction algorithms) in the case of linear self-adjoint operator A.

Method (5) with C = ∂I_{K} ( ∂I_{K} is the subdifferential of the indicator function for a convex set K) and with a positive definite, not self-adjoint, operator D has been studied in [10]. It has been shown that this method can be treated as a partial case of the splitting method for sum of two maximal monotone operators and its convergence has been proved. We also refer to [10] for bibliography of the papers concerning the splitting methods for variational inequalities.

The main purpose of this paper is to obtain the rate of convergence for iterative method (4) with block-diagonal operators. In applications for mesh variational inequalities this corresponds to joined using the non-overlapping domain decomposition methods and splitting iterative algorithms.

The rest of the article is organised as follows. In section 2 we derive a priori estimates for method (4) with nonlinear operator A. Then we prove some results on the rate of convergence and the optimal iterative parameters for the case of block-diagonal operators. A two-stage splitting iterative method with inner iterations to solve (7) in the case of linear operators A and B is investigated. One example of splitting iterative method for a finite element scheme constructed via domain decomposition technique for a variational inequality with nonlinear elliptic operator is studied in section 3. In section 4 we apply the theoretical results for solving a finite difference scheme for a continuous casting problem. Choice of the iterative parameters, implementation of the methods and their parallelisation are discussed.

2. A priori estimates for the splitting iterative method

We denote by z^{k} = u^{k} − u, where u^{k} is a k-th iteration of method (4) while u is the exact solution for problem (1). To shorten writing we also use the notations:

Au^{k} − Au ≡ A^{(k)}z^{k}, Bu^{k} − Bu ≡ B^{(k)}z^{k}.

Lemma 1. Let assumptions (2),(3) be fulfilled. Then the following error estimate is valid:

||D^{−1/2}B^{(k+1)}z^{k+1}|| ≤ 1/2||D^{−1/2}(2B^{(k)}z^{k} − z^{k} − DA^{(k)}z^{k})||+
1/2||D^{−1/2}(z^{k} − DA^{(k)}z^{k})||.  (8)

Proof. From (1),(4) we derive the inclusion:

D^{−1}(B^{(k+1)}z^{k+1} − B^{(k)}z^{k} + z^{k}) + C(Bu^{k+1} − Bu^{k} + u^{k}) − Cu ∋
D^{−1}z^{k} − A^{(k)}z^{k}.

Multiplying it by B^{(k+1)}z^{k+1} − B^{(k)}z^{k} + z^{k} and using the monotonicity of the operator C we have

||D^{−1/2}(B^{(k+1)}z^{k+1} − B^{(k)}z^{k} + z^{k})||^{2} ≤
(D^{−1}(B^{(k+1)}z^{k+1} − B^{(k)}z^{k} + z^{k}), z^{k} − DA^{(k)}z^{k}).

Obvious transformations of this inequality leads to

||D^{−1/2}B^{(k+1)}z^{k+1}||^{2} − (D^{−1}B^{(k+1)}z^{k+1}, 2B^{(k)}z^{k} − z^{k} − DA^{(k)}z^{k})+

(D^{−1}B^{(k)}z^{k}, B^{(k)}z^{k} − z^{k} − DA^{(k)}z^{k}) + (A^{(k)}z^{k}, z^{k}) ≤ 0.  (9)

As a consequence we get

||D^{−1/2}B^{(k+1)}z^{k+1}||^{2} − ||D^{−1/2}B^{(k+1)}z^{k+1}||||D^{−1/2}(2B^{(k)}z^{k} − z^{k}−
DA^{(k)}z^{k})|| + (D^{−1}B^{(k)}z^{k}, B^{(k)}z^{k} − z^{k} − DA^{(k)}z^{k}) + (A^{(k)}z^{k}, z^{k}) ≤ 0.

Because there is a quadratic trinomial of the ||D^{−1/2}B^{(k+1)}z^{k+1}|| in the lefthand side of the last inequality, then estimate (8) follows from it immediately.

Theorem 1. Let V = V_{1} × V_{2} × ... × V_{p}, where V_{i} are the finite-dimensional spaces with the inner products (., .)_{i} and the norms ||.||_{i} = (., .)_{i}^{1/2} and let A be the diagonal operator: A = diag(A_{1}, A_{2}, ..., A_{p}) with A_{i} : V_{i} → V_{i} satisfying for all i the following assumptions:

(A_{i}u_{i} − A_{i}v_{i}, u_{i} − v_{i})_{i} ≥ m_{i}||u_{i} − v_{i}||_{i}^{2}, m_{i} = const > 0,  (10)
||A_{i}u_{i} − A_{i}v_{i}||_{i}^{2} ≤ M_{i}(A_{i}u_{i} − A_{i}v_{i}, u_{i} − v_{i})_{i}, M_{i} = const > 0.  (11)

Let further D = diag(λ_{1}I_{1}, λ_{2}I_{2}, ..., λ_{p}I_{p}), I_{i} be the identity mappings in V_{i} and

either B = B_{1} ≡ I + DA or B = B_{2} ≡ 1/2(I + DA).

Then iterative method (4) converges for any λ_{i} > 0 and for the set of optimal parameters λ_{i} = 1/\sqrt{m_{i}M_{i}} the following estimate for the rate of convergence holds:

∥D^{−1/2}(I + DA^{(n)})z^{n}∥ ≤ q^{n/2}∥D^{−1/2}(I + DA^{(0)})z^{0}∥  (12)

with q = q_{1} = max_{1≤i≤p} \frac{1}{1 + \sqrt{η}_{i}} for B = B_{1} and with q = q_{2} = max_{1≤i≤p} \frac{1 − \sqrt{η_{i}}}{1 + \sqrt{η_{i}}} for

B = B_{2}, η_{i} = m_{i}/M_{i}.

Proof. Let B = 1/2(I +DA). Then it is easy to see that estimate (8) becomes:

∑_{i=1}^{p} λ_{i}^{−1}||z_{i}^{k+1} + λ_{i}A_{i}^{(k+1)}z^{k+1}||_{i}^{2} ≤ ∑_{i=1}^{p} λ_{i}^{−1}||z_{i}^{k} − λ_{i}A_{i}^{(k)}z^{k}||_{i}^{2}.

Using properties (10),(11) we get

||z_{i}^{k} + λ_{i}A_{i}^{(k)}z_{i}^{k}||_{i}^{2} ≤ (1/m_{i} + 2λ_{i} + M_{i}λ_{i}^{2})(A_{i}^{(k)}z_{i}^{k}, z_{i}^{k})_{i},

so,

∑_{i=1}^{p} λ_{i}^{−1}||z_{i}^{k+1} + λ_{i}A_{i}^{(k+1)}z_{i}^{k+1}||_{i}^{2} ≤ ∑_{i=1}^{p} λ_{i}^{−1}||z_{i}^{k} + λ_{i}A_{i}^{(k)}z_{i}^{k}||_{i}^{2} − 4 ∑_{i=1}^{p}(A_{i}^{(k)}z_{i}^{k}, z_{i}^{k})_{i}

≤ ∑_{i=1}^{p} \frac{(1/m_{i} − 2λ_{i} + M_{i}λ_{i}^{2})}{(1/m_{i} + 2λ_{i} + M_{i}λ_{i}^{2})}λ_{i}^{−1}||z_{i}^{k} + λ_{i}A_{i}^{(k)}z_{i}^{k}||_{i}^{2}.

The convergence of the iterative method and desired estimate are the straight consequences of the last inequality.

The proof for the case of B = I + DA is similar.

Remark 1. Let the diagonal operator A has linear Gateaux derivative A′(y) for all y ∈ V, A′(y) = diag(A_{1}(y), A_{2}(y), ..., A_{p}(y)), and

A_{i}(y) ∈ L(V_{i}; V_{i}), m_{i}I_{i} ≤ A_{i}(y) = A_{i}^{∗}(y) ≤ M_{i}I_{i} ∀i.

Then assumptions (10) and (11) are fulfilled.

Remark 2. Let B = I, D = R^{−1} with R = R^{∗} > 0 and

(Au − Av, u − v) ≥ m(R(u − v), u − v), m = const > 0,  (13)
||R^{−1/2}(Au − Av)||^{2} ≤ M(R(u − v), u − v), M = const > 0 (14)

for all u, v ∈ V. Then iterative method (5) converges for λ ∈ (0, 2m/M) and for λ = λ_{0} = m/M the following estimate is valid:

||R^{1/2}z^{n}|| ≤ q^{n/2}||R^{1/2}z^{0}|| ∀n, q = 1 − m^{2}/M.  (15)

In fact, from estimate (8) with B = I, D = λR^{−1} and using properties (13),(14) we deduce:

||R^{1/2}z^{k+1}||^{2} ≤ ||R^{1/2}z^{k} − λR^{−1/2}(Au^{k} − Au)||^{2} =
||R^{1/2}z^{k}||^{2} − 2λ(Au^{k} − Au, z^{k}) + λ^{2}||R^{−1/2}(Au^{k} − Au)||^{2} ≤
(1 − 2mλ + Mλ^{2})||R^{1/2}z^{k}||^{2},

which obviously implies the cited results.

Now let us study the convergence and derive the rate of convergence for two-stage iterative method with inner iterations (cf.[4]) resulting from an inexact inversion of the operator B in method (4). We analyse only the case of the linear operators A and B.

To do this we need one more error estimate for method (4), given in Lemma 2.

Lemma 2. Let A, B, and D be linear bounded, self-adjoint operators and let also they all commute. Then for the error z^{k} = u^{k} − u of iterative method (4) the following estimate is true:

(Sz^{k+1}, z^{k+1}) ≤ ||I − 2AS^{−1}|| · (Sz^{k}, z^{k}), S = B(D^{−1} + A). (16)

Proof. Let us further use for shortening the following notations: (u, v)_{d} = (D^{−1}u, v), Q = B − I − DA. Note, that Q commute with all operators A,B and D. We proceed as in Lemma 1 to obtain inequality (9) which can be rewritten as

||Bz^{k+1}||_{d}^{2} − (Bz^{k+1}, (B + Q)z^{k})_{d} + (Bz^{k}, Qz^{k})_{d} + (DAz^{k}, z^{k})_{d} ≤ 0

From this inequality we deduce:

(Bz^{k+1}, (B − Q)z^{k+1})_{d} − (Bz^{k}, (B − Q)z^{k})_{d} + 2(DAz^{k}, z^{k})_{d} ≤
−(Bz^{k+1}, (B + Q)z^{k+1})_{d} − (Bz^{k}, (B + Q)z^{k})_{d} + 2(Bz^{k+1}, (B + Q)z^{k})_{d}
= −(B(B + Q)(z^{k+1} − z^{k}), z^{k+1} − z^{k})_{d}.

Due to the positive definiteness of the self-adjoint operator B(B + Q) the right-hand side of this inequality is non-positive, so,

(Sz^{k+1}, z^{k+1}) ≤ (Sz^{k}, z^{k}) − 2(Az^{k}, z^{k}).

The last inequality leads to (16).

We consider iterative method (4) with the linear operators A and B and use an inner iterative method to solve the equation

Bv = F^{k}(= u^{k+1/2} − u^{k})  (17)

Let an m-th inner iteration v^{m} satisfies the equality

v^{m} = T_{m}v^{0} + Q_{m}F^{k}.  (18)

Obviously, the linear operators T_{m}, Q_{m} must be such that

v = T_{m}v + Q_{m}F^{k}

for the exact solution v of (17).

Choosing v^{0} = 0 we find

v^{m} = (I − T_{m})v.

It appears that the iterative method (4) with m inner iterations (18) for solving the equation (17) is equivalent to the iterative method (4) with operator B̃ = B(I − T_{m})^{−1} instead of B.

Theorem 2. Let V = V_{1} × V_{2} × ... × V_{p} and

\{ \begin{matrix}A = diag(A_{1}, A_{2}, ..., A_{p}), A_{i} ∈ L(V_{i}; V_{i}), \\ m_{i}I_{i} ≤ A_{i} = A_{i}^{∗} ≤ M_{i}I_{i} ∀i, \end{matrix}  (19)

D = diag(λ_{1}I_{1}, λ_{2}I_{2}, ..., λ_{p}I_{p}).  (20)

Let iterative method (4) with B = I + DA or B = 1/2(I + DA) and m inner iterations (18) for solving equation (17) be used. Let also the following assumptions be satisfied:

A = A^{∗} > 0, T_{m} = T_{m}^{∗}, AT_{m} = T_{m}A, ||T_{m}|| ≤ q < 1.  (21)

Then this two-stage iterative method converges for any λ_{i} > 0, i = 1, 2, ..., p and for λ_{i} = 1/\sqrt{m_{i}M_{i}} the following estimate holds:

\{ \begin{matrix}(Sz^{k+1}, z^{k+1}) ≤ ρ(Sz^{k}, z^{k}), \\ ρ = max_{1≤i≤p}(1 − 2(1 − q)\frac{η_{i}^{1/2}}{(1 + η_{i}^{1/2})^{2}} ), η_{i} = m_{i}/M_{i}. \end{matrix}  (22)

Proof. All the assumptions of Lemma 2 being fulfilled we can use estimate (16) which becomes:

(Sz^{k+1}, z^{k+1}) ≤ ||R|| · (Sz^{k}, z^{k})

with S = D^{−1}(I + DA)^{2}(I − T_{m})^{−1}, R = I − 2AD(I + DA)^{−2}(I − T_{m}).

The operator R is bounded and self-adjoint. Let r_{j} be its eigenvalues while a_{j}, d_{j} and t_{j} be the eigenvalues of A, D and T_{m}, respectively. Then r_{j} =

\frac{1 + 2t_{j}a_{j}d_{j} + (d_{j}a_{j})^{2}}{1 + 2a_{j}d_{j} + (d_{j}a_{j})^{2}} and r_{j} < 1, because t_{j} ≤ q < 1. It proves the convergence

of the two-stage iterative method for any choice of parameters λ_{i} > 0, i = 1, 2, ..., p.

To prove inequality (22) we fix an index j of eigenvalues. Let for this index they are the eigenvalues of i-th blocks of matrices, so, in particularly, d_{j} = λ_{i}. When choosing λ_{i} = 1/\sqrt{m_{i}M_{i}} we obtain for r_{j} the estimate

r_{j} = \frac{1 + 2t_{j}a_{j}λ_{i} + (λ_{i}a_{j})^{2}}{1 + 2a_{j}λ_{i} + (λ_{i}a_{j})^{2}} ≤ max_{m_{i}≤a≤M_{i}} \frac{1 + 2qaλ_{i} + (λ_{i}a)^{2}}{1 + 2aλ_{i} + (λ_{i}a)^{2}} =

(1 − 2(1 − q)\frac{η_{i}^{1/2}}{(1 + η_{i}^{1/2})^{2}} ).

This estimate implies (22).

Remark 3. To ensure assumptions (21) we can choose the conjugate gradient method or Chebyshev iterations as an inner iterative method.

Remark 4. In the case of problem (1) with linear non-selfadjoint operators or/and with noncommuting operators T_{m} and A we can only prove that twostage iterative method converges if the norm of operator T_{m} is sufficiently small, i.e. if we solve equation (17) with high accuracy. We prove this result supposing for simplicity that D = λI and using estimate (8) which is now

||Bz^{k+1}|| ≤ (1/2||(I − λA)B^{−1}|| + ||I − 1/2(I + λA)B^{−1}||)||Bz^{k}||. (23)

Let B = 1/2(I + λA) in iterative method (4) and we execute m inner iterations (18) for solving equation (17). Then the two-stage iterative method is equivalent to (4) with B = 1/2(I + λA)(I − T_{m})^{−1} and (23) reads as follows:

||Bz^{k+1}|| ≤ (||(I − λA)(I − T_{m})(I + λA)^{−1}||+
||I − (I + λA)(I − T_{m})(I + λA)^{−1}||)||Bz^{k}||.  (24)

Let us denote by S = (I − λA)(I + λA)^{−1} and by G_{m} = (I + λA)T_{m}(I + λA)^{−1}. Then the norm of operator in right-hand side of (24) is majorated by ρ = ||S||||I − G_{m}|| + ||G_{m}||. If we suppose that

||G_{m}|| < \frac{1 − ||S||}{1 + ||S||},  (25)

then ρ < 1 and the two-stage iterative method converges.

As ||G_{m}|| ≤ ||T_{m}||||I + λA||||(I + λA)^{−1}|| ≤ q · cond(I + λA) and ||S|| < 1 ∀λ > 0, then we get inequality (25) for sufficiently small q, in other words, after sufficiently many inner iterations.

3. Domain decomposition and splitting iterative method for a

finite element approximation of elliptic variational inequality

Let Ω be a polygonal domain in ℝ^{2} with the boundary ∂Ω, the functions k_{i}(x, t) ≡ k_{i}(x, t_{0}, t_{1}, t_{2}) be continuous and satisfy for all t, q, s ∈ ℝ^{3} and for all x ∈ \overline{Ω} the following assumptions:

∑_{i=0}^{2}(k_{i}(x, t) − k_{i}(x, q))(t_{i} − q_{i}) ≥ c_{0}||t − q||^{2},

∑_{i=0}^{2}(k_{i}(x, t) − k_{i}(x, q))s_{i} ≤ c_{1}( ∑_{i=0}^{2}(k_{i}(x, t) − k_{i}(x, q))(t_{i} − q_{i}))^{1/2}||s||.

Let also f(x) ∈ L_{2}(Ω) and ϕ(x, t) be continuous in x, t and convex in t function, |ϕ(x, t)| ≤ c · |t|^{r}, r < ∞, ∀x ∈ Ω̄, t ∈ ℝ^{3}.

We denote by V = {u ∈ H^{1}(Ω) : u(x) = 0 for a.a. x ∈ ∂_{1}Ω}, where ∂_{1}Ω is a part (possibly empty) of ∂Ω.

Let further

a(u, v) = ∫_{Ω}(∑_{i=1}^{2} k_{i}(x, u, ∇u)∂v/∂x_{i} + k_{0}(x, u, ∇u)v)dx,

f(v) = ∫_{Ω} f(x)v(x)dx, θ(u) = ∫_{Ω} ϕ(x, u(x))dx.

We consider the variational inequality

u ∈ V : a(u, v − u) + θ(v) − θ(u) ≥ f(v − u) ∀v ∈ V.  (26)

Let T_{h} be a regular triangulation of Ω̄ and diam δ = O(h) for any triangle δ ∈ T_{h}. Let Ω be divided into p non-overlapping subdomains Ω_{i} with the interfaces S_{ij} = \overline{Ω}_{i} ∩ \overline{Ω}_{j}. We suppose that all interfaces as well as \overline{∂_{1}Ω} consist of the sides of δ ∈ T_{h}.

We define several finite-dimensional linear spaces, namely, V_{h}^{0} = {u_{h} ∈ V ∩ C(Ω̄) : u_{h} ∈ P_{1} ∀δ ∈ T_{h}}, where P_{1} is the space of linear functions, and V_{h}^{i}, i = 1, 2, ..., p as the spaces of restrictions of all functions from V_{h}^{0} on the subdomains Ω_{i}. We also denote by V_{h} = V_{h}^{1} × V_{h}^{2} ⋯ V_{h}^{p}. It is easy to check that V_{h}^{0} is isomorphic to the subspace K_{h} of V_{h}, where K_{h} = {u_{h} = (u_{h}^{1}, u_{h}^{2}, ..., u_{h}^{p}) ∈ V_{h} : u_{h}^{i}(x) = u_{h}^{j}(x) for x ∈ S_{ij}, i, j = 1, 2, ..., p}.

Let Π_{h} be the local P_{1}-interpolant on the element δ ∈ T_{h} and for any continuous function v(x) and for any δ ∈ T_{h}

S_{δ}(v) = ∫_{δ} Π_{h}(v)dx.

To approximate variational inequality (26) by a mesh scheme we define the following forms and functionals:

a_{h}^{i}(u_{h}^{i}, v_{h}^{i}) = ∫_{Ω_{i}}(∑_{l=1}^{2} k_{l}(x, u_{h}^{i}, ∇u_{h}^{i})∂v_{h}^{i}/∂x_{l} + k_{0}(x, u_{h}^{i}, ∇u_{h}^{i})v_{h}^{i})dx,

f_{h}^{i}(v_{h}^{i}) = ∫_{Ω_{i}} f(x)v_{h}^{i}(x)dx, θ_{h}^{i}(u_{h}^{i}) = ∑_{δ⊂Ω_{i}} ∫_{δ_{h}} Π_{h}(ϕ(x, u_{h}^{i}(x)))dx.

Now by the mesh approximation of (26) we mean the problem:

\{ \begin{matrix}find u_{h} ∈ K_{h} such that ∀v_{h} ∈ K_{h} \\ ∑_{i=1}^{p} a_{h}^{i}(u_{h}^{i}, v_{h}^{i} − u_{h}^{i}) + ∑_{i=1}^{p} θ_{h}^{i}(v_{h}^{i}) − ∑_{i=1}^{p} θ_{h}^{i}(u_{h}^{i}) ≥ ∑_{i=1}^{p} f_{h}^{i}(v_{h}^{i} − u_{h}^{i}). \end{matrix} (27)

To rewrite this mesh scheme as an inclusion (1) we proceed as follows.

Let us put in the correspondence to the function u_{h}^{i} ∈ V_{h}^{i} the vector u^{i} ∈ ℝ^{N_{i}} of its nodal values for nodes from Ω̄_{i}\\overline{∂_{1}Ω} and denote this bijection by u^{i} ⇔ u_{h}^{i}. To u_{h} ∈ V_{h} corresponds the vector u ∈ ℝ^{N}, N = N_{1} + N_{2} + ⋯ + N_{p}.

Let the linear space of vectors u^{i} ∈ ℝ^{N_{i}} be equipped by the inner product

(u^{i}, v^{i})_{i} = ∫_{Ω_{i}} u_{h}^{i}(x)v_{h}^{i}(x)dx for u^{i} ⇔ u_{h}^{i}, v^{i} ⇔ v_{h}^{i}.

We denote this space by W^{i} and set W = W^{1} × W^{2} × ⋯ × W^{p}. Then the subspace K_{h} of V_{h} corresponds to a subspace of W, which we denote by K.

We define the nonlinear operators A_{i} : W^{i} → W^{i}, the functional θ^{i} : W^{i} → ℝ^{1} and the vectors F_{i} ∈ W^{i} as follows:

(A_{i}u^{i}, v^{i})_{i} = a_{h}^{i}(u_{h}^{i}, v_{h}^{i}), θ^{i}(u^{i}) = θ_{h}^{i}(u_{h}^{i}), (F^{i}, v^{i})_{i} = f_{h}^{i}(v_{h}^{i})
for u^{i} ⇔ u_{h}^{i}, v^{i} ⇔ v_{h}^{i}.

Note, that for every i = 1, …, p the function θ^{i} is a separable one, i.e. θ^{i}(u^{i}) =

∑_{l=1}^{N_{i}} θ_{l}^{i}(u_{l}^{i}) with the scalar convex functions θ_{l}^{i}.

Let further A = diag(A_{1}, A_{2}, ..., A_{p}) : W → W, and θ(u) = ∑_{i=1}^{p} θ^{i}(u^{i}) : W → ℝ^{1}, and F = (F^{1}, F^{2}, ..., F^{p}) ∈ W.

Below we denote by C = ∂θ(u) + ∂I_{K}, where I_{K} is the indicator function of the subspace K. As dom(∂θ) = ℝ^{N} then C is maximal monotone [8].

With these notations mesh scheme (27) becomes the partial case of problem (1).

Let now matrices R_{i}, i = 1, 2, ..., p be defined by the equalities

(R_{i}u^{i}, v^{i})_{i} = ∫_{Ω_{i}}(u_{h}^{i}(x)v_{h}^{i}(x) + ∇u_{h}^{i}(x)∇v_{h}^{i}(x))dx

and R = diag(R_{1}, R_{2}, ..., R_{p}). Then the assumptions on the ”coefficients” k_{j}(x, t) ensure the estimates:

(A_{i}u^{i} − A_{i}v^{i}, u^{i} − v^{i})_{i} ≥ c_{0}(R_{i}(u^{i} − v^{i}), u^{i} − v^{i})_{i},
(A_{i}u^{i} − A_{i}v^{i}, w^{i})_{i} ≤ c_{1}(A_{i}u^{i} − A_{i}v^{i}, u^{i} − v^{i})_{i}^{1/2}(R_{i}w^{i}, w^{i})_{i}^{1/2} ∀i.

As a consequence of the second inequality we have:

||A_{i}u^{i} − A_{i}v^{i}||_{i} ≤ c_{1}||R_{i}||_{i}^{1/2} · (A_{i}u^{i} − A_{i}v^{i}, u^{i} − v^{i})_{i}^{1/2} ∀i.

Now,

α_{i}(u^{i}, u^{i})_{i} = c_{i}h^{−2} ∫_{Ω_{i}}(u_{h}^{i})^{2}dx ≥ (R_{i}u^{i}, u^{i})_{i} ≥ ∫_{Ω_{i}}(u_{h}^{i})^{2}dx = (u^{i}, u^{i})_{i}

with α_{i} = O(h^{−2}). Summarising the received results we derive:

(A_{i}u^{i} − A_{i}v^{i}, u^{i} − v^{i})_{i} ≥ c_{0}||u^{i} − v^{i}||_{i}^{2},
||A_{i}u^{i} − A_{i}v^{i}||_{i}^{2} ≤ c_{1}^{2}α_{i}(A_{i}u^{i} − A_{i}v^{i}, u^{i} − v^{i})_{i}

for all i = 1, 2, ..., p. Thus, the assumptions of Theorem 1 are satisfied with m_{i} = O(1), M_{i} = O(h^{−2}). If we choose λ_{i} = O(h) in method (4) with either B = I + DA or B = 1/2(I + DA), D = diag(λ_{1}I_{1}, λ_{2}I_{2}, ..., λ_{p}I_{p}), where I_{i} is identity in V^{i}, then q_{1} = 1 − O(h), q_{2} = 1 − O(h) and the number of iterations to achieve accuracy ϵ is n(ϵ) = O(h^{−1} ln 1/ϵ).

Let us describe a possible implementation of iterative method (4) with B = I + DA, D = diag(λ_{1}I_{1}, λ_{2}I_{2}, ..., λ_{p}I_{p}), i.e. the method

D^{−1}(u^{k+1/2} − u^{k}) + Au^{k} + Cu^{k+1/2} ∋ F  (28)
(I_{i} + λ_{i}A_{i})(u^{i,k+1} − u^{i,k}) = u^{i,k+1/2} − u^{i,k}, i = 1, 2, ..., p,  (29)

where u^{k} = (u^{1,k}, u^{2,k}, ..., u^{p,k}).

Nonlinear equations (29) can be solved independently for i = 1, 2, ..., p.

Further, for the coordinates of u^{k+1/2}, corresponding to the internal nodes x ∈ Ω_{i}, i = 1, …, p, the operator C has the diagonal form and we get for these coordinates a system of non-coupled scalar (one-dimensional) inclusions. More precisely, let l be an index such that u_{l}^{i,k+1/2} ⇔ u_{h}^{k+1/2}(x) for a fixed node x ∈ Ω_{i}. Then u_{l}^{k+1/2} = u_{l}^{i,k+1/2} and (Cu^{k+1/2})_{l} = ∂θ_{l}^{i}(u_{l}^{i,k+1/2}). The corresponding to x ∈ Ω_{i} inclusion from system (28) becomes

\frac{1}{λ_{i}} u_{l}^{i,k+1/2} + ∂θ_{l}^{i}(u_{l}^{i,k+1/2}) ∋ ϕ_{l} ≡ \frac{1}{λ_{i}} u_{l}^{i,k} + (A_{i}u^{i,k})_{l}  (30)

Inclusion (30) is equivalent to the problem of minimisation of the scalar convex function

\frac{1}{2λ_{i}} (u_{l}^{i,k+1/2})^{2} + θ_{l}^{i}(u_{l}^{i,k+1/2}) − ϕ_{l}u_{l}^{i,k+1/2}

and can be solved by using a standard optimisation method. Moreover, in many applications (see below, as an example, a mesh scheme for continuous casting problem) ∂θ_{l}^{i} appears as a piecewise linear graph and inclusion (30) can be solved by the direct calculations.

For the nodes lying on the interfaces S_{ij} system (28) contains subsystems of two (if it is the interior node of the interface) or several (if it is a crosspoint of several interfaces) coupled equations. These subsystems can be also reformulated as the problems to minimise convex functions of one variable.

In fact, let for the simplicity a node x ∈ S_{ij} be an interior one. Then two coordinates of u^{k+1/2} correspond to it, let u_{l}^{i,k+1/2} and u_{s}^{j,k+1/2}. The subsystem of (28) for these two unknowns has the form

( \begin{matrix}\frac{1}{λ_{i}} & 0 \\ 0 & \frac{1}{λ_{j}} \end{matrix} ) ( \begin{matrix}u_{l}^{i,k+1/2} \\ u_{s}^{j,k+1/2} \end{matrix} ) + ( \begin{matrix}∂θ_{l}^{i} & 0 \\ 0 & ∂θ_{s}^{j} \end{matrix} ) ( \begin{matrix}u_{l}^{i,k+1/2} \\ u_{s}^{j,k+1/2} \end{matrix} ) +

∂I_{l,s}(u_{l}^{i,k+1/2}, u_{s}^{j,k+1/2}) = ( \begin{matrix}ϕ_{l,s}^{1} \\ ϕ_{l,s}^{2} \end{matrix} )  (31)

with known right-hand side ϕ_{l,s} and the subdifferential ∂I_{l,s} of the indicator function for the set K_{l,s} = {u = (u_{l}, u_{s}) : u_{l} = u_{s}}. Inclusion (31) can be rewritten as the minimisation problem for the convex function of two variables

\frac{1}{2λ_{i}} (u_{l}^{i,k+1/2})^{2}+θ_{l}^{i}(u_{l}^{i,k+1/2})+\frac{1}{2λ_{j}} (u_{s}^{j,k+1/2})^{2}+θ_{s}^{j}(u_{s}^{j,k+1/2})−ϕ_{l,s}^{1}u_{l}^{i,k+1/2}−ϕ_{l,s}^{2}u_{s}^{j,k+1/2}

with the constraint u_{l}^{i,k+1/2} = u_{s}^{j,k+1/2}. Obviously, it is transformed to scalar convex problem with one unknown u = u_{l}^{i,k+1/2} = u_{s}^{j,k+1/2} and can be handled as in the previous case.

4. Solution of a continuous casting problem

The continuous casting problem can be stated as follows. Let Ω be an open bounded domain in ℝ^{2} with the boundary Γ, T > 0, Q ≡ Ω×]0, T [ and Σ = Γ×]0, T [. We assume that the domain Ω is occupied by thermodynamically homogeneous and isotropic steel. We denote by H(x, t) the enthalpy related to unit mass and by u(x, t) the temperature for (x, t) ∈ Q.

Heat transfer within the two-dimensional domain is mathematically described by:

\frac{∂H}{∂t} + v(t)\frac{∂H}{∂x_{1}} = ∇ · (k(u)∇u) in Q,

where t is time, k(T) is thermal conductivity and v(t) is the actual casting speed, in x_{1} direction.

It is not restrictive to assume that k(u) ≡ 1, since otherwise we can substitute

u = ũ = ∫_{0}^{u} k(Θ) dΘ

and equation becomes:

\frac{∂H}{∂t} − Δu + v \frac{∂H}{∂x_{1}} = 0 in Q,

where Δ is Laplace operator.

We have also constitutive law

H = H(u) = ρ ∫_{0}^{u} c(Θ) dΘ + ρL (1 − f_{s}(u)) in Q,

where ρ is the density, c(u) is specific heat, L is latent heat and f_{s}(u) is solid fraction.

Graph H(u) is a increasing function ℝ → ℝ involving near vertical segments corresponding to the phase transition states, namely, for u ∈ [T_{L}, T_{S}] where 0 < T_{L} < T_{S} are melting and solidification temperatures, correspondingly.

We assume to know the initial distribution of the enthalpy (which determines also the initial temperature)

H(x, 0) = h_{0}(x) in Q.

We shall consider the case of heat flux imposed at the boundary

\frac{∂u}{∂n} + g(u) = 0 on Σ.

An example of the boundary condition, arising in continuous casting, is one of type of heat convection due to cooling water, namely: g(u) = h(u − u_{spray}) on Σ, where h is heat transfer coefficient and u_{spray} is cooling water temperature.

Let Ω = {0 < x_{1} < L, 0 < x_{2} < 1} be the rectangular domain with the boundary ∂Ω consisting of two parts:

Γ_{1} = {x ∈ ∂Ω : x_{1} = 0 ∨ x_{1} = L} and Γ_{2} = {x ∈ ∂Ω : x_{2} = 0 ∨ x_{2} = 1}.

We study the following boundary-value problem: find u = u(x, t) such that

(P)  \{ \begin{matrix}\frac{∂H(u)}{∂t} + v \frac{∂H(u)}{∂x_{1}} − Δu = 0  for x ∈ Ω, t > 0, \\ u = z(x_{2}, t) > 0 for x ∈ Γ_{1}, t > 0, \\ ∂u/∂n + au = g, a ≥ 0, g ≥ 0 for x ∈ Γ_{2}, t > 0, \\ u = u_{0}(x) > 0 for x ∈ Ω̄, t = 0. \end{matrix}

Here v = v(t) ≥ 0 is the speed of the slab being cast and n is outward unit normal vector.

The existence and uniqueness of the weak solution for problem (P) are proved in [11, 9].

To approximate problem (P) we rewrite it as the integral equality for fixed t > 0.

Let V = H^{1}(Ω), V^{0} = {u ∈ V : u(x) = 0 for x ∈ Γ_{1}}, V^{z} = {u ∈ V : u(x) = z for x ∈ Γ_{1}}. The solution of problem (P) for fixed t > 0 satisfies the following equality:

∫_{Ω}(∂H/∂t + v(t)∂H/∂x_{1})ηdx + ∫_{Ω} ∇u∇ηdx+
∫_{Γ_{2}} auηdΓ = ∫_{Γ_{2}} gηdΓ ∀η ∈ V^{0}, u(t) ∈ V^{z(t)}.

Let T_{h} be the triangulation of Ω in square elements of dimensions h_{1} × h_{2} and V_{h} = {u_{h}(x) ∈ H^{1}(Ω) : u_{h}(x) ∈ Q_{1} ∀δ ∈ T_{h}}, where Q_{1} is the space of bilinear functions. By Π_{h} we denote the local Q_{1}-interpolant. We also use the following notations: V_{h}^{0} = {u_{h}(x) ∈ V_{h} : u_{h}(x) = 0, ∀x ∈ Γ_{1}}, V_{h}^{z} = {u_{h}(x) ∈ V_{h} : u_{h}(x) = z_{h}, ∀x ∈ Γ_{1}} for the subsets of V_{h}. Here z_{h} is the V_{h} - interpolant of z on the boundary Γ_{1}. For any continuous function v(x) we put

S_{δ}(v) = ∫_{δ_{h}} Π_{h}(v)dx; S_{Ω}(v) = ∑_{δ∈T_{h}} S_{δ}(v).

S_{∂δ}(v) = ∫_{∂δ_{h}} Π_{h}(v)dx; S_{Γ_{2}}(v) = ∑_{∂δ∈T_{h}∩Γ̄_{2}} S_{∂δ}(v).

Let also ω_{τ} = {t_{k} = kτ, 0 ≤ k ≤ M, Mτ = T } be the uniform mesh in time on the segment [0, T].

To approximate the term (\frac{∂}{∂t} + v(t)\frac{∂}{∂x_{1}} ) H we use characteristics of this

first order differential operator (cf. eg. [2, 3]). Namely, if (x_{1}, x_{2}, t) is the mesh

point on the time level t we choose x̃_{1} = x_{1} −∫_{t−τ}^{t} v(ξ)dξ and approximate this

term as follows:

(\frac{∂}{∂t} + v(t)\frac{∂}{∂x_{1}} ) H ≈ \frac{1}{τ} (H(x_{1}, x_{2}, t) − H(x̃_{1}, x_{2}, t − τ)) .

If x̃_{1} < 0 then we put

H̃(x, t − τ) = H(0, x_{2}, t − τ) (= H(z(x_{2}, t − τ))).

In what follows we use the notation

dt̄H = \frac{1}{τ}(H(x, t) − H̃(x, t − τ))

for the difference quotient approximating the term (\frac{∂}{∂t} + v(t)\frac{∂}{∂x_{1}} ) H in each mesh point on time level t.

Then the finite difference scheme can be written as follows:

for all t ∈ ω_{τ}, t > 0, find u_{h} ∈ V_{h}^{z} such that
S_{Ω}(dt̄H_{h}η_{h}) + S_{Ω}(∇u_{h}∇η_{h}) + S_{Γ_{2}}(au_{h}η_{h}) = S_{Γ_{2}}(gη_{h}) ∀η_{h} ∈ V_{h}^{0}.  (32)

We rewrite the mesh scheme on the fixed time level t as nonlinear operator equation in Euclidian space of vectors. Let

ω̄ = {(x_{1}^{i}, x_{2}^{j}) = (ih, jh); 0 ≤ i ≤ N_{1}, 0 ≤ j ≤ N_{2}; hN_{1} = L, hN_{2} = 1},

ω = ω̄ \ {x_{1} = 0 ∨ x_{1} = L}, N = card ω_{.}

We put in the correspondence to the mesh function u_{h}(x) ∈ V_{h}^{0} the vector u ∈ ℝ^{N} of its nodal values: u_{h} ⇔ u.

Now we define the N × N matrices by the following relations:

for all V_{h}^{0} ∋ u_{h} ⇔ u ∈ ℝ^{N}, V_{h}^{0} ∋ η_{h} ⇔ η ∈ ℝ^{N}
(Au, η) = S_{Ω}(∇u_{h}∇η_{h}) + S_{Γ_{2}}(au_{h}η_{h}); (B_{0}u, η) = S_{Ω}(1/τu_{h}η_{h}).

By the similar way we define the vector f:

(f, η) = S_{Γ_{2}}(gη_{h}) + S_{Ω}(1/τH̃_{h}η_{h}).

Let now z̃_{h}(x) ∈ V_{h} be the function which is equal to z_{h} in Γ̄_{1} and 0 for all nodes in ω, then the vector f_{0} is defined by the equality:

(f_{0}, η) = S_{Ω}(∇z̃_{h}, ∇η_{h}) ∀η_{h} ∈ V_{h}^{0}.

Finally we get F = f − f_{0}.

In these notations the algebraic form for mesh scheme (32) at fixed time level can be written as follows:

Au + B_{0}H(u) = F.  (33)

Here A, B_{0} are symmetric and positive definite matrices (moreover B_{0} is diagonal one) and H(u) is vector with components (H(u))_{i} = H(u_{i}).

Thus, mesh scheme (32) at fixed time level is the partial case of (1) and we can use splitting algorithms of previous paragraphs for its solving.

Below several numerical results for the model problem are cited and analysed.

4.1. Solving the problem by using one processor. To validate the choice of the theoretical optimal parameter the following numerical example was considered. This numerical example has the same kind of boundary conditions as the real problem considering the continuous casting of steel.

Let Ω =]0, 1[×]0, 1[, with the boundary Γ divided in two parts such that Γ_{D} = {x ∈ ∂Ω : x_{1} = 0∨x_{1} = 1} and Γ_{N} = Γ\Γ_{D}, moreover let T = 1. Let us consider the case where the phase change temperature u_{SL} = 1 and the latent heat L = 1. Let the phase change interval be [u_{SL} − ε, u_{SL} + ε], ε = 0.01, and the velocity is v(t) = \frac{1}{5}.

Our numerical example is

\frac{∂H}{∂t} − ΔK + v(t)\frac{∂H}{∂x_{2}} = f(x; t)  on Ω,
u(x_{1}, x_{2}; t) = (x_{1} − \frac{1}{2})^{2} + (x_{2} − \frac{1}{2})^{2} − \frac{1}{2}e^{−4t} + 1  on Γ_{D},
\frac{∂u}{∂n} = 1  on Γ_{N},
u(x_{1}, x_{2}; 0) = (x_{1} − \frac{1}{2})^{2} + (x_{2} − \frac{1}{2})^{2} + \frac{1}{2}  on Ω,

where

H(u) = \{ \begin{matrix}2u & u < u_{SL} − ε, \\ (\frac{1+8ε}{2ε}) (u − 1) + \frac{5+4ε}{2} & u ∈ [u_{SL} − ε, u_{SL} + ε], \\ 6u − 3 & u > u_{SL} + ε, \end{matrix}

and

K(u) = \{ \begin{matrix}u & u < u_{SL} − ε, \\ \frac{3}{2}u − \frac{1−ε}{2} & u ∈ [u_{SL} − ε, u_{SL} + ε], \\ 2u − 1 & u > u_{SL} + ε. \end{matrix}

Furthermore

f(x; t) = \{ \begin{matrix}4e^{−4t} + \frac{1}{5}(4x_{2} − 2) − 4 & u < u_{SL}, \\ 12e^{−4t} + \frac{1}{5}(12x_{2} − 6) − 8 & u > u_{SL}. \end{matrix}

The exact solution of our problem is

u(x_{1}, x_{2}, t) = (x_{1} − \frac{1}{2})^{2} + (x_{2} − \frac{1}{2})^{2} − \frac{1}{2}e^{−4t} + 1.

We split the enthalpy function H(u) as follows:

H(u) = αu + H_{0}(u),

where α is the minimal slope of the enthalpy function. In our numerical example α = 2.

The algebraic form of is (33) with A = αI + τA_{0} where by A_{0} is noted the 5-point approximation of Laplacian operator, Cu = H_{0}(u) and F = τf + H̃(t − τ). Equation (33) is solved by using Douglas-Rachford splitting method (corresponds to choice B = I + λA):

λ^{−1}(u^{k+ \frac{1}{2}} − u^{k}) + Cu^{k+ \frac{1}{2}} ∋ F − Au^{k},  (34)
(I + λA)(u^{k+1} − u^{k}) = (u^{k+ \frac{1}{2}} − u^{k}).  (35)

Equation (35) is solved by using the conjugate gradient method.

The optimal choice of iterative parameter λ is λ = \frac{1}{\sqrt{mM}} with m = α +

τμ_{min}(A_{0}) and M = α + τμ_{max}(A_{0}).

In the numerical test the grid was 33 × 33 grid points in space and 64 time steps. In every time step the outer iterations was calculated till the L_{2}-norm of the residual is lower than the desired ϵ. In this case the ϵ = 10^{−3}.

λ Time [s] Maximum Average
0.05 7.41  41  22
0.075 6.98  28  15
0.1 6.24  22  12
λ_{0} 6.45  18  11
0.125 6.28  18  11
0.15 6.69  17  11
0.175 7.49  19  12
0.2 8.44  22  13
0.25 10.23  26  15
0.3 14.33  31  18

Table 1. The comparison of number of iterations for different λ.

In this case the theoretical optimal choice for the iterative parameter λ_{0} = 0.119587. From the table 1 we can see that the theoretical optimal parameter is very close to the practical optimal one. In this table Maximum means the maximum number of outer iterations in the whole time interval [0, T]. The Average means the average of outer iterations.

4.2. Solving the problem by using many processors. The problem was also solved by using the domain decomposition method. The domain Ω is divided into p non-overlapping subdomains by straight lines. In this case we have to solve the following kind of system of equations

D^{−1}(u^{k+ \frac{1}{2}} − u^{k}) + Cu^{k+ \frac{1}{2}} ∋ F − Au^{k}  (36)
(I_{i} + λ_{i}A_{i})(u^{k+1} − u^{k})_{i} = (u^{k+ \frac{1}{2}} − u^{k})_{i},  i = 1, .., p  (37)

where D = diag(λ_{1}I_{1}, ..., λ_{p}I_{p}).

Equation (37) is solved simultaneously in every subdomain by using conjugate gradient method.

We consider the same numerical example as in the case where we solve the problem by using only one processor in the whole domain Ω. The nonlinear equation is solved in the main processor and then the linear equation is solved in each subdomain simultaneously by using slave processors. We look for the case where the grid in space was 129 × 129 grid points and 256 time steps. The iterative parameter λ = 0.05 was the same in the every subdomain. It is close to the theoretical optimal value of the λ in the sequential case where it is λ = 0.06179.

The figure 1 shows the behaviour of the calculation time when the number of processor was changed.

processors Time [s] S  E
1 594.03 – –
2  417.87 1.42 0.711
4  237.53 2.5 0.625
6  190.96 3.11 0.518
8  148.61 3.99 0.499
10  129.89 4.57 0.457
12  115.00 5.17 0.43
14  105.70 5.62 0.40
16  100.43 5.92 0.369
18  90.98 6.53 0.364
20  90.753 6.55 0.326

Table 2. The behaviour of calculation time, speedup and efficiency when the grid is fixed and the number of processors is chanced.

Figure 1. The dependence of the calculation time on the number of processors.

We can see from figure 1 that the calculation time (solid line) first decrease quite rapidly when the number of processors is increased. After some number of processors the decreasing speed slows down. This is due to the fact that we use only one processor to solve the nonlinear equation and then send the quite big vector u^{k+ \frac{1}{2}} to the slave processors. When the number of processors increase the number of sending and receiving messages from processors also increase.

In the same picture is plotted the time consumption (dashed line) of solving nonlinear equation and sendings and receivings of the vector u^{k+ \frac{1}{2}} . From this picture we can see that when the number of processors is increased the time spend solving and sending of the vector u^{k+ \frac{1}{2}} becomes big comparing to the whole calculation time.

LARGE SPLITTING ITERATIVE METHODS 183

Another approach to solve the nonlinear equation is to solve it also simultaneously in slave processors. In this method we solve the nonlinear equation faster and don’t have to send the big vector u^{k+ \frac{1}{2}} to slave processors. On the other hand we have to send more smaller vectors from slave processors to another slave processors. These small vectors are in fact the values of the solution on the interfaces.

processors λ T_{old} [s] T_{new} [s]
2  0.05 417.87 359.07
4  0.05 237.53 233.73
8  0.05 148.61 147.47
12  0.05 115.0 133.79
14  0.05 105.7 118.99

Table 3. The comparison of calculation times when solving the nonlinear equation in main processor (T_{old}) and when the nonlinear equation is solved in slave processors (T_{new}).

We can see from the table 3 that when the number of processor is quite small the new technique is faster. When the number of processors increase the time consumption of the new method is bigger than in the old one. In this case also the time spend sending and receiving messages from processor to an other is a big part on the total calculation time.

References

[1] H. Brezis, Operateurs maximaux monotones et semigroups de contractions dans les espaces de Hilbert. (North–Holland Publishing Company, Amsterdam, 1973)

[2] Z. Chen, Numerical solutions of a two-phase continuous casting problem, Numerical Methods for Free Boundary Problem, P. Neittaanmaki, eds. International Series of Numerical Mathematics 99, Birkhauser, Basel, p. 103–121, 1991.

[3] J. Douglas Jr. and T.F. Russel, Numerical methods for convection-dominated diffusion problem based on combining the method of characteristic with finite element or finite difference procedures, SIAM J. Numer. Anal., V. 19, p. 871–885, 1982.

[4] E.G. D’yakonov, Optimisation in Solving Elliptic Problems, CRC Press, 562 p., 1995 [5] D. Gabay, Applications of the method of multipliers to variational inequalities // M.Fortin, R.Glowinski, eds. Augmented Lagrangian methods: applications to the numerical solution of boundary-value problems (North–Holland Publishing Company, Amsterdam 1983)

[6] P.L. Lions and B. Mercier, Splitting algorithms for the sum of two nonlinear operators, SIAM, J.Numer.Anal. V. 16, p. 964–979, 1979.

[7] A.V. Lapin and D.O. Solovyev, Splitting iterative methods for variational inequalities, Preprint No 783, Center of Calcul., Novosibirsk, 24 p., 1988

[8] R.T. Rockafellar, On the maximality of sums of nonlinear monotone operators, Trans.Amer.Math.Soc. V. 149, p. 75–88, 1970.

[9] J.F. Rodrigues and F. Yi, On a two-phase continuous casting Stefan problem with nonlinear flux, Euro. J. Appl. Math., V. 1, p. 259–278, 1990.

[10] P. Tseng, Further applications of splitting algorithm to decomposition in variational inequalities and convex programming, Math.Progr., V. 48, p. 249–263, 1991.

[11] F. Yi, An evolutionary continuous casting problem of two-phase and its periodic behaviour, J. Part. Diff. Eq., V. 2, p. 7–29, 1989.

184 E. LAITINEN, A.V. LAPIN, J.PIESKÄ

Department of Mathematical Sciences, University of Oulu, Oulu, Finland E-mail address: erkki.laitinen@oulu.fi, jpieska@cc.oulu.fi Department of Computational Mathematics and Cybernetics, Kazan State University, Kazan, Russia E-mail address: alapin@ksu.ru

Received June, 20, 2001