Parallel implementation of substance transport problems for restoration the salinity field based on schemes of high order of accuracy Alexander I. Sukhinov, Yulia V. Belova, Don State Technical Don State Technical University, Rostov-on-Don, Russia University, Rostov-on-Don, Russia sukhinov@gmail.com yuliapershina@mail.ru Alena A. Filina, Supercomputers and Neurocomputers Research Center, Taganrog, Russia j.a.s.s.y@mail.ru Abstract Paper covers the research of discrete analogs of convective and diffu- sion transfer operators of the fourth order of accuracy in the case of partial cell occupancy. According to the comparison of calculation re- sults of substance transport problem based on schemes of the second and fourth orders of accuracy, the accuracy was increased in 66.7 times for diffusion problem, and in 48.7 times for diffusion-convection prob- lem. A library of two-layer iterative methods was designed for solving two-dimensional diffusion-convection problem based on schemes of high order of accuracy. It has intended to solve the nine-diagonal grid equa- tions on a multiprocessor computer system. A mathematical algorithm was designed and numerically implemented for restoration the water salinity field based on hydrographic information (water salinity at sep- arate points or level isolines). The map of salinity of the Azov Sea was obtained using the proposed solution method. 1 Introduction 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 refining the spectral estimates of the preconditioned operator are presented in [Suk84]. Copyright 2019 for this paper by its authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0). In: S. Hölldobler, A. Malikov (eds.): Proceedings of the YSIP-3 Workshop, Stavropol and Arkhyz, Russian Federation, 17-09-2019–20-09-2019, published at http://ceur-ws.org 1 Often in applied problems, for example, in mathematical modeling of hydrodynamics [Suk11, Ka17, Rue05], heat and mass transfer [Suk18, Suk18’], geofiltration, population dynamics [Suk05], seismic exploration [Mur13] and other processes, it is necessary to solve the equations of convection-diffusion 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 modified 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 difference schemes are effective 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 filling 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 Problem Statement The substance transport problem can be represented by the diffusion-convection-reaction equation: 0 0 c0t + uc0x + vc0y = (µc0x )x + µc0y y + f with boundary conditions: c0n (x, y, t) = αn c + βn , where u,v are components of the velocity vector; µ is the turbulent exchange coefficient; 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 = {tn = nτ, xi = ihx , yj = jhy ; n = 0..Nt , i = 0..Nx , j = 0..Ny ; Nt τ = T, Nx hx = lx , Ny hy = ly } , where τ is a time step; hx , hy are spatial steps; Nt is an upper boundary on time; Nx , Ny are space boundaries. 0 Discrete analogues of convective uc0x and diffusive (µc0x )x transfer operators of the second order of approxima- tion error in the case of partially filled cells can be written as: ci+1, j − ci, j ci, j − ci−1, j (q0 )i, j uc0x ' (q1 )i, j ui+1/2, j + (q2 )i, j ui−1/2, j , (1) 2hx 2hx 0 ci+1, j − ci, j ci, j − ci−1, j (q0 )i, j (µc0x )x ' (q1 )i, j µi+1/2, j 2 − (q2 )i, j µi−1/2,j − (2) hx h2x αx ci, j + βx − (q1 )i, j − (q2 )i, j µi, j , hx where qi are coefficients, describing the fullness of control domain [Suk15]. 3 Schemes of High Order of Accuracy for Convective and Diffusive Transfer Op- erators 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 difference schemes: - the discrete analogue of the convective transport operator in absence of influence of domain boundary: ci+1, j − ci, j ci, j − ci−1, j uc0x ' ui+1/2, j + ui−1/2, j , (3) 2hx 2hx 2 - the discrete analogue of the diffusive transport operator in absence of influence of domain boundary: 0 ci+1, j − ci, j ci, j − ci−1, j (µc0x )x ' µi+1/2, j − µi−1/2,j , (4) h2x h2x The approximation error of expression (3) will take the following form: ci+1, j − ci, j ci, j − ci−1, j ui+1/2, j + ui−1/2, j = 2hx 2hx 0 00 000 0 00 0(ci, j ) (ui, j ) 2 ui, j (ci, j ) (ui, j ) (ci, j ) 2 h2x + hx + O h4x .  = ui, j (ci, j ) + hx + 4 6 4 Therefore, for approximation the convective transport operator uc0 by difference scheme of the fourth order of accuracy we have to approximate the operator uc0 − c0 u00 h2 /4 − uc000 h2 /6 − u0 c00 h2 /4 by the scheme of the second order of accuracy. The approximation of the convective transport operator uc0 by difference scheme of the fourth order of accuracy has the form: ui+1/2 (q1 )i+1    ui+1/2 (q1 )i (q0 )i L (c) = − (q1 )i ci+2 − − (q1 )i 2+ + (5) 12h (q0 )i+1 12h (q0 )i  (q2 )i+1    ui−1/2 (q1 )i  u i+1/2 (1) (2) ui+1/2 + (q2 )i + (q1 )i − + ki + ki ci+1 + − (q1 )i 2+ + 12h (q0 )i 2h 12h (q0 )i+1 (q1 )i−1    ui−1/2 ui−1/2 ui+1/2 (1) (2) + (q2 )i 2+ + (q2 )i − (q1 )i − ((q2 )i − (q1 )i ) ki + ((q2 )i + (q1 )i ) ki ci + 12h (q0 )i−1 2h 2h     ui+1/2 (q2 )i ui−1/2 (q2 )i u i−1/2 (2) (1) − − (q1 )i + (q2 )i 2+ + (q2 )i + ki − ki ci−1 − 12h (q0 )i 12h (q0 )i 2h ui−1/2 (q2 )i−1   − − (q2 )i ci−2 , 12h (q0 )i−1   (q ) (q ) (q ) (q ) where ki = (q10 )i (ui+1 − ui, ) − (q20 )i (ui − ui−1 ) / (8h), ki = (q10 )i ui+18h−ui + (q20 )i ui −u (1) (2) i−1 8h . i i i i The approximation error of expression (4) will take the following form: 0 2 ci+1, j − ci, j ci, j − ci−1, j  0  (IV ) hx µi+1/2, j − µi−1/2,j = µi, j (c i, j ) + µi, j (ci, j ) + h2x h2x 12 00 00 h2x 0 0 0 0 h2 x 000 0 h2 + (µi, j ) (ci, j ) x + O h4x .  + (µi, j ) (ci, j ) + (µi, j ) (ci, j ) 4 6 6 0 Therefore, for approximation the diffusive transport operator (µc0 ) by difference scheme of the fourth order 0 of accuracy we have to approximate the operator (µc0 ) − µc(IV ) h2 /12 − µ00 c00 h2 /4 − µ0 c000 h2 /6 − µ000 c0 h2 /6 by the scheme of the second order of accuracy. 0 The diffusive transport operator (µc0 ) by difference scheme of the fourth order of accuracy can be written as: (q0 )i L (c) = −Ai ci + B1,i ci+1 + B2,i ci−1 + B3,i ci+2 + B4,i ci−2 . (6) µ00i+1 − µ00i   µi+1/2 µi+1 (q1 )i µi−1 (q1 )i (3) B1,i = (q1 )i + (q 1 )i + 2 + (q 2 ) i − (q 1 ) k i i − (q1 )i , h2 12h2 (q0 )i 12h2 (q0 )i 12 µ00i − µ00i−1   µi−1/2 µi+1 (q2 )i µi−1 (q2 )i (3) B2,i = (q2 )i + (q 1 )i + (q2 )i + 2 − (q 2 ) k i i − (q2 )i , h2 12h2 (q0 )i 12h2 (q0 )i 12 µi+1 (q1 )i+1 µi−1 (q2 )i−1 B3,i = − (q1 )i , B4,i = − (q2 )i , 12h2 (q0 )i+1 12h2 (q0 )i−1 µi+1 (q2 )i+1   µi+1/2 µi−1/2 (3) Ai = (q1 )i + (q2 )i − ((q1 )i + (q2 )i ) ki + (q1 )i +2 + h2 h2 12h2 (q0 )i+1 3 (q1 )i−1 µ00 − µ00i−1 µ00 − µ00i   µi−1 µi−1 (q1 )i + (q2 )i + 2 − (q2 )i i − (q1 )i i+1 + (q2 )i + 12h2 (q0 )i−1 12 12 12h2 (q0 )i µi+1 (q2 )i µi+1 (q1 )i+1 µi−1 (q2 )i−1 + (q1 )i 2 − (q1 )i 2 − (q2 )i , 12h (q0 )i 12h (q0 )i+1 12h2 (q0 )i−1   (q ) −µi (q ) (q ) (q ) − (q20 )i µi −µ (3) where ki = (q10 )i µi+1 4h2 i−1 4h2 , µ00i = (q10 )i ci+1 − 2ci + (q20 )i ci−1 /h2 . i i i i 4 Comparison of Calculation Results of Substance Transport Problem Based on Schemes of the Second and Fourth Orders of Accuracy The field, describing the error of calculations obtained as the difference between the analytical and numerical solution of substance transport problem, is given in Fig. 1. The initial distribution was determined by the function:  sin (π (x − 10)) cos (π (y − 10)) , {x, y} ∈ D, D : {x ∈ [10, 20] , y ∈ [10, 20]} . C (x, y) = 0, {x, y} ∈/ D, 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 coefficient of turbulent exchange is 2 m2 /s. Figure 1: Computational accuracy of substances: overhand – for diffusion problem; below – for diffusion- convection problem (schemes of the second order of accuracy on the left, the fourth – on the right 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 diffusion problem, and in 48.7 times – for solution the diffusion problem-convection. 5 Parallel Implementation of Diffusion-Donvection Problem Solution A library of two-layer iterative methods for solution the nine-diagonal grid equations was designed for solution the two-dimensional diffusion-convection problem based on the schemes of high order of accuracy. This library 4 Figure 2: Values of substance concentration at initial and final time moments 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. Table 1: Dependences of amount iteration for SLAE solution by iteration methods from the time step Time step Number of iterations Jacobi Minimal Method of Zeidel Upper relax- MATM method correction steepest method ation method method descent 0.001 6 6 6 5 43 5 0.005 8 8 8 8 43 6 0.01 10 10 10 8 45 6 0.05 23 23 23 15 56 10 0.1 37 36 37 22 61 12 0.5 138 134 138 70 60 27 1 256 247 256 126 60 28 5 1138 1077 1138 558 131 50 10 2233 2110 2233 1073 246 72 50 10160 9523 10160 4774 1074 158 100 19966 18625 19966 9320 2096 218 500 99651 92789 99651 46383 10399 1281 1000 199295 185529 199295 92739 20781 4382 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 first 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 specified error, then the return to the calculation of the residual is performed. At calculation the value of the solution vector, only the first 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 5 Figure 3: Scheme for calculation the values of solution vector on the next time layer of processors. Thus, the problem of calculation (selection) the optimal amount of transferred data package occur. Values of acceleration and efficiency of parallel implementation of the software, designed to solve the two- dimensional diffusion-convection problem on the basis of high order accuracy schemes, are given in the Table 2. The grid equations were solved by the modified 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. Table 2: Number of Time, s. Acceleration Efficiency processors 1 1447.415 1 1 2 734.728 1.97 0.985 4 387.009 3.74 0.935 8 199.643 7.25 0.906 16 109.653 13.2 0.825 32 62.659 23.1 0.722 64 36.643 39.5 0.617 According to the table 2, the parallel algorithm of the modified alternating-triangular method can be applied to solve real problems, and the use of parallel technologies makes a significant contribution to reduce the calculation time. 6 Use the High Order Accuracy Schemes for Reconstruction the Salinity Field and Comparison of Interpolation Results With Other Algorithms 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 specified 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 calcu- lations of hydrodynamic processes, it is necessary to approximate the function of two variables describing the salinity field by more stable functions. Formulation the problem of calculation the salinity field. To determine the salinity function, we use the diffusion equation solution to which the Saint-Venant equation describing the transport of bottom materials is reduced [Sid17]. The solution of the diffusion problem for a long time intervals is reduced to the solution of the Laplace equation: ∆H = 0, (7) where H is a water salinity. 6 Figure 4: The original image of salinity isolines’ level in the Azov Sea This approach has a significant disadvantage due to the lack of smoothness at points where the salinity field values are specified. To resolve this problem, we can use the following equation: ∆2 H = 0. (8) The disadvantages of this approach include large outliers (deviation from the linear function). With the first 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: h2 2 ∆H − ∆ H = 0. (9) 12 Note that, the operator for the third problem can be written as a linear combination of operators for the first and second problems. The fundamental system of solutions for equation (7) is the following function: H1 (x) = 1, H2 (x) = x, (10) for equation (8): H1 (x) = 1, H2 (x) = x, H3 (x) = x2 , H4 (x) = x3 , (11) for equation (9): √ H1 (x) = 1, H2 (x) = x, H3 (x) = ch (kx) , H4 (x) = sh (kx) , k = 12/h. (12) In the first 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 field restoration. The proposed mathematical algorithm for determine the water salinity field was numerically implemented. The salinity isolines were obtained using the recognition algorithm (Fig. 5a) The salinity field 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 sufficient degree of smoothness at points of gluing functions and lower emissions compared to the cubic function used in the calculations. 7 Figure 5: a) Isolines of salinity of the Azov Sea; b) the result of interpolation by the proposed method Figure 6: Restored salinity field of the Azov Sea 7 Conclusion Schemes of high (fourth) order of accuracy for the convective and diffusive transfer operators, taking into account the filling of the cells, were constructed. A library of two-layer iterative methods was designed and implemented on MCS for solution the two-dimensional diffusion-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 diffusion problem, and in 48.7 times – for solution the diffusion problem-convection. The algorithms description of parallel implementation of iterative methods with preconditioners of triangular type and value of acceleration and efficiency of parallel variant of algorithm of the modified alternative triangular method is given. A mathematical algorithm was proposed to restore the water salinity field 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 sufficient 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 8 to calculate the three-dimensional velocity flow fields 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. References [Sam89] Samarskii A.A., Gulin A.V. Numerical methods. M.: Nauka, 1989. [Kon02] Konovalov A.N. On the theory of the alternating-triangular iterative method // Siberian mathematical journal. 2002. Vol. 43. No. 3. P. 552. [Suk84] Sukhinov A.I. Modified alternating-triangular method for heat conduction and filtration problems // Computing systems and algorithms. 1984. P. 52-59. [Suk11] Sukhinov A.I., Chistyakov A.E., Alekseenko E.V. Numerical realization of the three-dimensional model of hydrodynamics for shallow water basins on a high-performance system // Mathematical Models and Computer Simulations. 2011. 3 (5). P. 562-574. [Ka17] Ka C.O., Dembele J.-M., Cambier C., Stinckwich S., Lo M., Zucker J.-D. Deterministic convection- diffusion 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 difference scheme for the two-dimensional convection-diffusion 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. Su- percomputer modeling of hydrochemical condition of shallow waters in summer taking into account the influence 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 ba- sis 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 modified 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. [Iak05] 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 diffusion equations solved by schemes with weights // Mathematical Models and Computer Simulations. 2014. 6(3). P. 324-331. 9 [Pet13] Petrov I.B., Favorskaya A.V., Sannikov A.V., Kvasov I.. Grid-characteristic method using high order interpolation on tetrahedral hierarchical grids with multiple time step // Mathematical Models and Computer Simulations. 2013. Vol. 5. Issue 5. P. 409-415. [Lad09] Ladonkina M.E., Neklyudova O.A., Tishkin V.F., Chevanin V.S. About one choice of essentially non- oscillatory high occuracy order difference scheme for systems of conservation laws // Matem. Mod., 21:11 (2009), P. 19-32. [Suk15] Sukhinov A.I., Chistyakov A.E., Semenyakina A.A., Nikitina A.V. Parallel implementation of the objectives of the transport of substances and recovery of the bottom and on the top of the news on the basis of difference schemes of increased order of accuracy // Proceedings of the International scientific conference “Parallel computational technologies” (PCT-2015). 2015. P. 285-296. [Gus18] Gushchin V.A., Sukhinov A.I., Nikitina A.V., Chistyakov A.E., Semenyakina A.A. A Model of trans- port and transformation of biogenic elements in the coastal system and its numerical implementation // Computational Mathematics and Mathematical Physics. 2018. Vol. 58 (8). P. 1316-1333. [Sid17] 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. 10