Lobachevskii Journal of Mathematics Vol. 13, 2003, 67 – 80 \begin{matrix}c \\ ◯ \end{matrix}E. Laitinen, A.V. Lapin, and J. Pieskä

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

NUMERICAL EXPERIMENTS WITH MULTILEVEL SUBDOMAIN DECOMPOSITION METHOD

Abstract. In this paper we present a new numerical approach to solve the continuous casting problem. The main tool is to use so-called IPEC method and DDM similar to [6] with multilevel domain decomposition. On the subdomains we use the multidecomposition of the subdomains. The IPEC is used both in the whole calculation domain and inside the subdomains. Calculation algorithm is presented and numerically tested. Several conclusions are made and discussed.

1. Introduction

Theory of the so-called regional-additive schemes (splitting schemes with domain decomposition) for linear diffusion and convection-diffusion have been studied in [14],[15] and [16]. The stability have been proved and error estimates have been deduced. For the non-linear problems like our their technique is not available.

Several new finite-difference schemes for a nonlinear convection-diffusion problem are constructed and numerically studied in [6]. These schemes are constructed on the basis of non-overlapping domain decomposition and predictorcorrector approach. (Note that the term “predictor-corrector domain decomposition method” was introduced by Quarteroni and Valli [13]).

The paper of Lapin and Pieskä [6] was motivated by [2], [9], [10], [11], where TL^{1}, EP^{2} and EPIC^{3} methods have been studied and tested. The EPIC method was proved to be stable and scalable when solving on a big number of processors. In the paper [6] the scheme from [10], [11] was modified in such a way, that its implementation leads to IPEC^{4} method.

1time lagging

2explicit predictor

3explicit predictor-implicit corrector

4implicit predictor-explicit corrector

Key words and phrases. Stefan problem, continuous casting process, finite element method, predictor-corrector scheme, domain decomposition.

This work is supported by RFBR, project N 01-01-00068.

The main idea of these kind of algorithms is first to solve the problem in artificial boundaries predictor step). After the solution at the boundaries is known then it can be used as a Dirichlet-type boundary condition and the non-coupled subdomain problems can be solved in parallel. The last step of these methods is to correct the solution at the artificial boundaries (corrector step). The idea is similar to Schur’s complement type methods but in our approach the calculation of the artificial boundaries between subdomains are very easy and cheap to perform. There is also no need to construct a good preconditioner for the interface problems.

The advantages of predictor-corrector methods (IPEC or EPIC) is that we reduce the amount of information sending between processors. We need to send only once the subsolutions from slave processors to master processor. When we use Schwarz alternating methods with the overlapping subdomains, the number of sending and receiving is much more bigger. The numerical experiences show that the speedup of IPEC method is very good comparing with Schwarz methods and calculation times are roughly half of the times of Schwarz methods. Moreover, the linear speedup is reached in the numerical tests of [6].

The idea of multidecomposition method MDD is to use DDM with IPEC inside the subdomains. The subdomain is divided into smaller subdomains and then IPEC method is used to solve these smaller subproblems sequently.

The rest of the article is organized as follows. In the section 2 we present the continuous casting problem. The section 3 deals with the mesh approximation of the problem. The used domain decomposition and calculation algorithm is presented in the section 4. The section 5 is devoted to numerical testing of the new algorithm. Finally, in the section 6 the conclusions are made.

2. Problem statement

The continuous casting problem can be mathematically formulated as follows. Let the rectangular domain Ω ⊂ ℝ^{2}, Ω = (0, l_{1})×(0, l_{2}) be occupied by a thermodynamically homogeneous and isotropic steel. We denote by H̄(x, t) the enthalpy related to the unit mass and by T (x, t) the temperature for (x, t) ∈ Ω̄ × [0, t_{f}]. We have constitutive law

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

where ρ is the density, c(T) is the specific heat, L is the latent heat and f_{s}(T) is the solid fraction at temperature T of the form

f_{s}(T) = \begin{cases}1 for T < T_{S}, \\ \frac{T_{L} − T}{T_{L} − T_{S}} , for T_{S} < T < T_{L}, \\ 0 for T > T_{L}. \end{cases}

The graph H̄(T) is an increasing function ℝ → ℝ, involving near vertical segment, which corresponds to a phase transition state, namely, for T ∈ [T_{S}, T_{L}]. In the case of piecewise constant specific heat c(T), the enthalpy function is piecewise linear of the form

H̄(T) = \begin{cases}ᾱ_{1}T + β̄_{1} for T < T_{S}, \\ ᾱ_{2}T + β̄_{2}, for T_{S} < T < T_{L}, \\ ᾱ_{3}T + β̄_{3} for T > T_{L}. \end{cases}

Further by k(T) we denote the thermal conductivity coefficient, which is continuous and increasing in T.

A continuous casting process can be described by a boundary-value problem, formally written in the following pointwise form: find T (x, t) and H̄(x, t) such that

\frac{∂H̄}{∂t} + v \frac{∂H̄}{∂x_{2}} − ∇ · (k(T)∇T) = 0 for x ∈ Ω, t > 0,

T = z̄(x, t) for x ∈ Γ_{D}, t > 0,

k(T)\frac{∂T}{∂n} = g, for x ∈ Γ_{N}, t > 0,

H̄ = H̄_{0}(x) for x ∈ Ω̄, t = 0,

where v = const > 0 is a casting speed in x_{2}-direction, Γ_{D} ∪ Γ_{N} = ∂Ω is the boundary of the domain, below Γ_{D} = {x ∈ ∂Ω : x_{2} = 0 ∨ x_{2} = l_{2}}.

Using Kirchoff’s transformation u = K(T) = ∫_{0}^{T} k(ξ)dξ and the notation

H(u) = H̄(T) = H̄(K^{−1}(u)), we can rewrite the continuous casting problem as

(P) \{ \begin{matrix}\frac{∂H}{∂t} + v \frac{∂H}{∂x_{2}} − Δu = 0, for x ∈ Ω, t > 0, \\ u = z(x, t) for x ∈ Γ_{D}, t > 0, \\ \frac{∂u}{∂n} = g, for x ∈ Γ_{N}, t > 0, \\ H = H_{0}(x) for x ∈ Ω̄, t = 0, \end{matrix}

The existence and uniqueness of a weak solution for problem (P) are proved in [12].

In the case of piecewise constant specific heat c(u) the enthalpy function takes the form

(1)  H(u) = \begin{cases}α_{1}u + β_{1} for u < u_{S}, \\ α_{2}u + β_{2}, for u_{S} < u < u_{L}, \\ α_{3}u + β_{3} for u > u_{L}. \end{cases}

3. Mesh approximation of continuous casting problem

We approximate problem (P) by an implicit in time finite difference scheme and by a semi-implicit finite difference scheme, using for the approximation in space variables a finite element method with the quadrature rules.

Let Ξ_{h} be a partitioning of Ω in the rectangular elements δ of dimensions h_{1} × h_{2} and V_{h} = {u_{h}(x) ∈ H^{1}(Ω) : u_{h}(x) ∈ Q_{1} for all δ ∈ Ξ_{h}}, where Q_{1} is the space of bilinear functions. By Π_{h}v(x) we denote the V_{h}-interpolant of a continuous function v(x), i.e. Π_{h}v(x) ∈ V_{h} and coincides with v(x) in the mesh nodes (vertices of all δ ∈ Ξ_{h}). We also use an interpolation operator P_{h}, which is defined as follows: for any continuous function v(x) the function P_{h}v(x) is piecewise linear in x_{1}, piecewise constant in x_{2} and on δ = [x_{1}, x_{1} + h_{1}] × [x_{2}, x_{2} + h_{2}] it coincides with v(x) at (x_{1}, x_{2} + h_{2}) and (x_{1} + h_{1}, x_{2} + h_{2}).

Let further V_{h}^{0} = {u_{h}(x) ∈ V_{h} : u_{h}(x) = 0 for all x ∈ Γ_{D}}, V_{h}^{z} = {u_{h}(x) ∈ V_{h} : u_{h}(x) = z_{h} for all x ∈ Γ_{D}}. Here z_{h} is the bilinear interpolation of z on the boundary Γ_{D}. For any continuous function v(x) we define the quadrature formulas:

S_{δ}(v) = ∫_{δ} Π_{h}vdx, S_{∂δ}(v) = ∫_{∂δ} Π_{h}vdx, E_{δ}(v) = ∫_{δ} P_{h}vdx,

S_{Ω}(v) = ∑_{δ∈Ξ_{h}} S_{δ}v, S_{Γ_{2}}(v) = ∑_{∂δ∈Ξ_{h}∩Γ̄_{2}} S_{∂δ}(v), E_{Ω}(v) = ∑_{δ∈Ξ_{h}} E_{δ}(v).

Let also ω_{τ} = {t_{k} = kτ, 0 ≤ k ≤ M, Mτ = t_{f}} be an uniform mesh in time

on the segment [0, t_{f}] and ∂t̄H = \frac{1}{τ}(H(x, t) − H(x, t − τ)). Then the implicit in time finite difference scheme with up-wind approximation of the convective term v ∂H/∂x_{2} can be written as follows: for all t ∈ ω_{τ}, t > 0, find u_{h} ∈ V_{h}^{z} and H_{h} ∈ V_{h} such, that

(2) S_{Ω}(∂t̄H_{h}η_{h}) + E_{Ω}(v \frac{∂H_{h}}{∂x_{2}} η_{h}) + S_{Ω}(∇u_{h}∇η_{h}) = S_{Γ_{2}}(gη_{h}) for all η_{h} ∈ V_{h}^{0}.

When constructing the characteristic mesh scheme we approximate the term

(\frac{∂}{∂t} + v \frac{∂}{∂x_{2}} ) H by using the characteristics of the first order differential op-

erator (similar to [1], [3]). Namely, if (x_{1}, x_{2}, t) is the mesh point on the time

level t we choose x̃_{2} = x_{2} − ∫_{t−τ}^{t} v(ξ)dξ and approximate:

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

where we denote H̃(x, t − τ) = H(x_{1}, x̃_{2}, t − τ).

Near the boundary it can happen that x̃_{2} < 0. In that case we put H̃(x, t −

τ) = H(x_{1}, 0, t − τ). In what follows we use the notation dt̃H = \frac{1}{τ}(H(x, t) −

H̃(x, t − τ)) for the difference quotient in each mesh point on time level t.

Now, the characteristic finite difference scheme for problem (P) is: for all t ∈ ω_{τ}, t > 0, find u_{h} ∈ V_{h}^{z} and H_{h} ∈ V_{h} such that

(3)  S_{Ω}(dt̃H_{h}η_{h}) + S_{Ω}(∇u_{h}∇η_{h}) = S_{Γ_{2}}(gη_{h}) for all η_{h} ∈ V_{h}^{0}

Let N_{0} = card V_{h}^{0} and u ∈ ℝ^{N_{0}} be the vector of nodal values for u_{h} ∈ V_{h}^{0}. We use the writing u_{h} ⇔ u for this bijection. We define N_{0} × N_{0} matrices by the following relations: for all u, η ∈ ℝ^{N_{0}}, u ⇔ u_{h} ∈ V_{h}^{0} and η ⇔ η_{h} ∈ V_{h}^{0},

(Ãu, η) = S_{Ω}(∇u_{h}∇η_{h}), (Mu, η) = S_{Ω}(u_{h}η_{h}), (C̃u, η) = E_{Ω}(v \frac{∂u_{h}}{∂x_{2}} η_{h}),

A_{0} = M^{−1}Ã, C = M^{−1}C̃.

Let now z̃_{h}(x) ∈ V_{h} be the function, which is equal to z_{h} on Γ̄_{D} and 0 for all nodes in Ω ∪ Γ_{N}. Then a vector f is defined by the equality

(f, η) = S_{Γ_{2}}(gη_{h}) − S_{Ω}(∇z̃_{h}, ∇η_{h}) ∀η ∈ ℝ^{N_{0}}, η ⇔ η_{h} ∈ V_{h}^{0},

and we set F = M^{−1}f.

In these notations the algebraic form for implicit mesh scheme (2) at fixed time level is:

(4)  ∂t̄H + A_{0}u + CH = F,

while characteristic mesh scheme (3) becomes

(5)  dt̃H + A_{0}u = F,

It is easy to see, that A_{0} is the standard five-point finite difference approximation of Laplace operator, A_{0}u = −u_{x_{1}x̄_{1}} − u_{x_{2}x̄_{2}} for the internal mesh points with the notations u_{x_{1}} = h^{−1}(u(x_{1} +h_{1}, x_{2})−u(x_{1}, x_{2})), ux̄_{1} = h^{−1}(u(x_{1}, x_{2})− u(x_{1} − h_{1}, x_{2})), and similarly for u_{x_{2}} and ux̄_{2}. For more detailed writing of the explicit form for A_{0}u let us introduce several sets of the grid points. Namely, let ω̄ be the set of all grid points, i.e. vertices of finite elements, γ_{D} = Γ̄_{D}∩ω̄, γ_{N} = Γ_{N} ∩ ω̄, ω = Ω ∩ ω̄, γ_{N}^{−} = {x ∈ γ_{N} : x_{1} = 0}, γ_{N}^{+} = {x ∈ γ_{N} : x_{1} = l_{1}}.

Now A_{0} = A^{1} + A^{2} with

A^{1}u = \{ \begin{matrix}−u_{x_{1}x̄_{1}} for x ∈ ω, \\ −2h_{1}^{−1}u_{x_{1}} for x ∈ γ_{N}^{−}, \\ 2h_{1}^{−1}ux̄_{1} for x ∈ γ_{N}^{+}, \end{matrix}

A^{2}u = \{ \begin{matrix}−u_{x_{2}x̄_{2}} for x ∈ ω ∪ γ_{N}, \\ 0 for x ∈ γ_{D}. \end{matrix}

The term CH in the implicit scheme corresponds to an up-wind approximation of the nonlinear convective term, CH = Hx̄_{2} for x ∈ ω ∪ γ_{N}.

4. Domain decomposition by straight lines

In this section we present the IPEC algorithm [6]. We restrict our discussion to the case of decomposition by unidirect straight lines. More variations and possibilities of decomposition are discussed and tested in [6]. Methods presented here and in the following section can be easily implemented and used with more complicated decompositions of calculation domain. Moreover, the generalization to 3-D case is straightforward. In that case the 1-D lines are replaced by 2-D planes.

Let the domain Ω be decomposed into two subdomains Ω_{1} and Ω_{2} by a straight line S_{y} in x_{2}-direction, which is also a grid line. We denote by δ_{S_{y}} the characteristic function of this line, i.e., the mesh function δ_{S_{y}}(x) = 1 for x ∈ S_{y} ∩ ω̄, while δ_{S_{y}}(x) = 0 for other mesh points. Also, let ω̄_{k}, k = 1, 2 be the corresponding to the subdomains Ω̄_{k} sets of grid points, S_{y} being the common part of their boundaries.

Let A_{2}u = −δ_{S_{y}}u_{x_{1}x̄_{1}} and A_{1} = A_{0} − A_{2},

A_{1}u = \{ \begin{matrix}−(1 − δ_{S_{y}})u_{x_{1}x̄_{1}} − u_{x_{2}x̄_{2}} for x ∈ ω, \\ −2h_{1}^{−1}u_{x_{1}} − u_{x_{2}x̄_{2}} for x ∈ γ_{N}^{−}, \\ 2h_{1}^{−1}ux̄_{1} − u_{x_{2}x̄_{2}} for x ∈ γ_{N}^{+}. \end{matrix}

Now, instead of implicit scheme (4) we consider the following scheme on the time level t_{n+1} = (n + 1)τ:

(6)  \frac{1}{τ}(H^{n+1/2} − H^{n}) + A_{1}u^{n+1/2} + A_{2}u^{n} + CH^{n+1/2} = F,

(7)  \begin{matrix}\frac{δ_{S_{y}}}{τ} (H^{n+1} − H^{n}) + \frac{1 − δ_{S_{y}}}{τ} (H^{n+1} − H^{n+1/2}) & + \\ δ_{S_{y}}A_{1}u^{n+1/2} + A_{2}u^{n+1} + δ_{S_{y}}CH^{n+1} & = δ_{S_{y}}F. \end{matrix}

Similarly, characteristic scheme (5) is changed by

(8)  \frac{1}{τ}(H^{n+1/2} − H̃^{n}) + A_{1}u^{n+1/2} + A_{2}u^{n} = F,

(9) \frac{δ_{S_{y}}}{τ} (H^{n+1}−H̃^{n})+ \frac{1 − δ_{S_{y}}}{τ} (H^{n+1}−H^{n+1/2})+δ_{S_{y}}A_{1}u^{n+1/2}+A_{2}u^{n+1} = δ_{S_{y}}F.

Let us discuss the implementation of scheme (6),(7). In the points of S_{y} first equation of (6) has the form:

(10)  \frac{H^{n+1/2} − H^{n}}{τ} − u_{x_{2}x̄_{2}}^{n+1/2} − u_{x_{1}x̄_{1}}^{n} + Hn+1/2 ; x̄_{2} = F,

i.e. in the points of S_{y} we have one-dimensional problem (10), that we solve first. After, the rest part of the first equation in (6) is splitted into two noncoupled implicit schemes in the subdomains:

(11) \begin{cases}\frac{H^{n+1/2} − H^{n}}{τ} − u_{x_{1}x̄_{1}}^{n+1/2} − u_{x_{2}x̄_{2}}^{n+1/2} + Hn+1/2 ; x̄_{2} = F, for x ∈ ω_{1} ∪ ω_{2}, \\ \frac{H^{n+1/2} − H^{n}}{τ} − \frac{2}{h_{1}} u_{x_{1}}^{n+1/2} − u_{x_{2}x̄_{2}}^{n+1/2} + Hn+1/2 ; x̄_{2} = F, for x ∈ γ_{N}^{−}, \\ \frac{H^{n+1/2} − H^{n}}{τ} + \frac{2}{h_{1}} un+1/2 ; x̄_{1} − u_{x_{2}x̄_{2}}^{n+1/2} + Hn+1/2 ; x̄_{2} = F, for x ∈ γ_{N}^{+} \end{cases}

and these equations are accomplished by Dirichlet boundary conditions, given on γ_{D} and calculated from (10) on S_{y}.

Further, it is easy to check, that for the points x ∉ S_{y} the second equation (7) coincide with (11) and has the same Dirichlet boundary conditions on γ_{D} and on S_{y}, so, u^{n+1}(x) = u^{n+1/2}(x) for x ∉ S_{y} and it makes no sense to solve these equations. It remains only to solve the system of the equations, corresponding to x ∈ S_{y}:

\frac{H^{n+1} − H^{n}}{τ} − u_{x_{2}x̄_{2}}^{n+1/2} − u_{x_{1}x̄_{1}}^{n+1} + Hn+1/2 ; x̄_{2} = F.

As u^{n+1}(x) = u^{n+1/2}(x) for x ∉ S_{y}, this system becomes

(12) \begin{cases}\frac{H^{n+1} − H^{n}}{τ} + Hn+1/2 ; x̄_{2} + 2\frac{u^{n+1}(x_{1}, x_{2})}{h_{1}^{2}} − u_{x_{2}x̄_{2}}^{n+1/2} \\ −\frac{u^{n+1/2}(x_{1} − h_{1}, x_{2}) + u^{n+1/2}(x_{1} + h_{1}, x_{2})}{h_{1}^{2}} = F, x ∈ S_{y}, \end{cases}

i.e. we get the system of scalar equations for (u^{n+1}(x), H^{n+1}(x)), x ∈ S_{y}.

Thus, the algorithm for the implementation of (6),(7) consists of 3 steps: 1): Predictor step: solving one-dimensional problem (10); 2): Main step: concurrent solving subproblems (11); 3): Corrector step: solving the system of scalar equations (12).

With the slight modifications the implementation of scheme (8),(9) is similar. Namely, in the points of S_{y} first equation of (8) has the form:

(13)  \frac{H^{n+1/2} − H̃^{n}}{τ} − u_{x_{2}x̄_{2}}^{n+1/2} − u_{x_{1}x̄_{1}}^{n} = F,

i.e. in the points of S_{y} we have one-dimensional problem (13), that we solve first. After that the rest part of first equation in (8) is splitted into two noncoupled characteristic schemes in the subdomains:

(14)  \begin{cases}\frac{H^{n+1/2} − H̃^{n}}{τ} − u_{x_{1}x̄_{1}}^{n+1/2} − u_{x_{2}x̄_{2}}^{n+1/2} = F, for x ∈ ω_{1} ∪ ω_{2}, \\ \frac{H^{n+1/2} − H̃^{n}}{τ} − \frac{2}{h_{1}} u_{x_{1}}^{n+1/2} − u_{x_{2}x̄_{2}}^{n+1/2} = F, for x ∈ γ_{N}^{−}, \\ \frac{H^{n+1/2} − H̃^{n}}{τ} + \frac{2}{h_{1}} un+1/2 ; x̄_{1} − u_{x_{2}x̄_{2}}^{n+1/2} = F, for x ∈ γ_{N}^{+} \end{cases}

and these equations are accomplished by Dirichlet boundary conditions, given on γ_{D} and calculated from (13) on S_{y}.

Finally it remains to solve the system of the equations, corresponding to x ∈ S_{y}:

\frac{H^{n+1} − H̃^{n}}{τ} − u_{x_{2}x̄_{2}}^{n+1/2} − u_{x_{1}x̄_{1}}^{n+1} = F.

As u^{n+1}(x) = u^{n+1/2}(x) for x ∉ S_{y}, this system becomes

(15) \begin{cases}\frac{H^{n+1} − H̃^{n}}{τ} + 2\frac{u^{n+1}(x_{1}, x_{2})}{h_{1}^{2}} − u_{x_{2}x̄_{2}}^{n+1/2} \\ −\frac{u^{n+1/2}(x_{1} − h_{1}, x_{2}) + u^{n+1/2}(x_{1} + h_{1}, x_{2})}{h_{1}^{2}} = F, x ∈ S_{y}, \end{cases}

i.e. we get the system of scalar equations for (u^{n+1}(x), H^{n+1}(x)), x ∈ S_{y}.

Remark 1. Above we assumed that the calculation domain Ω is divided into two parts Ω_{1} and Ω_{2}. This is not restrictive and we refer to [6] for more detailed discussion of the decomposition into several subdomains, also with corner and cross points. Even the case of curvilinear decomposition is studied in [6] and found to be stable and accurate under the natural assumptions for the mesh parameters.

Remark 2. We consider the 2D-case. The proposed methods have natural extensions to the 3D-case. We notice that in predictor and corrector steps the one-dimensional problems corresponding to the mesh line S_{y} are replaced by the two-dimensional problems corresponding to a plane S_{xy}.

4.1. Multidecomposition method. The general idea of the multidecomposition is to divide the subdomain to smaller subdomains i.e. we use two-level decomposition of the calculation domain. The main reason is to decrease the algebraic size of the problems and thus make the calculation times much smaller. The division of the subdomains is presented in the figure 1. The idea of using multilevel decomposition is not new but the technique presented in this paper gives a good and effective method when solving large complex time dependent problems.

Figure 1. Used nonoverlapping domain decomposition and multidecomposition of the subdomains.

We use the notation Ω_{i} = ⋃_{j_{i}=1}^{p_{i}} Ω_{i,j_{i}}. The calculation algorithm for charac-

teristic mesh scheme (13)-(15) is presented below. The algorithm for implicit mesh scheme (10)-(12) is similar.

Algorithm 1. :

1.: On a time level n perform on the main processor the predictor step (13) on the artificial boundary S_{y}. 2.: Send the values of u^{n+1/2} and H^{n+1/2} on the boundary S_{y} to the slave processors. 3.: Concurrently on the slave processors perform the predictor step (13) on the artificial boundaries of the subdomains Ω_{i,j_{i}}, i = 1, 2, j_{1} = 1, ..., p_{1}, j_{2} = 1, ..., p_{2}. 4.: Concurrently on the slave processors perform sequentially the main step (14) on the subdomains Ω_{i,j_{i}}. 5.: Concurrently on the slave processors perform the corrector step (15) on the artificial boundaries of the subdomains Ω_{i,j_{i}}, i = 1, 2, j_{1} = 1, ..., p_{1}, j_{2} = 1, ..., p_{2}. 6.: Send the values of u^{n+1} and H^{n+1} from slave processors to the master processor.

7.: On the main processor perform the corrector step (15) on the artificial boundary S_{y}. 8.: Put n = n + 1, if the final time t_{f} is reached then STOP, else GOTO 1.

Remark 3. On the step 3. of the algorithm 1 we do not do the predictor step (13) on the artificial boundary S_{y}. Remark 4. On the steps 3.-5. we do the calculations concurrently. However, we do not synchronize our calculation in such a way that all processors move from step to another at the same time. The step 6. is synchronization point and it is performed at the same time.

The use of high number of subdomains inside the subdomain may increase the error dramatically. To overcome this feature we introduce so called smoothing steps to our method. Namely, our algorithm 1 is replaced by the following one. Algorithm 2. :

1.: On the time step n perform on the main processor the predictor step (13) on the artificial boundary S_{y}. 2.: Send the values of u^{n+1/2} and H^{n+1/2} on the boundary S_{y} to the slave processors. 3.: Concurrently on the slave processors perform the predictor step (13) on the artificial boundaries of the subdomains Ω_{i,j_{i}}, i = 1, 2, j_{1} = 1, ..., p_{1}, j_{2} = 1, ..., p_{2}. 4.: Concurrently on the slave processors perform sequentially the main step (14) on the subdomains Ω_{i,j_{i}}. 5.: Concurrently on the slave processors perform the corrector step (15) on the artificial boundaries of the subdomains Ω_{i,j_{i}}, i = 1, 2, j_{1} = 1, ..., p_{1}, j_{2} = 1, ..., p_{2}. 6.: On the slave processors perform the smoothing step i.e. few iterations of the MSOR-method (modified SOR-method) over the whole subdomain Ω_{i}. 7.: Send the subsolutions u^{n+1} and H^{n+1} from slave processors to the main processor. 8.: On the main processor perform the corrector step (15) on the artificial boundary S_{y}. 9.: On the master processor perform few iterations of the MSORmethod in the neighborhood of S_{y}. 10.: Put n = n + 1, if the final time t_{f} is reached then STOP, else GOTO 1.

Remark 5. In the algorithm 2 we use the smoothing on steps 6. and 9. According to our knowledge the good amount of MSOR-iterations is less than 10. This amount ensures the reducing of the error and do only slightly increase the calculation time. Remark 6. The step 9. in the algorithm 2 is performed in the neighborhood of the artificial boundary S_{y}. The good width of the smoothing area is about 10 grid lines i.e from artificial boundary 5 grid lines to both directions.

Remark 7. In both algorithms 1 and 2 we use the predictor step 3. and corrector step 5. which require the decomposition along and against the convection. There are also cross points. In the article of Lapin and Pieskä [6] are discussed more carefully the implementation of such kind of decomposition.

Remark 8. The use of only two subdomains and processors is not restrictive. The usage of more processors and subdomains is straightforward. The numerical results presented in the section 5 verify this.

4.2. Multidecomposition method with one processor. The multidecomposition algorithms 1 and 2 are very effective and extremely quick in the case of many processors. We used the idea of multidecomposition for the situation where we have only one processor. Now we do not have artificial boundaries which decouples problem into subproblems. Anyway, we still have boundaries between subdomains inside the whole domain. The multidecomposition algorithm for one processor with characteristic mesh scheme read as follows.

Algorithm 3. :

1.: On every time level n perform the predictor step (13) on the artificial boundaries of the subdomains Ω_{1,j}, j = 1, ..., p. 2.: Perform sequentially the main step (14) on the subdomains Ω_{1,j}. 3.: Perform the corrector step (15) on the artificial boundaries of the subdomains Ω_{1,j}, j = 1, ..., p. 4.: Perform the smoothing step i.e. few iterations of the MSOR-method over the whole calculation domain Ω. 5.: Put n = n + 1, if the final time t_{f} is reached then STOP, else GOTO 1.

5. Numerical experiments with one processor

Let Ω =]0, 1[×]0, 1[ with the boundary Γ divided in two parts such that Γ_{D} = {x ∈ ∂Ω : x_{2} = 0∨x_{2} = 1} and Γ_{N} = Γ\Γ_{D}, moreover let t_{f} = 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} − \frac{1}{2}e^{−4t} + \frac{5}{4}  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 Kirchoff’s temperature is according to it’s definition

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

and enthalpy

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

Furthermore, let the known right-hand side be

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

The stopping criterion of the calculations was the L_{2}-norm of residual: ∥r∥_{L_{2}(Ω)} < 10^{−4}. We divide the domain Ω into non-overlapping subdomains Ω_{i} presented in the figure 1.

We tested algorithm 3 in two different cases. First, we fixed the number of grid points and changed the number of inner subdomains. The calculation grid was 129 × 129 in space and we took 256 time steps. The calculation times for different number of inner subdomains are presented in the table 1.

♯ of inside subdomains Time [s]
1 × 1  112.9 s
2 × 2  76.2 s
3 × 3  66.1 s
4 × 4  56.9 s
8 × 8  39.9 s

Table 1. Calculation times in seconds for MDD when the grid size is fixed and number of inside subdomains are changed.

In the table 1 the 1 × 1 subdomain division means the calculation times for MSOR method for the whole calculation domain.

In the second test case we changed the calculation grid and number of time steps. We fixed the number of inner subdomains to be 4 × 4. For many calculation grid this decomposition is not the optimal one but it very clearly emphasizes the advantage of MDD. The results are in the table 2.

Grid  SEQ MDD
65 × 65 × 128 8.67 s 4.75 s
129 × 129 × 256 112.9 s 56.9 s
257 × 257 × 512 1425 s 687 s

Table 2. Calculation times in seconds for sequential MSOR (SEQ) and MDD when the calculation grid is changed.

In the figure 2 we measured the time which different algorithms (MSOR and MDD) spend at different time levels. In this case we had the grid size 129×129. Normally for this kind of mesh we would take 256 times steps. Now we liked to test how these methods behave when the solution goes to steadystate situation. We took 2560 time steps i.e. 10 times more than usually to see the behavior. The MDD is much better comparing to MSOR in the beginning of the calculations when we have big changes in our simulation.

When time increases the calculation time drop very rapidly and both methods spend equal amount of time for each time step. This is good result for MDD because it is very effective for the case when we have big changes. Near steady-state situation it is as good as MSOR.

Figure 2. Time spend on each time level for different methods.

6. Conclusions

The numerical examples show that the multidecomposition method (MDD) is very effective numerical method when solving continuous casting problem. The idea to divide the subdomains to smaller subdomains seems to be very good and profitable. The algebraic dimensions of the subproblems inside the subdomains are very small and thus they are very quick to solve. The dimension is in many cases so small that even direct solvers could be effective.

The numerical results for one processor seem to be very promising. We do not need to have big parallel computers to achieve advantages of parallel computer. Only few processors are enough and in some cases even only one processor is good. The table 2 very clearly shows the advantages of MDD.

NUMERICAL EXPERIMENTS 79

In that table the decomposition is clearly not optimal. The table 1 is good indicator to that. Anyhow, even with poor knowledge of the problem MDD method can be used. When the system is stable and big changes do not appear (like the number of grid points changes or time step changes) then MDD can be optimized good and numerical advantages comparing with MSOR method in the case of one processor and with additive Schwarz alternating method in the case of many processor becomes very dramatical.

References

[1] Z. Chen, Numerical solutions of a two-phase continuous casting problem, Numerical Methods for Free Boundary Problem, P. Neittaanmäki, eds. International Series of Numerical Mathematics 99., Birkhäuser, Basel., pp. 103-121, 1991.

[2] C. N. Dawson, Q. Du and T. F. Dupont, A finite difference domain decomposition algorithm for numerical solution of the heat equation, Mathematics of Computation, vol. 57, pp. 63-71, 1991.

[3] J. Jr. Douglas 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., pp. 871-885, 1982.

[4] E. Laitinen and A.V. Lapin, Semi-Implicit Mesh Scheme and Splitting Iterative Methods for the Solution of Continuous Casting Problem, Preprint, University of Oulu, Finland, 19p., 1999

[5] E. Laitinen, A.V Lapin and J. Pieskä, Splitting iterative methods and parallel solution of variational inequalities, Lobachevskii Journal of Mathematics, Vol. 8, p. 167-184, 2001. http://ljm.ksu.ru/vol8/latin.htm

[6] A.V Lapin and J. Pieskä, On the parallel domain decomposition algorithms for timedependent problems, Lobachevskii Journal of Mathematics, Vol. 10, p. 27-44, 2002. http://ljm.ksu.ru/vol10/lpp.htm

[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] J.M Ortega and W.C. Rheinboldt, Iterative solution of nonlinear equations in several variables, Academic Press Inc. Orlando, Florida 1970.

[9] W. Rivera, J. Zhu and D. Huddleston, An efficient parallel algorithm with application to computational fluid dynamics, To be appear in Computers and Mathematics with Applications.

[10] W. Rivera and J. Zhu, A scalable parallel domain decomposition algorithm for solving time dependent partial differential equations, Proceedings of the International Conference on Parallel and Distributed Processing Technology and Applications, edited by H. R. Arabnia, CSREA Press, Athens GA, pp. 240-246, 1999.

[11] W. Rivera, J. Zhu and D. Huddleston, An efficient parallel algorithm for solving unsteady nonlinear equations, Proceedings of the International Conference on Parallel Processing Workshops, edited by T. M. Pinkston, IEEE Computer Society, Los Alamitos, California, pp. 79-84, 2001.

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

[13] A. Quarteroni and A. Valli, Domain decomposition methods for partial differential equations, Clarendon Press Oxford, New York 1999.

[14] A.Samarskii, P.Vabischevich, Finite difference schemes with operator multipliers, Minsk, Institute of Mathematics of Belorussia, 442pp, 1998.

[15] A.Samarskii, P.Vabischevich, Factorized regional-additive schemes for convectiondiffusion problems (in Russian), Reports of Rusian Acad.Sciences (Mathematics), 1996, V.346, P/742-745.

80 E. LAITINEN, A.V. LAPIN, AND J. PIESKÄ

[16] P.Vabischevich, Parallel domain decomposition algorithms for time-dependent problems of mathematical physics , Advances in Numerical Methods and applications. Singapore: World scientific, 1994, P.293-299.

University of Oulu, Finland. Kazan State University, Russia. E-mail address: alapin@ksu.ru University of Oulu, Finland.

ReceivedOctober 15, 2003