Реализация на высокопроизводительных вычислительных системах математической модели ветровых течений в Керченском проливе Текст научной статьи по специальности «Математика»

МАТЕМАТИЧЕСКОЕ И ПРОГРАММНОЕ ОБЕСПЕЧЕНИЕ...

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

УДК 519.6:532.5

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

Представленная математическая модель описывается уравнениями движения несжимаемой вязкой жидкости. Задача решается конечно-разностными методами с использованием противопотоковых схем и численно реализована на высокопроизводительных вычислительных системах. Данная модель позволяет рассчитывать течения в случае изменения береговой линии, а ее трехмерность дает возможность вычислять скорости на различных горизонтах. Проведено сравнение эффективности параллельной реализации модели на различных вычислительных системах. Уравнения гидродинамики; ветровое течение; математическая модель; вычислительный эксперимент

Одним из основных этапов решения задач, связанных с моделированием гидрофизических процессов в водоемах, является расчет гидродинамических параметров течения. Определенный интерес представляют водоемы с большой неоднородностью поля глубин. Под большой неоднородностью глубин понимается наличие как глубоководной области, так и обширных мелководных районов (прибрежная зона, внутренние заливы и лиманы), глубина в которых соразмерна с перепадом уровня воды. В судоходной части Керченского пролива глубины составляют 8-10 м, а в южной его части до 19 м. В то же время в прибрежных районах, Таманском и Динском заливах глубина составляет 0,5^ м. Такое распределение глубин позволяет говорить о большой их неоднородности. В настоящее время при расчете параметров течений в водоемах с указанной топологией дна обычно используются специальные координатные системы (о-координаты), возможно использование криволинейных сеток или проведение предварительного преобразования исходной нерегулярной области в регулярную. Однако данные методики заметно усложняют как саму систему дифференциальных уравнений, так и ее численное решение.

В связи с этим для численного исследования течений в Керченском проливе использовалась двухслойная модель гидродинамики водоемов с большой неоднородностью глубин [1, 2]. Достоинство данного подхода заключается в том, что решение подобных задач происходит без предварительного преобразования расчетной области, используя конечно-разностные методы с применением равномерных прямоугольных сеток.

МАТЕМАТИЧЕСКАЯ МОДЕЛЬ

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

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

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

Движение воды в слое I задается уравнениями мелкой воды

Эи Эм Эи ^ эс

+^=-^+

Эу Эх

Г Э2. Э 2 Л (1)

д 2г

дх2 ду2

Т Т

+ _^----------— + р ( х, у),

Н Н хУ '

Контактная информация: (863)219-97-13

Эу Эу ду ЭС

+ и ------- + V ----- + Ом =- £ — +

дґ £ дх £ ду ’ *ду

Эх2

дС

+

д2

Л

ду2

+ ^1 -%! + Р (х, у),

Н Н л !

дґ

д(Нил.) + д(Ну£)

эу

дх

= 0.

(2)

(3)

В уравнениях (1)-(3)

= -1} НІ

мел-

1 с

= — [ vdz, Н = к + С; к = к(х, у) - глубина Н-н

ководного слоя; £ - перепад уровня воды; и/ = = иДх, у, ^, у^ = уДх, у, 0 - скорости в слое I; функции Fx(x, у) и Fy(x, у) описывают взаимодействие верхнего и нижнего слоев между собой; т^, т/у - проекции на оси ОХ и ОУ силы трения ветра о поверхность водоема; тЬх, тЬу -проекции на оси ОХ и ОУ силы трения жидкости о дно (или о глубоководный слой воды). Эти величины зависят от скорости ветра Wв = (ИХ; Жу] и течения WJ^ = (и/; иу/} и определяются так [1]:

т / =у|Wв К, Ть =р|Wг ^ ,

и£ + У£2

где = д/ Их2 + И2, |WГ| = ^

Коэффициент Р(х, у) трения верхнего слоя жидкости о дно (или о глубоководный слой) выбирался из диапазона 2• 10-4 ^5 • 10-4, а коэффициент у трения ветра о слой I из диапазона 10-5 +-10-6. Это делалось для лучшего совпадения расчетных значений с натурными данными.

Движение воды в слое II задается уравнениями движения несжимаемой вязкой жидкости

ди ди ди ди _ 1 др

------+ и--------+ V--------+ №---------------Ну =-----------— +

дґ дх ду дг р дх

( з2

2,Л

д2 и д и

дх2 ду2

+ -

д

ди

(4)

Эу Эу Эу Эу „ 1 др

------+ и-------+ у-------+ №---------+ І2и =----------------+

дґ

+у.

дх ^ д 2у

Эу Эг

р эу

____ ЭИ

дх2 Эу2

д ( Эу +—I п г — дг I г дг

ди + Эу + д^ о

дх ду дг

(5)

(6)

К системе уравнений (4)-(6) добавляется уравнение гидростатического давления

Р = &Р(С- г)+ Ра . (7)

В (4)-(7) и, у, ^ - компоненты скорости; х, у, z, t - пространственные переменные и время соответственно; уху = 15000^30000, уz = 0,2^0,4 -коэффициенты горизонтальной и вертикальной вязкости соответственно; р - плотность воды; g = 9,8 м/с2 - ускорение силы тяжести; ра - атмосферное давление.

Приведенные системы замыкаются соответствующими граничными условиями. Вдоль твердой границы ставятся условия скольжения или задаются скорости втекания или вытекания воды. На границе между слоями дО.1 ставится условие равенства скоростей

и\ = и , У = у .

|Эй; ^ 1Эй; *

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

Р, (х, у ) =

Н

ВЫЧИСЛИТЕЛЬНЫЕ

ЭКСПЕРИМЕНТЫ

Задача решается конечно-разностными методами. Алгоритм вычисления параметров течения воды на («+1)-м временном слое основан на том принципе, что каждое уравнение является «определяющим» для своего неизвестного. Все остальные переменные считаются известными и берутся с «-го слоя. При конечноразностной аппроксимации уравнений количества движения используются неявные «проти-вопотоковые» схемы. Перепад уровня воды и вертикальная компонента скорости определяются из разностных аналогов уравнений (3) и (6). Качество такого алгоритма оценивается хорошим совпадением расчетных значений с натурными данными, полученными на станциях I, II, III и IV (рис. 2). Более подробно эти сравнения приведены в [1].

Особенностью моделирования гидродинамики Керченского пролива является наличие косы Тузла. Тузлинский нанос - это песчаная полоса, полупогруженная в илы. Коса растет в результате типичной аккумуляции береговых наносов. В 1925 году во время сильного южного шторма в море прорвало косу вблизи Таманского берега, и она превратилась в остров, отделенный от мыса Тузла проливом, который достиг через некоторое время ширины нескольких километров. Средняя глубина пролива между мысом Тузла (Тамань) и надземной частью косы не превышает 0,3-0,6 м. На расстоянии 2,5 км от мыса имеется искусственная промоина глубиной порядка 1,3 м и шириной не более 30 м.

у

£

£

+

и

=

£

В 2003 г. со стороны Таманского полуострова на месте прежней косы была частично насыпана дамба.

Характерной ситуацией в Керченском проливе является смена течений с Азовоморского на Черноморское или наоборот. С помощью построенной математической модели была численно исследована возможная динамика изменения течения с Азовоморского на Черноморское. На первом этапе было рассчитано Азовоморское течение под действием северовосточного ветра силой 3 м/с. Затем под действием юго-юго-западного ветра силой 3,5 м/с проводился расчет течений на протяжении

7 часов модельного времени.

В течение первого часа направление течения изменилось, в основном, в примыкающих к проливу заливах. Затем Черноморское течение начинает формироваться в южной части пролива, но течения в протоках вокруг о. Тузла носят северный характер (рис. 2).

Через три часа после смены направления ветра Черноморское течение занимает почти всю южную часть пролива (рис. 3). Исключение составляет течение через створ между о. Тузла и Таманским полуостровом. Очень слабые течения наблюдаются в северной и Павловской узостях. К концу седьмого часа Черноморское течение полностью формируется.

С помощью построенной математической модели была численно исследована возможная картина течений в центральной части пролива при отсутствии или наличии дамбы. Создание дамбы вдоль косы Тузла в 2003 г. привело к изменению течений в центральной части Керченского пролива. К сожалению, за период с 2003 г. по настоящее время натурных данных по скоростям течений не имеется.

На рис. 4 и рис. 5 показаны поля скоростей при отсутствии или наличии дамбы в случае действия юго-западного ветра со скоростью

8 м/с в течение 8 часов. Так, например, при отсутствии дамбы примерно половина объема воды, проходящей через станцию III, нагоняется в Таманский залив, а другая половина движется вдоль северной стороны острова Тузла в сторону Азовского моря (рис. 4). При наличии дамбы расход воды через станцию III значительно сокращается, и вода в Таманский залив поступает, в основном, с северной стороны о. Тузла, двигаясь уже в юго-восточном направлении (рис. 5).

АНАЛИЗ ЧИСЛЕННОЙ РЕАЛИЗАЦИИ

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

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

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

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

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

• на рабочей станции QUAD с 4-ядерным процессором Q6600 с тактовой частотой 2.4 Ггц и оперативной памятью 4 Гбайт.

Рассматривалось три размера разностных сеток с шагами по горизонтали 60 м, 40 м, 20 м с максимальным числом неизвестных n = = 800000, n = 1700000 и n = 10300000 соответственно. При тестировании расчет проводился для 200 временных шагов.

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

Для самой грубой сетки (h = 60 м) линейное ускорение на этой платформе сохраняется для числа узлов np < 4. В то время как на кластере с однократной скоростью Infiniband ускорение начинает отклоняться от линейного уже для np < 2. И совершенно неоправданным представляется запуск нескольких MPI процессов на многоядерных вычислительных узлах (рабочая станция QUAD).

Рис. 2. Картина течения в Керченском проливе через 2 часа после смены ветра с северного на южное направление

Рис. 3. Картина течения в Керченском проливе через 3 часа после смены ветра с северного на южное направление

Рис. 5. Поле скоростей в Керченском проливе при наличии дамбы вдоль косы Тузла

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

Np Сетка h = 60 м; n = 800000

INFINI (с.) sp IBMX (сек.) sp QUAD (сек.) sp

1 1858 1,0 1677 1,0 1203 1,0

2 971 1,9 853 2,0 1008 1,2

3 714 2,6 570 2,9 1421 0,8

4 545 3,4 435 3,9 1648 0,7

5 470 4,0 357 4,7 -

6 405 4,6 313 5,4 -

7 362 5,1 273 6,1 -

Сетка h = 40 м; n = 1700000

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 881 4,3 -

6 903 4,3 758 5,1 -

7 802 4,9 681 5,6 -

Сетка h = 20 м; n = 10300000

1 mem - 15182 1,0 mem -

2 mem - 8097 1,9 mem -

3 mem - 5818 2,6 -

4 mem - 4503 3,4 -

5 mem - 3771 4,0 -

6 mem - 3304 4,6 -

7 mem - 3098 4,9 -

На сетке (h = 40 м) упомянутые закономерности ускорения модели на разных вычислительных системах сохраняются, однако коэффициент ускорения снижается.

На мелкой сетке (h = 20 м), когда число неизвестных становится очень большим, подобные задачи решаются только на кластере IBMX. При попытке решения этих задач на системах INFINI или QUAD возникают проблемы с оперативной памятью.

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

запускался только один счетный процесс. Использование второго ядра, так же, как и на системе QUAD приводило к деградации производительности узла.

ЗАКЛЮЧЕНИЕ

Сравнение результатов расчета с наблюдаемыми данными [1] показало, что представленная математическая модель достаточно адекватно описывает гидродинамику течений в Керченском проливе. Данная модель позволяет рассчитывать течения в случае изменения береговой линии, а ее трехмерность дает возможность вычислять скорости на различных горизонтах.

СПИСОК ЛИТЕРАТУРЫ

1. Моделирование гидрофизических процессов в водоемах с обширными районами мелководья / Л. А. Крукиер [и др.]. Ростов-на-Дону, 2009, Изд-во ЮФУ, 244 с.

2. Чикин А. Л., Дацюк В. Н., Дацюк О. В. Двухслойная математическая модель гидродинамики в водоемах с большой неоднородностью глубин // Научный сервис в сети Интернет: масштабируемость, параллельность, эффективность: труды Все-рос. суперкомпьютерной конф. (21-26 сентября 2009 г., г. Новороссийск). М.: Изд-во МГУ, 2009. С. 179182.

3. Расчет на МВС гидродинамики водоемов с большой неоднородностью глубин. / В. Н. Дацюк [и др.] // Суперкомпьютерные технологии: разработка, программирование, применение (СКТ-2010): Мат. Междунар. науч.-техн. конф. Т. 1. Таганрог, 2010.

ОБ АВТОРАХ

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

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

Чикин Алексей Львович, гл. науч. сотр. Южного научного центра РАН. Доктор физикоматематических наук.