CC BY
65
УДК 519.6
Геоинформатика
DOI: 10.14529/cmse160303
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ПРОЦЕССОВ ЭВТРОФИКАЦИИ В МЕЛКОВОДНЫХ ВОДОЕМАХ НА МНОГОПРОЦЕССОРНОЙ ВЫЧИСЛИТЕЛЬНОЙ
СИСТЕМЕ
© 2016 г. А.И. Сухинов1, А.В. Никитина2, А.Е. Чистяков2,
3 3 2
И.С. Семенов , А.А. Семенякина , Д.С. Хачунц
1
Донской государственный технический университет (344000 Ростов-на-Дону, пл. Гагарина, д. 1), 2Научно-исследовательский институт многопроцессорных вычислительных систем имени академика А.В. Каляева Южного федерального университета (347928 Таганрог, ул. Чехова, д. 2), ОАО «Научно-исследовательский центр суперЭВМ и нейрокомпьютеров» (347900 Таганрог, пер. Итальянский, д. 106) E-mail: sukhm@gmail.com, nikitina.vm@gmail.com, cheese_05@mail.ru, flanker555@yandex.ru, j.a.s.s.y@mail.ru, diana-hachunts@mail.ru Поступила в редакцию: 25.02.2016
Работа посвящена разработке методов решения модельной задачи эвтрофикации вод мелководного водоема, учитывающей движение водного потока, микротурбулентную диффузию, гравитационное оседание, пространственно-неравномерное распределение температуры и солености, а также загрязняющих биогенных веществ, кислорода, фито- и зоопланктона и др. В качестве объекта моделирования были выбраны мелководные водоемы — Азовское море и Таганрогский залив. Численное решение задачи основывается на градиентном методе вариационного типа — методе минимальных поправок. Ускорение и эффективность расчетов достигается при использовании многопроцессорной вычислительной системы Южного федерального университета. Решение поставленной задачи водной экологии позволит прогнозировать изменения качества вод мелководных водоемов, а также изучать механизмы формирования в них зон с пониженным содержанием кислорода.
Ключевые слова: математическая .модель, эвтрофикация, метод минимальных поправок, параллельные вычисления, Азовское море.
ОБРАЗЕЦ ЦИТИРОВАНИЯ
Сухинов А.И., Никитина А.В., Чистяков А.Е., Семенов И.С., Семенякина А.А., Хачунц Д.С. Математическое моделирование процессов эвтрофикации в мелководных водоемах на многопроцессорной вычислительной системе // Вестник ЮУрГУ. Серия: Вычислительная математика и информатика. 2016. Т. 5, № 3. С. 36-53. DOI: 10.14529/cmse160303.
Введение
Изменение природно-климатических условий и антропогенное воздействие приводят к эвтрофикации вод мелководных водоемов, приводящей к значительному росту популяций фитопланктона, многие виды которого являются токсичными, вызывают тяжелые заболевания у людей. Процессы эвтрофикации влияют на качество вод, состояние рыб-
Статья рекомендована к публикации программным комитетом Международной научной конференции «Параллельные вычислительные технологии - 2016».
ных запасов водоемов, в связи с этим разработка программно-алгоритмического инструментария имеет не только научное, но и важное народнохозяйственное значение, она позволит осуществлять предсказательное моделирование взаимосвязанных процессов эв-трофикации на основе натурных экспериментов, полученных в ходе экспедиционных работ в акватории мелководных водоемов.
Целью работы являлось построение и численная реализация математической модели эвтрофикации вод мелководного водоема, обладающей предсказательной ценностью. В соответствии с поставленной целью решены следующие задачи: изучена предложенная модель эвтрофикации вод мелководного водоема, получены достаточные условия существования и единственности ее решения; разработан и исследован дискретный аналог этой модели. В работе описана построенная библиотека двухслойных методов вариационного типа, предназначенных для решения сеточных уравнений с самосопряженными и несамосопряженными операторами, возникающими при дискретизации разработанной модели, на многопроцессорной вычислительной системе. Для равномерного распределения вычислительной работы между процессорами использовался к-теап алгоритм.
Предложенные модели и методы были применены при разработке программного комплекса, предназначенного для расчета значений концентраций планктона и загрязняющих примесей в Азовском море. Входными данными предложенной выше модели эвтрофикации вод мелководного водоема вида (1) - (3) является поле вектора скорости водной среды, рассчитанное по моделям Сухинова А.И., Чистякова А.Е. [3-6]. Статья организована следующим образом. В разделе 1 приводится разработанная математическая модель эвтрофикации вод мелководного водоема, а также ее исследование на непрерывном и дискретном уровне. Раздел 2 содержит описание разработанных методов вариационного типа для решения, возникающих в процессе дискретизации, сеточных уравнений. Раздел 3 посвящен описанию параллельного алгоритма реализации метода минимальных поправок. Раздел 4 содержит описание программного комплекса, численно реализующего модельную задачу эвтрофикации мелководного водоема. Раздел 5 посвящен описанию и анализу результатов численных экспериментов. В заключении приводятся итоги исследования.
1. Трехмерная математическая модель эвтрофикации мелководного водоема
Азовское море является внутренним морем расположено на востоке Европы. Его глубина не превышает 13,5 метров и является самым мелким морем в мире. Азовское море подвержено эвтрофикационным процессам в виду своей мелководности. Эвтрофи-кация (др. греч. еитроф^а — хорошее питание) — это процесс обогащения морей, озер и рек биогенными веществами, сопровождающийся повышением объемов растительности в водоемах. Эвтрофикация может возникать как результат антропогенного воздействия, так и естественного старения водоема. К химическим биогенным элементам, способствующим эвтрофикации водоема, относятся азот и фосфор. Морские течения, зависящие от рельефа водоема, оказывают влияние на биогенный режим.
При построении модели эвтрофикации вод (ЭВ) Азовского моря и Таганрогского залива использовались работы Матишова Г.Г., Ильичева В.Г., Якушева Е.В., Сухинова А.И. [1], Крукиера Л.А., посвященные математическому моделированию гидрохимических процессов. В данных моделях учтено, что при штилях и малых ветрах в придон-
полнении неравенств: max [m1,vi } - -0 max {|pz > о для всех i = 1,17 , где
ds_
P2 (Sj) Si + Wi ,1 * j ; Wi = LS_ — DS_ — P_S1 , LS2 = -t + div (SiU_)
DSz = —AS_ + d^ [ V, I ' Л0 = П2 (l / L2x + 1 / Ly + 1 / L) ; Lxf Ly, Lz — простран
ÖS,.
dz v dz
ственные максимальные размеры расчетной области; модельная задача эвтрофикации мелководного водоема ЭВ вида (1) - (3) имеет решение.
Теорема 2. Пусть S(xr y, z, t), y/1 e с2(цt) n c(ut), ju1 = const > о; и, v1(z) e с 1(G), 1 = 1, 17 . Тогда при выполнении неравенств: 2^i (l / l2x + 1 / L2y) + 2vi / Lz > p± для всех 1 = 1, 17 (где функции pi определяются источниками загрязняющего вещества (ЗВ)) модельная задача ЭВ мелководного водоема вида (1) - (3) имеет единственное решение. Для дискретизации модели (1) - (3) построим в области решения задачи связную сетку coh . hx, hy, hz — векторные параметры, характеризующие плотность расположения узлов:
Шь = {x = jhxr yk = khy, zi = lhz; j = 0^; k = 0^; l = 0^;
Nxhx = Lx' Nyhy = Ly' Nzhz = h } ' < = й X = К = T ^ = ,
где 1, j, k — индексы по направлениям x, y, z; Nx, N , Nz — количество узлов по координатным направлениям.
Расчетная область по пространственным направлениям x, y, z представляет собой объединение параллелепипедов. Дискретизируем модель (1) - (3) с помощью разностной
схемы с весами:
S(i)j,k,l - S(1)j,k,l
^+1+!,« + (1 - YS
j+y?k,i
2h
(i)j + 1,k,1 - — u,
?S(ni+' + (1 — Y)S(n I
^(i)j—1,k,1
j— 12,-m
2h
YSj+1,i + (1 — YS
j,k+12,i
2h,
(i)jk+1,1 — n
- — V -I /
—Ю,1
YS(i+1k—11 + (1 — Y)S(n I
J(i)j,k—1,1
2h
I W 1 / — w , x / \ J,k,1 +12 g (i) J,k,1 +12
YS_+ L,1+1 + (1 — Y)Sn
(i) j,k,1 +1
2h
— W w — w , .
\ J,k,1 — Ц g (i)j,k,1 —
' z
n +1 n
YSijk,i—1 + (1 — Y)S
(i) j,k,1 — 1
2h
h2 —i
h
У (4)
Н_+;+1Д,1 + (1 — Y)S(i)j+1,k,1 — ^j,^ — 2(1 — Y)Sni)j,k,1) — — (ys1))—w + (1 — Y)S(i)j—1,^,1)
2 Н_+1 k+1, 1 + (1 — Y)Sni)j, k+1 , 1 — 2YS(ni+j, k , 1 — 2(1 — y)S_)j, t, 1) — ^ (yS^+J,k—1, 1 + (1 — YS^j, k—1 , 1)
y hy
1 h
J(i)j,k + 1,1
YS(ni+j,k,1+1 + (1 — YS!) j,k,1+1 — YSij,k,1 — (1 — Y)Sni)j,k,1 Л
— (YSn+1
—i hi
К
Ys_+ j,k,1 + (1 — YSi) j,k,1 — Y^k!—1 — (1 — Y)S(
V
n \
(i)j,k,1 — 1
с1)j,k,1—
h
— W = 0,
1 < j < Ых - 1,1 < к < Ыу - 1,1 < 1 < - 1,1 = 1, 17.
К системе (4) добавим аппроксимированные начальные и граничные условия. Ис-
следуем устойчивость разностной схемы (4) для этого запишем ее в каноническом виде:
n
+
т
n
+
+
)
+
n
+
(i) j,k,1 +1
1
n
+
h
z
в „п + 1 + в „п + 1 + в „п + 1 + а ^ + 1 _ Д „п + 1 _
П1°(1) 1 + 1,к,1 П2°(1) ],к + 1,1 в3°(1 ) ],к,1 + 1 а1°(1 )],к,1 п4°(1) ]_1,к,1
_B5S(l)j,k _1,1 _ 1)1,к,1 _1 = _A2S(1 )],к,1 + B7S(1) 1 + 1,к,1 + БР^<1) ],к + 1,1 + ,_ч
п п п п - (5)
+Б^(1 )1,к,1 + 1 _ в^(1 )1_1,к,1 _ Б1^(1 )1,к _1,1 _ ) 1,к,1 _1 _ —1 ,
1 < 1 < ых _ 1, 1 < к < ыу _ 1, 1 < 1 < ыг _ 1, 0 < п < Nt _ 1,1 = 1, 17.
Достаточное условие монотонности и устойчивости модели (4) определяется на основе принципа максимума А. А. Самарского при следующих ограничениях:
щ < / и1с^), Щ < / vnILк),Щ < Ь / _ к )||с^) ,1 = ^ (6)
Исследование погрешности аппроксимации разработанной схемы с весами вида (4) показало, что при у = 0.5 (схема Кранка — Николсон) она имеет наибольший порядок
точности по временной и пространственным переменным: О (г2 + Щ|2), где
Щ = дЩ + Щ + Щ .
2. Метод решения сеточных уравнений
Дискретную модель (4) может быть представлена в виде операторного уравнения:
Аи = f (7)
с положительно определенным оператором А в гильбертовом пространстве н. Рассмотрим неявный двухслойный процесс вида
в Ук + 1 _ Ук + АУк = f, к = 0, 1, ... , ,
Г (8)
с неким начальным приближением у0 е н и оператором — предобуславливателем в [7]. Всякий двухслойный итерационный процесс (8), характеризуется операторами А и в, энергетическим пространством нв , в котором исследуется сходимость метода, и способом расчета итерационных параметров тк. К основным вопросам теории итерационных методов решения сеточных уравнений относится оптимальный выбор параметров тк [8-10].
Для вычисления параметров тк в методах вариационного типа не требуется априорная информация об операторах задачи (7) (кроме условий общего вида А = А* > 0, (вв_1А) = вв_1А ). Данные методы основаны на принципе: если задано значение ук, а
ук+1 находится по схеме (8), то параметр тк+1 рассчитывается из условия минимизации в нв энергетической нормы погрешности zk+1 = ук+1 _ и , где и — решение уравнения (7). Последовательность ук, рассчитанная на основе выражения (8), в которой значения тк вычисляются из приведенного выше условия, является минимизирующей для квадратичного функционала вида I (у) = (в (у _ и), у _ и). Данный функционал имеет ограничения снизу (в силу положительной определенности оператора в), и его значение достигает минимума, равного нулю при у = и. Выбор параметра тк+1 из приведенного условия гарантирует локальную минимизацию функционала I (у) при переходе от ук к ук+1. В случае задания оператора-предобуславливателя методом простой итерации (в = Е) переход от ук к ук+1 осуществляется согласно выражению
ук+1 = ук _ гк+1гк' гк = Аук _ f.
В случае положительно определенного самосопряженного оператора А переход от ук к ук+1 происходит в направлении _гк, которое противоположно направлению гради-
ента для функционала (A (у - u), у - u) в точке ук. В направлении антиградиента происходит наискорейшее убывание значения функционала. Параметр тк+1 в HD находится из условия минимизации нормы погрешности zk+1 = ук+1 - u [11].
Формула для вычисления итерационного параметра тк+1 находится в предположении невырожденности оператора A. Погрешности находятся из уравнений: zk = ук - u , к = 0,1, — Схема (8) с учетом ук = тк + u , примет вид: zk+1 = (е - ткв) zk, к = 0,1,..., z0 = у0 - u . С помощью замены zk = D~1/Zxk осуществляется переход к выражению, содержащему только один оператор:
Х+1 = sk+1Хк, sk = Е - ткС, с = d-1/2 (db~1a) d-1/2. (9)
С помощью равенства ||гк ||в = ||хк || у • ||в — норма в нВ, || • || — норма в н ) описанную задачу расчета параметра тк+1 можно описать следующим образом: рассчитывается значение параметр тк+1 из условия минимизации нормы хк+1 в пространстве н. Рас-
смотрим норму:
к-г = ((£Г-Г,"С) Хк-(Е-Тк"С) Хк)=2||ХкГ ;2Гк"(СХ- х2к ^ (СХк-СХк ) = (10)
= (СХк, СХк ) *Тк+1 - (СХк, Хк ) / (СХк, СХк )+ + ^Ц - (СХк, Хк ) / (СХк, СХк ) . Оператор с не вырожден, поскольку не вырожден оператор А. Таким образом, для любого хк имеем (Схк,Схк) > 0 , и минимальное значение нормы хк+1 достигается при
Тк +1 = (СХк , Хк ) / (СХк , СХк ) • (11)
Подставим (11) в (10) и получим:
||Хк + Л = Рк + 1 \\Хк\\ г (12)
где
Рк+1 = 1 - (СХк, Хк)2 / {(СХк, СХк) (Хк, Хк)} . (13)
Выражение (11) описывает оптимальные значения итерационного параметра тк+1.
1
Выражение (11) с учетом Хк = В 2zk запишется в виде:
Тк + 1 = (ВВZkХk ) / В) , к = 0, 1, ... . (14)
Принимая во внимание, что Azk = Аук - Аи = Аук - f = гк — невязка, а В _1гк = юк — поправка, выражение для расчета параметра тк+1 запишется в виде:
Тк+1 = (В®к, Zk ) / (В®к, ®к ) , к = 0,1, ... , (15)
а двухслойная итерационная схема (8) представима в явном виде:
Ук+1 = Ук - Тк+1®к' к = 0 1 (16)
то алгоритм для данного метода, запишется в виде:
1) для вектора ук вычисляется вектор невязки гк = Аук - f;
2) вычисляются значения вектора поправки из уравнения Вак = гк;
3) рассчитывается параметр тк+1 из выражения (15);
4) рассчитывается новое приближение ук+1 из выражения (16).
Рассмотрим двухслойные градиентные методы, которые использованы при решении модельной задачи эвтрофикации вод мелководного водоема. Каждый метод из данного класса имеет свою область применимости и определяется способом задания
оператора D. Оператор D выбирается таким образом, чтобы выражение (15) для итерационного параметра тк+1 содержало только известные величины.
Если оператор А самосопряжен и положительно определен в н, то для решения (8) можно использовать метод скорейшего спуска (МСС). Если оператор а невырожден и несамосопряжен, а оператор Б*А является положительно определенным, то применим метод минимальных невязок (ММН) [12]. С учетом Azk = Ахк - f = гк и А = А* и выражения (15) получим формулу для расчета итерационного параметра тк+1 согласно методу скорейшего спуска:
Т+1 = (Л, ®к) / (А®^ А®к), к = 0 1 ••• • (17)
Для случая Б = Е получим <пк = Б^1гк = гк и выражение для расчета тк+1 примет
вид:
Тк +1 = (Гк, Гк ) / (АГк, АГк ) , к = 0, 1, ... . (18)
При решении модельной задачи (1) - (3) методом скорейшего спуска по невязке параметр тк+1 рассчитывается по формуле (18), а метод скорейшего спуска по поправке реализуется с помощью расчетной формулы (17).
Метод минимальных поправок (ММП) можно применять для решения модельной задачи (1) - (3) в случае любого несамосопряженного невырожденного оператора а , кроме того требуется положительная определенность оператора Б1 А. Для метода минимальных невязок D = А* А.
Формула для итерационного параметра тк+1 в ММП имеет вид:
тк + 1 = (Аюк, гк) / (Аюк, Аюк ) , к = 0,1, ... . (19)
В случае (б = е ) требуется положительная определенность оператора А, а формула для тк+1 примет вид:
Тк+1 = (АГк , Гк ) / (АГк , АГк ) , к = 0, 1, ... . (20)
В таблице приведены результаты сравнения скоростей сходимости двухслойных градиентных методов вариационного типа, используемых для решения задачи эвтрофи-кации вод мелководного водоема вида (1) - (3).
В таблице используются следующие обозначения: метод минимальных невязок (ММН); метод минимальных поправок (ММП); метод скорейшего спуска по невязке (МССН); метод скорейшего спуска по поправке (МССП); р — номер итерации; п — номер временного слоя; 1р — количество итераций, необходимое для сходимости метода решения СЛАУ для уравнения, описывающего изменение концентрации р, р е ^, х) , где S = S5, X = S16.
Метод минимальных поправок (ММП) был выбран в качестве основного метода в виду его наибольшей скорости сходимости, согласно данным, приведенным в таблице.
Таблица
Сравнение скоростей сходимости методов вариационного типа
ММП
N. n p \ 1 2 3 4
1 Is = 48, Ix = 44 IS = 47, IX = 43 IS = 47, IX = 43 IS = 4 8, IX = 42
2 IS = 53, IX = 4 8 IS = 52, IX = 4 7 IS = 53, IX = 4 6 IS = 53, IX = 45
3 IS = 21, IX = 21 IS = 22, IX = 2 0 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 ^ = 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 ^ = 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 = 2 9 IS = 34, IX = 2 8
4 IS = 2 4, IX = 2 3 IX = 21, IX = 20
МССН
1 IS = 68, IX = 64 Is = 79, Ix = 67 IS = 7 8, IX = 66 IS = 7 8, IX = 65
2 IS = 86, IX = 74 IS = 8 7, IX = 7 3 IS = 86, IX = 73 IS = 8 6, IX = 72
3 IS = 53, IS = 4 9 IS = 52, IX = 45 IS = 51, IX = 42 IS = 49, IX = 41
4 IS = 40, IX = 38 IS = 36, IX = 35
3. Параллельный вариант метода решения сеточных уравнений
Для геометрического разбиения расчетной области с целью равномерной загрузки вычислителей (процессоров) МВС использовался метод k-means, который позволяет найти центры подобластей: Q = Q(3) на основе на минимизации функционала суммарной выборочной дисперсии функции, описывающей расположение элементов (расчетных узлов сетки). Пусть Xi — множество элементов i - ой подобласти, i е {1,..., m}, m —
количество подобластей. Q(3) = ^ -1-^ d2(x, ci) ^ min , где ci = ^ x — центр
i |Xi|xeXi xeXi
подобласти Xi, а d(x, ci) — расстояние между центром подобласти ci и элементом x в Евклидовой метрике. Метод k-means позволяет разбить исходную область на примерно равные подобласти.
Результат работы метода k-means для модельных двумерной и трехмерной областей представлен на рис. 2.
Рис. 2. Декомпозиция области на основе алгоритма k-means: для двумерной (на 9, 38, 150 подобластей); для трехмерной (на 6 и 10)
Опишем работу максиминного алгоритма. Начальное расположение центров подобластей задаются расчетными узлами сетки следующим образом:
1) первый центр подобласти задается первым элементом;
2) второй центр расположен на максимальном удалении от первого центра подобласти;
3) каждый следующий центр подобласти располагается из условия на максимального удаления от ближайшего центра.
Опишем алгоритм работы k-means.
1) Рассчитываются начальное расположение центров подобластей при помощи мак-симинного алгоритма.
2) Все элементы кластеризуются на m клеток Вороного согласно методу ближайшего соседа, т. е. для любого х е Хс, где Хс — подобласть, должно выполняться условие \\х - sll = min \\х - sj , где sc - центр области Хс.
3) Расчет новых центров осуществляется согласно выражению s(ck+1) = т—^А X х •
^ (к)
4) Выполняется проверка условия остановки алгоритма s{ck+1) = s{ck) для всех к = 1, ... г т. Если условие остановки алгоритма не выполняется, то переходим к пункту 2.
Элементы, лежащие на границах подобластей, обмениваются данными при параллельной реализации алгоритмов решения сеточных уравнений. Задача построения выпуклой оболочки решалась с помощью алгоритма Джарвиса. На основе данного алгоритма найдено расположение элементов, участвующих в обменах данными между вычислителями, и номера вычислителей, передающих и принимающих соответствующие данные [13, 14].
Метод сдваивания был использован для расчета итерационного параметра т при решении сеточных уравнений методом минимальных поправок. Следует отметить, что синхронизация разработанного алгоритма решения задачи (1) - (3) производится только при решении сеточных уравнений ММП на переходах между итерационными слоями.
4. Описание программного комплекса
Для решения модельной задачи эвтрофикации мелководного водоема (1) - (3) был создан комплекс прикладных программ, позволяющий производить расчеты концентраций загрязняющих веществ (ЗВ), фито- и зоопланктона в областях сложной формы (Азовское море и Таганрогский залив) [15].
С помощью созданного комплекса программ можно осуществлять:
• внедрение системы и совершенствование комплексного рыбохозяйственного мониторинга в водоемах (кормовых запасов и базы промысловых объектов наблюдение, прогноз и оценка состояния режима экосистем);
• разработку, согласование предложений и мероприятий по обеспечению оптимального режима, сохранения биоразнообразия промысловых ресурсов, экосистем мелководных водоемов;
• совершенствование методологии природоохранных исследований, разработка новых, апробация и внедрение перспективных методов изучения состояния водных экосистем и отдельных компонентов;
• разработку и совершенствование методов диагностики токсического воздействия ЗВ на гидробионты, в том числе ранней и дифференциальной диагностики токсикоза, а также поиск средств антидотной защиты водных экосистем;
• организацию и проведение исследований по выявлению тенденций и закономерностей изменения состояния водных экосистем под воздействием антропогенных факторов, разработка предложений и мероприятий по снижению и предупреждению таких воздействий;
• оценку ущербов рыбному хозяйству, наносимых разными видами хозяйственной деятельности, разработку предложений по предотвращению, уменьшению и адекватной компенсации ущербов.
Значения поля скоростей движения водной среды, рассчитанные в работах [16], относятся к входной информации для моделей гидробиологических процессов, описанных в первой главе. Для математического моделирования гидробиологических и гидродинамических процессов в трехмерной области сложной формы — Азовское море и Таганрогский залив использовались последовательно сгущающиеся прямоугольные сетки размерностями: 251 х 351 х 15, 502 х 702 х 30 , 1004 х 1404 х 60 ,....
Начальное распределение загрязняющих веществ и планктона было учтено в форме, соответствующей пространственно-временным масштабам моделируемых процессов. Разработанный алгоритм численного решения поставленной задачи позволяет свободно варьировать значениями соответствующих параметров, видами управляющих функций и граничные условиями [17]. Понимание механизма функционирования системы и знание ее основных характеристик позволили для преодоления трудностей при настройке программы использовать феноменологический подход. Эффективность такого подхода связана с тем, что адекватность описания поведение экологической системы часто определяется точностью не отдельных параметров, а их соотношений.
Калибровка и верификация разработанной модели эвтрофикации вод мелководного водоема проводились на основе экологических данных по Азовскому морю, полученных в ходе научно-исследовательских экспедиций, проводимых учеными ЮФУ, начиная с 2000 года. В ходе исследований акватории Азовского моря изучались: виды и концентрации основных загрязняющих воды Азовского моря веществ; пространственные рас-
пределения солености и температуры; кислородный режим; видовой состав фито- и зоопланктона; механизмы возникновения и развития заморов в центрально-восточной части водоема. Обработка экспедиционных данных заключалась в оцифровке, пересчете в стандартные единицы, классификации для использования в разных модельных задачах гидробиологии моря.
На рис. 3 приведена схема разработанного программного комплекса.
Рис. 3. Схема программного комплекса
Программный комплекс включает: блок управления, базы океанологических и метеорологических данных, системы интерфейсов, системы ввода — вывода и визуализации. Комплекс программ обладает удобным интерфейсом и обеспечивает ввод данных в диалоговом режиме. Построенный программный комплекс имеет универсальный характер и привязка его к условиям моделируемых районов и объектов осуществляется, как правило, на уровне входной информации. Для практического использования модулей требуется создание специальной информационной базы, содержащей сведения о параметрах, определяющих источники примесей, о климатических и физико-географических условиях исследуемых объектов.
При разработке программного комплекса использовался язык высокого уровня С++. Были разработаны две версии параллельного алгоритма решения модельной задачи эвтрофикации вод мелководного водоема. Первая предназначена для ЭВМ с операционной системой Windows, основана на технологии создания в этой системе дополнительных потоков. Вторая может использоваться для кластерных систем, и основана на технологии передачи сообщений (MPI).
5. Результаты численных экспериментов
При моделировании пространственно-неоднородных процессов эвтрофикации вод Азовского моря (1) - (3) учитывалась внешняя периодичность, приводящая к усложнению системы. При этом колебания плотности планктона могут быть столь велики, что не могут быть объяснены случайными флуктуациями, и визуальная картина такова, что сравнительно небольшие площади высокой плотности («облака», «пятна») разделены пространствами с различными плотностями, зачастую не фиксируемыми общепринятыми методами наблюдений. Особенно ярко это явление выражено в тех местах водоема, для которых характерна потребность в биогенных элементах. При моделировании процессов эвтрофикации учитывался вегетационный период фитопланктона.
Диффузные процессы сглаживают пространственное распределение и рассеивают «пятна». Одна из попыток объяснить парадокс стабильности «пятен» с помощью численных экспериментов заключалась в предположении об активном передвижении гетеротрофных организмов (зоопланктона и рыб) в направлении градиента «пищи», что обеспечивает закрепление пространственной неоднородности биогенных веществ в водной среде. Устойчивая пространственная неоднородность распределения может быть, например, связана с диффузионными процессами и наличием у фитопланктона механизма эктокринного регулирования, т.е. регулирования скорости роста водорослей посредством выделения в среду биологически активных метаболитов.
На рисунке 4. показаны результаты расчета концентрации загрязняющего биогенного вещества для модели (1) - (3) (начальное распределение полей течений при северном ветре). Приведенные ниже рисунки отражают влияние структур течений водного потока в Азовском море на распределение загрязняющего биогенного вещества и фитопланктона. Белым цветом выделены максимальные значения концентраций биогенного вещества (азота) и фитопланктона черным — минимальные.
Рис. 4. Распределение концентрации загрязняющих веществ в различные моменты времени
Результаты моделирования динамики фитопланктона в Азовском море представлены на рис. 5. (N — номер итерации, начальное распределение полей течений водного потока при северном направлении ветра).
Критерием проверки адекватности построенных моделей гидробиологии мелководного водоема служила оценка погрешности моделирования с одновременным учетом натурных данных по имеющимся п замерам, которая вычислялась по формуле:
8 = ! ^ паь , где Бк — значение концентрации загрязняющей
примеси, полученное с помощью натурных экспедиционных измерений; Бк — значение, рассчитанное с помощью модели (1) - (3).
Рис. 5. Распределение концентрации фитопланктона в различные моменты времени
Рассчитанные при различных ветровых ситуациях концентрации загрязняющих веществ и планктона принимались к рассмотрению, если относительная погрешность не превышала 30%.
Заключение
На основе экспедиционных исследований выполнена первичная верификация модели экосистемы Азовского моря. Реализована задача моделирования и прогноза состояния водной экосистемы Азовского моря в условиях антропогенного воздействия и всестороннего изучения уникального водного объекта, который в силу своей мелководности в значительной степени подвергается антропогенному влиянию.
С помощью модели ЭВ (1) - (3) могут быть описаны процессы ассимиляции ын4, нитратредукции (денитрификации), окисления н^, сульфатредукции, нитрифика-ции,окисления и восстановления марганца аммонификации, а также можно изучать механизм условий формирования заморов, прогнозировать изменения кислородного и биогенного режимов в мелководном водоеме. При построении модели были параметризованы процессы биогеохимических циклов химических элементов, ответственных за трансформацию аэробных условий в анаэробные [2].
Созданный программный комплекс объединяет математические модели и базы экологических данных, с помощью разработанного комплекса были изучены условия, при которых мелководные водоемы подвергаются эвтрофитрованию.
Отличительными особенностями разрабатываемых алгоритмов, реализующих поставленные гидробиологические модельные задачи, являются: высокая производительность, достоверность и точность получаемых результатов. Высокая производительность достигается за счет использования эффективных численных методов решения сеточных уравнений, реализованных на параллельных вычислительных системах в реальном и ускоренном масштабах времени. Достоверность достигается за счет учета определяющих физических факторов, таких как: стоки рек, испарение, сила Кориолиса, сложная геометрия береговой линии и дна водоема, турбулентный обмен, трение о дно и ветровые напряжения, динамическое перестроение расчетной области, а также за счет учета
отклонения значения поля давления от гидростатического приближения. Точность достигается применением подробных расчетных сеток, учитывающих степень «заполненности» расчетных ячеек, а также отсутствием неконсервативных диссипативных слагаемых и нефизичных источников (стоков), возникающих в результате конечно-разностных аппроксимаций.
Было проведено сравнение работы созданного программного комплекса, реализующего разработанные сценарии развития экологической обстановки в Азовском море с использованием численной реализации модельной задачи эвтрофикации вод мелководного водоема, с подобными работами в области математического моделирования гидробиологических процессов.
Анализ показал, что в результате разработки программного комплекса удалось повысить точность прогнозов изменения концентраций загрязняющих веществ, фито- и зоопланктона в мелководном водоеме на 10 - 15% в зависимости от решаемой модельной задачи водной экологии.
Работа выполнена при частичной финансовой поддержке Задания № 2014/174 в рамках базовой части государственного задания Минобрнауки России, а также при частичной финансовой поддержке РФФИ по проектам № 15-01-08619, № 15-07-08626, № 15-07-08408.
Литература
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 с.
6. Самарский А.А. Теория разностных схем. М.: Наука, 1989. 616 с.
7. Самарский А.А., Николаев Е.С. Методы решения сеточных уравнений. М.: Наука, 1978. 592 с.
8. Коновалов А.Н. К теории попеременно-треугольного итерационного метода // Сибирский математический журнал. 2002. 43:3. С. 552-572.
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. Никитина А.В., Пучкин М.В., Семенов И.С., Сухинов А.И., Угольницкий Г.А., Усов А.Б., Чистяков А.Е. Дифференциально-игровая модель предотвращения заморов в мелководных водоемах // Управление большими системами. М.: ИПУ РАН. 2015. Вып. 55. C. 343-361.
DOI: 10.14529/cmse160303
MATHEMATICAL MODELING OF EUTROPHICATION PROCESSES IN SHALLOW WATERS ON MULTIPROCESSOR COMPUTER SYSTEM
© 2016 A.I. Sukhinov1, A.V. Nikitina2, A.E. Chistyakov2, I.S. Semenov3, A.A. Semenyakina , D.S. Khachunts
1Don State Technical University (Gagarina Sq. 1, Taganrog, 344000 Russia), Kalyaev Scientific Research Institute of Multiprocessor Computer Systems at Southern Federal University (Chekhova 2, Taganrog, 347928 Russia), 3Scientific Research Center of Supercomputers and Neurocomputers (per. Italyansky 106,
Taganrog, 347900 Russia), E-mail: sukhin@gmail.com, nikitina.vm@gmail.com, cheese_05@mail.ru, flank-er555@yandex.ru, j.a.s.s.y@mail.ru, diana-hachunts@mail.ru Received: 25.02.2016
The paper covered the development of solution methods for model of eutrophication of shallow 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.
FOR CITATION
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. Bulletin of the South Ural State University. Series: Computational Mathematics and Software Engineering. 2016. vol. 5, no. 3. pp. 36-53. (in Russian) DOI: 10.14529/cmse160303.
References
1. Yakushev E.V., Sukhinov A.I., et al. Complex Oceanographic Research of the Sea of Azov in the 28-th Flight of the Research Vessel «Aquanaut». Okeanologiya [Oceanology]. 2003. vol. 43, no. 1. pp. 44-53. (in Russian)
2. Sukhinov A.I., Nikitina A.V., Chistyakov A.E. Simulation Scenario of Biological Rehabilitation of the Azov Sea. Matematicheskoe modelirovanie [Mathematical Modeling]. 2012. vol. 24, no. 9. pp. 3-21. (in Russian)
3. Sukhinov A.I., Chistyakov A.E., Alekseenko E.V. Numerical Implementation Three-Dimensional Model of Hydrodynamics for Shallow Water Basins on Superficialities System. Matematicheskoe modelirovanie [Mathematical Modeling]. 2011. vol. 23, no. 3. pp. 3-21. (in Russian)
4. Sukhinov A.I., Chistyakov A.E. Parallel Implementation of a Three-Dimensional Model of Hydrodynamics of Shallow Water Bodies on Superficialities System.
Vychislitel'nye metody i programmirovanie: Novye vychislitel'nye tekhnologii [Computational Methods and Programming: New Computing Technologies]. 2012. vol. 13. pp. 290-297. (in Russian)
5. Belotserkovsky O.M. Turbulentnost': novye podhody [Turbulence: New Approaches]. M.: Science, 2003. 286 p.
6. Samarskii A.A. Teoriya raznostnyh skhem [Theory of Difference Schemes]. M.: Science, 1989. 616 p.
7. Samarskii A.A., Nikolaev E.S. Metody resheniya setochnyh uravneniy [Methods of Solving Grid Equations]. M.: Science, 1978. 592 p.
8. Konovalov A.N. The Theory of the Alternating-Triangular Iterative Method. Sibirskiy matematicheskiy zhurnal [Siberian Mathematical Journal]. 2002. 43:3. pp. 552-572. (in Russian)
9. Sukhinov A.I., Chistyakov A.E. Adaptive Modified Alternating Triangular Iterative Method for Solving Grid Equations with Non-Selfadjoint Operator. Mathematical Models and Computer Simulations. 2012. vol. 4, no. 4. pp. 398-409. DOI: 10.1134/s2070048212040084.
10. Chistyakov A.E. Theoretical Estimation of Speed Up and Efficiency of Parallel Implementation of the Steepest Descent PTM. Izvestiya YuFU. Tekhnicheskie nauki [Izvestiya SFedU. Engineering Sciences]. 2010. no. 6(107). pp. 237-249. (in Russian)
11. Chistyakov A.E., Hachunts D.S., Nikitina A.V., Protsenko E.A., Kuz-netsova I.Yu. A Library of Parallel Iterative Methods Solvers of Linear Algebraic Equation for the Problem of Convection-Diffusion-Based Decomposition in One Spatial Direction. Sovremennye problemy nauki i obrazovaniya [Modern Problems of Science and Education]. 2015. no. 1. Available at: 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 YuFU. Tekhnicheskie nauki [Izvestiya SFedU. Engineering Sciences]. 2010. no. 6(107). pp. 113-116. (in Russian)
13. Nikitina A.V. Modelling of Biological Kinetics, Stabilizing the Ecological System of the Gulf of Taganrog. Izvestiya YuFU. Tekhnicheskie nauki [Izvestiya SFedU. Engineering Sciences]. 2009. no. 8(97). pp. 130-134. (in Russian)
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 YuFU. Tekhnicheskie nauki [Izvestiya SFedU. Engineering Sciences]. 2011. no. 8(121). pp. 6-13. (in Russian)
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. Vychislit-el'nye metody i programmirovanie: novye vychislitel'nye tekhnologii [Computational Methods and Programming: New Computing Technologies]. 2015. vol. 16. pp. 256-267. (in Russian)
16. Nikitina A.V., Rudneva T.V., Kamyshnikova T.V., Bokareva T.A., Duryagi-na V.V. To the Formation of Kill Zones in the Eastern Part of Azov Sea Coast. Sovremennye problemy nauki i obrazovaniya [Modern Problems of Science and Education]. 2015. no. 1. Available at: 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., Chistyakov A.E. Differential-Game Model of Prevention of Fish Kill in Shallow Ponds. Upravlenie bol'shimi sistemami [Managing Large Systems]. M.: IPU Russian Academy of Sciences. 2015. issue 55. pp. 343-361. (in Russian)
