INFORMATICA Informatica 1822-88440868-49520868-4952 Vilnius University INFOR481 10.15388/22-INFOR481 Research Article Parallel 3D ADI Scheme for Partially Dimension Reduced Heat Conduction Problem ČiegisRaimondasrc@vgtu.lt1

R. Čiegis has graduated from Vilnius University Faculty of Mathematics, in 1982. He received the PhD degree from the Institute of Mathematics of Byelorussian Academy of Science in 1985 and the degree of habil. doctor of mathematics from the Institute of Mathematics and Informatics, Vilnius, in 1993. He is a professor and the head of Mathematical Modeling department of Vilnius Gediminas Technical University. His research interests include numerical methods for solving nonlinear PDE, parallel numerical methods, mathematical modelling in nonlinear optics, porous media flows, technology, image processing, biotechnology.

 ČiegisRemigijusremigijus.ciegis@knf.vu.lt2

R. Čiegis has graduated from Kaunas Polytechnical Institute Faculty of Chemical Engineering in 1982 and Vilnius University Kaunas Faculty, in 1989. He received the PhD degree in economics from the Vilnius University in 1995 and the degree of habil. doctor of management from the Kaunas Vytautas Magnus University in 2002. He is a professor of Vilnius University. His research interests include environmental economics, sustainable economic development, macroeconomics, regional economics.

 SubočOlgaos@vgtu.lt1∗

O. Suboč has graduated from Vilnius University Faculty of Mathematics, in 1998. She received the PhD degree from the Institute of Mathematics and Informatics and Vytautas Magnus University, in 2002. She is an associated professor in Vilnius Gediminas Technical University. Her research interests include numerical methods for solving nonlinear PDE, mathematical modelling, technology, image processing, non-local conditions, biotechno-logy.

 Vilnius Gediminas Technical University, Saulėtekio av. 11, LT-10223 Vilnius, Lithuania Vilnius University Kaunas Faculty, Muitinės st. 8, LT-44280 Kaunas, Lithuania Corresponding author. 2022245202233347749711202142022 © 2022 Vilnius University2022 Open access article under the CC BY license.

This model describes the heat equation in 3D domains, and this problem is reduced to a hybrid dimension problem, keeping the initial dimension only in some parts and reducing it to one-dimensional equation within the domains in some distance from the base regions. Such mathematical models are typical in industrial installations such as pipelines. Our aim is to add two additional improvements into this methodology. First, the economical ADI type finite volume scheme is constructed to solve the non-classical heat conduction problem. Special interface conditions are defined between 3D and 1D parts. It is proved that the ADI scheme is unconditionally stable. Second, the parallel factorization algorithm is proposed to solve the obtained systems of discrete equations. Due to both modifications the run-time of computations is reduced essentially. Results of computational experiments confirm the theoretical error analysis and scalability estimates of the parallel algorithm.

 hybrid dimension model ADI scheme parallel factorization algorithm heat conduction piplines Introduction

An important trend in development of mathematical models, simulating various real world applications, deals with new type nonlocal boundary and conjugation conditions. The non-locality mechanism helps to simulate important physical processes more accurately. Still the existence and uniqueness of the solution, its stability is very sensitive to approximation of such nonlocal boundary conditions and require development and analysis of special discrete approximation techniques and modifications of implementation algorithms. As one important example we mention the simulation of enzyme-catalysed glucose oxidation and redox reactions over inhomogeneous surfaces (Skakauskas et al., 2018; Čiegis et al., 2018).

It is a standard situation in applications of the mathematical modelling technique that systems of nonstationary 3D differential equations should be solved. Even if state of the art high-order discrete approximations are constructed on adaptive meshes, the total complexity of the obtained computational algorithms is very large. It restricts our possibility to include such models into general optimization routines or to use them in real time decision making applications.

Recently, a new important technique has been proposed to solve such applied problems more efficiently. A model reduction technique is applied to avoid the effects of curse of dimensionality. The given 3D heat conduction problem is reduced to a hybrid dimension problem, keeping the initial dimension only in some parts of the domain and reducing it to one-dimensional equation within the domains in some distance from the base regions. Such mathematical models are typical in industrial installations such as pipelines.

Our aim is to add two additional improvements into this methodology. First, the economical ADI type finite volume scheme is constructed to solve the non-classical heat conduction problem. It is proved that the ADI scheme is unconditionally stable. Second, the parallel factorization algorithm is proposed to solve the obtained systems of discrete equations. Due to both modifications the run-time of computations is reduced essentially.

The method of partial dimension reduction was first introduced in Panasenko (1998) and then developed in Panasenko (2005). It is based on the asymptotic analysis of a solution of a partial differential equation set in a thin domain and on a projection of the variational model formulation on a subspace of functions having a form of the asymptotic expansion in all zones of regular behaviour of the solution. Then the dimension of the problem is reduced within the big part of the domain keeping the full dimension description only in small zones of a singular behaviour of the solution.

Applying this approach a model with nonlocal junctions of one dimensional and full dimensional equations is obtained. Note that nonlocal conjugation conditions make the main challenges in development of effective parallel numerical methods for this class of problems of hybrid dimension. Let us mention some papers where different techniques were used to discretize the obtained hybrid mathematical model. A finite volume setting was studied in Panasenko and Viallon (2015). The method of asymptotic partial decomposition of the domain was compared with several existing methods of dimension reduction in Amar and Givoli (2018) and its effectiveness with respect to the precision and time of the execution is shown. In a recent paper (Čiegis et al., 2020) ADI type finite volume scheme was considered for a 2D heat conduction problem.

However, the question about the most effective parallel numerical solvers for such problems of hybrid dimension is still an important topic. Non-classical conjugation conditions require to develop special solvers to get the discrete solutions. A straightforward way is to solve the obtained systems of linear equations by using general iterative solvers, e.g. well-known multigrid solvers. For multidimensional parabolic problems an alternative and most popular way is to use splitting methods (Hundsdorfer and Verwer, 2003; Samarskii, 2001). It is well-known that in many cases ADI type methods are the best selection.

Thus in this paper we have developed and analysed special efficient ADI schemes to solve the hybrid 3D partial dimension reduction models for a nonstationary heat conduction problem. Our second goal is to propose efficient parallel implementations of the constructed ADI schemes. All main topics required in the analysis of discrete schemes, i.e. the approximation, positivity and stability of the solution are investigated. The accuracy of the reduced dimension model and the scalability of the parallel algorithms are tested in a series of computational experiments.

The approximation error is investigated quite straightforwardly, since the regularity of the solution and dependence of it on the reduction parameter δ is known from results given in Amosov and Panasenko (2018), Panasenko (2005). A more complicated step is to investigate the stability of proposed parallel ADI schemes in the case of new non-local conjugation conditions.

The rest of the paper is organized in the following way. In Section 2 the problem is formulated. The nonstationary heat conduction equation is formulated in the 3D parallelepiped and the initial and boundary conditions are specified. In this section also, the partially reduced dimension model is constructed and an approximate solution is defined by considering a simplified domain when 3D small size subregions are coupled with a 1D line. The non-local conjugation conditions are defined at the truncations of the tube. In Section 3 the full problem in the 3D domain is solved by using the ADI stability correction scheme. It is proved that this scheme is unconditionally stable. Its implementation requires solving many linear systems with tridiagonal matrix and the classical factorization algorithm is used. The Wang parallel factorization algorithm also can be used for parallel computers. In Section 4 the modified ADI stability correction scheme is constructed for the reduced dimension heat conduction problem. It is proved that the unique solution of the obtained system of linear equations exists for each time level. The efficient factorization algorithm is defined, it is based on algorithms presented in Čiegis et al. (2020). A parallel version of the modified factorization algorithm for implementation of the ADI scheme approximating the reduced dimension heat conduction problem is constructed in Section 5. A theoretical scalability analysis of this parallel algorithm is done. In Section 6 results of computational experiments are presented and analysed. Some final conclusions are done in Section 7.

 Problem Formulation

Let us consider a parallelepiped Ω ⊂ R 3 , which has the following form Ω = ( 0 , X 1 ) × D, where D is a rectangle D = { ( x 2 , x 3 ) : 0 < x j < X j , j = 2 , 3 }.

We are interested to solve a linear heat equation in Ω × ( 0 , T ]: ∂ u ∂ t = ∑ j = 1 3 ∂ 2 u ∂ x j 2 + f ( x 1 , x 2 , x 3 , t ) , ( x 1 , x 2 , x 3 , t ) ∈ Q T = Ω × ( 0 , T ] , u ( 0 , x 2 , x 3 , t ) = g 1 ( x 2 , x 3 , t ) , u ( X 1 , x 2 , x 3 , t ) = g 2 ( x 2 , x 3 , t ) , ( x 2 , x 3 , t ) ∈ D ¯ × ( 0 , T ] , ∂ u ∂ η = 0 , ( x 1 , x 2 , x 3 , t ) ∈ ∂ Ω 1 × ( 0 , T ] , u ( x 1 , x 2 , x 3 , 0 ) = u 0 ( x 1 , x 2 , x 3 ) , ( x 1 , x 2 , x 3 ) ∈ Ω , where ∂ Ω 1 = ∂ Ω ∖ ( D × { x 1 = 0 } ∪ { x 1 = X 1 } ) denotes a part of the boundary, where the Neumann boundary condition is specified.

Next, following Amosov and Panasenko (2018), we consider a modified approximate problem. First, let S ( u ) S ( u ) = 1 | D | ∫ 0 X 2 ∫ 0 X 3 u ( x 1 , x 2 , x 3 , t ) d x 2 d x 3 denote the averaging operator. We assume that the initial condition u 0 and source function f satisfy the relations u 0 ( x 1 , x 2 , x 3 ) = S ( u 0 ) , f ( x 1 , x 2 , x 3 , t ) = S ( f ) , ( x 1 , t ) ∈ ( 0 , X 1 ) × ( 0 , T ] . It means that u 0 and f don’t depend on x 2 , x 3 within Ω.

Denote a reduced dimension domain Ω δ = { ( x 1 , x 2 , x 3 ) ∈ ( δ , X 1 − δ ) × D }. Function U is called an approximate solution to problem (1)–(4) if it satisfies the following problem (Amosov and Panasenko, 2018) ∂ U ∂ t = ∑ j = 1 3 ∂ 2 U ∂ x j 2 + f ( x 1 , t ) , ( x 1 , x 2 , x 3 , t ) ∈ ( Ω ∖ Ω δ ) × ( 0 , T ] , ∂ U ∂ t = ∂ 2 U ∂ x 1 2 + f ( x 1 , t ) , ( x 1 , x 2 , x 3 , t ) ∈ Ω δ × ( 0 , T ] , U ( 0 , x 2 , x 3 , t ) = g 1 ( x 2 , x 3 , t ) , U ( X 1 , x 2 , x 3 , t ) = g 2 ( x 2 , x 3 , t ) , ( x 2 , x 3 , t ) ∈ D ¯ × ( 0 , T ] , ∂ U ∂ η = 0 , ( x 1 , x 2 , x 3 , t ) ∈ ∂ Ω δ , 1 × ( 0 , T ] , U ( x 1 , x 2 , x 3 , 0 ) = u 0 ( x 1 , x 2 , x 3 ) , ( x 1 , x 2 , x 3 ) ∈ Ω , where ∂ Ω δ , 1 = ∂ Ω δ ∖ ( D × { x 1 = δ } ∪ { x 1 = X 1 − δ } ).

In Ω δ × ( 0 , T ] the solution U doesn’t depend on x 2 , x 3 , i.e. U ( x 1 , x 2 , x 3 , t ) = S ( U ) , ( x 1 , x 2 , x 3 , t ) ∈ Ω δ × ( 0 , T ] . Thus in the reduced domain it is sufficient to find 1D space function U ( x 1 , t ).

By analysing the weak form of the heat equation, it is shown in Amosov and Panasenko (2018) that the following conjugation conditions are valid at the truncation points of the space domain U | x 1 = δ − 0 = U | x 1 = δ + 0 , U | x 1 = X 1 − δ − 0 = U | x 1 = X 1 − δ + 0 , ∂ S ( U ) ∂ x 1 | x 1 = δ − 0 = ∂ U ∂ x 1 | x 1 = δ + 0 , ∂ U ∂ x 1 | x 1 = X 1 − δ − 0 = ∂ S ( U ) ∂ x 1 | x 1 = X 1 − δ + 0 . The conditions (10) are classical and mean that U is continuous at the truncation points. While the remaining two conditions (11) are nonlocal and they define the conservation of full fluxes along the separation planes.

 ADI Scheme for the 3D Differential Problem

The uniform spatial mesh Ω ¯ h = ω ¯ 1 × ω ¯ 2 × ω ¯ 3 is defined as ω ¯ k = { x k j : x k j = j h k , j = 0 , … , J k } , x k J k = X k , k = 1 , 2 , 3 . For simplicity of notations we also consider a uniform time mesh: ω ¯ t = { t n : t n = n τ , n = 0 , … , N } , t N = T . Let U i j k n be a numerical approximation to the exact solution u ( x 1 i , x 2 j , x 3 k , t n ) of problem (1)–(4) at the grid point ( x 1 i , x 2 j , x 3 k , t n ).

The following standard operators are defined for discrete functions: ∂ x 1 U i j k n : = U i j k n − U i − 1 , j , k n h 1 , ∂ x 2 U i j k n : = U i j k n − U i , j − 1 , k n h 2 , ∂ x 3 U i j k n : = U i j k n − U i , j , k − 1 n h 3 , A 1 h U i j k n : = − 1 h 1 ( ∂ x 1 U i + 1 , j , k n − ∂ x 1 U i j k n ) , 0 < i < J 1 A 2 h U i j k n : = − 2 h 2 ∂ x 2 U i 1 k n , j = 0 , − 1 h 2 ( ∂ x 2 U i , j + 1 , k n − ∂ x 2 U i j k n ) , 0 < j < J 2 , 2 h 2 ∂ x 2 U i J 2 k n , j = J 2 , A 3 h U i j k n : = − 2 h 3 ∂ x 3 U i j 1 n , k = 0 , − 1 h 3 ( ∂ x 3 U i , j , k + 1 n − ∂ x 3 U i j k n ) , 0 < k < J 3 , 2 h 3 ∂ x 3 U i j J 3 n , k = J 3 . The discrete approximation of diffusion operators is obtained by using the finite volume method.

Then the heat conduction problem (1)–(4) is approximated by the following alternating direction implicit (ADI) stability correction scheme (Hundsdorfer and Verwer, 2003) U i j k n + 1 4 − U i j k n τ + ∑ l = 1 3 A l h U i j k n = f i j k n + 1 2 , ( x 1 i , x 2 j , x 3 k ) ∈ ω 1 × ω ¯ 2 × ω ¯ 3 , U i j k n + l 4 − U i j k n + l − 1 4 τ + 1 2 A 5 − l h ( U i j k n + l 4 − U i j k n ) = 0 , l = 2 , 3 , 4 . At the first step the heat conduction problem (1)–(4) is approximated by the explicit Euler scheme. It is well known that this scheme is only conditionally stable. The next three steps of scheme (12) guarantee the unconditional stability of the discrete scheme and they increase the accuracy of approximation in time until the second order.

At all three steps 3D systems are split into a set of 1D problems along x 3 , x 2 and x 1 coordinates respectively, and each subproblem is solved by using the classical fast factorization algorithm.

If a solution of the problem (1)–(4) is sufficiently smooth, then the approximation error of ADI scheme (12) is O ( τ 2 + h 1 2 + h 2 2 + h 3 2 ).

For a self-completeness of the text and in order to explain the approximation property of the ADI scheme (12) we present a short proof of this result (for a detailed analysis see Hundsdorfer and Verwer, 2003). By adding all equations (12) the following discrete equation is obtained U i j k n + 1 − U i j k n τ + ∑ l = 2 4 1 2 A 5 − l h ( U i j k n + l 4 + U i j k n ) = f i j k n + 1 2 . In order to compute the approximation error of (12), we use the relations U i j k n ∼ u ( x 1 i , x 2 j , x 3 k , t n ) , U i j k n + l 4 ∼ u ( x 1 i , x 2 j , x 3 k , t n + 1 ) , l = 2 , 3 , 4 . Thus the approximation error of ADI scheme (12) with respect to time coordinate is the same as of the classical Crank-Nicolson method. The second order accuracy of the constructed approximations of space derivatives is shown, e.g. in Hundsdorfer and Verwer (2003), Samarskii (2001).  □

It is a straightforward task to show the result presented in the following lemma.

The discrete operators A 1 h , A 2 h and A 3 h are symmetric, A 1 h is positive definite and A 2 h , A 3 h are non-negative operators.

All operators A 1 h , A 2 h and A 3 h commute, thus there exists a complete system of functions { ϕ 1 i ( x 1 ) ϕ 2 j ( x 2 ) ϕ 3 k ( x 3 ) , 1 ⩽ i ⩽ J 1 − 1 , 0 ⩽ j ⩽ J 2 , 0 ⩽ k ⩽ J 3 } , where ϕ l m ( x l ) are eigenvectors of the operator A l h , l = 1 , 2 , 3: A l h ϕ l m = λ l m ϕ l m , l = 1 , 2 , 3 and λ l m are eigenvalues λ 1 i ⩾ λ 11 > 0 , λ l j ⩾ 0 , l = 2 , 3 . 0,\hspace{2em}{\lambda _{lj}}\geqslant 0,\hspace{1em}l=2,3.\]]]>

ADI scheme (12) is unconditionally stable.

For the given normal operators A l h , l = 1 , 2 , 3 the Fourier stability analysis can be used. Let us consider the solution of ADI scheme (12) in the case when boundary conditions g j = 0, j = 1 , 2 and f ( x , t ) = 0. The solution of (12) can be written as U i j k n = ∑ l = 1 J 1 − 1 ∑ m = 0 J 2 ∑ r = 0 J 3 c l m r n ϕ 1 l ( x 1 i ) ϕ 2 m ( x 2 j ) ϕ 3 r ( x 3 k ) . Substituting this formula into equations (12), we obtain for each mode the stability equations c l m r n + 1 = q l m r c l m r n , where the discrete stability function q l m r is defined as q l m r = ( 1 − τ 2 ( λ 1 l + λ 2 m + λ 3 r ) + τ 2 4 ( λ 1 l λ 2 m + λ 1 l λ 3 r + λ 2 m λ 3 r ) + τ 3 8 λ 1 l λ 2 m λ 3 r ) ( 1 + 0.5 τ λ 1 l ) ( 1 + 0.5 τ λ 2 m ) ( 1 + 0.5 τ λ 3 r ) . Since eigenvalues satisfy estimates (14), then | q l m r | < 1 and the ADI scheme (12) is unconditionally stable in the L 2 norm.  □

It is easy to see that the discrete stability function approximates the continuous stability function of the differential equation (1) with the second order accuracy.

 The ADI Scheme for the Partially Dimension Reduced Problem

Let K 1 and K 2 define the truncation points of the domain x 1 K 1 = δ, x 1 K 2 = X 1 − δ. Then the spatial mesh ω 1 is splitted into three parts: ω 11 = { x 1 i : x 1 i = i h 1 , i = 1 , … , K 1 − 1 } , ω ¯ 12 = { x 1 i : x 1 i = i h 1 , i = K 1 , … , K 2 } , ω 13 = { x 1 i : x 1 i = i h 1 , i = K 2 + 1 , … , J 1 − 1 } . We note that more general meshes are used in Viallon (2013); Panasenko and Viallon (2015), where some atypical cells along the interface are constructed, and that lead to non-admissible meshes. In this paper we use the finite volume approach. For more general non-uniform meshes it can be recommended to use a general framework of finite element and finite difference schemes for construction of discrete approximations on non-matching grids (Eymard et al., 2000; Faille, 1992). A numerical comparison of these methods for 1D-2D discrete schemes is done in Viallon (2013).

The heat conduction problem (5)–(11) is approximated by the modified version of the ADI scheme U i j k n + 1 4 − U i j k n τ + ∑ l = 1 3 A l h U i j k n = f i j k n + 1 2 , ( x 1 i , x 2 j , x 3 k ) ∈ ( ω 11 ∪ ω 13 ) × ω ¯ 2 × ω ¯ 3 , U i 00 n + 1 4 − U i 00 n τ + A 1 h U i 00 n = f i 00 n + 1 2 , x 1 i ∈ ω ¯ 12 ∖ { x 1 K 1 , x 1 K 2 } , U K 1 00 n + 1 4 − U K 1 00 n τ + 1 h 1 2 ( − S h ( U K 1 − 1 n ) + 2 U K 1 00 n − U K 1 + 1 , 0 , 0 n ) = f K 1 00 n + 1 2 , U K 2 00 n + 1 4 − U K 2 00 n τ + 1 h 1 2 ( − S h ( U K 2 + 1 n ) + 2 U K 2 00 n − U K 2 − 1 , 0 , 0 n ) = f K 2 00 n + 1 2 , U i j k n + l 4 − U i j k n + l − 1 4 τ + 1 2 A 5 − l h ( U i j k n + l 4 − U i j k n ) = 0 , l = 2 , 3 , ( x 1 i , x 2 j , x 3 k ) ∈ ( ω 11 ∪ ω 13 ) × ω ¯ 2 × ω ¯ 3 , U i j k n + 1 − U i j k n + 3 4 τ + 1 2 A 1 h ( U i j k n + 1 − U i j k n ) = 0 , ( x 1 i , x 2 j , x 3 k ) ∈ ( ω 11 ∪ ω 13 ) × ω ¯ 2 × ω ¯ 3 , U i 00 n + 1 − U i 00 n + 3 4 τ + 1 2 A 1 h ( U i 00 n + 1 − U i 00 n ) = 0 , x 1 i ∈ ω ¯ 12 ∖ { x 1 K 1 , x 1 K 2 } , U K 1 00 n + 1 − U K 1 00 n + 3 4 τ + 1 2 h 1 2 ( − S h ( U K 1 − 1 n + 1 ) + 2 U K 1 00 n + 1 − U K 1 + 1 , 0 , 0 n + 1 + S h ( U K 1 − 1 n ) − 2 U K 1 00 n + U K 1 + 1 , 0 , 0 n ) = 0 , U K 2 00 n + 1 − U K 2 00 n + 3 4 τ + 1 2 h 1 2 ( − S h ( U K 2 + 1 n + 1 ) + 2 U K 2 00 n + 1 − U K 2 − 1 , 0 , 0 n + 1 + S h ( U K 2 + 1 n ) − 2 U K 2 00 n + U K 2 − 1 , 0 , 0 n ) = 0 . Here S h denotes the discrete averaging operator: S h ( U i n ) = 1 | D | ∑ j = 0 J 2 ∑ k = 0 J 3 d 2 j d 3 k U i j k n h 2 h 3 , where d l m = 1, if 0 < m < J l and d l m = 0.5, if m = 0 , J l , l = 2 , 3. Notation of the argument U i n indicates that the two dimensional average of discrete function U n is computed at the grid point x 1 i .

Note that equations (19), (20) approximate the nonlocal flux conjugation conditions (11).

We define two discrete meshes Ω h , R D = ( ω ¯ 2 × ω ¯ 3 × ( ω 11 ∪ ω 13 ) ) ∪ ω ¯ 12 and Ω ¯ h , R D = Ω h , R D ∪ ( ω ¯ 2 × ω ¯ 3 × ( x 10 ∪ x 1 J 1 ) ). Let us consider discrete functions U i j k = U ( x 1 i , x 2 j , x 3 k ) defined on the spatial mesh Ω ¯ h , R D . We denote the set of such vectors by D h .

Let us define three operators for U ∈ D h : A 2 h U = A 2 h U , ( x 1 i , x 2 j , x 3 i ) ∈ ( ω 11 ∪ ω 13 ) × ω ¯ 2 × ω ¯ 3 , 0 , x 1 i ∈ ω ¯ 12 , A 3 h U = A 3 h U , ( x 1 i , x 2 j , x 3 i ) ∈ ( ω 11 ∪ ω 13 ) × ω ¯ 2 × ω ¯ 3 , 0 , x 1 i ∈ ω ¯ 12 , A 1 h U = A 1 h U , ( x 1 i , x 2 j , x 3 i ) ∈ ( ω 11 ∪ ω 13 ) × ω ¯ 2 × ω ¯ 3 , A 1 h U i 00 , x 1 i ∈ ω ¯ 12 ∖ { x 1 K 1 , x 1 K 2 } , 1 h 1 2 ( − S h ( U K 1 − 1 ) + 2 U K 1 00 − U K 1 + 1 , 0 , 0 ) , x 1 = K 1 , 1 h 1 2 ( − S h ( U K 2 + 1 ) + 2 U K 2 00 − U K 2 − 1 , 0 , 0 ) , x 1 = K 2 .

Then we can write the ADI scheme as U n + 1 4 − U n τ + ( A 1 h + A 2 h + A 3 h ) U n = f n + 1 2 , U n + l 4 − U n + l − 1 4 τ + 1 2 A 5 − l h ( U n + l 4 − U n ) = 0 , l = 2 , 3 , 4 .

Next, our goal is to prove that the discrete operators A 1 h , A 2 h and A 3 h are symmetric.

For U , V ∈ D h , such that U 0 j k = U J 1 j k = 0, V 0 j k = V J 1 j k = 0, ( x 2 , x 3 ) ∈ ω ¯ 2 × ω ¯ 3 the formulas ( U , V ) = ∑ j = 0 J 2 ∑ k = 0 J 3 d 2 j d 3 k ( ∑ i = 1 K 1 − 1 U i j k V j k + ∑ i = K 2 + 1 K − 1 U i j k V i j k ) h 1 h 2 h 3 + | D | ∑ i = K 1 K 2 U i 00 V i 00 h 1 , ‖ U ‖ = ( U , U ) 1 / 2 define a scalar product and a norm in this vector space. We remind that d l m = 1, if 0 < m < J l and d l m = 0.5, if m = 0 , J l , l = 2 , 3.

The discrete operators A 1 h , A 2 h and A 3 h are symmetric, A 1 h is positive definite and A 2 h , A 3 h are non-negative operators.

First, we investigate the operator A 2 h . Applying the summation by part formula, we get ( A 2 h U , V ) = ∑ j = 0 J 2 ∑ k = 0 J 3 d 2 j d 3 k ( ∑ i = 1 K 1 − 1 A 2 h U i j k V i j k + ∑ i = K 2 + 1 K − 1 A 2 h U i j k V i j k ) h 1 h 2 h 3 + ∑ j = 1 J 2 ∑ k = 0 J 3 d 3 k ( ∑ k = 1 K 1 − 1 ∂ x 2 U i j k ∂ x 2 V i j k + ∑ k = K 2 + 1 K − 1 ∂ x 2 U i j k ∂ x 2 V i j k ) h 1 h 2 h 3 = ( U , A 2 h V ) , ( A 2 h U , U ) ⩾ 0 . In a similar way we get the estimates ( A 3 h U , V ) = ∑ j = 0 J 2 ∑ k = 0 J 3 d 2 j d 3 k ( ∑ i = 1 K 1 − 1 A 3 h U i j k V i j k + ∑ i = K 2 + 1 K − 1 A 3 h U i j k V i j k ) h 1 h 2 h 3 + ∑ j = 0 J 2 ∑ k = 1 J 3 d 2 k ( ∑ k = 1 K 1 − 1 ∂ x 3 U i j k ∂ x 3 V i j k + ∑ k = K 2 + 1 K − 1 ∂ x 3 U i j k ∂ x 3 V i j k ) h 1 h 2 h 3 = ( U , A 3 h V ) , ( A 3 h U , U ) ⩾ 0 . It follows from the obtained estimates that A 2 h and A 3 h are symmetric and non-negative definite operators.

Next, we investigate the operator A 1 h . Applying the summation by parts formula and taking into account the boundary conditions for vectors U, V, and the nonlocal conjugation conditions (19), (20) we get ( A 1 h U , V ) = ∑ j = 0 J 2 ∑ k = 0 J 3 d 2 j d 3 k ( ∑ i = 1 K 1 − 1 A 1 h U i j k V i j k + ∑ i = K 2 + 1 K − 1 A 1 h U i j k V i j k ) h 1 h 2 h 3 + | D | ( ∑ i = K 1 + 1 K 2 − 1 A 1 h U i 00 V i 00 h 1 + 1 h 1 ( − S h ( U K 1 − 1 ) + 2 U K 1 00 − U K 1 + 1 , 0 , 0 ) V K 1 00 + 1 h 1 ( − S h ( U K 2 + 1 ) + 2 U K 2 00 − U K 2 − 1 , 0 , 0 ) V K 2 00 ) = ∑ j = 0 J 2 ∑ k = 0 J 3 d 2 j d 3 k ( ∑ i = 1 K 1 ∂ x 1 U i j k ∂ x 1 V i j k + ∑ i = K 2 + 1 K ∂ x 1 U i j k ∂ x 1 V i j k ) h 1 h 2 h 3 + | D | ∑ i = K 1 + 1 K 2 ∂ x 1 U i 00 ∂ x 1 V i 00 h 1 = ( U , A 1 h V ) , ( U , A 1 h U ) > 0 . 0.\end{aligned}\]]]> Here we use the notation U K 1 j k = U K 1 00 , U K 2 j k = U K 2 00 for ( x 2 j , x 3 k ) ∈ ω ¯ 2 × ω ¯ 3 .

It follows from the obtained estimates that A 1 h is a symmetric and positive definite operator.  □

Due to nonlocal conjugation conditions (19), (20) the classical factorization algorithm should be modified in order to solve 1D subproblems (17)–(20).

The unique solution of the linear system of equations (15)–(20) exists and it can be computed by using the efficient factorization algorithm.

The given proof also defines the constructive algorithm to solve (17)–(20). Since the first part (15) of the discrete ADI scheme defines an explicit algorithm and 1D subproblems (16) in the second part of the scheme are efficiently solved by using the classical factorization algorithm, then it is sufficient to consider in detail the subproblem (17)–(20). For each ( x 2 j , x 3 k ) ∈ ω ¯ 2 × ω ¯ 3 the solution is factorized separately in subdomains ω 11 , ω 12 and ω 13 . Let us write equations of the system (17)–(20) as − a i j k U i − 1 , j , k n + 1 + c i j k U i j k n + 1 − b i j k U i + 1 , j , k n + 1 = d i j k , a i j k , b i j k , c i j k ⩾ 0 , c i j k ⩾ a i j k + b i j k .

1. Domain ω 11 . The solution is presented in the following form: U i j k n + 1 = α i j k U i + 1 , j , k n + 1 + γ i j k , 0 ⩽ i < K 1 , α 0 j k = 0 , γ 0 j k = g 1 ( x 2 j , x 3 k , t n + 1 ) , α i j k = b i j k c i j k − a i j k α i − 1 , j , k , γ i j k = d i j k + a i j k γ i − 1 , j , k c i j k − a i j k α i − 1 , j , k . By induction it can be proved that the estimates 0 ⩽ α i j k ⩽ 1 are valid.

2. Domain ω 12 . The solution is presented in the following form: U i 00 n + 1 = α i 00 U K 1 00 n + 1 + β i 00 U K 2 00 n + 1 + γ i 00 , K 1 < i < K 2 . This factorization is done in two steps. First, the solution is written in the form U i 00 n + 1 = α ˜ i 00 U K 1 00 n + 1 + β ˜ i 00 U i + 1 , 0 , 0 n + 1 + γ ˜ i 00 , K 1 < i < K 2 , α ˜ K 1 + 1 , 0 , 0 = a K 1 + 1 , 0 , 0 c K 1 + 1 , 0 , 0 , β ˜ K 1 + 1 , 0 , 0 = b K 1 + 1 , 0 , 0 c K 1 + 1 , 0 , 0 , γ ˜ K 1 + 1 , 0 , 0 = d K 1 + 1 , 0 , 0 c K 1 + 1 , 0 , 0 , α ˜ i 00 = a i 00 c i 00 − a i 00 β ˜ i − 1 , 0 , 0 α ˜ i − 1 , 0 , 0 , β ˜ i 00 = b i 00 c i 00 − a i 00 β ˜ i − 1 , 0 , 0 , γ ˜ 0 k = a i 00 γ ˜ i − 1 , 0 , 0 + d i 00 c i 00 − a i 00 β ˜ i − 1 , 0 , 0 . Let us assume that 0 ⩽ α ˜ i − 1 , 0 , 0 , β ˜ i − 1 , 0 , 0 ⩽ 1 and α ˜ k − 1 , 0 , 0 + β ˜ k − 1 , 0 , 0 ⩽ 1. Then it follows from the given formulas that 0 ⩽ β ˜ i 00 ⩽ 1 and α ˜ i 00 ⩾ 0. It remains to prove that α ˜ i 00 ⩽ 1 and α ˜ i 00 + β ˜ i 00 ⩽ 1. It follows from simple inequalities a i 00 α ˜ i − 1 , 0 , 0 + b i 00 + a i 00 β ˜ i − 1 , 0 , 0 = a i 00 ( α ˜ i − 1 , 0 , 0 + β ˜ i − 1 , 0 , 0 ) + b i 00 ⩽ a i 00 + b i 00 ⩽ c i 00 , that a i 00 α ˜ i − 1 , 00 + b i 00 ⩽ c i 00 − a i 00 β ˜ i − 1 , 0 , 0 . Thus we get the estimate 0 ⩽ α ˜ i 00 + β ˜ i 00 ⩽ 1 . The proof by induction is completed.

Next, we compute coefficients α i 00 , β i 00 and γ i 00 in formula (25) α K 2 − 1 , 0 , 0 = α ˜ K 2 − 1 , 0 , 0 , β K 2 − 1 , 0 , 0 = β ˜ K 2 − 1 , 0 , 0 , α i 00 = α ˜ i 00 + β ˜ i 00 α i + 1 , 0 , 0 , β i 00 = β ˜ i 00 β i + 1 , 0 , 0 , γ i 00 = γ ˜ i 00 + β ˜ i 00 γ i + 1 , 0 , 0 , i = K 2 − 2 , … , K 1 + 1 . From estimates for α ˜ 0 k , β ˜ 0 k and (29) it follows that 0 ⩽ α i 00 , β i 00 ⩽ 1 , 0 ⩽ α i 00 + β i 00 ⩽ 1 .

3. Domain ω 13 . For each ( x 2 j , x 3 k ) ∈ ω 2 × ω 3 the solution is presented in the following form: U i j k n + 1 = β i j k U i − 1 , j , k n + 1 + γ i j k , K 2 < i ⩽ J 1 , β J 1 j k = 0 , γ J 1 j k = g 2 ( x 2 j , x 3 k , t n + 1 ) , β i j k = a i j k c i j k − b i j k β i + 1 , j , k , γ i j k = d i j k + b i j k γ i + 1 , j , k c i j k − b i j k β i + 1 , j , k . By induction it can be proved that the estimates 0 ⩽ β i j k ⩽ 1 are valid.

Substituting (23)–(30) into equations (19) and (20) we get a linear system of two equations to find U K 1 00 n + 1 , U K 2 00 n + 1 : A 11 U K 1 00 n + 1 + A 12 U K 2 00 n + 1 = B 1 , A 21 U K 1 00 n + 1 + A 22 U K 2 00 n + 1 = B 2 , where coefficients are defined by A 11 = 1 τ + 1 2 h 1 2 ( 2 − α K 1 + 1 , 0 , 0 − S h ( α K 1 − 1 ) ) , A 12 = − 1 2 h 1 2 β K 1 + 1 , 0 , 0 , A 21 = − 1 2 h 1 2 α K 2 − 1 , 0 , 0 , A 22 = 1 τ + 1 2 h 1 2 ( 2 − β K 2 − 1 , 0 , 0 − S h ( β K 2 + 1 ) ) , B 1 = U K 1 00 n + 3 4 τ − 1 2 h 1 2 ( S h ( U K 1 − 1 n ) − 2 U K 1 00 n + U K 1 + 1 , 0 , 0 n + γ K 1 + 1 , 0 , 0 − S h ( γ K 1 − 1 ) ) , B 2 = U K 2 00 n + 3 4 τ − 1 2 h 1 2 ( S h ( U K 2 + 1 n ) − 2 U K 2 00 n + U K 2 − 1 , 0 , 0 n + γ K 2 − 1 , 0 , 0 − S h ( γ K 2 + 1 ) ) . From the estimates given above we have that 2 − α K 1 + 1 , 0 , 0 − S h ( α K 1 − 1 ) ⩾ β K 1 + 1 , 0 , 0 , 2 − β K 2 − 1 , 0 , 0 − S h ( β K 2 + 1 ) ⩾ α K 2 − 1 , 0 , 0 , thus the determinant of the matrix of system (31) is positive and the unique solution U K 1 00 n + 1 , U K 2 00 n + 1 exists. Then the backward factorization step is applied and U i j k n + 1 are computed.  □

It is noted in Hundsdorfer and Verwer (2003) that the stability analysis of such ADI schemes for three dimensional cases is far from easy and even the linear commuting case needs a closer attention. Since in our case operators A 1 h , A 2 h , A 3 h don’t commute, we can’t apply the spectral stability analysis.

Here we restrict to writing the stability matrix R = ( I + τ 2 A 1 h ) − 1 ( I + τ 2 A 2 h ) − 1 ( I + τ 2 A 3 h ) − 1 ( I − τ 2 ( A 1 h + A 2 h + A 3 h ) ) . It follows from Lemma 4, that ‖ I + τ 2 A j h ‖ ⩽ 1 , j = 1 , 2 , 3 . Then it is sufficient to assume that the discrete problem is such that the stability estimate ‖ R n ‖ ⩽ C , n ⩾ 1 is satisfied. The validity of this assumption should be tested for each case of the discrete problem experimentally.

 Parallel Algorithm

In this section a parallel version of the ADI algorithm (15)–(20) is presented. We note that development of efficient parallel algorithms for implementation of ADI type solvers is quite a non-trivial task. First the factorization algorithm for 3D problems severely depends on data movement in the memory of computers (from and to in different levels of the memory). Only few arithmetical operations are done with each small size portion of data, therefore cost of data movement is very important issue (see, Guo and Lu, 2016; Otero et al., 2020). Second, nonlocality of some parts of the discrete scheme lead to modifications of the standard factorization algorithm and these changes should be resolved efficiently by the parallel algorithm (see Čiegis et al., 2014, where a parallel ADI algorithm is constructed to solve multidimensional parabolic problems with nonlocal boundary conditions). Third, different approaches of parallelization can be used for different architectures of parallel computers. At least three main classes can be considered, including GPU processors (Imankulov et al., 2021; Xue and Feng, 2018; Otero et al., 2020), shared memory processors (multicore processors) and general global memory parallel machines, when implementation of parallel algorithms can be done, e.g. by using MPI library (Čiegis et al., 2014). In this paper we restrict to the construction of parallel algorithms for general distributed memory computers. It is obvious that these algorithms can be efficiently used also for multicore computers and distributed clusters of multicore nodes.

First, in defining a parallel version of the ADI algorithm for the reduced dimension model we restrict to a a simple modification of the standard sequential factorization algorithm. Only two parts are used to decompose the domain Ω h in each dimension. Therefore the maximum number of parallel processes is restricted to eight processes.

We start by defining a parallel algorithm for the solution of a tridiagonal system of linear equations: c 0 U 0 − b 0 U 1 = F 0 , − a j U j − 1 + c j U j − b j U j + 1 = F j , j = 1 , … , J − 1 , − a J U J − 1 + c J U J = F J . Let us assume that two processes are used to solve it in parallel. System (33) is partitioned among both processes. The first process is responsible for the first ( J + 1 ) / 2 equations 0 ⩽ j < ( J + 1 ) / 2 and the second process solves the remaining equations. Both processes implement the classical sequential factorization algorithm in parallel but in different directions. The first process applies this algorithm from left-to-right: U j = α j U j + 1 + β j , j = 0 , … , ( J + 1 ) / 2 − 1 while the second process computes in the opposite direction: U j = α j U j − 1 + β j , j = J , … , ( J + 1 ) / 2 . Then both processes exchange coefficients ( α J + 1 2 − 1 , β J + 1 2 − 1 ), ( α J + 1 2 , β J + 1 2 ), solve 2 × 2 dimension systems and find solutions U J + 1 2 − 1 , U J + 1 2 . The second backward part of the factorization algorithm is again implemented in parallel.

Now we are ready to describe the main decomposition details of the parallel version of ADI algorithm (15)–(20).

1. The case of p = 2 processes. The discrete grid ω ¯ 1 is splitted into two parts ω ¯ 1 = ω 1 , p 0 ∪ ω 1 , p 1 , where ω 1 , p 0 = ω 11 ∪ ω ¯ 12 , ω 1 , p 1 = ω 13 . The first process computes the solution on the discrete grid ( ω 11 × ω ¯ 2 × ω ¯ 3 ) ∪ ω ¯ 12 , and the second process computes the solution on the grid ω 13 × ω ¯ 2 × ω ¯ 3 .

Before computations of the explicit step (21), the first process sends the value U K 2 00 n to the second process, and the second process computes the value of S ( U K 2 + 1 n ) and sends it to the first process. Then both processes compute their parts of U n + 1 4 in parallel.

The l = 2 , 3 steps of the ADI scheme (22) don’t require any communication between processes and can be done in parallel.

The last step l = 4 of the ADI scheme (22) is implemented in the following way. The first process computes the factorization coefficients of (23) and (25), and the second process computes coefficients of (30). These computations can be done in parallel. Then the second process computes the values of S ( β K 2 + 1 n ), S ( γ K 2 + 1 n ) and sends them to the first process. The latter solves system (31) and sends the value U K 2 00 n + 1 to the second process. Then both process compute their parts of U n + 1 in parallel.

A simple but still accurate complexity model of the parallel ADI algorithm can be constructed T 2 = c J 1 J 2 J 3 2 + 2 ( α + 2 β ) , where T 2 is computation time of the one time step of the parallel algorithm, α is data communication start-up time, β defines time required to send one element.

We note that this model takes into account only main computational and data communication costs. In real computations the values of these parameters can vary. One more important detail, i.e. automatic memory caching is also not estimated in this formula. Thus efficiency estimates of such models are only asymptotically valuable, but exactly such information is sufficient in most cases.

Let us define the speed-up S p and efficiency E p of the parallel algorithm S p = T 0 T p , E p = S p p , where p is the number of processes, T 0 is the computation time of the sequential algorithm.

Since T 0 = c J 1 J 2 J 3 and additional data communication costs are constant, then even for moderate size grids ω 1 × ω 2 × ω 3 it follows from (34) that S 2 ≈ 2 , E 2 ≈ 1 , i.e. the parallel algorithm is scalable.

2. The case of p = 4 processes. In addition to splitting the discrete grid ω ¯ 1 , the grid ω ¯ 2 is also divided into two parts ω ¯ 2 = ω 2 , p 0 ∪ ω 2 , p 1 , where ω 2 , p 0 = { x 2 j , 0 ⩽ j ⩽ J 2 2 } , ω 2 , p 1 = { x 2 j , J 2 2 + 1 ⩽ j ⩽ J 2 } . The first process computes the solution on the grid ( ω 11 × ω ¯ 2 , p 0 × ω ¯ 3 ) ∪ ω ¯ 12 , the second process computes the solution on the grid ω 13 × ω ¯ 2 , p 0 × ω ¯ 3 , the third process computes the solution on the grid ( ω 11 × ω ¯ 2 , p 1 × ω ¯ 3 ) ∪ ω ¯ 12 , and the fourth process computes the solution on the grid ω 13 × ω ¯ 2 , p 1 × ω ¯ 3 .

Before parallel computations of the explicit step (21), the first process sends the value U K 2 00 n to the second and fourth process, these two processes compute their part of S ( U K 2 + 1 n ) and send the obtained values to the first process. In a similar way the first process sends the value U K 1 00 n to the third process, which computes its part of S ( U K 1 + 1 n ) and sends the obtained value to the first process.

The second part of data communication is done among two pairs of processes: the first and third processes exchange the values of U i , J 2 / 2 , k n and U i , J 2 / 2 + 1 , k n , 0 < i < K 1 , 0 ⩽ k ⩽ J 3 , and the second and fourth processes exchange the values of U i , J 2 / 2 , k n and U i , J 2 / 2 + 1 , k n , K 2 < i < J 1 , 0 ⩽ k ⩽ J 3 . These communications are done in parallel. Then all four processes compute their parts of U n + 1 4 in parallel.

The l = 2 step of the ADI scheme (22) don’t require any communication between processes and can be done in parallel.

The l = 3 step of the ADI scheme (22) is implemented by using the parallel factorization algorithm, which is described above for the system of linear equations (33). The data communication part between the given pairs of the processes is equivalent to the one defined for the explicit part of the ADI algorithm (21). This step again can be done in parallel.

The last step l = 4 of the ADI scheme (22) is implemented as for p = 2 processes with one small modification. The second and fourth process compute the values of S ( β K 2 + 1 n ), S ( γ K 2 + 1 n ) and the third process computes the values of S ( α K 1 − 1 n ), S ( γ K 1 − 1 n ). Then they send these values to the first process. It solves system (31) and sends the values U K 1 00 n + 1 , U K 2 00 n + 1 to the remaining three process. Then all four process compute their parts of U n + 1 in parallel.

The complexity model of the parallel ADI algorithm for p = 4 processes is given by T 4 = c J 1 J 2 J 3 4 + 4 α + J 1 J 3 2 β .

It follows from this theoretical model, that now the data communication costs depend linearly on the size of J 1 J 3 . Thus only for sufficiently large J 2 these costs can be neglected and estimates S 4 ≈ 4 , E 4 ≈ 1 are valid. We get that the parallel algorithm is scalable for sufficiently large J 2 and this specific size depends on both parameters – computational complexity c and data communication costs β.

3. The case of p = 8 processes. In addition to splitting discrete grids ω ¯ 1 , ω ¯ 2 the discrete grid ω ¯ 3 is also divided into two parts ω ¯ 3 = ω 3 , p 0 ∪ ω 3 , p 1 , where ω 3 , p 0 = { x 3 k , 0 ⩽ k ⩽ J 3 2 } , ω 3 , p 1 = { x 3 k , J 3 2 + 1 ⩽ k ⩽ J 3 } . The parallel algorithm is modified to include data distribution steps in the x 3 dimension, all communications are defined by using analogous formulas given for p = 4 case in the x 2 dimension.

The complexity model of the parallel ADI algorithm for p = 8 processes is given by T 8 = c J 1 J 2 J 3 8 + 4 α + ( J 1 J 3 4 + J 1 J 2 4 ) β . It follows from this model that for sufficiently large grid sizes, when J 1 J 2 J 3 ≫ J 1 ( J 2 + J 3 ), the estimates S 8 ≈ 8 , E 8 ≈ 1 are valid and the parallel algorithm is scalable.

Let us consider the case when the number of processors/cores is larger than 8. In order to use all p processors, we modify the constructed parallel algorithm. We assume that p is an even number and factorize it as p = 2 × p 2 × p 3 . The decomposition of the discrete grid Ω h is done in two steps. First, the ω ¯ 1 grid is divided into two parts following the algorithm defined for p ⩽ 8 processes. Thus the implementation of the parallel algorithm in x 1 dimension remains unchanged. Next, we divide the grids ω ¯ j into p j , j = 2 , 3 parts. If p j > 22$]]>, then in this dimension the proposed parallel factorization algorithm is changed to a more expensive but general algorithm. We recommend to use Wang’s or Cyclic reduction parallel factorization algorithms (Kumar et al., 1994; Imankulov et al., 2021).

 Computational Experiments

Consider the problem (1)–(4) in the domain Ω with geometry parameters X 1 = 6, X 2 = 1, X 3 = 1. We solve this problem till T = 1 with given initial, boundary conditions and the zero source function u 0 ( x 1 , x 2 , x 3 ) = 0 , g 1 ( x 2 , x 3 , t ) = 5 t e − 10 ( ( x 2 − 0.5 ) 2 + ( x 3 − 0.5 ) 2 ) , g 2 ( x 2 , x 3 , t ) = 2 t e − 6 ( ( x 2 − 0.5 ) 2 + ( x 3 − 0.5 ) 2 ) , f ( x 1 , x 2 , x 3 , t ) = 0 .

The accuracy of time integration of ADI schemes. First, we have tested the time integration accuracy of the ADI solver for the full dimension model. A uniform space grid Ω h with J 1 = 120, J 2 = 20, J 3 = 20 is used. Table 1 gives for a sequence of decreasing time step widths τ the errors e ( τ ) and the experimental convergence rates ρ ( τ ) of the discrete solution for ADI scheme (12) in the maximum norm: e ( τ ) = max ( x 1 i , x 2 j , x 3 k ) ∈ Ω h | U i j k N − U ( x 1 i , x 2 j , x 3 k , T ) | , ρ ( τ ) = log 2 ( e ( 2 τ ) / e ( τ ) ) , where the reference solution U ( x 1 i , x 2 j , x 3 k , T ) is computed by using the very small time step τ = 6.25 × 10 − 5 . Then the integration error introduced by the ADI scheme can be measured accurately.

Errors e ( τ ) and experimental convergence rates ρ ( τ ) for the discrete solution of ADI scheme (12) for a sequence of time steps τ.

τ e ( τ ) ρ ( τ )
0.002 3.0352 × 10 − 5 2.016
0.001 7.5272 × 10 − 6 2.012
0.0005 1.8549 × 10 − 6 2.021
0.00025 4.411 × 10 − 7 2.072

It follows from the presented results, that the accuracy of ADI scheme agrees well with the theoretical prediction.

The accuracy of the reduced dimension model. Next we investigate the accuracy of the reduced dimension model (5)–(11). For the space discretization a uniform grid Ω h with J 1 = 300, J 2 = 50, J 3 = 50 is used, and integration in time is done with τ = 0.001.

The numerical tests have been performed on the computer with Intel® Xeon® processor E5-2670, 8GB RAM, and the algorithms have been implemented using C++ language.

We note that CPU time for computing the full model solution by using the ADI scheme (12) is 55.9 seconds.

Table 2 gives for a sequence of reduction parameter δ errors e ( δ ) e ( δ ) = max ( x 1 i , x 2 j , x 3 k ) ∈ Ω h h | U i j k N − U i j k N ( δ ) | of the reduced dimension model (5)–(11) discrete solution in the maximum norm, where U i j k N is the solution of the ADI scheme (12) and U i j k N ( δ ) is the solution of the reduced dimension ADI scheme (15)–(20). CPU times are also presented for the scheme (15)–(20).

Errors e ( δ ) of the discrete solution of the reduced dimension model (5)–(11) for a sequence of truncation parameter δ.

δ = 1 δ = 1.5 δ = 2 δ = 2.8
e ( δ ) 4.304 × 10 − 3 1.748 × 10 − 4 7.201 × 10 − 6 4.522 × 10 − 8
CPU 19.4 29.1 40.8 60.2

It follows from the presented results, that starting from δ = 1 the accuracy of the reduced dimension model is sufficient for most real world applications, while CPU time for the reduced dimension model is reduced up to 3 times.

The efficiency of the parallel ADI scheme (15)–(20). In this section we present results of the parallel scalability tests. All parallel numerical tests in this work were performed on the computer cluster “HPC Vanagas” at the High Performance Computing Laboratory of Vilnius Gediminas technical university. We have used up to 8 cores of Intel® Xeon® processor E5-2630 with 20 cores (2.20 GHz) and 32 GB DDR4 of RAM. The parallel algorithms have been implemented using MPI library.

Our goal is to investigate the efficiency of the proposed parallel version of the ADI scheme (15)–(20). We have solved two problems with different sizes of the discrete meshes 450 × 75 × 75 and 600 × 100 × 100. The CPU time T p , speed-up S p and efficiency E p coefficients are presented in Table 3.

Scalability analysis of the parallel algorithm. The speed-up S p and efficiency E p coefficients for two problems of dimension 450 × 75 × 75 and 600 × 100 × 100.

T p ( 450 ) S p ( 450 ) E p ( 450 ) T p ( 600 ) S p ( 600 ) E p ( 600 )
p = 1 77.4 – – 254.1 – –
p = 4 17.5 4.42 1.12 57.79 4.40 1.10
p = 8 9.85 7.86 0.97 29.2 8.70 1.08

Two conclusions follow from the presented results of computational experiments. First, the constructed parallel algorithm scales well even for problems of moderate sizes.

Second, MPI library enables efficient usage of the allocated local memory for each process, thus the efficiency of the parallel algorithm E p over 1 is achieved.

As a future task, it is interesting to consider the scalability of the modified parallel algorithm from Remark 1, when the number of processes is larger than 8.

 Conclusions

In this paper we have applied the new efficient approach how to solve the heat equation in 3D domains. This problem is reduced to a hybrid dimension problem, keeping the initial dimension only in some parts and reducing it to one-dimensional equation within the domains in some distance from the base regions. Such mathematical models are typical in industrial installations such as pipelines. We added two additional improvements into this methodology. First, the economical ADI type finite volume scheme is constructed to solve the non-classical heat conduction problem. A finite volume method is used to approximate space differential operators with non classical conjugation conditions between the 3D and 1D parts. It is proved that the ADI scheme is unconditionally stable.

Second, the parallel factorization algorithm is proposed to solve the obtained systems of discrete equations. The ADI scheme leads to non-iterative implementation algorithm and a set of one dimensional linear systems are solved by using the parallel versions of the factorization algorithm for the parallel solving tridiagonal equations systems in each direction. An efficient modification of the basic factorization algorithm is developed to resolve non-local conjugation conditions.

The presented results of numerical experiments confirm both main theoretical conclusions. Due to both modifications the run-time of computations is reduced essentially. First, the hybrid reduced dimension models can be used to simulate heat conduction models for a broad set of domains and CPU time of the serial algorithm is reduced up to factor 3. Second, the developed parallel version of the ADI algorithm enables to achieve effective acceleration of computations and it scales well even for discrete problems of moderate sizes.

 References Amar, H., Givoli, D. (2018). Mixed-dimensional modeling of time-dependent wave problems using the Panasenko construction. Journal of Theoretical and Computational Acoustics, 26(3), 1850034. https://doi.org/10.1142/S2591728518500342. Amosov, A., Panasenko, G. (2018). Partial dimension reduction for the heat equation in a domain containing thin tubes. Mathematical Methods in the Applied Sciences, 41(18), 9529–9545. https://doi.org/10.1002/mma.5311. Čiegis, R., Katauskis, P., Skakauskas, V. (2018). Numerical methods for fractional diffusion. Nonlinear Analysis: Modelling and Control, 23(4), 234–250. https://doi.org/10.15388/NA.2018.2.6. Čiegis, R., Suboč, O., Bugajev, A. (2014). Parallel numerical algorithms for three-dimensional parabolic and pseudoparabolic problems with different boundary conditions. Nonlinear Analysis: Modelling and Control, 19(3), 382–395. https://doi.org/10.15388/NA.2014.3.5. Čiegis, R., Panasenko, G., Pileckas, K., Šumskas, V. (2020). ADI scheme for partially dimension reduced heat conduction models. Computers & Mathematics with Applications, 80(5), 1275–1286. https://doi.org/10.1016/j.camwa.2020.06.012. Eymard, R., Gallouët, T., Herbin, R. (2000). Finite volume methods. Handbook of Numerical Analysis, 7, 713–1018. https://doi.org/10.15388/NA.2018.2.6. Faille, I. (1992). A control volume method to solve an elliptic equation on a two-dimensional irregular mesh. Computer Methods in Applied Mechanics and Engineering, 100(2), 275–290. https://doi.org/10.1016/0045-7825(92)90186-N. Guo, G., Lu, S. (2016). Unconditional stability of alternating difference schemes with intrinsic parallelism for two-dimensional fourth-order diffusion equation. Computers & Mathematics with Applications, 71, 1944–1959. Hundsdorfer, W., Verwer, J. (2003). Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations. Springer Series in Computational Mathematics, Vol. 33. Springer, Berlin, Heidelberg, New York, Tokyo. https://doi.org/10.1007/978-3-662-09017-6. Imankulov, T., Daribayev, B., Mukhambetzhanov, S. (2021). Comparative analysis of parallel algorithms for solving oil recovery problem using CUDA and OpenCL. International Journal of Nonlinear Analysis and Applications, 12(1), 351–364. Kumar, V., Grama, A., Gupta, A., Karypis, G. (1994). Introduction to Parallel Computing: Design and Analysis of Algorithms. Benjamin/Cummings, Redwood City, CA. Otero, B., Rojas, O., Moya, F., Castillo, J.E. (2020). Alternating direction implicit time integrations for finite difference acoustic wave propagation: parallelization and convergence. Computers and Fluids, 205, 104584. Panasenko, G. (1998). Method of asymptotic partial decomposition of domain. Mathematical Models and Methods in Applied Sciences, 8(1), 139–156. https://doi.org/10.1142/S021820259800007X. Panasenko, G. (2005). Multi-scale Modeling for Structures and Composites. Springer Netherlands, Dordrecht. https://doi.org/10.1007/1-4020-2982-9. Panasenko, G., Viallon, M.C. (2015). Finite volume implementation of the method of asymptotic partial domain decomposition for the heat equation on a thin structure. Russian Journal of Mathematical Physics, 22(2), 237–263. https://doi.org/10.1134/S1061920815020107. Samarskii, A.A. (2001). The Theory of Difference Schemes. Marcel Dekker, New York. https://doi.org/10.1007/1-4020-2982-9. Skakauskas, V., Katauskis, P., Čiegis, R. (2018). Modelling of the NO + CO reaction over inhomogeneous surfaces. Journal of Mathematical Chemistry, 56(9), 2626–2642. https://doi.org/10.1007/s10910-018-0908-3. Viallon, M.C. (2013). Error estimate for a 1D–2D finite volume scheme. Comparison with a standard scheme on a 2D non-admissible mesh. Comptes Rendus Mathematique, 351(1), 47–51. https://doi.org/10.1016/j.crma.2013.01.011. Xue, G., Feng, H. (2018). A new parallel algorithm for solving parabolic equations. Advances in Difference Equations. 174. https://doi.org/10.1186/s13662-018-1617-8.