Two-Dimensional-One-Dimensional Alternating Direction Schemes for Coastal Systems Convection-Diffusion Problems

The initial boundary value problem for the 3D convection-diffusion equation corresponding to the mathematical model of suspended matter transport in coastal marine systems and extended shallow water bodies is considered. Convective and diffusive transport operators in horizontal and vertical directions for this type of problem have significantly different physical and spectral properties. In connection with the above, a two-dimensional–one-dimensional splitting scheme has been built—a three-dimensional analog of the Peaceman–Rachford alternating direction scheme, which is suitable for the operational suspension spread prediction in coastal systems. The paper has proved the theorem of stability solution with respect to the initial data and functions of the right side, in the case of time-independent operators in special energy norms determined by one of the splitting scheme operators. The accuracy has been investigated, which, as in the case of the Peaceman–Rachford scheme, with the special definition of boundary conditions on a fractional time step, is the value of the second order in dependency of time and spatial steps. The use of this approach makes it possible to obtain parallel algorithms for solving grid convection-diffusion equations which are economical in the sense of total time of problem solution on multiprocessor systems, which includes time for arithmetic operations realization and the one required to carry of information exchange between processors.


Introduction
In this article the splitting scheme for the 3D convection-diffusion equation has been considered as an analog of the Peaceman-Rachford alternating direction scheme [1][2][3]. This scheme assumes splitting the original problem into a set of two-dimensional convectiondiffusion problems along horizontal directions and one-dimensional problems along vertical coordinates. Despite great attention to the possibilities of using such schemes, at present, we do not have rigorous mathematical results on their stability and accuracy for convection-diffusion problems. Until now, conclusions about the stability, accuracy and performance of these schemes have been drawn not on the theoretical study basis, but on the basis of the analysis of the computational experiments on the approximate solution of some model problems. In this work, an attempt has been made to fill this gap.
The choice made by the authors in favor of these splitting schemes is based on the following physical motivation for convection-diffusion problems in coastal marine systems and shallow water bodies [4][5][6].
Turbulent diffusion coefficients in horizontal directions Ox 1 , Ox 2 : k 1 → x , k 2 → x very weakly depend on spatial coordinates and they can be considered constant, assuming that: The influence of the terms of convective transfer along the horizontal coordinates, as a rule, significantly exceeds the contribution of the terms describing turbulent exchange in the same directions to the formation of the distributions of suspended matter or solutes [7][8][9]. In the vertical coordinate direction, turbulent exchange in presence of wind conditions in coastal systems largely determines the concentration of suspended matter. The turbulent exchange coefficient in the vertical direction k 3 → x significantly depends on → x (x 1 , x 2 , x 3 ), especially from x 3 , where x 3 is the vertical coordinate. The typical values of the velocities of water medium movement in the vertical direction differ significantly. It has been changed from 10 −4 m/s to 10 −2 m/s. When the problem is discretized, additional asymmetry arises due to the fact that along the horizontal directions the typical number of grid nodes is 10 3 − 10 4 , when in the vertical direction this number is several tens, no more than 10 2 . The foregoing motivates the two-component splitting of the three-dimensional convectiondiffusion operator, which at each time step assumes the splitting of the three-dimensional problem into two-dimensional in horizontal directions and one-dimensional in the vertical direction, which at each time step may be symbolically represented as: In previous expression C (α) and D (α) denote operators of convective and diffusion transportations respectively in directions of Ox α axes, α = 1, 2, 3. In this article, the theorem has been proved for the stability of different solutions with respect to the initial data and functions of the right side in the case of time-independent operators in a special energy norm determined by one of the operators of the splitting scheme. The accuracy has been investigated, which, as in the case of the Peaceman-Rachford scheme [10,11], with the special setting of boundary conditions on a fractional time step, is the value of the second order of accuracy for the time and spatial steps.

Problem Statement
Let us consider the formulation of the initial-boundary value problem for the spatial three-dimensional convection-diffusion equation in a parallelepiped G [12,13]: Applied to the coastal hydrodynamics problems function w → x , t can be the concentration of pollutants.
Assume that there is a classical equation w → x , t of the initial-boundary value problem (2)-(4), satisfying the accepted smoothness conditions [14], as well as the input data: the ∂G and also the velocity vector of the medium To construct and study the two-dimensional-one-dimensional Peasman-Reckford scheme, it is necessary to strengthen the requirements for the smoothness of the solution function mentioned in [15]. Namely, assume that G × [0 ≤ t ≤ T] there are bounded derivatives of 4-5 orders of the solution function, that is, the following inequalities hold: where M α , α = 1, 2, . . . , 6 are positive constants independent of x α , α = 1, 2, 3 and t.
Operator A is the sum of three one-dimensional convection-diffusion operators: where C (α) 0 is the operator of convective transfer along the coordinate directions in symmetric form: Note also that each of the one-dimensional convection-diffusion operators presented in Formula (5) is positive defined, because: Difference analogs of these differential operators will inherit this property.

Alternating Direction Schemes Construction
In order to simplify the investigation, assume that the operators A (α) , α = 1, 2, 3 do not depend on the time variable.
With respect to the right-hand side of Equation (2), assume: The splitting form (8) of the function f → x , t , responsible for the source function is dictated by the physical features of two-dimensional and one-dimensional problems. With x , t the smoothness conditions are preserved as Let's build a temporary grid: ω τ = {t n = (n + α/2)τ, α ∈ {0, 1}; n = 0, 1, . . . , N t ; N t τ ≡ τ}, at each step of which consider the chain of interconnected initial-boundary value problems: With boundary conditions, which will be specified below.
Let's build a uniform spatial mesh ω h = ω 1 × ω 2 × ω 3 , where: For each of the grids ω α , α = 1, 2, 3 select the set of boundary nodes γ α = ω α \ω α , where ω α contains the inner nodes of the grid in the direction Ox α , On the grid ω τ × ω h the chain of problems was approximated using the balance relations method, bearing in mind the analog of the Peaceman-Rachford alternating direction schemes, which has the accuracy O |h| 2 + τ 2 , |h| 2 h are different analogs of the corresponding differential operators and y n+α/2 , α = 1, 2 are difference analogs of the functions w (α) problems (9)−(13): where the function µ → x , t n+1/2 is defined similarly to how it is done for the Peaceman-Rachford scheme in the case of a two-dimensional problem.
It can be shown that: where: It should be noted that problem (14)−(17) can be written in the operator form: where the use of the initial condition and boundary conditions are not separately detailed for brevity; E is the unit (identity) operator. The numerical implementation of difference equations (19) can be performed by an adaptive alternating triangular iterative method for grid equations with a non-self-adjoint operator [17].
For the numerical implementation of the adaptive alternating-triangular iterative method, an essential constraint is imposed on the grid Peclet number Pe h , which is formulated for the operators A (1) h and A (2) h : This limitation is met with a margin for discrete models of coastal marine systems and extended shallow water bodies, which can be seen from the estimates: Difference Equation (20) are numerically implemented by the Thomas algorithm for three-point difference equations [18].

Investigating Scheme Stability
Let us formulate a theorem about the stability of scheme (14)-(18) with respect to the initial condition and the right-hand side function.
where H h -Hilbert space of grid functions. Then, for the different solution of problem (14)-(18) on an arbitrary(n + 1)-in the n-th time step, the following inequality takes place-an estimate of stability independent of initial data and right-hand side functions: Proof of Theorem 1. Let us first prove the auxiliary inequality Consider first the scalar product on the left-hand side of Inequality (22): Similarly, the scalar product on the right-hand side of relation (22) is transformed to the form: Comparing S 1 , S 2 and considering that: Arrive at the inequality: which means the validity of Inequality (22).
Let's go on to the main proof. From Equality (19) we have the estimate: On the other hand, from Equality (20) the inequality is obtained: Inequalities (25) and (26) imply the estimate: Use Inequality (22) h , v ≡ y n and get: Substituting estimate (28) into the right-hand side of Inequality (27), we arrive at an Inequality in the form of: Since this Inequality (29) is valid for an arbitrary n, n = 0, 1, . . . , N t − 1, then we have a chain of inequalities: Inequalities (29) and (30) imply the inequality: The condition y 0 = u 0 , x ∈ ω h leads to Inequality (21), as required.

Remark 1.
Formal splitting of three-dimensional operator A = 3 ∑ α=1 A (α) the original differential problem for three one-dimensional convection-diffusion operators and the construction on this basis of a scheme of alternating directions, which can be symbolically represented as leads, as a rule, to a computationally unstable scheme.

Remark 2.
The transition to additive schemes of total approximation, which can be symbolically represented by relation (32), allows obtaining stable schemes. However, in the case of parallel implementation on multiprocessor systems, after each fractional time step (solving the three-point difference problem), it is necessary to organize exchanges of blocks of three-dimensional grid information. This process is similar to the parallel implementation of the transposition of a three-dimensional matrix based on the united operational memories of processors (and/or cores) before the solution of individual three-point problems and is an extremely time-consuming procedure [19].

Investigation of the Two-Dimensional-One-Dimensional Alternating Direction Schemes Accuracy
The study of the scheme accuracy (14)- (18) involves the estimation of errors at integer and fractional time steps.
The errors of solutions need the estimation, the approximation errors of the difference equations (14) and (16) were estimated. At the same time, it has been a motivation to set the boundary condition (18). It is clear that these procedures are of a standard manner. However, in order not to refer the reader to not too accessible sources, we present these investigations here. Substitute instead of functions y n , y n+1/2 , y n+1 respectively functions: We get: Add up these equalities term by term, we arrive at the equalities: w n+1 − w n 0.5τ In what follows, we will assume f n ≡ 0.5 f n . Subtracting Equality (34) from Equality (33), we have: In fact, the choice w n+1/2 is, in a sense, the best, because in this case the right side and ψ n,τ equal to zero. In addition, expression (36) remains correct on the lateral boundary Γ ver , Γ ver ≡ ω 3 × (γ 1 ∪ γ 2 ) of a parallelepiped G, and therefore, the boundary condition has the form (18). Now it is easy to estimate the errors z n ≡ y n − w n иz n+1/2 ≡ y n+1/2 − w n+1/2 . The tasks for them are written in an obvious way: where: Using Equality (36) from Equation (39), we arrive at the following estimate: Let us take into account that when discretizing diffusion and convective terms written in symmetric form using central differences, their approximation error with respect to grid steps is the value O |h| 2 , that is A Let us complete the estimation of the accuracy (errors) in some special (energy) norm · , operator-defined A h > 0. Note that problems (37), (38) in structure coincide with problems (19), (20), for which, in connection with the study of stability, estimate (21) was obtained. Let us take into account only that z 0 = 0 (the initial condition has been defined «exactly») and we obtain the estimation, bearing in mind Inequality (41): Remark 3. If, for some reason, it is more preferable to obtain an estimate of the stability and accuracy in the energy space norm, defined operator E + τ 2 A , then in the construction of the difference scheme (14)-(18) the order of using the operators A h and A (3) h in different schemes should be reversed. Moreover, all the previous reasoning will be similar.

Discussion
A wide class of heat and mass transfer problems in liquid and gaseous media are described by convection-diffusion equations. The two-dimensional-one-dimensional scheme of alternating directions considered in the article takes into account the peculiarities of convection and diffusion processes, in particular, suspended particles and salts in coastal marine systems and can be used to predict the spread of pollutants in shallow seas like the Azov and White Seas of Russia and in other water bodies.
For the numerical solution of multidimensional problems of mathematical physics, the splitting method or the method of fractional steps has become widespread since the fifties of the XX century [20][21][22][23]. Currently, the list of works related to this topic includes a lot of thousands of units and continues to grow [24][25][26][27][28].
An important factor in motivating the development of this method is the widespread use of multiprocessor systems and parallel computations, the introduction of which led to the need to expand and, in some sense, rethink the concept of economical algorithms. In most cases, the user is interested in minimizing the total time expenditure of solving the problem, the main part of which is the cost of performing arithmetic, logical operations and information exchange operations between the processor units. It seems expedient to consider such a parallel algorithm as economical, which, among others, has the minimum total time expenditures. It turned out that the one-dimensional splitting schemes, which imply the replacement of the original multidimensional problem by a chain of one-dimensional problems, being economical in terms of time spent on performing computations, in some important cases are not economical in sense of time for performing information exchanges between processors.
In this article, using the example of a three-dimensional convection-diffusion problem, one of the families of splitting schemes is considered, which implies replacing the original problem with a chain of two-dimensional and one-dimensional problems. When implemented on multiprocessor systems with fixed connections between processorsfor example-line type of communications between neighboring processor units, then two-dimensional-one dimensional schemes, in this case, will require less time to perform information exchanges in comparison to one-dimensional splitting schemes. The time expenditure at each time step in the case of using a one-dimensional splitting scheme for information exchange is proportional to the value (N 1 · N 2 · N 3 · p) · t e , where p-is the number processors, t e -is the time required for one number (machine word) transfer between neighboring processor units. For the two-dimensional -one dimensional splitting schemes this estimation is proportional to value (N 1 · N 2 · N 3 /p) · t e · N it (ε), where N it (ε) is the number of iterations for the numerical realization of two-dimensional difference convection-diffusion tasks with the required precision ε. For a large enough number of p (for massively parallel supercomputers) in the case of using advanced iteration methods, the proposed two-dimensional-one-dimensional splitting scheme is more preferable in comparison to the one-dimensional one. Detailed analysis is required for special investigation and may be the subject of another article.

Conclusions
For the numerical realization of the initial-boundary value problem corresponding to the 3D model of suspended matter convection-diffusion in an aqueous medium, the authors proposed an approach based on two-dimensional-one-dimensional splitting schemes taking into account the specific features of convection-diffusion processes in the coastal system. An analog of such schemes is the Peaceman-Rachford alternating direction scheme. The use of two-dimensional-one-dimensional schemes gives the possibility to obtain low-cost algorithms in sense of total time for numerical realization on massively parallel computer systems. In the norm of the Hilbert space of grid functions, an estimate of stability is obtained with respect to the initial data and the right-hand side function. The study of accuracy (errors) was carried out for the proposed difference scheme in a certain special (energy) norm determined by the operator A (3) h and has been estimated the rate of convergence, which is the value of O |h| 2 + τ 2 .
The study of the problem of choosing a parallel algorithm that is economical in terms of total time costs is a rather complicated and time-consuming procedure and deserves a separate article, which the authors are currently working on.
Finally, the following text has been added: System study of this problem requires taking into account the following factors: -The nodes number of the spatial grid; -Parallel computing system architecture, including p the number of calculators (processors), the structure of information transmission channels between processors and their performance, in particular, the characteristic time of the exchange of the number t exc and characteristic latency time t lat , execution time of one arithmetic operation t ar the method used for solving grid equations and its parallel analog, etc.
To solve two-dimensional grid equations in a two-dimensional-one-dimensional scheme, the authors previously used the adaptive iterative alternating-triangular method developed by them for problems with a non-self-adjoint operator. The articles [17,29] are devoted to the study of its characteristics.
It was found that for the number of computers from 16 to 128, this algorithm has an efficiency of more than 50% when using a distributed memory system installed at the Southern Federal University, Russia (SFedU).
To solve one-dimensional three-point problems in the case of two-dimensional-onedimensional schemes, the Thomas algorithm is used and no additional information exchanges between the processors are required for this. A completely different situation develops if one-dimensional splitting schemes are used, which, after each fractional time step, lead to the need to transpose a set of two-dimensional arrays based on the union of processor memories-to rotate them relative to the coordinate axes.
The authors for two multiprocessor systems: SFedU and Sirius University of Science and Technology with the characteristics are indicated in Table 1. Experiments were carried out to solve the model problem of convection-diffusion on a grid N x 1 × N x 2 × N x 3 , where N x 1 = N x 2 = 1024, N x 3 = 128 for a different number p processors p = 8, 32, 128.
The time required to perform arithmetic operations t ar and exchange operations t exc information between processors when switching to a new integer time layer, depending on the number of processors, is indicated in Tables 2 and 3.  A preliminary conclusion can be drawn that for p ≥ 32, N ≥ 10 8 , where N = N 1 · N 2 · N 3 are one-dimensional splitting schemes on these multiprocessor systems have low efficiency, since the overwhelming time expenditure-more than 99.5%-is the time spent on exchanging information between processors, and here two-dimensional onedimensional schemes are preferable (best) in terms of total time costs.