Моделирование экстремального наводнения в дельте Дона на многопроцессорных вычислительных системах

МОДЕЛИРОВАНИЕ ЭКСТРЕМАЛЬНОГО НАВОДНЕНИЯ В ДЕЛЬТЕ ДОНА НА МНОГОПРОЦЕССОРНЫХ ВЫЧИСЛИТЕЛЬНЫХ СИСТЕМАХ1

В.Н. Дацюк, Л.А. Крукиер, А.Л. Чикин, Л.Г. Чикина

По материалам ежедневных гидрометеорологических наблюдений на береговой базе Южного научного центра РАН в период с 20 марта по 26 марта 2013 года, проведено восстановление картины аномального затопления дельты Дона. Расчеты уровня воды и объема поступающей в дельту Дона воды проводились с помощью двухслойной математической модели. Расчеты проводились на многопроцессорных вычислительных системах, установленных в Южном федеральном университете.

Ключевые слова: математическая модель, нагон воды, многопроцессорная вычислительная система.

Введение

23-24 марта 2013 г. в результате сильного штормового нагона воды из акватории Азовского моря значительная территория в дельте Дона чрезвычайно быстро была затоплена водой, в 20 населенных пунктах более 2,4 тыс. домовладений и свыше 5,3 тыс. человек оказались пострадавшими от наводнения. Материальные потери для населения и экономики региона оцениваются суммой более 500 млн. руб. Ветровая ситуация в период с 20 марта по 26 марта 2013 года складывалась таким образом, что ветер существенно менял направление и скорость: первую половину 23 марта дул ветер восточных направлений 3-11 м/с, к 19 часам направление ветров разворачивается на юг, усиливаясь до 5-12 м/с, после полуночи 24 марта ветра приняли западное направление со скоростью до 15 м/с. К 14-17 часам 24 марта порывы ветра достигали 20-22 м/с, что привело к очень быстрому затоплению прибрежных территорий, строений и объектов в населенных пунктах. Уровень воды поднимался со скоростью 11-12 см/час. К 18 часам 24 марта наблюдался максимальный уровень воды (280 см) по данным самописца, установленном на посту в п. Донской.

Ситуацию с затоплением усугубил повышенный (до 400 м3/с) попуск воды через гидроузел Цимлянского водохранилища в период с 21 по 25 марта, при этом дошедшая до низовьев Дона за двое суток волна половодья могла быть причиной дополнительного подпора вод, поступающих из Таганрогского залива.

По метеорологическим наблюдениям, полученным с береговой научноэкспедиционной базы Южного научного центра РАН на период с 20 марта по 26 марта 2013 года было промоделировано развитие нагонного явления в дельте Дона.

1 Статья рекомендована к публикации программным комитетом Международной научной конференции "Параллельные вычислительные технологии 2014".

1. Постановка задачи

При моделировании использовалась двухслойная математическая модель гидродинамики водоема, достаточно адекватно описывающая ветровые течения в Азовском море в целом, и в Таганрогском заливе в частности [3]. Данные водоемы имеют как относительно глубоководные районы, так и обширные районы мелководья.

Суть нашего подхода состоит в разбиении исходной трехмерной области расчета ^ (водной толще водоема) ограниченной сверху акваториальной, а снизу донной поверхностями на две — глубоководную и мелководную части. Для декомпозиции пространственной области моделирования ^ проведем горизонтальную секущую плоскость Р, отстоящую от невозмущенной поверхности водоема Р0 на глубине, равной максимальной глубине мелководья (рис. 1). Таким образом плоскость Р разделила исходную область на две подобласти: верхний слой Ш (слой I) -все мелководье и верхняя часть глубоководного слоя, и глубоководный слой Ш (слой II). Предполагается, что эффект осушения из-за сгона воды может присутствовать только в мелководных районах.

Границы расчетной области О могут быть твердыми дО^ (донная поверхность, переходящая в береговую линию), участками втекания или вытекания воды дО^ , свободной поверхностью дО £ .

Считается, что на движение воды в слое I влияет ветер, а движение в слое II инициируется как градиентами давления, так и движением слоя I.

Считаем, что слой I достаточно мелкий (значения возможных возмущений уровня воды и глубины слоя близки), а и и V не зависят от z. Движение воды в слое I описывается уравнениями мелкой воды:

Рис. 1. Вертикальный разрез водоема с обширным мелководьем

’*

дх ду

Здесь н = И + д ; И = И(х, у) - глубина мелководного слоя; и5 = и5 (х, у, /), V5 = V,

скорости в слое I; функции

описывают взаимодействие верхнего и

нижнего слоев между собой; С, = С, (х, у, t) - возмущение уровня воды; О -коэффициент

Кориолиса; т

sx^ sy

проекции на оси ОХ и OY силы трения ветра о поверхность водо-

ема; тЬх ,тЬу - проекции на оси ОХ и OY силы трения жидкости о дно (или нижний слой воды). Эти величины зависят от скорости ветра ЖБ ^¥х ,Жу } и течения Жт {и, У5 } и определяются так:

Т = г\Щ%, т = в\&т\жт,

где в ( х, у) - коэффициент трения верхнего слоя жидкости о дно (или о глубоководный слой); у - коэффициент трения ветра о слой I.

Движение воды в глубоководном слое II описывается системой, состоящей из уравнений количества движения, уравнения неразрывности среды и уравнения гидростатического давления:

ди ди ди ди ^ 1 др

------+ и-------+ V-------+ w---------V =------------------+

дt дх

ду

дz

( ^<2 ^2

д и д и

+

р дх ^ дх ду

ди дz

+

д + — дz

Ву Ву Ву Ву

1 др

В1 дх ду

дz

Р дУ

ху

д у д у

---2 +----2

удх ду ,

г Вул

V 2 в

+

д +— Вz

/

ди ду дм _ — + — + — = 0.

дх ду дz

Р = 8Р (С - ^) + Ра.

Здесь и = и(х, у, 7, ^), V = у(х, у, 7, ^), и> = н^х, у, 7, ^) - компоненты вектора скорости; р(х, у, 7,1) - давление; х, у, z, t - пространственные переменные и время соответственно; Ра = Ра(хУ) - атмосферное давление; Vху, Vг (z) - коэффициенты горизонтальной и вертикальной вязкости соответственно; р - плотность воды; % = 9,8 м/с2 - ускорение силы тяжести.

Граничные условия на твердой границе 30.т задаются условиями скольжения: дУ

дп

= 0

дПт

где Уп - нормальная составляющая вектора скорости, \т - касательная составляющая вектора скорости. В местах втекания или вытекания воды 30к задаются соответству-

ющие значения скоростей и|дп = и1, у|

= V , и \ = и ,, V \ = V,.

1ап я !’ * ш я ^ ’ * ш я

На границе между слоями сЮ; ставится условие равенства скоростей

и

= и , V = V

ш, * ’ I ш, *

Функции ^ (х, у) и Р, (х, У) , описывающие взаимодействие I и II слоя, задаются следующим образом:

F (x, У) = f, F.. (x, y) = f.

В качестве начальных данных можно задавать известное распределение скоростей и уровня воды

I о I 0 1 о I о I 0 ¡Л /-0

u = U , и = U , V = V , v I = v , W = w , Z = Z

11=0 ’ slt=0 s’ 11=0 ’ slt=0 s’ 11=0 ’ ^ It=0 ^

или считать эти значения нулевыми.

Уравнения модели решаются конечно-разностными методами. Алгоритм вычисления параметров течения воды на (п+1)-м временном слое основан на том принципе, что каждое уравнение является «определяющим» для своего неизвестного. Все остальные переменные считаются известными и берутся с n-го слоя. При конечно-разностной аппроксимации уравнений количества движения используются неявные "противопотоко-вые" схемы. Перепад уровня воды и вертикальная компонента скорости определяются из разностных аналогов уравнений дифференциальных уравнений.

Для решения систем линейных алгебраических уравнений, возникающих при дискретизации исходных дифференциальных уравнений, использовалась библиотека параллельных подпрограмм Aztec. В этой библиотеке реализован набор итерационных методов Крылова для решения систем линейных алгебраических уравнений с разреженными матрицами. Распараллеливание выполнено в парадигме «распределенной памяти» с использованием коммуникационной библиотеки MPI. Такой подход не только ускоряет вычисления, но значительно снижает требования к объему оперативной памяти на вычислительных узлах. При увеличении размера сетки при дискретизации решаемой системы уравнений достаточно просто увеличить количество вычислительных узлов для решения задачи. Особенностью итерационных методов является многократное использование операции умножения матрицы на вектор. Профилирование программы показало, что на эту процедуру приходится более 50% всех вычислительных затрат [1]. Поэтому оптимизации параллельной версии этой подпрограммы для минимизации объемов пересылаемых данных было уделено большое внимание. Тем не менее, из-за большого удельного веса коммуникационных операций при работе с разреженными матрицами трудно рассчитывать на высокую степень масштабируемости этой процедуры и соответственно всей программы в целом. Кроме того, прямое тестирование подтвердило, что использование этих методов предъявляет повышенные требования к коммуникационной среде.

На задаче расчета течений было протестировано время счета на трех вычислительных системах [2]:

1) на кластере IBMX с 2-х ядерными процессорами Intel Xeon 5160 и коммуникационной сетью DDR Infiniband;

2) на кластере INFINI с процессорами Intel Pentium 4 3.4 Ггц и коммуникационной сетью SDR Infiniband;

3) на рабочей станции QUAD с 4-х ядерным процессором Q6600 с тактовой частотой 2.4 Ггц и оперативной памятью 4Гбайт. Использовалась разностная сетка с максимальным числом неизвестных n=1700000. При тестировании расчет проводился для 200 временных шагов.

Как и следовало ожидать, наибольшее ускорение (Sp) прослеживается на кластере IBMX, с коммуникационной сетью DDR Infiniband (см. таблицу). Линейное ускорение на этой платформе сохраняется для числа узлов. В то время как, на кластере с одно-

кратной скоростью Infiniband ускорение начинает отклоняться от линейного уже для n < 4

р . И совершенно неоправданным представляется запуск нескольких MPI процессов на многоядерных вычислительных узлах (рабочая станция QUAD).

Сравнение производительности различных вычислительных платформ на задаче расчета ветровых течений в Керченском проливе.

Np INFIN I (с) Sp IBM X (с) Sp QUA D (с) Sp

1 3896 1,0 3828 1,0 2740 1,0

2 2132 1,8 2015 1,9 2154 1,3

3 1542 2,5 1376 2,8 2865 0,95

4 1212 3,2 1066 3,6 2920 0,93

5 1029 3,8 8 OO 1—1 4,3 -

6 903 4,3 758 5,1 -

7 802 4,9 6 OO 1—1 5,6 -

Численное исследование показало, что при использовании технологии MPI на вычислительной системе QUAD с 4-х ядерным процессором производительность немного возрастает при подключении 2-го ядра, а затем падает при подключении последующих. В то время, как на системах с распределенной памятью INFINI и IBMX наблюдается хороший рост ускорения при увеличении числа вычислительных узлов. Заметим, что на вычислительном кластере IBMX с 2-ядерным процессором на каждом узле запускался только один счетный процесс. Использование второго ядра, на малом числе узлов приводило к небольшому ускорению работы программы, а при увеличении числа узлов так же, как и на системе QUAD к замедлению работы программы.

2. Анализ результатов расчета

Качество алгоритма оценивается совпадением расчетных значений с натурными данными, полученными на метеопосту в п. Донской (рис. 2).

Расчет проводился с шагом по времени dt = 5 минут, и данные снимались через каждые 3 часа. Видно, что модель достаточно хорошо описывает развитие гидрологической ситуации в данном районе.

Для расчета объема поступающей в дельту Дона воды (куб.м/с) из Таганрогского залива было рассмотрено поперечное сечение восточной части залива длиной около 30 км. Считалось, что на затопление расходуется только нагонная составляющая перепада уровня, которая с учетом скорости течения давала объемный расход поступающей воды.

На рис. 3 приведен график изменения объема поступающей из Таганрогского залива воды в дельту Дона.

Видно, что максимальный объем воды поступал в ночь с 23 на 24 марта и составил около 10 тыс. куб. м в секунду. Затем поступление стало резко уменьшаться, но все равно скорость поступления воды оставалось положительной, что вело к продолжению процесса затопления. И только к вечеру 24 марта поступление воды прекратилось, и во-

да стала убывать. Вечером 25 марта расход убывающей в Таганрогский залив воды составил примерно 5 тыс. куб.м/с.

Рис. 2. Изменение относительного уровня воды на период с 20 по 26 марта 2013 года на

посту «Кагальник»

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

а) начало нагона при юго-западном ветре после штилевой погоды и б) начало нагона после сгонного явления, то есть после действия восточного ветра. Сгон моделировался при скорости восточного ветра 9 м/с, который затем изменялся на ветер юго-западного направления скоростью 15 м/с.

В случае предварительного сгона значительно увеличивается скорость нагонного течения, а соответственно, и объем поступающей в дельту воды. Хотя в обоих случаях затопление дельты длился 22-23 часа, что соответствует динамике повышения уровня 2324 марта, изменения уровня воды существенно различаются. Если нагон начинается при штилевой погоде, то уровень воды поднимается от нулевой отметки, соответствующей естественной глубине водоема (рис. 4), и дальше монотонно возрастает по мере поступления излишков воды. Если перед нагоном наблюдался сгон воды, то сначала уровень отрицателен (что характерно для сгонного явления), но затем быстро увеличивается и, в конце концов, становится выше нагона после штиля. Этот результат показывает, что при резкой смене восточного ветра на западный затопление происходит быстрее и оно сильнее по масштабу, чем при постоянно действующим западном ветре.

Наибольшая

Рис. 3. Изменение объемного расхода воды на период с 20 по 26 марта 2013 года

Рис. 4. Поведение скорости течения и уровня воды при различных начальных условиях

нагона

Заключение

Сравнение результатов расчета с наблюденными данными показало, что представленная математическая модель достаточно адекватно описывает гидродинамику течений в северо-восточной части Таганрогского залива, а также сгонно-нагонные явления. Данная модель позволяет рассчитывать ветровые течения и сгоны-нагоны в случае изменения береговой линии, что позволит в дальнейшем прогнозировать степень опасности затопления дельты Дона.

Литература

1. Дацюк В.Н., Крукиер Л.А., Чикин А.Л.. Реализация на высокопроизводительных вычислительных системах математической модели ветровых течений в Керченском проливе //Вестник УГАТУ, Т.15, № 5 (45), 2011. — С.155-160.

2. Букатов А.А., Дацюк В.Н., Дацюк О.В., Крукиер Л.А. О реальной производительности вычислительных кластеров с многоядерными узлами // Известия вузов. Северо-Кавказский регион. Естественные науки. Спецвыпуск «Теоретические и прикладные проблемы математики, экономики, образования». 2011. — С. 9—15.

3. Chikin A.L. A technique for evaluating flow parameters in water bodies with a highly heterogeneous depth. Water Resources, 2005, Vol. 32, No. 1, January. — P. 50-55.

Дацюк Виктор Николаевич, к.ф.-м.н., начальник отдела, Южно-Российский региональный центр информатизаци, г. Ростов-на-Дону, Россия, dvn@sfedu.ru

Крукиер Лев Абрамович, д.ф.-м.н., профессор, директор Южно-Российский региональный центр информатизаци, г. Ростов-на-Дону, Россия, krukier@sfedu.ru

Чикин Алексей Львович, д.ф.-м.н., г.н.с., Институт аридных зон Южного научного центра РАН, г. Ростов-на-Дону, Россия, chikin@sfedu.ru

Чикина Любовь Григорьевна, д.ф.-м.н., профессор, Южный федеральный университет, г. Ростов-на-Дону, Россия, lchikina@sfedu.ru

Поступила в редакцию 17 февраля 2014 г.

Bulletin of the South Ural State University Series “Computational Mathematics and Software Engineering”

2014, vol. 3, no. 1, pp. 80-88

MODELING OF EXTREME FLOODS IN THE DELTA OF DON RIVER ON THE MULTIPROCESSOR COMPUTER SYSTEMS

V.N. Datsuk, Southern Federal University (Rostov-on-Don, Russian Federation),

L.A. Krukier, Southern Federal University (Rostov-on-Don, Russian Federation),

A.L. Chikin, Institute of Arid Zones Southern Scientific Centre RAS (Rostov-on-Don, Russian Federation),

L.G. Chikina, Southern Federal University (Rostov-on-Don, Russian Federation),

According to the materials of daily meteorological observations in the period from March 20 to March 26, 2013 the restoration of abnormal flooding Don delta were carried out. For calculations of water level and volume of incoming water into the delta of the Don the two-layer mathematical model were carried out. The calculations were carried on multiprocessor computer systems installed in the Southern Federal University.

Keywords: Mathematical model, surge, multiprocessor computer system.

References

1. Datsuk V.N., Krukier L.A., Chikin A.L. realizatsiya na VVS matematicheskoy modeli vetrovih techeniy v Kerchenskom prolive [Implementation on High Performance Computing Systems mathematical model of the wind currents in the Kerch Strait]. Vestnik UGATU [Vestnik Ufa State Aviation Technical University]. 2011. Vol. 15, No. 5. P. 155160

2. Bukatov A.A., Datsuk V.N., Datsuk O.V. O real’noy proizvoditelnosty vichislitelnih klasterov s mnogoyadernimi uzlamy [About the actual performance computing clusters of multicore nodes]. Izvestiya vuzov. Severo-Kavkazskiy region. Estestvennie nauky. Spetsvipusk “Teoreticheskie I prikladnie problemi matematiki, ekonomiki, obrazovaniya” [Izvestiya Vuzov. North Caucasus region. Natural sciences. Special issue of "Theoretical and applied problems of mathematics, economics, education"]. 2011. P. 9-15.

3. Chikin A.L. A technique for evaluating flow parameters in water bodies with a highly heterogeneous depth. Water Resources, 2005, Vol. 32, No. 1, January. - P. 50-55.

Received 17 February 2014