Mathematical modeling of eutrophication processes in Azov Sea on supercomputers Текст научной статьи по специальности «Математика»

UDC 519.6

Mathematical modeling of eutrophication processes in Azov Sea on supercomputers*

A. V. Nikitina1, A. A. Semenyakina2**

1 Southern Federal University, Rostov-on-Don, Russian Federation

2 Supercomputers and Neurocomputers Research Center, Taganrog, Russian Federation

Introduction. The paper is devoted to the research of effective numerical methods for solving eutrophication problem of shallow waters taking into account water environment, spatially-nonuniform distribution of temperature and salinity, microturbulent diffusion, gravitational sedimentation, and spreading of biogenic pollution, oxygen, phyto- and zooplankton, etc. The simulation objects are shallow waters - the Azov Sea and Taganrog Bay.

Materials and Methods. The mathematical model of eutrophication of shallow waters was developed. Parallel implementation was performed for computationally laborious convection-diffusion problems, taking into account the architecture and parameters of supercomputers based on the decomposition methods of grid domains. We determined that the maximum acceleration was 43 times on 128 computational nodes. We developed two algorithms including the algorithm based on the k-means method for data distribution between processors in parallel implementation. Due to the using these algorithms, the efficiency of algorithm for solving the problem is increased on 15% compared to the algorithm of the standard partition of computational domain.

Results. New mathematical models and software complex were developed for mathematical modeling of eutrophication processes in shallow waters. The concentrations of pollutants and plankton calculated for different wind situations were taken into consideration, if the relative error did not exceed 30%.

Due to expedition researches the primary verification of the model of ecosystem of the Azov Sea was performed. The problem of modeling and forecasting the state of water ecosystems of the Azov Sea in conditions of anthropogenic influence and comprehensive research of the unique water object was implemented. Because of the water object is shallow, it's more affected by anthropogenic influence.

The software complex, combining mathematical models and databases, was designed. Using this complex we researched conditions which are contributed by the eutrophication processes in shallow waters.

Discussion and Conclusions. Due to the solving the water ecology problem we can forecast different scenario of changing the water quality in shallow waters, and to investigate the mechanisms of formation of zones with low oxygen content.

Keywords: mathematical modeling, eutrophication, minimum corrections method, parallel computing, supercomputer.

Introduction. Researching of shallow waters is important and significant in mathematical modeling of aquatic ecosystems. The example of such waters is the Azov Sea, the great shallow water. Such water areas experience a tremendous anthropogenic influence, but not many of them are

* This paper was partially supported by the grant No.17-11-01286 of the Russian Science Foundation.

** E-mail: nikitina.vm@gmail.com , j.a.s.s.y@mail.ru

unique fish-producing ecological systems. The biogenic matters enter the shallow waters with the river flows which causing the growth of the algae - «water bloom».

The results of satellite monitoring of the Earth are used in this paper to control the quality modeling of processes of hydrodynamics and biological kinetics [1, 2]. The satellite monitoring data of the Azov Sea, obtained by SRC «Planeta», are given in Fig.1 [3]. The analysis of satellite data reveals the water areas of suffocations.

Colorized images of Azov Sea

(spectral channels: 0.620-0.670 mem, 0.54S-0.565 mem, 0.459-0.479 mem)

Fig. 1. The wide areas of the «water bloom» in the Azov Sea

The 3D spatially heterogeneous mathematical model was performed for the reconstruction of the «water bloom». This process caused the suffocation in the South-Eastern part of the Azov Sea in July 2013. The information about the wind velocity and direction in the Temryuk Bay in July 2016, provided the meteorological station in Kerch city (WMO_ID 33983) and shown in Fig.2, was used for this model as input data.

CB c CCB j \ C3 3C3

/\ / \ r^

\ ccbZ-j \ /

CC3 CC3 C \ /

11.07.15 13.07.13 15.07.13 17.07.13 19.07.13 21.07.13 23 12.07.13 14.07.13 16.07.13 18.07.13 20.07.13 22.07.13

Fig. 2. Wind velocity and direction, July 16, 2013 Water temperature in the computational domain for the simulated time interval is shown in

Fig.3.

Fig. 3. Water temperature, July 16, 2013

We used classes, described by the Institute for nature protection and reserves Ministry of ecology of Russia and given in Table 1, for analyzing the water quality of the Azov Sea. In Table 1: 1 - trophicity, class quality, the nature of saprobity, options; 2 - oligotrophy, very clean, xenosaprobic (I); 3 - esotropia, clean, batalhas-fineness (II); 4 - esotropia, moderately polluted, al-fareros-fineness (III); 5 - eutropia, polluted, betamatch-fineness (IV); 6 - polytrope, dirty, alphabets-fineness (V); 7 - hypereutrophy, very dirty, polycarobonate (VI); *) - watercourses; **) - stagnant waters; ***) - possible the higher values.

Table 1

Parameters and alphabet of classes of the state estimation of aquatic ecosystems

1 2 3 4 5 6 7

Water transparency by Secchi disk, m, *) **) 3 6 0,7-3 4 0,5-0,7 4 0,3-0,5 2 0,1-0,3 1 0,05-0,1 0,5

Specific conductivity of water, mks/sm *) **) >400 150 700 250 700 300***) 1100 500 1300 1000 1600 1000

The saprobity phytoplankton index according to Watanabe 85-100 70-85 50-70 30-50 15-30 0-15

Concentration, Cl «a», mkg/l, **) 3 8 8 15 30 60

Phytoplankton biomass, mg/l 0-0,1 0,1-0,5 0,5-1,0 1,0-5,0 5,0-50,0 50-100

Gross daily phytoplankton production , gO/sq.m 0-1,5 1,5-3 3-4,5 4,5-7,5 7,5-10,5 10,5-12

The Shannon index, H. 0-4 1-4.5 0-5 0-5 1,5-4,5 0-4

Range of variation, H. H min. - H max. 0-1,5 3-4 1-2 4-4,5 0-2 4,5-5 0-2 4,5-5 1,5-2 4-4,5 0-1,5 2-4

Materials and methods. Hydrodynamic Mathematical Model. The initial equation of hydrodynamics of shallow waters is the next: - motion equations (The Navier-Stokes motion equations):

u't + uu'x + vu'y + wu'2 = -1 p'x + + ^y )y + {vwz )Z +

- p'x + (K )x +([iu'y )y +(v< )Z + 2Q(vsin 6 - wcos 6), v + uvl. + W, , + wv; =--p'y + ((<)x +([iv'y)y +(vv'2)Z -2Qusin6 , (1)

w[ + uw'x + wy + ww'z = --p; +(^wX )x +(^wy ^ +(vwZ ); + 2Qucos 6 + g{p0 / p- continuity equation was written for the case of variable density:

p't + (pu )x + (pv)y + (pw); = 0, (2)

where u = {u,v, w} - velocity vector components; p is an excess pressure above the undisturbed fluid hydrostatic pressure; p is density; Q is Earth's angular velocity; 0 is an angle between the angular velocity vector and the vertical vector; p, v are horizontal and vertical components of turbulent exchange coefficient.

We consider the equation system (1), (2) with the following boundary conditions:

- at the entrance ( the mouth of Don and Kuban rivers):

u(x,y,z,t) = u(t),v(x,y,z,t) = v(t), p'n(x,y,z,t) = 0,un(x,y,z,t) = 0,

- the lateral boundary (beach and bottom):

pXU' )n (x У, z,t) = -T x( t), pv^( v ' )n (x У, z,t) = -T y( t), Un(x,y,z,t) = 0, p'n(x,y,z,t) = 0,

- the upper boundary:

pmTu' )n(xУ,z,t) = -tx(t), pv(v n)n(xУ,z,t) = -ty(t), w(x,y,t) = -ra-p'n/pg, p'n(x,y,t) = 0, (3)

- at the output (Kerch Strait): p'n (x, y, z, t) = 0, u'n(x,y,z,t) = 0,

where ra is a liquid evaporation intensity; tx ,t are tangential stress components (Van-Dorn law);

pv is suspension density.

Tangential stress components for free surface are in the form:

T x =paCp ( w )wxH , T y = p aC p ( A KM ,

where w is a wind velocity vector relative to the water; pa is atmosphere density,

f0.0088; x < 6.6 m/s C„ (x) = i - non-dimensional coefficient.

pW [0.0026; x > 6.6 m/s

Tangential stress components for bottom are in the form:

tx = pcp (ulVM, t y =Pcp (ul)vu .

We can define the coefficient of the vertical turbulent exchange with inhomogeneous depth on the basis of the measured velocity pulsation:

(4)

v = CfA'i \ + [*] ,

s 2\ [a*J [a* J

where A is a grid scale; Cs is non-dimensional empirical constant, defined on the basis of attenuation process calculation of homogeneous isotropic turbulence.

Grid method was used for solving the problem (1) - (3) [4]. The approximation of equations by time variable was performed on the basis of splitting schemes into physical processes [5 - 9] in the form of the pressure correction method.

Three-dimensional mathematical model of eutrofication processes of shallow water. Computational domain G (Fig. 4) is a closed area, limited by the undisturbed water surface £0, bottom £H = £H(x,y), and the cylindrical surface, the undisturbed surface o for 0 < t < T0. Z = Z0 oo - the sectionally smooth boundary of the domain G [10 - 14].

Fig. 4. Diagram of the computational domain G

The developed model is described by equations:

n, dSt dS Slt + u—- + v—- + (w - w ,

' dx dy g dz

(w - = ^ ASi

d f dSt ^

dz

(5)

where is the concentration of i-th impurity, index i indicates the substance type, i = 1,17: 1 is hydrogen sulfide (H2 S ) ; 2 is elemental sulfur (S ) ; 3 are sulfates (SO4 ); 4 are thiosulfates (and sulfites); 5 is the total organic nitrogen

(N )

; 6 is ammonium (nH4 ) (ammonia nitrogen); 7 are nitrites (NO2 ) ; 8 are nitrates (NO3 ) ; 9 is dissolved manganese (DOMn); 10 is suspended manganese (POMn); 11 is dissolved oxygen (O2); 12 are silicates ( SiO3 is the metasilicate; SiO4 is the orthosilicate); 13 are phosphates (PO4); 14 is ferrum (Fe2+); 15 is silica (H2SiO3 is meta-silicic; H2SiO4 is orthosilicic); 16 is phytoplankton; 17 is zooplankton; u = (u,v,w) is the water flow velocity; w is the gravitational sedimentation velocity of i-th component (in a suspended state); ^i,vi is the coefficient of turbulent diffusion in horizontal and vertical directions correspondingly; yi is a chemical-biological source or component describing the aggregation if the corresponding component is suspension.

We consider the system (5) with the following boundary conditions:

-S

Si = 0 on a,ifUn < 0; -i = 0 on a, if Un > 0;

on (6)

S'l z = 0 on E0 ;S',Z =-S;S; on S^ where s, is the absorption coefficient of the i-th component by the bottom material.

We has to add the following initial conditions to (5):

Si|t=0 = Si0 (x,y,z),i = I-6- (7)

Water flow velocity fields, calculated according to the model (1) - (3), are used as input data for the model (1) - (3). The discretization of models (1) - (4), (5) - (7) was performed on the basis of the high-resolution schemes which are described in [15].

Using the eutrofication model (5) - (6) the processes of ammonification, nitrification, nitra-tereductase (denitrification), assimilation NH4, oxidation H2S, sulfate reduction, oxidation and recovery of manganese can be described, and we can research the formation of zones of reduced oxygen and forecast changes of oxygen and nutrient regimes in shallow waters.

The processes of biogeochemical cycles of chemical elements, due to which the transfer of aerobic conditions in anaerobic is performed, were parameterized in developing the model.

The three-dimensional eutrophication model of the Azov Sea (5) - (6) was investigated, and sufficient conditions for the existence and uniqueness were obtained and formulated as theorems.

Theorem 1. Let Si(x,y,z,t), v e C2mt) oC^t), where U,t = G x (0 < t < T0) ;

^ = const > 0; U = {u,v,w- wg),Vj(z) e Cl(G); Sio e C(G), i = 1,17. Thus, when the following inequalities are held: max{\xi,vi}--1 max\pi|}>0 for all i = 1,17, where v = pi{Sj)S ,

G 0 G

as, ^ _ _ aT SS

v, az

i * J; Vi = LS, - DS, - pS, LSi =-S + div^U,), DSt = ^ AS, +

ot -z __

21 2 2 2 |

X0 = {1 /Lx +1 /L.y +1 /Lz); Lx,Ly,Lz - the spatially maximum dimensions of the computational domain; the model problem of eutrofication of shallow water in Eq. (5) - (7) has a solution.

Theorem 2. Let Si(x,y,z,t), vi e C2(Цt) oC^t), ^ = const > 0; U, v(z) e Cl(G),

- • • • • i 2 2 1 2

i = 1,17. Thus, when the following inequalities are held: 2^ {1 / L2 +1 / Ly )+ 2^/L2 > ^ for all

i = 1,17 (where are functions, defined by Sources of pollutant) the model problem of eutrofication of shallow water in Eq. (5) - (7) has a unique solution.

We constructed the coherent grid &h for discretization the model (5) - (6). hx-hy-hz are vector parameters, characterizing the density of nodes:

®h ={xj = Jhx,yk = khy,zl = lhz; J = 0,Nx;k = 0,Ny;l = 0,NZ;

Nxhx = Lx,Nyhy = Ly,Nzhz = Lz }®hx =®h X®T , =ltn = n^,n = 0- Nt }

where i, J, k are indexes in x,y, z directions; Nx, N^, Nz is the number of nodes in the coordinate directions.

The computational domain in spatial directions x, y, z is represented by the combination of

parallelepipeds. We performed the discretization of model (5) - (7) using the difference schemes with weights:

1 < j < Nx -1, 1 < k < N -1, 1 < l < Nz -1, i = 1, 17.

We added the approximated initial and boundary conditions to the system (8). We investigated the stability of difference schemes. The difference scheme (8) in the canonical form:

D ora+1 . D n»+l D n»+l A on+1 D fl»+l

B1S(i)j+1,k,l + B2S(i)j,k+1,l + B3S(i)j,k,l+1 + AS(i)-B4S(i)j-1,k,l -

D Çtn+1 D on+1 _ i f>» D On i D On ,

- B5S(i)j,k-1, / - B6S(i)j, k,/-1 S(i)j,k,l + B7S(i)j+1, k ,/ + B8S(i)j, k+1, / +

+ B9S"i)j,k,/+1 - B10S(1)j-1,k,/ - B11S(i)j,k-1,/ - B12S"i)j,k,/-1 ,.,

1 < j < Nx -1,1 < k < Ny -1,1 < / < Nz -1, 0 < n < Nt -1, i = 117.

Sufficient condition for the stability and monotonicity of the model (8) is determined based on the maximum principle of A.A. Samarsky with the following limitations:

K <

^ /U"ILA <lh /V"IL>K <ll2v< -*"«)

,i = 1,17.

C)' y 1 C)' z 1 )' The analysis of the accuracy of approximation of the scheme with weights in the form (8)

showed that it has the highest order of accuracy on temporary and spatial variables: o(t 2

where |h| + h, + hZ, , in the case of y = 0.5 (the Crank-Nicholson scheme).

The initial data for proposed model of eutrofication of shallow water in the form (5) - (7) is the velocity vector field of the water flow calculated by model (1)-(4).

Method for solving grid equations. The grid equations, obtained in the finite-difference approximations of tasks (1) - (3), (5) - (8), can be presented in the matrix form [16]:

Ax = f, (9)

where A is a linear, positive definite operator (A > 0). We use the implicit iterative process for solving problem (8):

xm+1 — xm

B-+ Axm = f. (10)

^ m+1

In Eq. (10) m is the number of iteration, t > 0 is an iterative parameter, and B is an invertible operator (a stabilizer). The inverting of the operator B in Eq. (10) should be significantly easier than the directly inverting of the original operator A in Eq. (9). We construct B using the additive representation of operator A0, i.e., the symmetric part of the operator A:

Ao = R + R2, R = R,, (11)

* *

where A = A0 + A1, A0 = A* , A = ~A .

The operator-stabilizer can be written as follows:

B = (D + oRj )D~l(D + OR2 ),D = D* > 0, Q> 0, (12)

where D is some, generally diagonal, operator.

Relations (11), (12) define the modified alternating triangular method (MATM) for solving the problems if the operators Rx ,R2 are defined and methods of determining the parameters tm+1, ®

and the operator D are specified.

The algorithm of the adaptive modified alternating triangular method of minimal corrections for calculating the grid equations with nonself-adjoint operators is in the form:

~(Dwm,wm)

rm = Axm — f , B(<»m)wm = rm, Sm =

\p~lR2wm,R2w™) , (13)

(Ao wm,wmJ _{B~lA1wm,A1wm )

(bA0wm,A0wm\Bwm,wm) , m "(b-1 A0wm,A0wm)

1 -.

v 2 k

sm km

v(l + km )

1 + k (1 - s

ra ■

Cm+1 0m'

a wm,wm ) (s- a0 wm,a wm ) '

m+1

x

x

m+1

co„

where rm is the residual vector, wm is the correction vector, the diagonal part of the operator A is used as the operator D . [17, 18].

The estimation of convergence rate of this method is in the form:

v -1 * , ,

p< —-, v -v(V 1 + k +

v +1

(VT+Ï+4k )2,k-(s-A^vd,

V ' (S-1 A0 am,A0 ©m )

where V is the condition number of the operator C0, C0 = B 172A0B

Parallel implementation. We describe the parallel algorithms, which are used for solving problems (1) - (3), (5) - (7), with different types of domain decomposition.

Algorithm 1. Each processor receive its computational domain after the partition of the initial computational domain into two coordinate directions, as shown in Fig. 5. The adjacent domains overlap by two layers of nodes in the perpendicular direction to the plane of the partition [19].

The residual vector and it uniform norm are calculated after that as each processor will receive the information for its part of the domain. Then, each processor determines the maximum element in module of the residual vector and transmits its value to all remaining calculators. Now receiving the maximum element on each processor is enough to calculate the uniform norm of the residual vector [20].

-1 / 2

Fig. 5. Domain decomposition

The parallel algorithm for calculating the correction vector is in the form:

(D + amRi )D_1 (D + amR2 )wm = rm, where R is the lower-triangular matrix, and R2 is the upper-triangular matrix. We should solve consistently the next two equations for calculating the correction vector:

(D + »mRi )ym = rm, (D + »mR2 )wm = Dym.

At first, the vector ym is calculated, and the calculation is started in the lower left corner. Then, the correction vector wm is calculated from the upper right corner. The calculation scheme of the vector ym is given in Fig. 6 (the transferring elements after the calculation of two layers by the first processor is presented).

In the first step of calculating the first processor operate on with the top layer. Then the transfer of overlapping elements is occurred to the adjacent processors. In the next step the first processor operate on with the second layer, and its neighbors - the first. The transfer of elements after calculating two layers by the first processor is given in Fig. 6. In the scheme for the calculation

of the vector ym only the first processor does not require additional information and can

independently operate on with its part of the domain. Other processors wait the results from the previous processor, while it transfers the calculated values of the grid functions to the grid nodes, located in the preceding positions of this line. The process continues until all the layers will be calculated. Similarly, we can solve the systems of linear algebraic equations (SLAE) with the upper-triangular matrix for calculating the correction vector. Further, the scalar products are calculated (12), and the transition is proceeded to the next iteration layer.

Fig. 6. Scheme of calculation the vector ym

We constructed the theoretical estimate of the time. It's required to perform the MATM step for SLAE with seven-diagonal matrix with using decomposition in two spatial directions on a cluster of distributed calculations. All computational domain is distributed among processors (n is the total number of processors, n = nx • ny, nx > ny), i.e. each of them received the domain by the size

N/n, N = NxNyNZ, where Nx,Ny,NZ is the number of nodes in the spatial directions; t0 is an execution time of one arithmetic operation; tx is response times (latency); tn is the time, required to transfer the floating point numbers.

14 N

For computing the residual vector the - arithmetical operations should be performed by

n

each processor, and 2(nx -1) transfers at NyNZ variables and 2(ny -1) transfers at NxNZ variables should be performed. Thus, the total computational time for residual vector is following:

14 N

-1o + 2(nx - 1)tnNyNz + 2(ny - 1)tnNxNz.

n

We defined the computational time of the vector wm .

N Ny

—x—y is the number of elements in one layer;

nx ny

NxNy

9-— • t0 is the computational time by one processor of each layer;

n

1,2,...,ny -1,ny,ny,...,ny,ny -1,...,2 is the sequence of starting processor,

nx-ny+1

9(n

nx + ny -

xNxNv

2)-— 10 is the latency time of starting the last processor (without the time for

n

transferring data);

N

9—10 is the computational time by the last processor; n

(nx -1 ){tnNy + txny)Nz + (ny -1 ){tnNx + txnx )Nz is the total transferring time at computing the vector ym .

We obtained that the computational time for vector wm which equals to the following: 18Nt0 /n + 18(nx + ny - 2)NxNyt0 7 n + + 2Nz {(nx - 1){tnNy + txny ) + (ny - 1){tnNx + txnx )) 16Nt0 / n - the computational time for calculating scalar products of fast descent. Time for the transition to the next iteration was equal to 2Nt^ / n.

Thus, the computation time for calculating on n processors is following: 50Nt0/n + 18{nx + ny - 2)NxNyt0 /n +

+ 2Nz {(nx - 1^{2tnNy + txny )+ (ny - 1){2t„Nx + txnx ))

Let's consider the case nx = n =*Jn ; so, the computational time is following: tN = 50N • t0 -1){36NxNyt0 /n + 4Nz{tn{Nx + Ny )+)).

The obtained theoretical estimations [21] of acceleration and efficiency of the parallel algorithm 1:

S(1) - -

50 Ntn

tN 50 N • tn +{4n - 1J36-"

n

10

n

t0 + 4 Nz (tn (Nx + Ny )+ tx4n )|

or

S( 1 ) - '

n

, (r A 36 4n

1 + y n -1)-+-

v h 50N 50t 0

f f

11

— + —

NN

+ -

tyfn

\\

E(1) - ■

S.

( 1 )

V v x y

1

NxNy x y ))

n , i r A 36 4n( f

1 + y n -1 )-+-

v h 50NZ 50t0

1 1 ^ tx4n ^

-+- +—-

,NX Nv NxNv

VV x y J x y yy

t„

t

n

We considered the case of the problem solution with the rectangular domain. The domain has a complex shape in the case of real water. At the same time the real acceleration is less than its theoretical estimation. The dependence of the acceleration, obtained in the theoretical estimates, can be used as the upper estimate of the acceleration for parallel implementation of the MATM algorithm by the domain decomposition in two spatial directions.

We describe the domain decomposition in two spatial directions with using the k-means algorithm.

Algorithm 2. The k-means method was used for geometric partition of the computational domain for the uniform loading of MCS calculators (processors). This method is based on the minimization of the functional of the total variance of the element scatter (nodes of the computational grid) relative to the gravity center of subdomains: Q = Q(3). Let Xi is the set of computational grid nodes, included in the i-th subdomain,ie-[l,...,m}, m is the given number of subdomains.

q( 3 ) = £-L- £ d2 (x,ci) ^ min, where ci = -¡-^ £ x is the center of the subdomain Xt, and

i |Xi|xeXi |X,|xeXi

d(x,ci ) is the distance between the calculated node and the center of the grid subdomain in the Euclidean metric. The k-means method converges only when all subdomain will be approximately equal.

The algorithm of k-means method.

1) The initial centers of subdomains are selected with using maximum algorithm.

2) All calculated nodes are divided into m Voronoi's cells by the method of the nearest neighbor, i.e. the current calculated grid node x e Xc, where Xc is a subdomain, which is chosen according

to the condition ||x -sc || = min\x -s^, where the sc is the center of the subdomain Xc.

1<i<m

3) New centers are calculated by the formula: sf+1 ) = ^ £ x .

X[k) xeX(k)

4) The condition of the stop is checked s(ck+l) = s(ck), k = 1,...,m. If the condition of the stop is not

performed, then the transition proceeds to the item 2 of the algorithm.

The result of the k-means method for model domains is given in Fig. 7 (arrows indicate the exchanges between subdomains). All points in the boundary of each subdomains are required to data exchange in the computational process. The Jarvis's algorithm was used for this aim (the task of constructing the convex hull). The list of the neighboring subdomains for each subdomain was created, and an algorithm was developed for data transfer between subdomains.

Fig. 7. Domain decomposition

Theoretical estimates of the acceleration and efficiency of the algorithm 2 were obtained similarly to the corresponding estimates of the algorithm 1:

S (2)='

n-X

, (r J 36 4n

1 + y n -1)-+-

v h 50N 50t0

f f

\

E( 2 ; = ■

S.

( 2 ;

i

V^ IN + NT y

X

+

tx4n

\\ •

x y yy

n

, /r ,1 36 4n

1 + y n -1)-+-

v h 50N 50tn

11

-+ —

yNx Nv V V x y y

+ -

t4n

W'

x y yy

where x is the ratio of the number of computational nodes to the total number of nodes (computational and fictitious).

Parallel algorithms of the adaptive alternating triangular method was implemented on supercomputer with the peak performance 18.8 TFlops. The system includes 8 computational racks. The computational field is designed on the basis of the HP BladeSystem c-class infrastructure with integrated communication modules, power and cooling systems. The computational nodes are 128 single-type 16-core HP ProLiant BL685c Blade-servers, each of which has the four 4-core AMD Op-teron 8356 2.3 GHz processors and the operative memory in the volume of 32 GB. The total number of cores in the complex - 2048, the total amount of RAM - 4 TB.

The comparison of the developed parallel algorithms 1 and 2 for the solution (1) - (3), (5) -(7) was performed. The results are given in Table 2.

Table 2

Comparison of acceleration and efficiency of algorithms

" '(1) S' S(i) S(1) t(2) E' E(2) E E(2)

1 7.490639 1.0 1.0 6.072899 1.0 1.0

2 4.151767 1.653577 1.804205 3.121229 1.181126 1.945675

4 2.549591 3.256077 2.937976 1.810628 2.325769 3.354028

8 1.450203 6.317738 5.165234 0.996729 4.512670 6.092825

16 0.88242 11.928279 8.488745 0.619345 8.520199 9.805356

32 0.458085 21.482173 16.352072 0.317173 15.344409 19.146924

64 0.265781 35.954877 28.18350 0.183929 25.682055 33.017611

128 0.171535 54.617841 43.668283 0.116936 39.012744 51.933099

In Table 2: n is a number of processors; t(k),S(k),E(k) are the calculation time, acceleration and efficiency of the k -th algorithm; S'(k), E'(k) are the theoretical estimations acceleration and efficiency of the k -th algorithm, k = {1,2}.

According to Table 1, we can conclude that the developed algorithms based on the decomposition method in two spatial directions and k-means method can be effectively used for solving hydrodynamics problems in the case the sufficiently large number of computational nodes.

The graphs of accelerations of algorithm 1 and 2 for solving the WB problem (5) - (7), obtained theoretically and practically, are given in Fig. 8.

Fig. 8. Graphs of accelerations for the developed parallel algorithms: 1 - the theoretical estimation of acceleration algorithm 1; 2 - the acceleration of the algorithm 2, obtained practically; 3 - the acceleration of the algorithm 1, obtained practically; 4 - the theoretical estimations acceleration of algorithm 2

The estimation was used for comparison the efficiency values of the algorithm 1 and 2, obtained practically:

Thus, the use of the algorithm 2 and k-means method for the WB problem (5) - (6) increased the efficiency on 15%.;

Results. The external frequency, leading to complication of the system, is taken into account in modeling spatially heterogeneous euthophication water processes of the Azov Sea (1) - (3). In this case the fluctuations of plankton density can be so much that it cannot be explained by the random fluctuations. Relatively small areas of high density ("slicks", "clouds") are separated by zones with low densities, sometimes not fixed by standard observational methods. Especially clearly this phenomenon is expressed in water areas which are characterized by the necessity for nutrient elements. Vegetative period of phytoplankton were taken into account in modeling eutrophication processes.

Diffusion processes occur in the direction of smoothing the spatial distribution and dispersion of "slicks". One of the attempts to explain the paradox of stability "slicks" with the help of numerical experiments is to assume the active movement of heterotrophic organisms (zooplankton and fish) in the direction of the gradient "food" that provides consolidation of spatial heterogeneity of nutrients in the aquatic environment. Sustainable heterogeneity of the spatial distribution can be, for example, due to diffusion processes and the presence of phytoplankton mechanism actoring regulation, i.e. regulation the rate of growth through selection in the environment of the biologically active metabolites.

Results of modeling the concentration of pollution (total organic nitrogen) for the eutrophication model of the Azov Sea (the initial distribution of water flow fields with the northern wind) are given in Fig 9. The influence of water flow structures in the Azov Sea on the distribution of biogenic pollution and phytoplankton are shown in the figures below.

The value 8 in the Eq. (13) was calculated by:

Fig. 9. Distribution of pollution concentration in different time interval

Results of modeling the phytoplankton dynamics in the Azov Sea are given in Fig. 10. (N is a number of iteration).

Maximum concentration of biogenic pollution is assigned by white color, minimum concentration of phytoplankton is assigned by the black color.

The verification criterion of the developed models (1) - (3), (5) - (7) was an estimate of the error modeling taking into account the available field data measurements at the same time, calculated according to the formula: 8= tt {sknat - Sk y / /¿S£2nat , where Sknat is the value of the

V k=1 V k=1

harmful algae concentration, obtained through field measurements; S^, is the value of the harmful algae concentration, calculated by the model (1) - (3).

<0f I

0f\ #1

ÄjTj ¿fl

Fig. 10. Distribution of phytoplankton concentration in different time interval

The concentrations of pollutants and plankton calculated for different wind situations were taken into consideration, if the relative error did not exceed 30%.

Discussion and Conclusions. Due to expedition researches the primary verification of the model of ecosystem of the Azov Sea was performed. The problem of modeling and forecasting the state of water ecosystems of the Azov Sea in conditions of anthropogenic influence and comprehensive research of the unique water object was implemented. Because of the water object is shallow, it's more affected by anthropogenic influence.

The software complex, combining mathematical models and databases, was designed. Using this complex we researched conditions which are contributed the eutrophication processes in shallow waters.

The distinctive features of the developed algorithms, implementing hydrobiological model problems, are the following: high performance, reliability and accuracy of the results. High performance is achieved by using efficient numerical methods for solving grid equations, aimed for use on parallel computer systems in real and accelerated time intervals. The accuracy is achieved by taking into account the important physical factors, such as: the Coriolis force, turbulent exchange, the complex geometry of bottom and coastline, evaporation, river flows, a dynamic rebuild of the computational region, wind stress and friction on the bottom, and also taking into account the deviation of the pressure field from the hydrostatic approximation. The accuracy is achieved by using detailed computational grids, taking into account the degree of "fullness" of computational cells, and the absence of nonconservative dissipative terms and revision sources arising from finite difference approximations.

The comparison of the developed software complex that implements the designed scenarios for the development of ecological situation in the Azov Sea using the numerical realization of model plankton evolution problems of the biological kinetics use with the similar works in the mathematical modeling of hydro-biological processes.

The analysis showed that due to using the developed software system we increased the accuracy of forecasts of changes in concentrations of pollutants, phyto- and zooplankton in the Azov Sea on 10 - 15% depending on the model problem of water ecology.

References

1. Samarsky, A.A., Nikolaev, E.S. Methods of solving grid equations. Moscow, Science, 1978, pp. 588.

2. Sukhinov, A.I., Chistyakov, A.E. Adaptive modified alternating triangular iterative method for solving grid equations with non-selfadjoint operator. Mathematical modeling, 2012, vol. 24, no. 1, pp. 3-20.

3. State Research Center "Planeta", available at: http://planet.iitp.ru/english/index_eng.htm.

4. Samarsky, A.A. Theory of difference schemes. Moscow, Nauka, 1989.

5. Konovalov, A.N. The method of steepest descent with adaptive alternately-triangular pre-amplification. Differential equations, 2004, vol. 40, no. 7, pp. 953.

6. Konovalov, A.N. The theory of alternating-triangular iterative method. Siberian mathematical journal, 2002, vol. 43, no. 3, pp. 552.

7. Sukhinov, A.I., Chistyakov, A.E., Shishenya, A.V. Error estimate of the solution of the diffusion equation on the basis of the schemes with weights. Mathematical modeling, 2013, vol. 25, no. 11, pp. 53-64.

8. Chetverushkin, B., Gasilov, V., Iakobovski, M., Polyakov, S., Kartasheva, E., Boldarev, A., Abalakin, I., Minkin, A. Unstructured mesh processing in parallel CFD project GIMM. Parallel Computational Fluid Dynamics, Amsterdam, Elsevier, 2005, pp. 501-508.

9. Petrov, I.B., Favorsky, A.V., Sannikov, A.V., Kvasov, I.E. Grid-characteristic method using high order interpolation on tetrahedral hierarchical meshes with a multiple time step. Mathematical modeling, 2013, vol. 25, no. 2, pp. 42-52.

10. Sukhinov, A.I., Chistyakov, A.E., Semenyakina, A.A., Nikitina, A.V. Parallel realization of the tasks of the transport of substances and recovery of the bottom surface on the basis of schemes of high order of accuracy. Computational methods and programming: new computing technology, 2015, vol. 16, no. 2, pp. 256-267.

11. Chistyakov, A.E., Hachunts, D.S., Nikitina, A.V., Protsenko, E.A., Kuznetsova, I.Y. Parallel Library of iterative methods of the SLAE solvers for problem of convection -diffusion-based decomposition in one spatial direction. Modern problems of science and education, 2015, no. 1-1, pp. 1786.

12. Sukhinov, A.I., Nikitina, A.V., Semenyakina, A.A., Protsenko, E.A. Complex programs and algorithms to calculate sediment transport and multi-component suspensions on a multiprocessor computer system. Engineering journal of Don. 2015, vol. 38, no. 4 (38), pp. 52.

13. Nikitina, A.V., Abramenko, Y.A., Chistyakov, A.E. Mathematical modeling of oil spill in shallow waters. Informatics, computer science and engineering education, 2015, № 3 (23), pp. 49-55.

14. Chistyakov, A.E., Nikitina, A.V., Ougolnitsky, G.A., Puchkin, V.M., Semenov, I.S., Sukhinov, A.I., Usov, A.B. A differential game model of preventing fish kills in shallow waterbodies. Game Theory and Applications, 2015, vol. 17, pp. 37-48.

15. Sukhinov, A.I., Nikitina, A.V., Semenyakina, A.A., Chistyakov, A.E., A set of models, explicit regularized schemes of high order of accuracy and programs for predictive modeling of consequences of emergency oil spill. Parallel computational technologies (PCT ' 2016), proceedings of the international scientific conference, 2016, pp. 308-319.

16. Nikitina, A.V., Semenyakina, A.A., Chistyakov, A.E. Parallel implementation of the tasks of diffusion-convection-based schemes of high order of accuracy. Vestnik of computer and information technology, 2016, no. 7 (145), pp. 3-8.

17. Sukhinov, A.I., Chistyakov, A.E., Semenyakina, A.A., Nikitina, A.V. Numerical modeling of an ecological condition of the Sea of Azov with application of schemes of the raised accuracy order on the multiprocessor computing system. Computer researches and modeling, 2016, vol. 8, no. 1, pp. 151-168.

18. Sukhinov, A.I., Nikitina, A.V., Semenyakina, A.A., Chistyakov, A.E. Complex of models, explicit regularized schemes of high-order of accuracy and applications for predictive modeling of aftermath of emergency oil spill. CEUR Workshop Proceedings. 10th Annual International Sci-

entific Conference on Parallel Computing Technologies, PCT 2016, Arkhangelsk, Russian Federation, 2016, vol. 1576, pp. 308-319.

19. Sukhinov, A.I., Nikitina, A.V., Chistyakov, A.E., Semenov, I.S., Semenyakina, A.A., Khachunts, D.S. Mathematical modeling of eutrophication processes in shallow waters on multiprocessor computer system. CEUR Workshop Proceedings. 10th Annual International Scientific Conference on Parallel Computing Technologies, PCT 2016, Arkhangelsk; Russian Federation, 2016, vol. 1576, pp. 320-333.

20. Nikitina, A.V., Sukhinov, A.I., Ougolnitsky, G.A., Usov, A.B., Chistyakov, A.E., Puchkin, M.V., Semenov, I.S. Optimal management of sustainable development at the biological rehabilitation of the Azov Sea. Mathematical modeling, 2016, vol. 28, no. 7, pp. 96-106.

21. Nikitina, A.V., Sukhinov, A.I., Ougolnitsky, G.A., Usov, A.B., Chistyakov, A.E., Puchkin, M.V., Semenov, I.S. Optimal control of sustainable development in the biological rehabilitation of the Azov Sea. Mathematical Models and Computer Simulations, 2017, no. 9 (1), pp. 101 -107.

Authors:

Nikitina Alla Valerevna, Associate Professor, Department of intellectual and multiprocessor systems, South state university (SFU) (Bolshaya Sadovaya street, 105/42, Rostov-on-Don, Russian Federation, 344006) Dr.Sci. (Eng.), Associate Professor

Semenyakina Alena Aleksandrovna, Junior researcher, Supercomputers and Neurocomputers Research Center (Italyansky lane, 106, Taganrog, Rostov region, Russian Federation, 347900)

УДК 519.63

Математическое моделирование процессов эвтрофикации в Азовском море на супер-ЭВМ*

А. В. Никитина1, А. А. Семенякина2**

1 Южный федеральный университет, г. Ростов-на-Дону, Российская Федерация

2 Научно-исследовательский центр супер-ЭВМ и нейрокомпьютеров, г. Таганрог, Российская Федерация

Введение. Статья посвящена исследованию эффективных численных методов решения задачи эвтрофикации вод мелководного водоема с учетом множества факторов, таких, как движение водного потока, пространственно-неравномерное распределение температуры и солености, микротурбулентная диффузия, гравитационное оседание, а также распространение загрязняющих биогенных веществ, кислорода, фито- и зоопланктона и др. Объектом моделирования выступают мелководные водоемы - Азовское море и Таганрогский залив. Материалы и методы. Разработана математическая модель эвтрофикации мелководных водоемов. Параллельная реализация выполнена на основе методов декомпозиции сеточных об-

* Работа выполнена при поддержке РНФ, проект № 17-11-01286

** E-mail: nikitina.vm@gmail.com , i.a.s.s.y@mail.ru

ластей для вычислительно трудоемких задач диффузии-конвекции, учитывающих архитектуру и параметры супер-ЭВМ. Было установлено, что максимальное ускорение достигалось на 128 вычислительных узлах и составило 43 раза. При реализации параллельного алгоритма решения задачи на МВС для распределения данных между процессорами были разработаны два алгоритма, в том числе алгоритм на основе метода к-шват, применение которого позволило повысить эффективность алгоритма решения поставленной задачи на 15% по сравнению с алгоритмом, основанном на стандартном разбиении расчетной области. Результаты исследования. Разработаны новые математические модели и программное обеспечение для математического моделирования процессов эвтрофикации мелководных водоемов. Рассчитанные при различных ветровых ситуациях концентрации загрязняющих веществ и планктона принимались к рассмотрению, если относительная погрешность не превышала 30%.

С помощью экспедиционных исследований проведена первичная верификация модели экосистемы Азовского моря. Реализована задача моделирования и прогноза состояния водной экосистемы Азовского моря в условиях антропогенного воздействия и всестороннего изучения уникального водного объекта, который в силу мелководности в большей степени подвержен антропогенному влиянию.

Создан программный комплекс, объединяющий разработанные математические модели и базы данных, с помощью которого изучены условия, при которых мелководные водоемы подвергаются эвтрофитрованию.

Обсуждение и заключения. Решение поставленной задачи водной экологии позволит прогнозировать возможные сценарии изменения качества вод мелководных водоемов, а также изучать механизмы формирования в них зон с пониженным содержанием кислорода.

Ключевые слова: математическое моделирование, эвтрофикация, метод минимальных поправок, параллельные вычисления, супер-ЭВМ.

Авторы:

Никитина Алла Валерьевна, доцент кафедры «Интеллектуальные и многопроцессорные системы» Южный федеральный университет (ЮФУ) (РФ, г. Ростов-на-Дону, ул. Б. Садовая, 105/42), доктор технических наук, доцент

Семенякина Алёна Александровна, младший научный сотрудник Научно-исследовательского центра супер-ЭВМ и нейрокомпьютеров (РФ, г. Таганрог, пер. Итальянский, д. 106)