Paper covers the research of discrete analogs of convective and di usion transfer operators of the fourth order of accuracy in the case of partial cell occupancy. According to the comparison of calculation results of substance transport problem based on schemes of the second and fourth orders of accuracy, the accuracy was increased in 66.7 times for di usion problem, and in 48.7 times for di usion-convection problem. A library of two-layer iterative methods was designed for solving two-dimensional di usion-convection problem based on schemes of high order of accuracy. It has intended to solve the nine-diagonal grid equations on a multiprocessor computer system. A mathematical algorithm was designed and numerically implemented for restoration the water salinity eld based on hydrographic information (water salinity at separate points or level isolines). The map of salinity of the Azov Sea was obtained using the proposed solution method.
One of the main problems of computational mathematics is the problem of solving systems of linear algebraic equations. Direct and iterative methods are used to obtain an approximate solution of systems of equations. One of the most successful method among the two-layer iterative methods is the alternately triangular method (ATM) proposed by A.A. Samarsky [Sam89]. Later, the academician A.N. Konovalov developed an adaptive version of ATM [Kon02]. The technique for increasing the convergence rate of ATM with a priori information by re ning the spectral estimates of the preconditioned operator are presented in [Suk84].
Often in applied problems, for example, in mathematical modeling of hydrodynamics [Suk11, Ka17, Rue05], heat and mass transfer [Suk18, Suk18'], geo ltration, population dynamics [Suk05], seismic exploration [Mur13] and other processes, it is necessary to solve the equations of convection-di usion type. In the case of implicit schemes and schemes with weights, such problems lead to linear algebraic equations with a non-self-adjoint operator. One of approaches to solving such problems is the Gaussian symmetry method. The disadvantage of this method is the squared increase of the condition number of operator, which leads to a decrease the convergence rate of iterative methods for solving grid equations. This fact contributed to the creation by the author's team the version of the modi ed iterative alternating triangular method of minimum corrections for solving grid equations with non-self-adjoint operator [Suk12'].
The use of rapidly converging iterative method as well as the use of parallal computations [Che06, Iak05] and the choice of di erence schemes are e ective ways to reduce the running time of the algorithm.
To increase the time step, we can use schemes with the optimal value of weight parameter [Suk14]. In addition, we can use the splitting computational grids, but it leads to increasing the calculation time. To increase the accuracy of calculations, it is possible to use schemes of higher order of accuracy [Pet13, Lad09] and schemes that take into account the lling of cells [Suk15]. In the second case, the accuracy is increased due to a better approximation of the boundary of computational domain.
Within the framework of this research, a library of iterative methods was designed to solve grid equations with self-adjoint and non-self-adjoint operators, arising in the solution of applied problems, schemes of high order of accuracy, taking into account the fullness of cells on a multiprocessor computing system. 2
The substance transport problem can be represented by the di usion-convection-reaction equation: with boundary conditions: c0t + uc0x + vc0y = ( c0x)0x +
c0y 0y + f c0n (x; y; t) = nc + n; where u,v are components of the velocity vector; is the turbulent exchange coe cient; f is a function, describing the intensity and distribution of sources.
We introduced a uniform grid for numerical implementation of the discrete mathematical model: wh = ftn = n ; xi = ihx; yj = jhy; n = 0::Nt; i = 0::Nx; j = 0::Ny; Nt = T; Nxhx = lx; Nyhy = lyg ; where is a time step; hx, hyare spatial steps; Nt is an upper boundary on time; Nx, Ny are space boundaries.
Discrete analogues of convective uc0x and di usive ( c0x)0xtransfer operators of the second order of approximation error in the case of partially lled cells can be written as:
(q0)i; j uc0x ' (q1)i; j ui+1=2; j (q0)i; j ( c0x)0x ' (q1)i; j i+1=2; j (q1)i; j ci; j + (q2)i; j ui 1=2; j
ci; j (q2)i; j i 1=2;j
ci; j xci; j + x ; hx
ci 1; j ; 2hx
ci 1; j h2 x where qi are coe cients, describing the fullness of control domain [Suk15]. 3
Expressions (1)-(2) can be considered in the case if (q1)i; j = (q2)i; j = 1. To increase the approximation order of equations (1)-(2), it's necessary to research the following di erence schemes: - the discrete analogue of the convective transport operator in absence of in uence of domain boundary: uc0x ' ui+1=2; j ci; j + ui 1=2; j ci; j
ci 1; j ; 2hx - the discrete analogue of the di usive transport operator in absence of in uence of domain boundary: The approximation error of expression (3) will take the following form:
(q1)i) ki(1) + ((q2)i + (q1)i) ki(2) ci+ (q1)i 2 + ui+1=2
2h (q2)i (q0)i ui 1=2 (q2)i 1 (q2)i 12h (q0)i 1 ci 2; where ki(1) = ((qq10))ii (ui+1 ui;) ((qq20))ii (ui
The approximation error of expression (4) will take the following form: ui 1) = (8h), ki(2) = ((qq10))ii ui+1 ui + ((qq20))ii ui ui 1 .
8h 8h + ( i; j)00 (ci; j)00 h42x + ( i; j)0 (ci; j)000 h62x + ( i; j)000 (ci; j)0 h62x + O h4x : Therefore, for approximation the di usive transport operator ( c0)0 by di erence scheme of the fourth order c(IV )h2=12 00c00h2=4 0c000h2=6 000c0h2=6 by of accuracy we have to approximate the operator ( c0)0 the scheme of the second order of accuracy.
The di usive transport operator ( c0)0 by di erence scheme of the fourth order of accuracy can be written as: (q0)i L (c) = Aici + B1;ici+1 + B2;ici 1 + B3;ici+2 + B4;ici 2: (6) B1;i = (q1)i ih+21=2 + (q1)i 12h2 i+1 (q1)i + 2 + (q2)i 12h2 (q0)i
i 1 (q1)i (q0)i i+1 (q2)i + (q2)i 12h2 B2;i = (q2)i ih21=2 + (q1)i 12h2 (q0)i i 1 (q2)i + 2 (q0)i (q1)i ki(3) (q2)i ki(3) (q1)i (q2)i 0i0+1 00 i 12 12 00 i ; 0i0 1 ; B3;i =
i+1 (q1)i+1 ; B4;i = (q1)i 12h2 (q0)i+1 (4) (5) i 1 + (q2)i 12h2 12 00 i 1 (q1)i + i + (q2)i 12h2 (q0)i i+1 (q2)i + (q1)i 12h2 (q0)i The eld, describing the error of calculations obtained as the di erence between the analytical and numerical solution of substance transport problem, is given in Fig. 1. The initial distribution was determined by the function:
C (x; y) = sin ( (x 10)) cos ( (y 0; fx; yg 2= D;
10)) ; fx; yg 2 D; D : fx 2 [10; 20] ; y 2 [10; 20]g :
Simulation was performed on the grid by dimension of 100x100 computational nodes. Simulation parameters: the dimensions of the computational domain lx=100 m, ly=100 m, and the time step is ht=0.001 s; the time period is 100 s; the horizontal component is 4 m/s, vertical { 3 m/s; the coe cient of turbulent exchange is 2 m2/s.
According to the comparison of results of numerical experiments based on schemes of the second and fourth orders of accuracy (see Fig. 1), the accuracy was increased in 66.7 times for solution the di usion problem, and in 48.7 times { for solution the di usion problem-convection. 5
A library of two-layer iterative methods for solution the nine-diagonal grid equations was designed for solution the two-dimensional di usion-convection problem based on the schemes of high order of accuracy. This library for solution the systems of linear algebraic equations (SLAE) include the following methods: the Jacobi method; the method of minimal corrections; the method of steepest descent; the Seidel method; the method of upper relaxation; the adaptive MATM of variation type.
Dependences of the number of iterations, required to solve the model problem on the time variable step, are given in the Table 1.
The idea of parallel algorithm of iterative methods with preconditioners of triangular type [Suk12'] (Zeidel method, upper relaxation method, alternative triangular method) on a system with distributed memory is as follows: at the rst step, each processor receives a subdomain, obtained by partition of the source domain into parts in one or more coordinate directions with an intersection of two nodes in each direction. Then, the SLAE solution with the upper-triangular operator is carried out, as a result of which the vector of solutions is calculated at the next iteration. The order of traversal of grid nodes in calculations and data exchanges in the case of decomposition in one spatial direction are shown in Fig. 3 and denoted by arrows. On the next step the residual vector and its uniform norm (the maximum modulo element) is calculated. In this case, each processor determines the maximum modulo element of the residual vector and transfers its value to all other processors. After data exchanges, processors calculate the maximum element in which the norm of the residual vector will be stored. If the norm of the residual vector is greater than the speci ed error, then the return to the calculation of the residual is performed.
At calculation the value of the solution vector, only the rst processor does not require additional information and can process its part of the region independently of other calculators; other processors are waiting for data transfer from the previous one.
Data transfer for one element is not optimal, because there're time costs associated with the organization of transfers. It can be minimized by increasing the size of data package; but it increases the delay time of the start of processors. Thus, the problem of calculation (selection) the optimal amount of transferred data package occur.
Values of acceleration and e ciency of parallel implementation of the software, designed to solve the twodimensional di usion-convection problem on the basis of high order accuracy schemes, are given in the Table 2. The grid equations were solved by the modi ed alternating-triangular method. The computational grid consist of 2000x2000 nodes. Parallel implementation of the developed algorithms was based on Message Passing Interface (MPI) technologies. The peak performance of the multiprocessor computer system (MCS) is 18.8 TFlops. As computing nodes 128 one-type16-core HP ProLiant BL685c Blade-servers were used, each of which is equipped with four 4-core processors AMD Opteron 8356 2.3 GHz and 32GB RAM.
According to the table 2, the parallel algorithm of the modi ed alternating-triangular method can be applied to solve real problems, and the use of parallel technologies makes a signi cant contribution to reduce the calculation time. 6
One of the urgent problems that arise at mathematical modeling of hydrodynamics of shallow waters [Suk18'] is the problem of hydrographic information processing. Typically, the salinity is speci ed at separate points or level isolines (see Fig. 4).
Using these maps for construction the computational grids is undesirable because of the error of calculations related to the "coarse" setting geometry of computational domain. Thus, for increasing the accuracy of calculations of hydrodynamic processes, it is necessary to approximate the function of two variables describing the salinity eld by more stable functions.
Formulation the problem of calculation the salinity eld. To determine the salinity function, we use the di usion equation solution to which the Saint-Venant equation describing the transport of bottom materials is reduced [Sid17]. The solution of the di usion problem for a long time intervals is reduced to the solution of the Laplace equation: (7) where H is a water salinity.
H = 0;
This approach has a signi cant disadvantage due to the lack of smoothness at points where the salinity eld values are speci ed. To resolve this problem, we can use the following equation:
The disadvantages of this approach include large outliers (deviation from the linear function). With the rst two approaches, we can get functions that do not have a direction, but each approach has disadvantages.
To determine a smooth salinity function, we can also apply the equation solution used to obtain schemes of the high order of accuracy for the Laplace equation:
2H = 0: H
In the rst case, the interpolation is performed by segments of lines passing through neighboring points; in the second case, the interpolation is based on cubic splines; in the third case { on function splines (12). The algorithm for one-dimensional interpolation based on the function (12) is described below, and the proposed approaches are compared.
Results of salinity eld restoration. The proposed mathematical algorithm for determine the water salinity eld was numerically implemented. The salinity isolines were obtained using the recognition algorithm (Fig. 5a) The salinity eld was obtained using the describing above interpolation algorithm in a rectangle (Fig. 5b). The map of salinity of the Azov Sea was obtained by applying the boundaries of the region (Fig. 6).
Note that the proposed algorithm has a su cient degree of smoothness at points of gluing functions and lower emissions compared to the cubic function used in the calculations.
Schemes of high (fourth) order of accuracy for the convective and di usive transfer operators, taking into account the lling of the cells, were constructed. A library of two-layer iterative methods was designed and implemented on MCS for solution the two-dimensional di usion-convection problem based on the schemes of high order of accuracy.
The comparison of calculation results of substance transport problem on the basis of schemes of the second and fourth orders of accuracy was performed. According to the comparison of results of numerical experiments, the accuracy was increased in 66.7 times for solution the di usion problem, and in 48.7 times { for solution the di usion problem-convection. The algorithms description of parallel implementation of iterative methods with preconditioners of triangular type and value of acceleration and e ciency of parallel variant of algorithm of the modi ed alternative triangular method is given. A mathematical algorithm was proposed to restore the water salinity eld on the basis of hydrographic information (water salinity at separate points or level isolines), and its numerical implementation was performed. The map of the salinity of the Azov Sea was obtained and based on the proposed method for solving the problem. The developed algorithm has a su cient degree of smoothness at points of gluing functions and lower emissions in the one-dimensional case compared to the cubic function used in calculations. Note that the proposed schemes were also used for development a software package designed to calculate the three-dimensional velocity ow elds in shallow waters [Suk11]. In the future, the developed schemes will be software implemented for calculation the biological kinetic problems [Gus18] and transport of bottom materials [Sid17]. 7.0.1
Acknowledgements This study was supported in part by task No. 2.6905.2017/BP within the basic part of the state task of the Ministry of Education and Science. [Ka17]
Ka C.O., Dembele J.-M., Cambier C., Stinckwich S., Lo M., Zucker J.-D. Deterministic convectiondi usion approach for modeling cell motion and spatial organization: Experi-mentation on avascular tumor growth // IEEE International Conference on Bioinformatics and Biomedicine. 2017. P. 556-560. [Rue05] Ruether N., Singh J.M., Olsen N.R.B., Atkinson E. 3-D computation of sediment transport at water intakes // Proceedings of the Institution of Civil Engineers: Water Management. 2005. Vol. 158 (1). P. 1-8. [Suk18] Sukhinov A.I., Belova Yu.V., Chistyakov A.E. The di erence scheme for the two-dimensional convection-di usion problem for large peclet numbers // MATEC Web of Conferences. 2018. Vol. 226.
No. 04030. [Suk18'] Sukhinov A.I., Chistyakov A.E., Nikitina A.V., Belova Y.V., Sumbaev V.V., Semenyakina A.A. Supercomputer modeling of hydrochemical condition of shallow waters in summer taking into account the in uence of the environment // Communications in Computer and Information Science. 2018. Vol. 910.
P. 336-351. [Suk05] Sukhinov, A.I., Sukhinov A.A. Reconstruction of 2001 ecological disaster in the Azov sea on the basis of precise hydrophysics models. Parallel Computational Fluid Dynamics 2004: Multidisciplinary Applications. 2005. P. 231-238. [Mur13] Muratov M.V., Petrov I.B. Calculation of wave responses from systems of subvertical macrofractures using the grid-characteristic method // Matem. Mod. 2013. Vol. 25. No. 3. P. 89-104. [Suk12'] Sukhinov A.I., Chistyakov A.E. Adaptive modi ed alternating triangular iterative method for solving grid equations with a non-self-adjoint operator// Mathematical Models and Computer Simulations. 2012. Vol. 4. No. 4. P. 398-409. [Che06] Chetverushkin B., Gasilov V., Iakobovskij M., Polyakov S., Kartasheva E., Boldarev A., Abalakin I., Minkin A. Unstructured mesh processing in parallel CFD project GIMM // Parallel Computational Fluid Dynamics 2005. 2006. P. 501-508.
Iakobovskij M.V. Incremental algorithm of graphs decomposition // Vestnik of Lobachevsky University of Nizhni Novgorod. Seria: Mathematical modeling and optimal control. 2005. No. 1. P. 243. [Suk14] Sukhinov A.I., Chistyakov A.E., Shishenya A.V. Error estimate for di usion equations solved by schemes with weights // Mathematical Models and Computer Simulations. 2014. 6(3). P. 324-331.
Sidoryakina V.V., Sukhinov A.I. Well-posedness analysis and numerical implementation of a linearized two-dimensional bottom sediment transport problem // Computational Mathematics and Mathematical Physics. 2017. 57 (6). P. 978-994.