CC BY
38UDC 51-77 10.23947/2587-8999-2019-2-2-83-103
Mathematical modeling of nonlinear effects in dynamic of interacting plankton and fish populations of Azov Sea*
A. I. Sukhinov1, V. V. Sidoryakina2, A. V. Nikitina3,4, A. E. Chistyakov1, A. A. Filina4, V. N. Litvinov1,5
:Don State Technical University, Rostov-on-Don, Russian Federation
2Taganrog institute named after A.P. Chekhov (branch) OF FSBIU «Rostov State University (RINH)», Taganrog, Russian Federation
3Southern Federal University, Rostov-on-Don, Russian Federation
Supercomputers and Neurocomputers Research Center, Taganrog, Russian Federation
5Azov-Black Sea Engineering Institute of Don State Agrarian University, Zernograd, Russian Federation
Paper covers the research of nonlinear effects in population dynamics of the pelengas commercial fish of the Azov Sea at low and high size taking into account the Allee effect, competition for resources, taxis, catching, spatial distribution of biogenic matter and detritus based on a multi-species model of plankton and fish interaction. Discrete analogue of developed model problem of water ecology, included in a software complex, were calculated using schemes of second order of accuracy taking into account the partial filling of computational cells. The system of grid equations of large dimension, arising at discretization, has been solved on the basis of a two-layer variational type method - the minimum corrections method having the maximum convergence rate. Effective parallel algorithms were developed for numerical implementation of biological kinetics problem and oriented on NVIDIA Tesla K80 graphics accelerator with the data storage format modification. Due to it, the reproduction processes of biogeocenose populations have been analyzed in real and accelerated time.
Keywords: Allee effect, taxis, mathematical model of population interaction, biogydrocenosis, Azov Sea, pelengas, detritus, parallel algorithm, modified data storage format, graphics accelerator, software.
Introduction. Among the fundamental problems of ecology and hydrobiology are researching the dynamics of biological communities, the circulation of matter and energy balance, the identification of cause-and-effect relationships and mechanisms of individual processes in ecological systems. Shallow waters, such as the Azov Sea and Taganrog Bay, are particularly strongly affected by anthropogenic pressure, associated with the construction of hydrotechnical constructions, various types of pollution. It leade to the emergence of environmental processes that violate the conditions of sustainable development of water.
* E-mail: sukhinov@gmail.com. cvv9@mail.ru. nikitina.vm@gmail.com. cheese 05@mail.ru. j.a.s.s.y@mail.ru. litvinovvn@rambler.ru
The reported study was funded by RFBR according to the research project № 19-01-00701
The experimental research of biological kinetics processes is based on large scale experiments with the marine ecosystem and is not only extremely expensive, but also dangerous in terms of consequences. Therefore, the creation of mathematical and software-algorithmic tools for predictive modeling of significant spatial hydrobiological processes and using data from relatively inexpensive experiments that do not have a negative impact on the ecosystem is relevant for today.
The development of concepts for researching of aquatic ecosystems involved many scientists such as Lotka A.J. [1], Volterra V. [2], Svirezhev Yu.M., Logofet D O. [3], Abakumov V.A. [4], Konstantinov A.S., Vinberg G.G. [5], Odum H.T. [6], Rosenberg G.S [7]. Methods for assessing the results of hydrobiological monitoring were developed by Makrushin A.V., Bakanov A.I., Mordukhay-Boltovskoi D.D., Tansley A.G., Rozenzweig M.L., MacArthur R.H. [8], Luck R.F. [9], Berryman A.A. [10]. Nonlinear effects of interaction of biological populations in water were researched by Berdnikov S.V., Tyutyunov Yu.V. [11], Stephens P.A. [12], Sutherland W.J., Arditi R. [13], Zhou S.-R., Liu Y.-F., Wang G. [14].
Ichthyological researches of fish population dynamics are practically important and theoretically complex areas of mathematical modeling. According to the Azov Research Institute of fisheries (AzNIIRKH), commercial fishing of pelengas, bream and walleye was prohibited in the Azov Sea waters in August 2018. There are only three species of commercial fish out of 40 commercial fish: the seal, the goby and the hamsa. Extinction is due to sea pollution and uncontrolled fishing.
According to the existing regulations, the time for decisions and eliminating contingencies of technogenic or natural characters should be from several hours to two to three days. Therefore, the time for precision the ecological situation of coastal systems in the event of emergencies is limited.
Several researches in mathematical modeling of processes of hydrophysics and biological kinetics are devoted to the parallel implementation of problems of this class. Although the conditions for the development of catastrophic and unfavorable phenomena in shallow waters, it is necessary to forecast the development of such phenomena and make decisions within tens of minutes - units of hours. It, in turn, requires the modeling of hydrobiological processes on multiprocessor computer systems on the accelerated time mode.
Universal modeling packages of hydrobiological processes (FlowVision, FLUENT, GAS DYNAMICS TOOL, PHOENICS, Star-CD, etc.) are focused on multiprocessor systems, but the versatility of these packages is to use the limited number of models, algorithms and methods to variety of different cases. Programs, aimed for solving the particular problems, have the potential to address these challenges more effectively. Most of the known specialized software (ADAM, CAL3QHC, Chensi, TASCflow, ISC-3, PANACHE, REMSAD, UAM-IV, ЭКОЛОГ, ПРИЗМА, VITECON), designed for calculation the pollution spread, is focused on single-processor systems. Only separate modules of specialized software (for example, ECOSIM и MAQSIP), oriented on multiprocessor systems, are parallelized. This fact does not allow achieving high efficiency computing in some cases.
These methods use simplified models of hydro-biological processes for water objects with the slightly varying salinity, and, in most cases, with the varying depth; do not provide the operational forecasts of the environmental situation of shallow waters after the disaster; cannot be replicated to other aquatic ecosystems, as focused on the description of the biogeochemical cycles and species composition of biological plankton populations and their interactions in single water objects.
To obtain operational forecasts of ecological situation of shallow water, it is necessary to design software, oriented on high-performance computing systems that allows simulating hydrobiological processes in water and their possible consequences in limited time. For complex geoinformational analysis of spatial-temporal processes and phenomena the expedition data [15], the data of the Unified State System of Information on the Situation in the World Oce an ("ESIMO") portal [16] were used for modeling hydrobiological processes in the Azov Sea.
2. Problem statement. Creating and using models of ecosystems is an integral part of a comprehensive research of ecological systems to assess its condition, to study the eutrophication processes, the transport and transformation of pollutants, reproduction processes, including fish productivity of water, solution of problems of rational use of its resources, preservation of water quality and model predictions further development of the ecosystem under the influence anthropogenic and climatic factors.
According to the literature [17, 18], biological purification of water can be carried out by stocking it by valuable pelengas commercial fish (Mugil soiuy Basilewsky). Its introduction into the water was recommended in the late 1960s to restore the fish production of the Azov Sea and improve its sanitary regime. Pelengas is a detritophagus; it feeds on decomposed organic residues that are part of detritus and partially settle to the sea bottom and seyston (planktonic organisms), periphyton. Pelengas withstands high salinity and chemical pollution of water, low content of dissolved oxygen (up to 1.4 ml/l), which is detrimental for most Azov fish, and has a high fertility.
2.1. Modeling the dynamics of commercial fish. A multi-species model of interaction between plankton and pelengas commercial fish was developed to research nonlinear effects in the dynamics of biogydrocenosis of a shallow water (Azov Sea). It takes into account the spatially inhomogeneous distribution of Coscinodiscus diatoms, Copepoda zooplankton, and based on the models of Ricker W.E. [19], Rosenzweig M.L., MacArthur R.H., Ginzburg V.L., Holling K.S., Arditi R., Tyutyunov Yu.V. [20], Latun V.S., Kosyan R.D., Pryakhin Yu.V. [21]:
SS S ( SS ^
f + div ( USi ) = +— V1~i = g (SJ. S3 )~S1S!S2 "4s! -O-^A.
S z
dt y • i i dz
¥2= g2{SVS2)-^2S2~Ô2S2> ¥i=rAS^-gî{S^Sî) + B(^-Sî) + f^ (!)
¥a = \S, - g4 (S4,S5 ) + ^S2 -A4S4, ^ = g5 (S4,S5 )-X5S5 S5S5,
where S is the concentration of i-th component, i =1,5 ; ^ is a chemical-biological source (runoff) or a term describing aggregation (conglutination-unconglutination) if the corresponding component is a suspension, i is a type of substance, i =1,5 : 1 is the Coscinodiscus phytoplankton concentration (X), 2 is the Copepoda zooplankton (Z), 3 is the nutrient (S), 4 is the detritus (D), 5 is the pelengas (P); u is the velocity field of water flow; U - u + u0; is the rate of convective mass transfer; u0i is
the deposition rate of i -th substance under the gravity, i e 1,4; g are trophic functions for substances i e 1,5. Let assume that g (S, S ) = /aSS, g (S, S2 ) = /2^1 SS, g3 (S > S ) = aSS, g4 ( S4 > S5 ) = P4S4S5 5 g5 (S1,S4,S5) = (r5^4S4 +^5a1S1 )S5, where a3 is the consumption ratio of S by phytoplankton; /,r5 are transfer coefficients of trophic functions; r is a fraction of S in the
biomass of X; \ is the coefficient taking into account the mortality and metabolism of X ; Sx is the loss of X by eating out of Z ; X2 ,AS are elimination (mortality) rates of Z, P respectively; S2 is the loss of Z by eating out of fish; S5 is the loss of P by eating out of fish and catching; S} is the maximum possible concentration of S; f = f (t, x, y, z) is the function of pollution source; B -specific rate of receipt of S; X4 is the decomposition ratio of detritus; pA is the consumption rate of organic residues of P ; ox is the loss ratio of X as a result of consuming it by P ; is the transfer coefficient of concentration growth of P due to X ; ^, vi are diffusion coefficients in horizontal and vertical directions of substances i e 1,5.
The computational domain G (Azov Sea) is a closed area bounded by the undisturbed water surface E0, the bottom EH = EH (x, y) and the cylindrical surface a for 0 < t < T0. Z = Z0 is a piecewise smooth border of G, n is the surface normal vector; Un is a normal component of the water velocity vector to the Z.
Let's define initial conditions
S*|t=0 = So(x, y, z), i =1,5 (2)
and boundary conditions:
8S dS dS
St = 0 on a if Un < 0; = 0 on a if Un > 0; = p(S,) on Eo; 1sL = -eiSi on ,
on oz oz
(3)
where sl ,£3 are non-negative constants; sl ,s3 ,s5 take into account the lowering of plankton and fish to the bottom and their flooding; £ ,£ take into account the absorption of nutrients and detritus by bottom sediments; <p is the given function.
2.2. Research the nonlinear effects of biogeocenose dynamics of shallow water. One of the most characteristic properties of biological objects is a wide range of nonlinear effects, within which the system functions normally [24]. The most important property of water ecosystems is the high velocity of biotic cycle. This is facilitated by both the mobility of water environment, mixing, conditions favorable for bacterial decomposition, dissolution and distribution of mineralization products, and intensive exchange due to the small size of producers (phytoplankton) and enzymes (zooplankton, fish).
Research the influence of taxis effect on spatial distribution of commercial fish populations. The taxis of pelagic (living in the water column) fish was taken into account explicitly at modeling, since it determines the movement of fish crowding in water area at any scale. Accounting for taxis, especially fish taxis, is a difficult problem, because the movement activity of organisms is directly related to the spatial structure and temporal variability of several characteristics of the ecotone (transition band between easily physiologically distinguishable communities), which include dependent variable models. In addition, gregarious effects affect the taxis of fish.
In mathematical biology, the taxis is considered in two ways. In one case, the weighted gradient of the stimulus is defined in accordance with the velocity of the active relative movement of
the organism (fast taxis), in the other - the acceleration (slow taxis). The proposed model simulates the spatial-temporal variability of the taxis rate of pelengas and specific energy equivalents of plankton, fish, nutrient, detritus biomass.
The direction and velocity of the collective taxis are established in the process of random wanderings and purposeful movement of individuals. The effect of such wanderings on the collective taxis is analogous to the effect of turbulent velocity pulsations on large-scale fluid flow. It can be taken into account as a kind of diffusion in the equation for the pelengas feed taxis.
We added the following equation to the system (1) - (3):
0 + div ( U5u5 ) = HU AU5 +-°{vu "T1! - auu5 + Kgr^dS1 + k4 gradS4 ,
ot oz \ oz ^
where U5 = u + u5 is the rate of convective pelengas transfer; u5 is the velocity of fish movement relative to the water; kx, k4 are taxis coefficients; /uu, vu are coefficients of horizontal and vertical diffusion component of taxis velocity; au is coefficient of inertial motion of fish.
As pelengas is a detritofag, he eats an excessive amount of deposits of organic silts (the remnants of phyto- and zooplankton). It causes the lack deficit of oxygen in water and leads to constant suffocation. Therefore, we add the following equation to the system (1) - (3) to account the distribution of oxygen concentration in water:
S
C = C0 - ko-4, C(x, y, z, t) = 0(x, y, z, t), V(x, y, z, t) e DG, h0
where C0 is the oxygen concentration in the absence of organic impurities; k0 is the amount of oxygen required to oxidize one conventional unit of detritus; h is constant value, the thickness of the bottom layer in which the detritus is located; H is the maximum water depth, h0 << H; Q is the given function.
Research the Allee effect to prevent collapse of commercial fish population. Complex regulatory processes, among which positive, ensure structural and functional stability of ecosystems and negative feedbacks are of great importance. Such fundamental process in ecosystems as the recycling of nutrients in a simple form can be represented in the form of positive feedback. It is carried out in "microbial loops", which can be considered as real "biogenic traps", through which up to 5060% of the net production of phytoplankton passes. For rational use, population prediction, identification of limiting factors, protection and restoration of populations of shallow water biohydrocenosis, it is necessary to take into account such important factor as the Allee effect [28]. The Allee effect is to combine biological individuals into groups (aggregation). This, on the one hand, increases competition between groups for food resources and living space, but on the other hand -increases the ability of survival of the group of individuals as a whole. According to Allee's principle, there is an optimal group size and optimal population density for each species. Both overpopulation and underdevelopment (or lack of aggregation) have an adverse impact on the ecological situation of the water.
The Allee effect is based on two fundamental laws of ecology about limiting factors: by Liebig and Shelford. According to the Liebig law, the limiter is the limiting factor, which is in the minimum; the Shelford law - in a maximum. All this is reduced to Shelford's law of tolerance, according to
which the prosperity of the organism is limited to zones of maximum and minimum of certain environmental factors. Each species is characterized by its tolerance - the ability to tolerate deviations from optimal environmental factors; a similar situation is manifested in higher-order communities -populations. For this reason, any population aspires to dynamic balance, i.e. to homeostasis at which achievement each population has strictly certain genetic, phenotypic, sex-age and other structure. It cannot consist of fewer or more individuals than is necessary to ensure the stable implementation of this structure and the stability of the population to environmental factors. A change in the population of a species is the result of an imbalance between the biotic potential of the water (the addition of individuals) and the resistance of the environment (the death of individuals).
Table 1. Functional forms of the productivity of ecosystems
(Stock-Recruitment (SR) models)
Name Expression Source
Ricker SR model n{P) = a-P-e bP, P - spawning stock; a is related to productivity (recruits per stock unit at small stock size) and b to density dependence. (a, b > 0) Ricker W.E., 1954
Beverton-Holt and Myers SR model = a ■ pj+b ■ Py a is the maximum recruitment (asymptotically) and b is the stock level needed to produce the half of maximum recruitment a/2 (a, b > 0) Beverton R.J.H. and Holt S.J., 1957; Myers R.A., 1995
Saila-Lorda, Cabral SR model // (i5) = a ■ Pc ■ e~bp, a,b,c are given parameters Saila S., 1988; Cabral R.B., 2010
Eppley, Gin SR models ju{T) = j exp(ruT), ju0 is the specific rate of phytoplankton growth at 0oC, ru is a constant Eppley R.W., 1972;Gin K.Y.H et al., 1998
Sidko SR model (p(I) = (bI + c2) - al - c, a, b, c are given parameters Sidko F Y. et. al., 1975
Belyanin SR model %(I) = %(1 - exp (-I / b)) / (l- aI), %, a, b are given parameters Belyanin V.N. et al.,1980
Dombrovskiy SR model j(%) = Uoexp )( -U% + U2 , T}, JU, (7%, JU, U are given parameters Dombrovsky Yu.A., 1990
Fennel SR Model j(T,I) = ao exp(aT)(I /IoPt)exp(1 -I /IoPt), «o is the growth rate of phytoplankton at 0oC (0,8), a is an empirical coefficient Fennell D.A., 1999
Logistic dependence taking into account the Allee effect R (SSB) = L(l + exp(-k(SSB - SSB ))) is the mass of fish capable for producing the offspring in an environment with a limited capacity L, k; SSB0 are parameters of the degree of the Allee effect (at low population densities). Vasilyeva N.A., Vladimirov A.A., Winter A.M., 2017
The influence of effect on the pelengas commercial fish population was research in this paper. It has been found that intensive fishing can lead to population collapse, causing serious and often
irreversible damage to the ecosystem. Table 1 shows the main functional forms used at mathematical modeling of the biogeocenose evolution of water taking into account limiting biotic and abiotic factors (C is the salinity, T is the temperature, I is the illumination) [29].
For description the chemical-biological sources we added observation models [7-9]. Various dependences in observational models are used to account the effects of C, T and I on the plankton
productivity function in water (Fig. 1): a) ¿u(C) = p0 exp -1(C- Copt) / ac j - jufi + ^2
values:
C = 12,
opt 9
^C = 10,
A = 0.001,
A = 0.1,
parameter A = 0.01; b)
A ( I, T, C ) = «0 exp ( aT ) ( I /1 op, ) exp (1 -1 / Ipt ) ^0 exp - {(C - Copt ) / } -Vf + %
a0 = 0.8,
a = 0.063, I = It = 86, Copt = 12, ^c = 15, V = 0.001, V = 0.1, V2 = 0.1, I = 10, Iopt = 86.
a)
b)
Fig. 1. Functional dependences of observational models.
In [24], nonlinear effects arising in the model of the age structure of commercial fish populations were analyzed. To research the weak or strong Allee effect, it is necessary to analyze the intensity of commercial fish catch and select its threshold value, at which collapse is impossible (see Fig. 2). Intensive commercial fishing near the coast of Canada led to the collapse of the cod population in 1992 [11] (Fig. 3).
Fig. 2. Allee effect in the ichthyological model
Fig. 3. Collapse of commercial fish population
In recent decades, a member of the cephalic family, an acclimatizant, pelengas has made the greatest progress in the ichthyofauna of the Azov Sea basin. After pelengas became a commercial fish of the Azov Sea, the first maximum catch of it was fixed by Ukraine in 1996 and amounted to 1031 tons; in 1997 - 2600 tons. The increase in commercial catch of pelengas in the Azov Sea was fixed until 2000. The second maximum was recorded in 2006. Over the next three years, there was a decrease in commercial stocks of pelengas almost twice. In the period from 2011 to 2013, there was a clear trend of decrease in commercial catch of pelengas in the Azov Sea by both Ukraine and Russia. In comparison with 2006, there was a decrease in catches in Russia by 8 times in 2013, in Ukraine -in 6 times. In 2014, the catch of pelengas by Russia in the Azov Sea was amounted to 215.5 tons. In 2002, the catch of pelengas was 80.4 tons in the Black Sea-Crimean region; it decreased significantly to 0.2 tons in 2017. This decrease is due not only to the existing natural pattern of population fluctuation, but also to anthropogenic factors. In particular, the reproduction conditions in the main spawning grounds of the Azov Sea - Milk estuary and Eastern Sivash - were violated. The increased fishing load on this facility since 2000, especially in winter, has also had a significant impact. The demographic Allee effect is due to the properties of the population itself and always manifests itself when its population is low. This effect can occur due to the characteristics of the food chain, as well as under the influence of external factors, including salinity, temperature and illumination. 3. Approximation of the biological kinetics problem
Each equation of the system (1) - (3) can be represented as the diffusion-convection-reaction equation in the two-dimensional case:
c't + uc'x + vc'y = (juc'x + (juc'y) ^ + f (5)
with boundary conditions:
c'„ (x y, t) = anc+pn,
where u, v are water velocity components; p is the turbulent exchange coefficient; f is the function, describing the intensity and distribution of sources.
For numerical implementation of the discrete mathematical model we used a uniform grid:
={tn = nx,x, = ihx,y = jhy; n = 0Nt,' = 0NX,j = 0Ny; Ntx = T,Nxhx = lx,Nyhy = ly}, where x is the time step; hx, h are spatial steps; Nt is the upper time boundary; Nx, N are spatial boundaries.
Discrete analogs of convective uc'x and diffusive (juc'x)' operators of the second order of accuracy in the case of partially filled cells can be written as:
(<7o< □ (<7l ",+1/x, C 2/; C + (<?2w,-i/2.y C 2/;C >
(),, j (Mc'x )'x □ (q ),, j M+1/2, j c'+1,j, 2 C'J - (42 ),, j M-1/2,j C'J h - |( ?1 )i, j - (42 ),, j | Pi, j axC'Jh +Px , (6)
hx hx hx
where ql, l e{0,1,2} are coefficients, describing the «fullness» of computational domains. Each
equation of the system (1) - (3) after linearization was represented as the convection-diffusion equation. Schemes of high-order of accuracy taking into account the partialling of computational cells
were used for model discretization [13]. The adaptive modified alternately triangular variational method was used to solve the obtained grid equations [13, 14, 15].
4. Parallel implementation. For numerical implementation of the developed model of biological kinetics (1) - (3), a parallel algorithm oriented on a multiprocessor computer system (MCS) was developed [30]. NVIDIA Tesla K80 graphics accelerator was used to improve the efficiency of the algorithm of numerical implementation of the set multi-species model of interacting populations of biogydrocenosis.
We describe parallel algorithms with various types of domain decomposition for solving problems (1) - (3) on MCS. Parallel algorithms for the MATM were implemented on MCS of the Southern Federal University (SFU). MCS technical parameters: the peak performance is 18.8 TFlops; 8 computational racks; the computational field of MCS is based on the HP BladeSystem c-class infrastructure with integrated communication modules, power supply and cooling systems; 512 single-type 16-core HP ProLiant BL685c Blade servers are used as computational nodes, each of which is equipped with four 4-core AMD Opteron 8356 2.3 GHz processors and 32 GB RAM; the total number of computational nodes is 2048; the total amount of RAM is 4 TB.
The k-means method was used for geometric partition of computational domain for the uniform loading of MCS calculators (processors) []. It based on the minimization of the functional of the total variance of the element scatter (nodes of computational grid) relative to the gravity center of subdomains. The result of the k-means method for model two-dimensional and three-dimensional regions is shown in Fig. 4.
a) b) c)
Fig. 4. Domain decomposition for the k-means method into 150 (two-dimensional domain) (a), into 10 subdomains (three- dimensional domain) (b); data transfer between neighboring subdomains (c).
For data exchange in the computing process, it is necessary to define all points on the border of each subdomain. For this, the Jarvis algorithm (the problem of constructing a convex hull) was used. The list of neighboring subdomains for each subdomain was formed, and the algorithm of data transfer between subdomains was developed (Fig. 6c, arrows indicate exchanges between subdomains).
Theoretical estimates of acceleration and efficiency of the developed parallel algorithm:
E = S'/p = X
1 +
(VP - 1)'i
36 50 N
4 p 50tn
( f 'n
V V
N.
N.
y J
NN„
\\
' JJ
where % is the ratio of the number of computing nodes to the total number of nodes (computing and
fictitious); p is the total number of processors; t0 is the execution of an arithmetic operation; tx is
the response times (latency); Nx, N , Nz are the number of nodes in the spatial directions.
Results of parallel implementation the proposed algorithm for solution the problem (1) - (3) were compared and presented in the Table 2, where t, S, E are operating time, acceleration and
efficiency of the algorithm; St, Et are theoretical estimates of the acceleration and efficiency of the
algorithm.
Table 2. Comparison of acceleration and efficiency of algorithms
p t, 5 S' S E' E
1 6.073 1.0 1.0 1 1
2 3.121 1.181 1.946 0.59 0.973
4 1.811 2.326 3.354 0.582 0.839
8 0.997 4.513 6.093 0.654 0.762
16 0.620 8.520 9.805 0.533 0.613
32 0.317 15.344 19.147 0.48 0.598
64 0.184 25.682 33.018 0.401 0.516
128 0.117 39.013 51.933 0.305 0.406
According to the Table 4, the developed algorithm based on the k-means method can be effectively used for solving hydrodynamic problems at sufficiently large number of computing nodes.
Dependence of acceleration on the number of processors is given in Fig. 5: a continuous curve is a theoretical dependence; a broken curve is an experimental one.
Fig. 5. Dependence of acceleration on the number of processors
The estimation for comparing the efficiency values of algorithms (Ei is efficiency of the standard algorithm; E2 is efficiency of the algorithm on the basis of k-means) has the form:
S =
E(2)k E(1)k) / JEE(2)k . (1)
y k=1 V k=1
We obtained that the efficiency is increased on 10 - 20% using of algorithm on the basis of k-means method for solving problem in the form (1) - (3) compared to the standard algorithm.
For numerical implementation of proposed interrelated mathematical model of biological kinetics, we developed parallel algorithms which will be adapted for hybrid computer systems using the NVIDIA CUDA architecture. The NVIDIA Tesla K80 computing accelerator has the high computing performance and supports all modern both the closed (CUDA) and open (OpenCL, DirectCompute) technologies. The NVIDIA Tesla K80 specifications: the GPU frequency of 560 MHz, the GDDR5 video memory of 24 GB, the video memory frequency of 5000 MHz, the video memory bus digit capacity is equaled to 768 bits. The NVIDIA CUDA platform characteristics: Windows 10 (x64) operating system, CUDA Toolkit v10.0.130, Intel Core i5-6600 3.3 GHz processor, DDR4 of RAM 32 GB, the NVIDIA GeForce GTX 750 Ti video card of 2GB, 640 CUDA cores.
Using the GPU with the CUDA technology is required to address the effective resource distribution at solving the system of linear algebraic equations (SLAE). The dependence of the SLAE solution time on the matrix dimension and the number of nonzero diagonals was obtained for implementation the corresponding algorithm (see Fig. 6).
■ 0.00-10.00 «10.00-20.00 ■ 20.00-30.00 30.0040.00 ■ 40.00-50.00 ■ 50.00-60.00 " 60.00-70.00 « 70.00-80.00 «80.00-90 00
Fig. 6. Dependence of SLAE solution time on matrix dimension and the number of nonzero diagonals
Due to it, in particular, we can choose the grid size and to determine the time for SLAE solution based on the amount of nonzero matrix diagonals.
Modified storage format of the sparse matrix with a repeating sequence of elements. The
solution of biological kinetics problem (1) - (3) by the finite difference method (FDM) on uniform grids leads to the necessary operate with sparse matrices, elements of which are a repeating sequence for internal nodes. In the case of high-dimensional problems, this leads to inefficient memory consumption [31]. Using the CSR (Compressed Sparse Rows) matrix storage format avoids the need to store their null elements. However, all nonzero elements, including many repeating, are stored in the corresponding array. This disadvantage is not critical at using computing systems with shared memory. However, this can adversely affect performance at data transferring between nodes in heterogeneous and distributed computing systems. We performed the modification of the CSR format to improve the efficiency of data storage with a repeating sequence of CSR1S elements for modeling
continuous biological kinetics processes by the finite difference method. In this case, it is enough to change them in an array that preserves a repeating sequence to change the differential operator, instead of repeatedly finding and replacing values of non-zero elements in an array.
The memory capacity for the CSR format was estimated as:
Pcsr = NnzBnz + (Nnz+R + 1)Bidx,
for the CSR1S format:
Pcsr1s = (Nnz - Nseq(Rseq + 1))Bnz+(Nnz - Rseq(Nseq+1) + R + 1)Bidx, where R is the number of matrix rows; Rseq is the number of matrix rows that contain a repeating sequence of elements; Nnz is number of non-zero matrix elements; Nseq is the number of elements in a repeating sequence; Bnz is the memory capacity to store one non-zero element; Bidx is the memory capacity to store one index.
Pcsns = Bnz[Nnz(ki + 1)- Nseq{kikrR + krR + 1) - kt(krR -R- 1)], where kr=^,k
R Bnz
Efficient function libraries [Cuda docs] were developed for solution of grid equations at discretization of model problem (1) - (3) in CSR format on GPUs using CUDA technology. The developed algorithm for solving the problem uses a CSR1S modified data storage format with further conversion to CSR format to solve the resulting SLAE on a graphics accelerator using NVIDIA CUDA technology. This raises the problem to develop an algorithm of matrix conversion from CSR1S format to CSR format in minimal time. A number of computational experiments with fivefold repetition and fixation of the average value of the calculation time was performed (Table. 3 - 5). An experimental research of the dependence of processing time of the matrix conversion algorithm from CSR1S to CSR format on the number of elements of the repeating sequence Nseq and the ratio of matrix rows, containing the sequence, to the total number of matrix rows for the sequence (Table 3), parallel using the TPL library (Table 4) and parallel using the NVIDIA CUDA platform (Table 5) is performed. Nseq values ranged from 3 to 19 in increments of 2. kr was changed in the range from 0 to 1 in increments of 0.1. The matrix dimension is 106.
Table 3. The runtime of matrix conversion (sequential algorithm), ms
kr Nseq
3 5 7 9 11 13 15 17 19
0 174 207 278 313 384 403 459 611 564
0.1 179 220 269 316 407 398 470 509 662
0.2 172 218 285 385 381 449 510 553 694
0.3 166 225 288 382 392 423 487 572 659
0.4 186 239 318 378 430 466 512 634 672
0.5 186 234 327 384 433 468 567 600 728
0.6 190 238 345 431 470 487 565 614 837
0.7 194 243 376 397 467 614 617 714 836
0.8 206 258 356 436 498 564 628 690 861
0.9 206 276 366 438 507 576 647 742 860
1 222 275 406 494 537 633 708 760 857
According to the analysis of the Table 3-5, the algorithm using the parallel computing on the CPU is more efficient than the serial algorithm and the algorithm using NVIDIA CUDA.
The graphs of the dependence of the matrix conversion time from CSR1S format to CSR format by sequential (a) and parallel (b) algorithms from fa and Nseq are shown in Fig. 7.
Table 4. The runtime of matrix conversion (parallel algorithm), ms
kr Neq
3 5 7 9 11 13 15 17 19
0 114 120 127 156 165 211 253 242 346
0.1 103 105 119 154 177 182 245 225 364
0.2 101 101 135 227 167 190 213 229 349
0.3 96 106 138 143 200 205 261 225 378
0.4 93 97 135 144 165 195 200 241 287
0.5 97 123 132 183 184 184 206 239 285
0.6 96 121 133 149 177 206 205 228 281
0.7 96 123 141 184 179 217 228 229 296
0.8 98 116 130 166 217 247 272 294 281
0.9 103 110 172 160 167 192 220 226 290
1 101 132 135 167 194 216 233 275 313
Table 5
The runtime of matrix conversion (using NVIDIA CUDA technology), ms
kr Nseq
3 5 7 9 11 13 15 17 19
0 320 359 442 504 576 644 760 781 849
0.1 301 332 420 525 537 631 645 800 817
0.2 287 359 440 455 537 623 647 748 796
0.3 273 330 390 443 527 582 612 711 764
0.4 278 304 384 435 492 585 594 629 764
0.5 279 317 356 431 488 529 636 678 782
0.6 271 303 352 416 463 495 549 616 683
0.7 270 301 373 399 431 479 517 606 615
0.8 256 291 364 386 415 458 504 542 546
0.9 264 309 325 402 443 472 509 501 561
1 264 285 320 365 413 448 476 510 586
Fig. 7. The runtime of matrix conversion from CSR1S to CSR format
Graphs of the efficiency of algorithms are given in Fig. 8: sequential implementation (Fig. 8a), parallel implementation using the TPL library (Fig. 8b) and parallel implementation using the NVIDIA CUDA platform (Fig. 8c).
i-fiK), Nseq - 19
900 800 700 600 500 400 300 200 100 0
Y = -315.27X +863.36 -1
y = 297. 9lx -+ 599.23 R2 = 0.8874
—»—en da
Линейная (cuda^
a
Ö J
0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
y=-0 .02x + 0 .8329
R: г = 0,9267
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 O.S 0.9 fc, o.e.
il
13
15
17
19
N,
seq
c) d)
Fig. 8. Graphs of algorithm efficiency at sequential and parallel implementation
Analysis of the graphs in Fig. 8 (a-c) shows that the algorithm using NVIDIA CUDA technology is more efficient at Nseq > 7. In this case, the point of equal efficiency decreases, starting from kr, equal to 0.7.
Regression equation y = -0.02x + 0.8329 (Fig. 8d) with a determination factor of 0.9276, describes the boundary of equal time-consuming sequential algorithm and the algorithm using NVIDIA CUDA. Therefore, we can calculate the minimum value of kr, above which the second algorithm will be more efficient than the first, substituting the value of Nseq in it. It is necessary to implement an improved iterative method that directly operates with data in the CSR1S format. Therefore, we can perform the scheme CSR1S CUDA.
Analysis of the CUDA architecture characteristics showed the algorithms for numerical implementation of the developed mathematical model of hydrobiological processes could be used for creation high-performance information systems.
5. Description of software complex. The software complex (SC) «Azov3D» was developed for solution of biological kinetics problem with implementation on MCS and graphic accelerator. The SC includes the following modules: control module, oceanological and meteorological databases, application library for solving hydrobiology grid problems, integration with various geoinformation systems (GIS), Global Resource Database (GRID) for geotagging and access to satellite data collection systems, NCEP/NCAR Reanalysis database (Fig. 9).
Fig. 9. Scheme of software complex
a)
b)
Fig. 10. Distribution of the S3 (a) and S (b) concentrations
The use of GIS provides additional possibilities for more qualitative and complex spatial analysis, and solutions based on it are more accurate. The problem was solved on a high-performance computer system [23] capable for performing a large amount of complex calculations and huge data processing in limited time.
Calculation results of pollution biogenic substance concentration (S3) for the model problem of biogeocenose evolution of shallow water (1) - (3) is given in Fig. 10a, and changes in phytoplankton concentration (S) - in Fig. 10b (initial distribution of water flow fields for north wind).
Using the developed SC, we researched the mechanism of suffocation areas in shallow waters [32]. Simulation results of possible scenarios for the development of the Azov Sea ecosystem (changes in the concentration of pelengas commercial fish) are given in Fig. 11, 12 (the initial distribution of water flow fields for the north wind direction). The maximum values of pelengas and detritus concentrations indicate by white color
Fig. 11. Distribution of detritus concentration, time Fig ,2 Distribution offish concentrations, time mterval T=62 days. mterval t=76 days.
Using the results of numerical experiment we can analyze the possible scenario of waters stocking the of the Azov Sea by pelengas involved in a bottom-land reclamation areas of detritus accumulation. According to the presented results, there is a decrease in the detritus concentration at the value of the time interval, starting from 61 days [33]. This means a decrease in the concentration of bottom sediments in the Central-Eastern part of the Azov Sea. It ultimately will lead to a decrease of suffocation areas and improve the water quality in this water basin.
With the help of developed special software, we can explore questions about acclimatization of pelengas on the environmental regime of a shallow water, to assess the specificity of the water conditions. This will avoid unforeseen negative technological impacts and plan to increase the production of this fish. The analysis of obtained results allowed concluding the properties of the detritus-pelengas mathematical model and the possibilities of water quality management of shallow waters, such as the Azov Sea, using methods of mathematical modeling.
Conclusion. Analysis of nonlinear effects in the dynamics of commercial fish populations at low and high numbers, taking into account the Allee effect, competition for resources, taxis, catch,
spatial distribution of nutrients and detritus was performed on the basis of a multi-species model of interaction between plankton and pelengas commercial fish. Due to it, we researched influence of competition for resources, features of the food chain, biotic and abiotic factors on the reproduction processes in water. Discretization of the developed model problem of water ecology as the part of software was performed on the basis of schemes of second order of accuracy taking into account the partial fulness of the computational cells. It allowed to significantly reducing the solution error in the computational domain of complex shape. The comparative analysis of two layer methods of variational type was performed. MCM was used as the main method of solving the system of grid equations arising in the discretization in view of its highest convergence rate. Effective parallel algorithms were developed for numerical implementation of biological kinetics problem and oriented on multiprocessor computer system and NVIDIA Tesla K80 graphics accelerator with the modification of data storage format. Due to it, the reproduction processes of biogeocenose populations have been analyzed in real and accelerated time.
References
1. Lotka A.J. Contribution to the energetics of evolution // Proc. Natl. Acad. Sci, 1922. No. 8. pp. 147-150.
2. Volterra V. Variations and fluctuations of the number of individuals in animal species living together // Rapp. P. - V. Reun. Cons. Int. Explor. Mer. 1928. Vol. 3. pp. 3-51.
3. Logofet D.O., Lesnaya E.V. The mathematics of Markov models: what Markov chains can really predict in forest successions // Ecological Modelling, 2000. Vol. 126, P. 285-298.
4. Abakumov A.I., Israel'skii Yu.G. Stabilizing role of the fish population structure in terms of fisheries impacts in a random environment // Computer Research and Modeling. 2017, Vol. 9, No.4, P. 609-620 (In Russian).
5. Vinberg, G.G. Biological productivity of waters // Ecology. 1983. No. 3. P. 3-12. (In Russian)
6. Odum, H.T. System Ecology // New York: Wiley. - 1983. - 644 p.
7. Rosenberg G.S., Zinchenko T.D. Stability of hydroecosystems: a review of the problem / / Arid ecosystems. 2014. No. 4 (61). P. 11-23 (In Russian).
8. Rosenzweig M.L., MacArthur R.H., 1963. Graphical representation and stability conditions of predator-prey interactions // American Naturalist. Vol. 97. No. 895. P. 217-223. https://doi.org/10.1086/282272
9. Luck R.F., 1990. Evaluation of natural enemies for biological control: a behavioral approach // Trends in Ecology and Evolution. Vol. 5. No. 4. P. 196-199. https://doi.org/10.1016/0169-5347(90)90210-5
10. Berryman A.A., 1999. The theoretical foundations of biological control / Hawkins B.A., Cornell H.V., eds. Theoretical Approaches to Biological Control. Cambridge: Cambridge Univ. Press. P. 3-21.
11. Perevarukha A.Yu. Chaotic regimes in the models of the formation theory of the replenishment of populations // Non-linear world. 2009. Vol. 7. No. 12. P. 925-932 (In Russian).
12. Stephens P.A., Sutherland W.J., 1999. Consequences of the Allee effect for behaviour, ecology and conservation // Trends in Ecology and Evolution. Vol. 14. No. 10. P. 401 -405. (In English) https://doi.org/10.1016/S0169-5347(99)01684-5
13. Arditi R., Ginzburg L.R., 1989. Coupling in predator-prey dynamics: ratio-dependence // Journal of Theoretical Biology. Vol. 139. No. 3. P. 311-326. https://doi.org/10.1016/S0022-5193(89)80211-5
14. Zhou S.-R., Liu Y.-F., Wang G., 2005. The stability of predator-prey systems subject to the Allee effects // Theoretical Population Biology. Vol. 67. P. 23-31. https://doi.org/10.1016Zj.tpb.2004.06.007
15. Sukhinov, A., Nikitina, A., Chistyakov, A., Sumbaev, V., Abramov, M., Semenyakina, A. Predictive modeling of suffocation in shallow waters on a multiprocessor computer system (2017) Lecture Notes in Computer Science, pp. 172-180. DOI: 10.1007/978-3-319-62932-2_16.
16. Esimo. http://portal.esimo.ru/portal
17. Matishov G.G., Ilyichev V.G. On optimal exploitation of water resources. The concept of internal prices // Reports of the Academy of Sciences. 2006. Vol. 406. No. 2. P. 249-251 (In Russian)
18. Menshutkin V. V., Rukhovets L. A., Filatov N. N. Modeling of freshwater lakes ecosystems (review). 2. Models of freshwater lakes ecosystems // Water Resources. 2014. Vol. 41. No. 1. P. 24-38 (In Russian).
19. Ricker W.E. 1954. Stock and recruitment. Journal of the Fisheries Research Board of Canada. Vol. 11. P. 559-623.
20. Tyutyunov Yu.V., Titova L.I., Senina I.N., 2017. Prey-taxis destabilizes homogeneous stationary state in spatial Gause-Kolmogorov-type model for predator-prey system // Ecol. Complexity. Vol. 31. P. 170-180. https://doi.org/10.1016/j.ecocom.2017.07.001
21. Pryakhin, Yu.V. Features of pelengas accounting in the Azov Sea //: Proceedings of the 7th All-Russian conference "Problems of fishing forecasting" - Murmansk: PINRO Publishing house. - 1998. - P. 155-156 (In Russian).
22. Sukhinov, A.I., Chistyakov, A.E., Ugol'nitskii, G.A., Usov, A.B., Nikitina, A.V., Puchkin, M.V., Semenov, I.S. Game-theoretic regulations for control mechanisms of sustainable development for shallow water ecosystems (2017) Automation and Remote Control, 78 (6), pp. 10591071. DOI: 10.1134/S0005117917060078
23. Alexander I. Sukhinov, Alexander E. Chistyakov, Alla V. Nikitina, Yulia V. Belova, Vladimir V. Sumbaev, and Alena A. Semenyakina. Supercomputer Modeling of Hydrochemical Condition of Shallow Waters in Summer Taking into Account the Influence of the Environment // Communications in Computer and Information Science. P. 336-351. https://doi.org/10.1007/978-3-319-99673-8.
24. Vasilieva N.A., Vladimirov A.A., Winter A.M. Consideration of nonlinear effects in the model of the age structure of the Atlantic cod population // Proceedings of the Fifth National scientific conference with international participation "Mathematical modeling in ecology", October 16-20, 2017. - Puschino, Institute of Physico-Chemical and Biological Problems of the Russian Academy of Sciences, 2017. p.42-44 (In Russian).
25. Marchuk G.I., Sarkisyan A.S. Mathematical modeling of ocean circulation. M.: Nauka, 1988. 297 p. (In Russian)
26. Sukhinov, A., Chistyakov, A. Isayev, A. Nikitina, A., Sumbaev, V., Semenyakina, A.: Complex of Models, High-Resolution Schemes and Programs for the Predictive Modeling of Suffocation in Shallow Waters. In: Communications in Computer and Information Science, vol. 753, pp. 169-185 (2017) DOI: 10.1007/978-3-319-67035-5_13.
27. Konovalov, A.N. The theory of alternating-triangular iterative method // Siberian mathematical journal, 2002. 43:3. C. 552-572. (In Russian)
28. Allee W.C., 1931. Animal Aggregations: A Study in General Sociology. Chicago: Chicago Univ. Press. 431 p.
29. Beverton, R.J.H., Holt S.J. 1957. On the dynamics of exploited fish populations, Fisheries Investigations (Series 2), vol. 19. United Kingdom Ministry of Agriculture and Fisheries, 533 pp.
30. Gergel V.P. high-Performance computing for multiprocessor multicore systems. - M.: Publishing house of Moscow University, 2010. - 544 p. (In Russian)
31. Voevodin V.V., Voevodin Vl.B. Parallel computing. - SPB.: BHV-Petersburg, 2002. - 608 p. (In Russian).
32. Gushchin, V.A., Sukhinov, A.I., Nikitina, A.V., Chistyakov, A.E., Semenyakina, A.A.: A Model of Transport and Transformation of Biogenic Elements in the Coastal System and Its Numerical Implementation. In: Computational Mathematics and Mathematical Physics, vol. 58, No. 8, pp. 1316-1333 (2018).
33. Alla Nikitina, Ludmila Kravchenko, Ilya Semenov, Yuliya Belova, and Alena Semenyakina. Modeling of production and destruction processes in coastal systems on a supercomputer // MATEC Web of Conferences 22, 04025 (2018), Volume 226 (2018), DTS-2018. DOI: https://doi.org/10.1051/matecconf/201822604025.
Authors:
Sukhinov Alexander Ivanovich, Don State Technical University (1st Gagarin Square, Rostov-on-Don, Russian Federation), Doctor of Science in Physics and Maths, Professor
Sidoryakina Valentina Vladimirovna, Taganrog Institute of A.P. Chekhov (branch) RSUE (Initiative Street, Taganrog, Russian Federation), Candidate of Science in Physics and Maths, Associate professor
Nikitina Alla Valeryevna, Southern Federal University (Chekhov str., 2, Taganrog, Russian Federation), Doctor of Technical Science, Associate professor
Chistyakov Alexander Evgenievich, Don State Technical University (1st Gagarin Square, Rostov-on-Don, Russian Federation), Doctor of Science in Physics and Maths, Associate professor
Filina Alena Aleksandrovna, Supercomputers and Neurocomputers Research Center, Co Ltd. (Italyansky lane, 106, Taganrog, Russian Federation), Candidate of Technical Science, Researcher
Litvinov Vladimir Nikolaevich, Azov-Black Sea Engineering Institute of FSBEI HE Don State Agrarian University (21 Lenin st., Zernograd, Russian Federation), Candidate of Technical Sciences, Head of the Thermal Engineering and Information Management Systems department
УДК 51-77 10.23947/2587-8999-2019-2-2-83-103
Математическое моделирование нелинейных эффектов в динамике взаимодействия планктонных и рыбных популяций Азовского моря^
A. И. Сухинов1, В. В. Сидорякина2, А. В. Никитина3'4, А. Е. Чистяков1, А. А. Филина4,
B. Н. Литвинов1'5
1 Донской государственный технический университет, г. Ростов-на-Дону, Российская Федерация
2 Таганрогский институт имени А.П. Чехова (филиал) ФГБОУ ВО «РГЭУ (РИНХ)», г. Таганрог, Российская Федерация
3 Южный федеральный университет, г. Ростов-на-Дону, Российская Федерация
4 НИЦ супер-ЭВМ и нейрокомпьютеров, г. Таганрог, Российская Федерация
5 Азово-Черноморский инженерный институт ФГБОУ ВО Донской ГАУ, г. Зерноград, Российская Федерация
Статья посвящена исследованию нелинейных эффектов в динамике популяций промысловой рыбы пеленгас Азовского моря в малых и больших размерах с учетом эффекта Аллее, конкуренции за ресурсы, отлова, пространственного распределения биогенного вещества и детрита на основе мультивидовой модели взаимодействия планктона и рыбы. Дискретный аналог разработанной модельной задачи водной экологии, входящий в программный комплекс, рассчитывался по схемам второго порядка точности с учетом частичного заполнения вычислительных ячеек. Система сеточных уравнений большой размерности, возникающих при дискретизации, была решена на основе двухслойного метода вариационного типа - метода минимальных поправок, имеющего максимальную скорость сходимости. Разработаны эффективные параллельные алгоритмы для численной реализации задачи биологической кинетики, ориентированные на графический ускоритель NVIDIA Tesla K80 с модификацией формата хранения данных. В связи с этим процессы воспроизводства популяций биогеоценоза анализировались в реальном и ускоренном времени.
Ключевые слова: Эффект Allee, математическая модель взаимодействия популяций, биогидроценоз, Азовское море, пеленгас, детрит, параллельный алгоритм, модифицированный формат хранения данных, графический ускоритель, программное обеспечение.
Авторы:
Сухинов Александр Иванович, Донской государственный технический университет (344000 Ростов-на-Дону, пл. Гагарина, д. 1), доктор физико-математических наук, профессор
Сидорякина Валентина Владимировна, Таганрогский институт им. А.П. Чехова (филиал) РГЭУ (РИНЭ) (347936 Таганрог, улица Инициативная, д. 48), кандидат физико-математических наук, доцент
^Исследование выполнено при финансовой поддержке РФФИ в рамках научного проекта № 19-01-00701 E-mail: sukhinov@gmail.com. cvv9@mail.ru. nikitina.vm@gmail.com, cheese 05@mail.ru, i.a.s.s.v@mail.ru, litvinovvn@rambler.ru
Никитина Алла Валерьевна, Южный Федеральный Университет (347928, Россия, Ростовская область, г. Таганрог, пер. Некрасовский, 44), доктор технических наук, профессор.
Чистяков Александр Евгеньевич, Донской государственный технический университет (344000 Ростов-на-Дону, пл. Гагарина, д. 1), доктор физико-математических наук, доцент
Филина Алёна Александровна, Научно-исследовательский центр супер-ЭВМ и нейрокомпьютеров (347900, Россия, Ростовская область, г. Таганрог, пер. Итальянский, дом 106), кандидат технических наук, научный сотрудник.
Литвинов Владимир Николаевич, Азово-Черноморский инженерный институт ФГБОУ ВО Донской ГАУ (РФ, г. Зерноград, ул. Ленина, 21), Кандидат технических наук, доцент.
