Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt Математическое моделирование процессов эвтрофикации в мелководных водоемах на многопроцессорной вычислительной системе* А.И. Сухинов1, А.В. Никитина2, А.Е. Чистяков2, И.С. Семенов2, А.А. Семенякина2, Д.С. Хачунц2 Донской государственный технический университет1, Научно-исследовательский институт многопроцессорных вычислительных систем им. академика А. В. Каляева Южного федерального университета2 Работа посвящена разработке методов решения модельной задачи эвтрофикации вод мелководного водоема, учитывающей движение водного потока, микротурбулентную диффузию, гравитационное оседание, пространственно-неравномерное распределение температуры и солености, а также загрязняющих биогенных веществ, кислорода, фито- и зоопланктона и др. В качестве объекта моделирования были выбраны мелко- водные водоемы – Азовское море и Таганрогский залив. Численное решение задачи основывается на градиентном методе вариационного типа – методе минимальных по- правок. Ускорение и эффективность расчетов достигается при использовании много- процессорной вычислительной системы Южного федерального университета. Реше- ние поставленной задачи водной экологии позволит прогнозировать изменения каче- ства вод мелководных водоемов, а также изучать механизмы формирования в них зон с пониженным содержанием кислорода. Ключевые слова: математическая модель, эвтрофикация, метод минимальных попра- вок, параллельные вычисления, Азовское море. 1. Введение Азовское море – внутреннее море на востоке Европы. Это самое мелкое море в мире, его глубина не превышает 13,5 метров. В силу своей мелководности оно наиболее подвержено эв- трофикационным процессам. Известно, что эвтрофикация (др. греч. εὐτροφία – хорошее питание) – это обогащение рек, озeр и морей биогенами, сопровождающееся повышением продуктивности растительности в водоeмах. Эвтрофикация может быть как результатом естественного старения водоeма, так и антропогенных воздействий. К основным химическим биогенным элементам, спо- собствующим эвтрофикации водоема, можно отнести фосфор и азот. Рельеф и течения моря находятся в сильной взаимосвязи друг с другом и оказывают влияние на биогенный режим. 2. Трехмерная математическая модель эвтрофикации мелководного во- доема При построении модели эвтрофикации вод (ЭВ) Азовского моря и Таганрогского залива ис- пользовались работы Матишова Г.Г., Ильичева В.Г., Якушева Е.В., Сухинова А.И. [1], Кру- киера Л.А., посвященные моделированию гидрохимических процессов. Учитывался тот факт, что при штилях и близких к ним ветровых ситуациях возникают анаэробные условия в придон- ных слоях Азовского моря. Восстановление поверхностного водонасыщенного ила влечет за со- бой высвобождение в раствор (кроме сероводорода) сульфатов, двухвалентного марганца и же- леза, органических соединений, аммония, силикатов и фосфатов. Расчетная область G (рис. 1) представляет собой замкнутый бассейн, ограниченный невозмущенной поверхностью водоема * Работа выполнена при частичной финансовой поддержке Задания № 2014/174 в рамках базовой части государственного задания Минобрнауки России, а также при частичной финансовой поддержке РФФИ по проектам № 15-01-08619, № 15-07-08626, № 15-07-08408. 320 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt 0 , дном H  H(x, y) и цилиндрической поверхностью  для 0  t  T0.   0  H  – кусочно-гладкая граница области G . x 0 y  H z Рис. 1. Схема расчетной области G Модель будет иметь вид: Si Si S   Si  Si,t  u v  w  w gi  i  i Si  i    i, (1) x y z z  z  где Si – концентрация i -й примеси, i  1,17 , индекс i указывает на вид субстанции, i  1,17 : 1 – сероводород H 2S  ; 2 – элементная сера S  ; 3 – сульфаты SO4  ; 4 – тиосуль- фаты (и сульфиты); 5 – общий органический азот N  ; 6 – аммоний NH4  (аммонийный азот); 7 – нитриты NO2  ; 8 – нитраты NO3  ; 9 – растворенный марганец DOMn  ; 10 – взвешенный марганец POMn  ; 11 – растворенный кислород O2  ; 12 – силикаты ( SiO3 – метасиликат; SiO4 – ортосиликат); 13 – фосфаты PO4  ; 14 – железо Fe2  15 – кремнекислота ( H2SiO3 – ме- такремневая; H2SiO4 – ортокремневая); 16 – фитопланктон; 17 – зоопланктон; u  u, v,w  – скорость водного потока; w g – скорость гравитационного осаждения i-й компоненты, если она i находится во взвешенном состоянии; i,i – коэффициенты турбулентного обмена соответ- ственно по горизонтальному и вертикальному направлениям;  i – химико-биологический ис- точник (сток) или член, описывающий агрегирование (слипание-разлипание), если соответству- ющая компонента является взвесью. К системе (1) добавим граничные условия: Si Si  0 на , если U n  0;  0 на , если U n  0; Si,z  0 на 0; Si,z  iSi на H , 2 n где i – коэффициент поглощения i -й компоненты донными отложениями. Необходимо также добавить начальные условия: Si t 0  Si0(x, y, z), i  1,17. (3) С помощью модели ЭВ (1) – (3) могут быть описаны процессы аммонификации, нитрифика- ции, нитратредукции (денитрификации), ассимиляции NH4 , окисления H2S , сульфатредукции, окисления и восстановления марганца, а также можно изучать механизм условий формирования заморов, прогнозировать изменения кислородного и биогенного режимов в мелководном водо- еме. При построении модели были параметризованы процессы биогеохимических циклов хими- ческих элементов, ответственных за трансформацию аэробных условий в анаэробные [2]. Проведено исследование трехмерной модели ЭВ Азовского моря вида (1) – (3), получены достаточные условия существования и единственности ее решения, сформулированные в виде теорем. Теорема 1. Пусть Si(x, y,z,t),  i  C 2(Цt)  C(Цt), где Цt  G  (0  t  T0); i  const  0 ; U  u, v, w  w gi  , i(z)  C 1(G); Sio  C(G), i  1,17 . Тогда при выпол- нении неравенств: max i, i   1 max  pi   0 для всех i  1,17 , где  i  pi Sj  Si   i 0 G G Si   Si  ,i  j ;  i  LSi  DSi  piSi , LSi   div Si Ui  , DSi  i Si  i , t z  z  321 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt 0   2 1 / L2x  1 / L2y  1 / L2z  ; Lx, Ly , Lz – пространственные максимальные размеры рас- четной области; модельная задача эвтрофикации мелководного водоема ЭВ вида (1) – (3) имеет решение. Теорема 2. Пусть Si(x, y,z,t),  i  C 2(Цt)  C(Цt), i  const  0 ; U ,  i(z)  C 1(G), i  1,17 . Тогда при выполнении неравенств: 2i 1 / L2x  1 / L2y   2vi / L2z  i для всех i  1,17 (где функции i определяются источниками загрязняющего вещества (ЗВ)) модель- ная задача ЭВ мелководного водоема вида (1) – (3) имеет единственное решение. Для дискрети- зации модели (1) – (3) построим в области решения задачи связную сетку h . hx, hy , hz – вектор- ные параметры, характеризующие плотность расположения узлов: h  xj  jhx, y k  khy ,zl  lhz; j  0, N x; k  0, N y;l  0, N z; N xhx  Lx, N y hy  Ly , N zhz  Lz  , h  h   ,   tn  n, n  0, Nt ,   где i, j, k – индексы по направлениям x, y,z ; N x, N y , N z – количество узлов по координатным направлениям. Расчетная область по пространственным направлениям x, y,z представляет собой объеди- нение параллелепипедов. Дискретизируем модель (1) – (3) с помощью разностной схемы с ве- сами: Sni j1,k,l  Sni j,k,l   Sni j11,k,l  (1   )Sni j 1,k,l  Sni j11,k,l  (1   )Sni j 1,k,l   ujn 1 ,k,l  ujn 1 ,k,l     2 2hx 2 2hx      Sni j1,k 1,l  (1   )Sni j,k 1,l  Sni j1,k 1,l  (1   )Sni j,k 1,l    vjn,k  1 ,l  vjn,k  1 ,l    2 2hy 2 2hy     Si j,k,l 1  (1   )Si j,k,l 1   n 1  n   w jn,k,l  1  w gn i j,k,l  1   2 2 2hz   Sni j1,k,l 1  (1   )Sni j,k,l 1    w j,k,l  1  w g i j,k,l  1 n 2 n  2 2hz     (4)       i2  Sni j1 1,k,l  (1   )Sni j 1,k,l  2 Sni j1,k,l  2(1   )Sni j,k,l  i2  Sni j11,k,l  (1   )Sni j1,k,l  hx hx        2i  Sni j1,k  1,l  (1   )Sni j,k  1,l  2 Sni j1,k,l  2(1   )Sni j,k,l  2i  Sni j1,k 1,l  (1   )Sni j,k 1,l  hy hy  1  n  Sni j1,k,l 1  (1   )Sni j,k,l 1   Sni j1,k,l  (1   )Sni j,k,l    i j,k,l  1   hz  2 hz   1  n  Sni j1,k,l  (1   )Sni j,k,l   Sni j1,k,l 1  (1   )Sni j,k,l 1    i j,k,l  1    i  0, hz  2 hz   1  j  N x  1,1  k  N y  1,1  l  N z  1,i  1,17. К системе (4) добавим аппроксимированные начальные и граничные условия. Исследуем устойчивость разностной схемы (4) для этого запишем ее в каноническом виде: B1Sni j1 1,k,l  B2Sni j1,k  1,l  B3Sni j1,k,l  1  A1Sni j1,k,l  B4Sni j11,k,l  B5Sni j1,k 1,l  B6Sni j1,k,l 1  A2Sni j,k,l  B7Sni j 1,k,l  B8Sni j,k 1,l  (5) B9Sni j,k,l  1  B10Sni j 1,k,l  B11Sni j,k 1,l  B12Sni j,k,l 1   i , 1  j  N x  1,1  k  N y  1,1  l  N z  1,0  n  Nt  1, i  1,17. Достаточное условие устойчивости и монотонности модели (4) определяется на основе принципа максимума А. А. Самарского при следующих ограничениях: hx  2i / u n C  , hy  2i / v n C  , hz  2i / w n  w gn  ,i  1,17. (6)     h h i C h  322 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt Исследование погрешности аппроксимации разработанной схемы с весами вида (4) пока- зало, что при   0.5 (схема Кранка-Николсон) она имеет наибольший порядок точности по временной и пространственным переменным: O  2  h  2  , где h  h  h  h . 2 x 2 y 2 z Входными данными предложенной выше модели эвтрофикации вод мелководного водоема вида (1) – (3) является поле вектора скорости водного потока, рассчитанное по моделям Сухи- нова А.И., Чистякова А.Е. [3–6]. 3. Метод решения сеточных уравнений Представим в виде операторного уравнения дискретную модель вида (4): Au  f (7) с невырожденным оператором A , заданным в вещественном гильбертовом пространстве H . Рас- смотрим неявную двухслойную итерационную схему вида yk 1  yk B  Ayk  f, k  0,1, (8) k с произвольным начальным приближением y0  H и невырожденным оператором B [7]. Любой двухслойный итерационный метод, построенный на основе схемы (7), характеризуется операто- рами A и B , энергетическим пространством H D , в котором доказывается сходимость метода, и набором итерационных параметров  k . Основным вопросом теории итерационных методов яв- ляется вопрос об оптимальном выборе параметра  k [8–10]. В двухслойных итерационных методах вариационного типа для вычисления параметров  k не требуется никакой априорной информации об операторах схемы (8) (кроме условий общего вида A  A*  0, DB 1A   DB 1A и т.д.), и построение этих методов основано на следующем * принципе: если задано приближение y k , а yk 1 находится по схеме (8), то итерационный пара- метр  k 1 выбирается из условия минимума в H D нормы погрешности zk 1  yk 1  u , где u – решение уравнения (7). Последовательность y k , построенная по формуле (8), в которой пара- метры  k выбираются из указанного выше условия, является минимизирующей для квадратич- ного функционала вида I  y   D  y  u , y  u . Этот функционал (в силу положительной определенности оператора D ) ограничен снизу, достигает минимума, равного нулю, при реше- нии уравнения (7), т.е. при y  u . Выбор параметра  k 1 из указанного условия обеспечивает локальную минимизацию функционала I  y  при переходе от y k к yk 1 , т.е. за один итераци- онный шаг. В случае явной схемы  B  E  переход от y k к yk 1 осуществляется по формуле yk 1  yk   k 1rk, rk  Ayk  f . Для самосопряженного положительно определенного оператора A переход от y k к yk 1 происходит по направлению rk , которое совпадает с направлением антиградиента для функци- онала  A  y  u , y  u в точке y k . По направлению антиградиента происходит наибольшее убывание значения функционала. Параметр  k 1 будем находить из условия минимума в H D нормы погрешности zk 1  yk 1  u [11]. Формула для вычисления итерационного параметра  k 1 находилась в предположении, что оператор A невырожден. Выпишем сначала уравнение для погрешности: zk  yk  u , k  0,1, . Подставляя yk   k  u в схему (8), получим zk  1  E   kB 1A  zk , 1  k  0,1, , z0  y0  u . Замена zk  D 2xk позволяет перейти к уравнению, содержащему только один оператор: 1 1 xk  1  Sk  1xk , Sk  E   kC, C  D 2 DB 1A  D 2.   (9) 323 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt Поставленную выше задачу о выборе параметра  k 1 можно сформулировать, используя ра- венство zk D  xk   D  норма в H D ,   норма в H  . Следующим образом: выбрать па- раметр  k 1 из условия минимума нормы xk 1 в пространстве H . Вычислим норму:    E   k  1C  xk, E   k  1C  xk   xk  2 k  1 Cxk, xk    k2 1 Cxk,Cxk   2 2 xk  1 (10)  Cxk, Cxk   k  1  Cxk, xk  / Cxk, Cxk    xk  Cxk, xk  / Cxk, Cxk  . 2 2 2 Так как оператор A не вырожден, то не вырожден и оператор C . Поэтому для любого xk имеем Cxk,Cxk   0 и минимум нормы xk  1 достигается при  k 1  Cxk, xk  / Cxk,Cxk  . (11) Подставляя (11) в (10), получим: xk  1  k  1 xk , (12) где k2 1  1  Cxk, xk  / Cxk, Cxk   xk, xk  . 2 (13) Формула (11) определяет оптимальное значение итерационного параметра  k 1 . Подставляя 1  xk  D 2zk в (11) получим:  k  1  DB 1Azk, zk xk  / DB 1Azk, B 1Azk  , k  0,1, , (14) Учитывая, что Azk  Ayk  Au  Ayk  f  rk – невязка, а B rk  k – поправка, фор- 1 мулу для параметра  k 1 можно записать в следующем виде:  k 1  Dk,zk  / Dk, k  , k  0,1, , (15) а итерационная схема (8) – в виде явной формулы для вычисления yk 1 : yk  1  yk   k  1k , k  0,1,... . (16) Алгоритм, реализующий построенный метод, можно кратко описать следующим образом: 1) по заданному y k вычисляется невязка rk  Ayk  f; 2) решается уравнение для поправки Bk  rk; 3) по формуле (15) вычисляется параметр  k 1. 4) по формуле (16) находится новое приближение yk 1. Рассмотрим частные случаи двухслойных градиентных методов, которые мы будем исполь- зовать для решения модельной задачи эвтрофикации вод мелководного водоема. Каждый кон- кретный метод определяется выбором оператора D и имеет свою область применимости. Опера- тор D будем выбирать так, чтобы в формулу (15) для итерационного параметра  k 1 входили только известные в процессе итераций величины. Если оператор A самосопряжен и положительно определен в H , то для решения (8) можно использовать метод скорейшего спуска (МСС). Если оператор A несамосопряженный и невыро- жденный, а оператор B *A положительно определен, то можно использовать метод минимальных невязок (ММН) [12]. Пусть оператор A самосопряжен и положительно определен в H . Для метода скорейшего спуска (МСС) D  A . Оператор B должен быть положительно определен в H . Учитывая соот- ношения Azk  Axk  f  rk и A  A * из (15) получим формулу для итерационного пара- метра  k 1 в неявном методе скорейшего спуска:  k 1  rk, k  /  Ak, Ak  , k  0,1, . (17) Для случая явной двухслойной схемы (11) B  E получим k  B 1rk  rk и формула для  k 1 примет вид:  k 1  rk,rk  /  Ark, Ark  , k  0,1, . (18) При решении модельной задачи (1) – (3) методом скорейшего спуска по невязке параметр  k 1 рассчитывается по формуле (18), а метод скорейшего спуска по поправке реализуется с по- мощью расчетной формулы (17). 324 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt Метод минимальных поправок (ММП) можно применять для решения модельной задачи (1) – (3) в случае любого несамосопряженного невырожденного оператора A , кроме того требуется положительная определенность оператора B * A . Для метода минимальных невязок D  A *A . Формула для итерационного параметра  k 1 в методе минимальных поправок имеет вид:  k 1  (Ak, rk)/ (Ak, Ak), k  0,1, . (19) В случае явной схемы (8) B  E  требуется положительная определенность оператора A, а формула для  k 1 примет вид:  k 1   Ark,rk  /  Ark, Ark  , k  0,1, . (20) Был проведен анализ эффективности описанных выше градиентных методов вариационного типа на основе численного решения модельной задачи эвтрофикации вод мелководного водоема вида (1) – (3). В таблице 1 приведены результаты сравнения скоростей сходимости двухслойных градиентных методов вариационного типа, используемых для решения задачи эвтрофикации вод мелководного водоема вида (1) – (3): минимальных невязок (ММН), минимальных поправок (ММП), скорейшего спуска по невязке (МССН), скорейшего спуска по поправке (МССП). Таблица 1. Сравнение скоростей сходимости методов вариационного типа ММП n 1 2 3 4 p 1 IS  48,IX  44 IS  47,IX  43 IS  47,IX  43 IS  48, IX  42 2 IS  53, IX  48 IS  52, IX  47 IS  53, IX  46 IS  53, IX  45 3 IS  21, IX  21 IS  22, IX  20 IS  22, IX  19 IS  21, IX  19 4 IS  17,IX  15 IS  14, IX  13 ММН 1 IS  83, IX  65 IS  77, IX  64 IS  76, IX  63 IS  76, IX  62 2 IS  87,IX  73 IS  87,IX  73 IS  86, IX  72 IS  85, IX  71 3 IS  51, IX  51 IS  50, IX  46 IS  49, IX  43 IS  48,IX  40 4 IS  40, IX  38 IS  38, IX  35 МССП 1 IS  69, IX  62 IS  67,IX  61 IS  67,IX  60 IS  67,IX  59 2 IS  71, IX  68 IS  74, IX  67 IS  74, IX  67 IS  74, IX  66 3 IS  35, IX  32 IS  35, IX  30 IS  35, IX  29 IS  34, IX  28 4 IS  24, IX  23 IX  21, IX  20 МССН 1 IS  68,IX  64 IS  79, IX  67 IS  78, IX  66 IS  78, IX  65 2 IS  86, IX  74 IS  87,IX  73 IS  86, IX  73 IS  86, IX  72 3 IS  53, IS  49 IS  52, IX  45 IS  51, IX  42 IS  49, IX  41 4 IS  40, IX  38 IS  36, IX  35 В таблице 1 используются следующие обозначения: p – номер итерации; n – номер времен- ного слоя; I – количество итераций, необходимое для сходимости метода решения СЛАУ для уравнения, описывающего изменение концентрации  ,   S, X , где S  S5, X  S16 . Метод минимальных поправок (ММП) был выбран в качестве основного метода в виду его наибольшей скорости сходимости, согласно данным, приведенным в таблице 1. 325 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt 4. Параллельный вариант метода решения сеточных уравнений Для геометрического разбиения расчетной области с целью равномерной загрузки вычисли- телей (процессоров) МВС использовался метод k-means, основанный на минимизации функцио- нала суммарной выборочной дисперсии разброса элементов (узлов расчетной сетки) относи- тельно центра тяжести подобластей: Q  Q(3). Пусть Xi – множество расчетных узлов сетки, входящих в i  ю подобласть, i  1,..., m , m  заданное количество подобластей. 1 1 Q(3)    d 2(x,ci)  min , где ci   x – центр подобласти Xi , а d(x,ci)– рассто- i Xi xXi Xi xXi яние между расчетным узлом сетки x и центром подобласти ci в Евклидовой метрике. Метод k- means является сходящимся только тогда, когда все подобласти будут примерно равны. Результат работы метода k-means для модельных двумерной и трехмерной областей пред- ставлен на рис. 2. Рис. 2. Результат работы алгоритма k-means для разбиения модельной двумерной области на 9, 38, 150 подобластей; для трехмерной на 6 и 10 Опишем работу максиминного алгоритма. В качестве центров подобластей выбирает расчет- ные узлы сетки следующим образом: 1) первый центр – первый расчетный узел области; 2) второй центр находится в расчетном узле сетки, расположенном на максимальном расстоя- нии от первого центра; 3) если количество подобластей больше 3-х, то каждый следующий центр находится на макси- мальном удалении от ближайшего центра. Опишем алгоритм работы k-means. 1) Выбираются начальные центры подобластей при помощи максиминного алгоритма. 2) Все расчетные узлы разбиваются на m клеток Вороного по методу ближайшего соседа, т. е. текущий расчетный узел сетки x  Xc , где Xc – подобласть выбирается из условия x  sc  min x  si , где sc – центр области Xc . 1 i  m 1 3) Рассчитываются новые центры по формуле s(ck  1)   x. X(ik) xX(ik) 4) Проверяется условие остановки s(ck  1)  s(ck) для всех k  1,..., m . Если условие остановки не выполняется, то осуществляется переход на пункт 2 алгоритма. Для организации обмена данными в вычислительном процессе требуется найти все точки, лежащие на границе каждой подобласти. Для этой цели использовался алгоритм Джарвиса (за- дача построения выпуклой оболочки). Был сформирован список соседних подобластей для каж- дой подобласти и разработан алгоритм пересылки данных между подобластями [13,14]. 326 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt При решении СЛАУ методом минимальных поправок для расчета итерационного параметра  использовался метод сдваивания. Синхронизация алгоритма решения задачи (1) – (3) требуется только в ММП при переходе на следующую итерацию. 5. Описание программного комплекса Для решения модельной задачи эвтрофикации мелководного водоема (1) – (3) был создан комплекс прикладных программ, позволяющий производить расчеты концентраций ЗВ, фито- и зоопланктона в областях сложной формы (Азовское море и Таганрогский залив) [15]. С помощью созданного комплекса программ можно осуществлять:  совершенствование и внедрение системы комплексного рыбохозяйственного мониторинга в водоемах (наблюдение, оценка и прогноз состояния режима экосистем, кормовой базы и за- пасов промысловых объектов);  разработку, согласование предложений и мероприятий по обеспечению оптимального ре- жима, сохранения биоразнообразия промысловых ресурсов, экосистем мелководных водое- мов;  совершенствование методологии природоохранных исследований, разработка новых, апро- бация и внедрение перспективных методов изучения состояния водных экосистем и отдель- ных компонентов;  разработку и совершенствование методов диагностики токсического воздействия ЗВ на гид- робионты, в том числе ранней и дифференциальной диагностики токсикоза, а также поиск средств антидотной защиты водных экосистем;  организацию и проведение исследований по выявлению тенденций и закономерностей изме- нения состояния водных экосистем под воздействием антропогенных факторов, разработка предложений и мероприятий по снижению и предупреждению таких воздействий;  оценку ущербов рыбному хозяйству, наносимых разными видами хозяйственной деятельно- сти, разработку предложений по предотвращению, уменьшению и адекватной компенсации ущербов. Поля скоростей водного потока, рассчитанные в работах [16], относятся к входным данным для моделей гидробиологических процессов, описанных в первой главе. Для математического моделирования гидробиологических и гидродинамических процессов в трехмерной области сложной формы – Азовское море и Таганрогский залив использовались последовательно сгуща- ющиеся прямоугольные сетки размерностями: 251  351  15 , 502  702  30 , 1004  1404  60 , . Начальное распределение загрязняющих веществ и планктона было учтено в форме, соот- ветствующей пространственно-временным масштабам моделируемых процессов. Реализован- ный алгоритм численного решения позволяет свободно варьировать граничные условия, вид управляющих функций и значения соответствующих параметров [17]. Понимание механизма функционирования системы и знание еe основных характеристик позволили для преодоления трудностей при настройке программы использовать феноменологический подход. Эффектив- ность такого подхода достаточно высока ещe и потому, что поведение системы часто определя- ется точностью не отдельных параметров, а их соотношений. Калибровка и верификация разработанной модели эвтрофикации вод мелководного водоема проводились на основе экологических данных по Азовскому морю, полученных в ходе научно- исследовательских экспедиций, проводимых учеными ЮФУ, начиная с 2000 года. В ходе иссле- дований экватории Азовского моря изучались: виды и концентрации основных загрязняющих воды Азовского моря веществ; пространственные распределения солености и температуры; кис- лородный режим; видовой состав фито- и зоопланктона; механизмы возникновения и развития заморов в центрально-восточной части водоема. Обработка экспедиционных данных заключалась в оцифровке, пересчете в стандартные еди- ницы, классификации для использования в разных модельных задачах гидробиологии моря. На рис. 4 приведена схема разработанного программного комплекса. 327 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt Рис. 4. Схема программного комплекса Программный комплекс включает: блок управления, базы океанологических и метеорологи- ческих данных, системы интерфейсов, системы ввода – вывода и визуализации. Программный продукт обладает удобным интерфейсом, что позволяет вводить нужную информацию в диало- говом режиме. Построенный программный комплекс имеет универсальный характер и привязка его к условиям конкретных объектов и районов осуществляется, как правило, на уровне входной информации. Это значит, что для практического использования модулей требуется создание спе- циальной информационной базы, содержащей сведения о физико-географических и климатиче- ских условиях исследуемых объектов, о параметрах, определяющих источники примесей. При разработке программного комплекса использовался язык высокого уровня С++. Были разработаны две версии параллельного алгоритма решения модельной задачи эвтрофикации вод мелководного водоема. Первая предназначена для ЭВМ с операционной системой Windows, ос- нована на технологии создания в этой системе дополнительных потоков. Вторая может исполь- зоваться для кластерных систем и основана на технологии передачи сообщений (MPI). 6. Результаты численных экспериментов При моделировании пространственно-неоднородных процессов эвтрофикация вод Азов- ского моря (1) – (3) учитывалась внешняя периодичность, приводящая к усложнению системы. При этом колебания плотности планктона могут быть столь велики, что не могут быть объяснены случайными флуктуациями, и визуальная картина такова, что сравнительно небольшие площади высокой плотности ("пятна", "облака") разделены пространствами с низкими плотностями, ино- гда не фиксируемыми стандартными методами наблюдений. Особенно ярко это явление выра- жено в тех местах водоема, для которых характерна потребность в биогенных элементах. При моделировании процессов эвтрофикации учитывался вегетационный период фитопланктона. Диффузные процессы действуют в направлении сглаживания пространственного распреде- ления и рассеивания "пятен". Одна из попыток объяснения парадокса стабильности "пятен" с помощью численных экспериментов заключается в предположении об активном передвижении гетеротрофных организмов (зоопланктона и рыб) в направлении градиента "пищи", что обеспе- чивает закрепление пространственной неоднородности биогенных веществ в водной среде. Устойчивая неоднородность пространственного распределения может быть, например, обуслов- лена диффузионными процессами и наличием у фитопланктона механизма эктокринного регу- лирования, т.е. регулирования скорости роста посредством выделения в среду биологически ак- тивных метаболитов. На рисунке 5. показаны результаты расчета концентрации загрязняющего биогенного веще- ства для модели динамики вредоносной водоросли (начальное распределение полей течений при 328 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt северном ветре). Приведенные ниже рисунки отражают влияние структур течений водного по- тока в Азовском море на распределение загрязняющего биогенного вещества и фитопланктона. Рис. 5. Распределение концентрации загрязняющих веществ в различные моменты времени Результаты моделирования динамики фитопланктона в Азовском море представлены на рис. 6. ( N – номер итерации, начальное распределение полей течений водного потока при северном направлении ветра). Рис. 6. Распределение концентрации фитопланктона в различные моменты времени Белым цветом выделены максимальные значения концентраций биогенного вещества (азота) и фитопланктона черным – минимальные. Критерием проверки адекватности построенных моделей гидробиологии мелководного во- доема служила оценка погрешности моделирования с одновременным учетом натурных данных по имеющимся n замерам, которая вычислялась по формуле:  S  n n S 2   k nat  Sk 2 k nat , где Sk nat – значение концентрации загрязняющей при- k 1 k 1 меси, полученное с помощью натурных экспедиционных измерений; Sk – значение, рассчитан- ное с помощью модели (1) – (3). Рассчитанные при различных ветровых ситуациях концентрации загрязняющих веществ и планктона принимались к рассмотрению, если относительная погрешность не превышала 30%. 329 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt 7. Заключение С помощью экспедиционных исследований проведена первичная верификация модели эко- системы Азовского моря. Реализована задача моделирования и прогноза состояния водной эко- системы Азовского моря в условиях антропогенного воздействия и всестороннего изучения уни- кального водного объекта, который в силу мелководности в большей степени подвержен антро- погенному влиянию. Создан программный комплекс, объединяющий математические модели и базы данных, с помощью которого изучены условия, при которых мелководные водоемы под- вергаются эвтрофитрованию. Отличительными особенностями разрабатываемых алгоритмов, реализующих поставленные гидробиологические модельные задачи, являются: высокая производительность, достоверность и точность получаемых результатов. Высокая производительность достигается за счет использо- вания эффективных численных методов решения сеточных уравнений, ориентированных для применения на параллельных вычислительных системах в реальном и ускоренном масштабах времени. Достоверность достигается за счет учета определяющих физических факторов, таких как: сила Кориолиса, турбулентный обмен, сложная геометрия дна и береговой линии, испаре- ние, стоки рек, динамическое перестроение расчетной области, ветровые напряжения и трение о дно, а также за счет учета отклонения значения поля давления от гидростатического приближе- ния. Точность достигается применением подробных расчетных сеток, учитывающих степень «за- полненности» расчетных ячеек, а также отсутствием неконсервативных диссипативных слагае- мых и нефизичных источников (стоков), возникающих в результате конечно-разностных аппрок- симаций. Было проведено сравнение работы созданного программного комплекса, реализующего раз- работанные сценарии развития экологической обстановки в Азовском море с использованием численной реализации модельной задачи эвтрофикации вод мелководного водоема, с подобными работами в области математического моделирования гидробиологических процессов. Анализ показал, что в результате разработки программного комплекса удалось повысить точность прогнозов изменения концентраций загрязняющих веществ, фито- и зоопланктона в мелководном водоеме на 10 – 15% в зависимости от решаемой модельной задачи водной эколо- гии. Литература 1. Якушев Е.В., Сухинов А.И. и др. Комплексные океанологические исследования Азовского моря в 28-м рейсе научно-исследовательского судна «Акванавт» // Океанология. 2003. Т. 43, №1. С. 44–53. 2. Сухинов А.И., Никитина А.В., Чистяков А.Е. Моделирование сценария биологической реа- билитации Азовского моря // Математическое моделирование. 2012. Т. 24, №9. С. 3–21. 3. Сухинов А.И., Чистяков А.Е., Алексеенко Е.В. Численная реализация трехмерной модели гидродинамики для мелководных водоемов на супервычислительной системе // Математиче- ское моделирование. 2011. Т. 23, № 3. С.3–21. 4. Сухинов А.И., Чистяков А.Е. Параллельная реализация трехмерной модели гидродинамики мелководных водоемов на супервычислительной системе // Вычислительные методы и про- граммирование: Новые вычислительные технологии. 2012. Т. 13. С. 290–297. 5. Белоцерковский О. М. Турбулентность: новые подходы. М.: Наука, 2003. 286 c. 6. Самарский А.А. Теория разностных схем. М.: Наука, 1989. 616 c. 7. Самарский А.А., Николаев Е.С. Методы решения сеточных уравнений. М.: Наука, 1978. 592 c. 8. Коновалов А.Н. К теории попеременно-треугольного итерационного метода // Сибирский ма- тематический журнал. 2002. 43:3. С. 552–572. 330 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt 9. Сухинов А.И., Чистяков А.Е. Адаптивный модифицированный попеременно-треугольный итерационный метод для решения сеточных уравнений с несамосопряженным оператором // Математическое моделирование. 2012. Т. 24, №1. С. 3–20. 10. Чистяков А.Е. Теоретические оценки ускорения и эффективности параллельной реализации ПТМ скорейшего спуска // Известия ЮФУ. Технические науки. 2010. №6(107). С 237–249. 11. Чистяков А.Е., Хачунц Д.С., Никитина А.В., Проценко Е.А., Кузнецова И.Ю. Библиотека па- раллельных итерационных методов решателей СЛАУ для задачи конвекции-диффузии на ос- нове декомпозиции по одному пространственному направлению // Современные проблемы науки и образования. 2015. № 1. URL: http://www.science-education.ru /121-19510 (дата обра- щения 4.06.2015). 12. Никитина А.В. Численное решение задачи динамики токсичных водорослей в Таганрогском заливе // Известия ЮФУ. Технические науки. 2010. №6(107). С 113–116. 13. Никитина А.В. Модели биологической кинетики, стабилизирующие экологическую систему таганрогского залива // Известия ЮФУ. Технические науки. 2009. № 8 (97). C. 130–134. 14. Сухинов А.И., Чистяков А.Е., Бондаренко Ю.С. Оценка погрешности решения уравнения диффузии на основе схем с весами // Известия ЮФУ. Технические науки. 2011. №8 (121). С. 6–13. 15. Сухинов А.И., Чистяков А.Е., Семенякина А.А., Никитина А.В. Параллельная реализация задач транспорта веществ и восстановления донной поверхности на основе схем повышен- ного порядка точности // Вычислительные методы и программирование: новые вычисли- тельные технологии. 2015. Т. 16. C. 256–267. 16. Никитина А.В., Руднева Т.В., Камышникова Т.В., Бокарева Т.А., Дурягина В.В. К вопросу о формировании заморных зон в восточной части Азовского моря // Современные проблемы науки и образования. 2015. №1. URL: http://www.science-education.ru/121-19509 (дата обра- щения 4.06.2015). 17. Никитина А.В., Пучкин М.В., Семенов И.С., Сухинов А.И., Угольницкий Г.А., Усов А.Б., Чистяков А.Е. Дифференциально-игровая модель предотвращения заморов в мелководных водоемах // Управление большими системами, Выпуск 55. М.: ИПУ РАН. 2015. C. 343–361. 331 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt Mathematical modeling of eutrophication processes in shallow waters on multiprocessor computer system A.I. Sukhinov1, A.V. Nikitina2, A.E. Chistyakov2, I.S. Semenov2, A.A. Semenyakina2, D.S. Khachunts2 Don State Technical University1, Kalyaev Scientific Research Institute of multiprocessor computing systems of South Federal University2 The paper covered the development of solution methods for model of eutrophication of shal- low water in view the movement of the water flow, microturbulent diffusion, gravitational settling, spatially uneven distribution of temperature and salinity, and pollutants of nutrients, oxygen, phyto- and zooplankton, etc. Shallow waters of the Azov Sea and Taganrog Bay were selected as the object of simulation. Numerical solution of the problem is based on gradient method of variation type – the method of minimal corrections. The acceleration and efficiency of calculations is achieved with using multiprocessor computer system of Southern Federal University. The solution of the problem of water ecology will allow to predict changes of water quality of shallow reservoirs, and to study the mechanisms of formation of zones with low oxygen content. Keywords: mathematical model, eutrophication, minimum bet, wok, parallel computing, Azov Sea. References 1. Yakushev E.V., Sukhinov A.I. Complex Oceanographic research of the Sea of Azov in the 28-th flight of the research vessel "Aquanaut" // Oceanology. 2003. Vol. 43, No. 1. P. 44–53. 2. Sukhinov A.I., Nikitina A.V., Chistyakov A.E. Simulation scenario of biological rehabilitation of the Azov Sea // Mathematical modeling. 2012. Vol. 24, No. 9. P. 3–21. 3. Sukhinov A.I., Chistyakov A.E., Alekseenko E.V. Numerical implementation three-dimensional model of hydrodynamics for shallow water basins on superficialities system // Mathematical mod- eling. 2011. Vol. 23, No. 3. P. 3–21. 4. Sukhinov A.I., Chistyakov A.E. Parallel implementation of a three-dimensional model of hydrody- namics of shallow water bodies on superficialities system // Computational methods and program- ming: New computing technologies. 2012. T. 13. P. 290–297. 5. Belotserkovsky O.M. Turbulence: new approaches. M.: Science, 2003. 286 p. 6. Samarskii A.A. Theory of difference schemes. M.: Science, 1989. 616 p. 7. Samarskii A.A., Nikolaev E.S. Methods of solving grid equations. M.: Science, 1978. 592 p. 8. Konovalov A.N. The theory of the alternating-triangular iterative method // Siberian mathematical journal. 2002. 43:3. P. 552–572. 9. 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. P. 3–20. 10. Chistyakov A.E. Theoretical estimation of speed up and efficiency of parallel implementation of the steepest descent PTM // Izvestiya SFEDU. Technical Sciences. 2010. No. 6(107). P. 237–249. 332 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt 11. Chistyakov A.E., Hachunts D.S., Nikitina A.V., Protsenko E.A., Kuznetsova I.Yu. A library of parallel iterative methods solvers of linear algebraic equation for the problem of convection-diffu- sion-based decomposition in one spatial direction // Modern problems of science and education. 2015. No. 1. URL: http://www.science-education.ru /121-19510 (accessed: 4.06.2015). 12. Nikitina A.V. Numerical solution of problems of dynamics of toxic algae in the Taganrog Bay // Izvestiya SFedU. Engineering sciences. 2010. No. 6(107). P. 113–116. 13. Nikitina A.V. modelling of biological kinetics, stabilizing the ecological system of the Gulf of Ta- ganrog // Izvestiya SFedU. Engineering sciences. 2009. No. 8 (97). P. 130–134. 14. Sukhinov A.I., Chistyakov A.E., Bondarenko Yu.S. Error estimate of the solution of the diffusion equation on the basis of the schemes with weights // Izvestiya SFedU. Engineering sciences. 2011. No. 8 (121). P. 6–13. 15. Sukhinov A.I., Chistyakov A.E., Semenikhina A.A., Nikitina A.V. Parallel realization of the tasks of the transport of substances and the recovery of the bottom surface on the basis of the schemes of high order of accuracy // Computational methods and programming: new computing technolo- gies. 2015. T. 16. P. 256–267. 16. Nikitina A.V., Rudneva T.V., Kamyshnikova T.V., Bokareva T.A., Duryagina V.V. To the for- mation of kill zones in the Eastern part of Azov Sea coast // Modern problems of science and edu- cation. 2015. No. 1. URL: http://www.science-education.ru/121-19509 (accessed: 4.06.2015). 17. Nikitina A.V., Puchkin M.V., Semenov I.S., Sukhinov A.I., Ougolnitsky G.A., Usov A.B., Chisty- akov A.E. Differential-game model of prevention of fish kill in shallow ponds // Managing large systems, Issue 55. M.: IPU Russian Academy of Sciences. 2015. P. 343–361. 333