CC BY
72
УДК 519.63
ПАРАЛЛЕЛЬНЫЕ РЕАЛИЗАЦИИ МЕТОДА КОНЕЧНЫХ ЭЛЕМЕНТОВ ДЛЯ КРАЕВОЙ ЗАДАЧИ ДЛЯ УРАВНЕНИЙ МЕЛКОЙ ВОДЫ
Е.Д. Карепова, В.В. Шайдуров, М. С. Вдовенко
PARALLEL IMPLEMENTATION OF THE FEM FOR BOUNDARY-VALUE PROBLEM FOR THE SHALLOW WATER EQUATIONS
E.D. Karepova, V.V. Shaidurov, M.S. Vdovenko
Проведено исследование эффективности двух параллельных реализаций алгоритма численного решения краевой задачи для уравнений мелкой воды, выполненных с помощью библиотеки MPI для языка Си. Представлены результаты численных экспериментов на модельной сетке и неструктурированной сетке для акватории Охотского моря. Приведены сравнительные результаты ускорения вычислений в зависимости от количества процессов, способа реализации коммуникаций, способа декомпозиции вычислительной области.
Ключевые слова: модели мелкой воды, декомпозиция согла-
сованной триангуляции, метод конечных элементов, блокирующие/неблокирующие передачи, ускорение вычислений
The autors carry out the efficiency research of two parallel realizations of the algorithm of the numerical solution of a boundary-value problem for the shallow water equations solved with the help of the MPI library for C. They set out the results of the numerical experiments on the model grid and on the unstructured grid for the Sea of Okhotsk water area. The authors provide the comparative outputs of the computations speedup in relation to the number of processes, the mode of communication realization, the mode of decomposing of the calculating area.
Keywords: shallow water model, decomposition of the coordinated triangulation, FEM, blocking/nonblocking transmissions, computational speedup
Введение
Модели мелкой воды хорошо описывают большой круг природных явлений, таких как крупномасштабные поверхностные волны, возникающие в морях и океанах, цунами, приливные течения, поверхностный и русловой сток, гравитационные колебания поверхности океанов [1, 2]. В работах [2, 3, 4] рассмотрено численное моделирование поверхностных волн в больших акваториях с учетом сферичности Земли и ускорения Кориолиса на основе уравнений мелкой воды. В работе [2] для дифференциальной постановки задачи выведены полезные априорные оценки, обеспечивающие устойчивость решения и однозначную разрешимость задачи. В [3, 4] для этой же задачи построен метод конечных элементов, для которого получены необходимые априорные оценки. Там же приведены результаты численных экспериментов на модельных сетках, для акваторий Охотского моря и Мирового океана.
В настоящей работе проведено исследование эффективности двух параллельных реализаций алгоритма численного решения краевой задачи для уравнений мелкой воды, выполненных с помощью библиотеки МР1 для языка Си. Первая из них основана на декомпозиции области с перекрытием, ширина которого диктуется шаблоном дискретного аналога. Вторая параллельная реализация явно использует проведение вычислений в методе конечных элементов по треугольникам с декомпозицией без использования теневых граней.
Численные эксперименты проводились на модельной сетке и сетке для акватории Охотского моря. Приведены сравнительные результаты ускорения вычислений в зависимости от количества процессов, способа реализации коммуникаций (блокирующие, неблокирующие передачи), способа декомпозиции вычислительной области.
Статья рекомендована к печати программным комитетом международной научной конференции «Параллельные вычислительные технологии 2009» (http://agora.guru.ru/pavt).
1. Моделирование поверхностных волн в водоемах методом конечных элементов
Для стандартной сферической системы координат (г, А, в) с началом в центре земного шара будем использовать вместо угла в географическую широту (р = тг — 0, так что
О < ср < 7г. Через Л обозначена географическая долгота 0 < Л < 2тг. Полагаем всюду г = Д#, где Яе — радиус Земли, который считается постоянным.
Рассмотрим задачу в следующей постановке. Пусть О, — заданная область на сфере с границей Г=Гх и Г2, где - часть границы, проходящая вдоль берега, а Г2 = Г \ Гх
— часть границы, проходящая по морю. Обозначим через т\ и — характеристические функции соответствующих участков границы. Для простоты считается, что точки <р = О и (р = 7Г (полюса) не входят в Г2. Относительно неизвестных функций и = г^(£, Л, <р), V = ?;(£, А, (р) и £ = £(£, А, у?) запишем в х (О, Т) уравнения баланса импульсов и уравнение неразрывности [2]:
где га, v — компоненты вектора скорости U по осям А и ip соответственно; £ — отклонение свободной поверхности от невозмущенного уровня; Н(\,ср) > 0 — глубина водоема в точке (А,<£>); функция Rf = r*|U|/iif учитывает силу трения о дно, г* — коэффициент трения; I = —2ш cos <р — параметр Кориолиса; т = 1/(Re siny?); п = 1/Re; д — ускорение силы тяжести; /1 = /i(£,A,<£), /2 = /2(t,A,<^) и /3 = fs(t,X,(p) — заданные функции внешних воздействий.
Граничные условия рассмотрим в следующем виде [2]:
где ип — и • п, п = (щ, —П2) — вектор внешней нормали к границе; /З Є [0,1] — заданный
т
параметр, (і = А, ср) — заданная граничная функция, определенная на границе Г2. Зададим также начальные условия
-=Ш + т9--Я,и + /1
= lv + mg — — RfU + /1,
(1)
Нип + Рт2^ДН$ = m2 у/дНd на Г х (0,Т),
(2)
п
и(0,\,<р) = и0(\,<р), v(0,\,<p) = v0(\,<p), £(0,А,<р) = £о(А,¥>).
(3)
Для дискретизации по времени разобьем временной отрезок [О, Г] на К интервалов:
О = <о < *1 < ■ • • < Ък = Т с шагом т = Т/К. Аппроксимируем производные по времени левыми разностями и рассмотрим систему (1)-(2) на временном интервале (<*,£*+1):
и — 1ь — тд ^ = /1 + ^ик в О,
(-^ + V + 1и — пд ^ = /2 + ^-ук в П, (4)
(&<я->+&(£*•))=Л+^ в "•
Нип + /Зт2УДНС = гп2л/дН(1 на Г, к = 0,1,..., К — 1, (5)
где /(£ь А, ф) = /*, А, ф) = /к+1 = /. Индекс к + 1 в разностных выражениях далее
будем опускать там, где это не вызывает двойного толкования. Донное трение задается в виде Л} = г*\Х]к\/Н. Заметим, что граничные условия (5) являются естественными для задачи (4).
Для построения метода Бубнова - Галеркина рассмотрим некоторую согласованную три-ангуляцию Т = {с*}|Й области Л, состоящую из невырожденных треугольников с прямолинейными сторонами в координатах А и и содержащую область Л. Согласованность означает, что каждое ребро треугольника либо является граничным и принадлежит лишь одному треугольнику, либо является общим для двух соседних треугольников, внутренность которых попарно не пересекается. Сетка в общем случае может быть неструктурированной (рис. 1). Пусть — множество узлов (т.е. вершин треугольных элементов) общим числом Нпй, ^/г — множество Внутренних уЗЛОВ.
Для каждого узла х3 Е введем базисную функцию Ф^(А, <р), которая равна единице в равна нулю во всех остальных узлах и линейна в каждом треугольнике. Обозначим линейную оболочку этих функций через = 8рап{Ф^}^*.
Для действительных вектор-функций Фь = Фд = (й/1,#/г,£/г) € Н/ДОд) =
рассмотрим [2, 3, 4] дискретный аналог скалярного произведения
(Ф/г? Ф/г)/г = ^ ^ ^ ^ 81*1 (^г?) уНц^и^йц + 'ОцУц) + • (6)
*=1 ,7=0
Здесь через 5^ обозначена площадь г-го треугольного элемента, вершины которого занумерованы через 0, 1, 2, следовательно = /(А~ значение функции в j-й вершине г-го элемента.
В терминах скалярного произведения (6) сформулируем метод Бубнова - Галеркина [3, 41 следующим образом: для фиксированного момента времени найти вектор-функцию
Ипй №П(1 №П(1
ЛА,^)=Х>“ Ук(\,ср) = У£а]Ъ,(\,<р), £А(А,р) = 2>«Ф,-(А,¥>)
3—1 3—1 *=1
такую, что тождество
ал(Фл, Wfc) = /л№) + Ън{йн, ТУЛ) (7)
справедливо V W^l = (гу^,ш^,гу|) Е Н^. Для аппроксимации интегралов при получении билинейной а/1(*,*) и линейных /Л(-)? ^(0 форм использовалась формула трапеций и ее двумерный аналог.
Занумеровав узлы 0,^ от 1 до ЛГП(*, задачу (7) можно записать в векторно-
В работе [3] показан второй порядок сходимости решений в норме, порождаемой скалярным произведением (6) на равномерной сетке.
Для решения системы (8) использовался итерационный метод Якоби, который обладает хорошим параллелизмом, а диагональное преобладание для его сходимости легко обеспечивается выбором шага по времени т. В векторно-матричной форме метод Якоби запишется в следующем виде:
Здесь у — номер итерации, индексы (& + 1) для шага по времени опущены, однако компоненты глобальной матрицы жесткости и вектора правой части зависят от времени и должны пересчитываться в начале каждого временного шага.
Отметим некоторые особенности реализации, диктуемые методом конечных элементов. Глобальная матрица жесткости А зависит от времени и должна пересчитываться на каждом временном шаге. Однако для реализации метода Якоби на конечных элементах не требуется явного хранения глобальной матрицы жесткости. В программе насчитываются только элементы локальных матриц жесткости (причем только их диагональные элементы зависят от времени и перевычисляются на каждом временном шаге). Вычисление невязки АФМ —¥ в (9) производится по треугольникам с использованием элементов локальных матриц.
Сходимость метода определялась малостью разности Ф^ для двух соседних итераций в дискоретном аналоге равномерной нормы:
2. Параллельный алгоритм
2.1. Декомпозиция области
Метод конечных элементов для задачи (7) дает семиточечный шаблон, поэтому в методе
Используя явный параллелизм по данным, исходную расчетную область можно разбить на несколько пересекающихся только по границе разбиения подобластей. Расчеты в каждой подобласти можно проиводить независимо друг от друга в рамках итерации. Однако, после каждой итерации Якоби необходимо проводить согласование данных на границах подобластей. Имеет место, по крайней мере, два варианта разбиения.
1. Декомпозиция без перекрытий. Исходная область разрезается на взаимно неперекры-вающиеся полосы (рис. 3). На границах каждой подобласти невязка насчитывается только по части треугольников, принадлежащих области вляния точки, лежащей в подобласти. При обмене данными после каждой итерации Якоби требуется дополнительное суммирование для значений невязки в граничных в подобласти точках.
2. Декомпозиция с теневыми гранями. Исходная область включает взаимно перекрывающиеся подобласти, определяемые шаблоном (см. рис. 2). В этом случае невязка
матричной форме: для фиксированного момента времени £^+1 найти вектор
у&+1 = (г41?..., и^п^У 1,..., 1>лгп<р £г, • • • ? €мпЛ), удовлетворяющий системе линейных
алгебраических уравнений
^/с+1-у&+1 __ р&+1
(8)
(9)
тах |Ф^+1^ — Ф^| < е.
(10)
Якоби для вычисления значения Ф^+1^ на (и + 1)-й итерации требуются значения в семи точках Ф^ и-й итерации (рис. 2).
(у+1)-ая итерация Якоби
(2)
V -ая итерация Якоби
Рис. 1. Фрагмент триангуляции
Рис. 2. Зависимость в итерациях Якоби
в граничных точках для подобласти ьго процесса насчитываются в подобласти соседнего процесса. С учетом семиточечного шаблона и согласованности триангуляции достаточно перекрытия областей в два слоя расчетных точек (рис. 4). При обмене данными после каждой итерации Якоби дополнительная обработка не требуется, однако для каждого процесса необходимо хранить больше данных (теневые грани шириной в две расчетные точки).
Рис. 3. Декомпозиция без перекрытий Рис. 4. Декомпозиция с теневыми гранями
Таким образом, первый способ декомпозиции области более экономичен по памяти, прост в программировании, однако предполагает дополнительное суммирование на каждой итерации Якоби. Можно предположить, что его преимущества проявятся только на больших сетках с протяженными границами подобластей. Кроме того, очевидно его достоинство для неструктурированных сеток, когда границы подобластей не являются последовательным множеством точек.
2.2. Ускорение и масштабируемость
В рамках выбранной схемы распределения данных все процессы осуществляют одни и те же вычисления, но только над разными подобластями. Структура обменов также является однородной, за исключением первого и последнего процессов.
На рис. 5 приведен псевдокод параллельной программы на одном шаге по времени для случая декомпозиции области без перекрытий. При использовании декомпозиции с теневыми гранями на Этапе 3 необходимо сначала вычислить новое приближение решения, а затем произвести обмены без дополнительного суммирования.
Реализация параллельной программы осуществлялась на языке программирования Си с применением функций библиотеки передачи сообщений МР1.
Выполним теоретические оценки возможного ускорения каждого из параллельных алгоритмов, следуя [7].
К начало цикла по времени
Этап 1: Подготовка данных для временного шага
Локальные вычисления в своей двумерной подобласти:
I Расчет трения, внешних воздействий, коэффициентов локальных I матриц жесткости и правых частей, зависящих от времени;
Этап 2: Подготовка диагонали глобальной матрицы жесткости I Вычисление диагонали глобальной матрицы жесткости;
Обмен с соседями граничными значениями для диагонали глобальной матрицы жесткости
Этап 3: Определение решения для текущего шага по времени Локальные вычисления в своей двумерной подобласти:
I Пока не достигнута точность метода Якоби выполнять I { (начало цикла по итерациям метода Якоби)
I { (начало цикла по треугольникам)
I Вычисление невязки АФ-Р;
I } (конец цикла по треугольникам)
Обмен с соседями граничными значениями невязки.
Для декомпозиции без перекрытия подобластей -суммирование по границам подобластей
Локальные вычисления в своей двумерной подобласти:
I Вычисление нового приближения решения,
I Вычисление локальной ошибки в норме С_Ь;
Глобальная операция МР1_11е<1исе (... МР1_МАХ ...) для вычисления нормы ошибки
Локальные вычисления в своей двумерной подобласти:
I } (конец цикла по итерациям метода Якоби)
1} (конец цикла по времени)
Рис. 5. Псевдокод фрагмента программы Предположим, что все операции, в том числе и обмены, имеют одинаковое время исполнения £. Пусть ИШ1 - общее количество точек сетки расчетной области, к - количество
операций, выполняемых в одной расчетной точке, в - количество шагов по времени. Тогда
общий объем вычислений Vcalc В алгоритме определяется соотношением Vcalc = skNndj а время выполнения алгоритма на одном процессоре соответственно можно оценить следующим образом: Т\ ~ skNndt. Тогда средняя степень параллелизма г данного алгоритма согласно [7] равна
skNnd АТ
г-------г = Nnd-
sk
Потенциальное ускорение алгоритма оценивается как отношение времени вычисления
-^1
на одном процессоре Т\ к времени вычислений на р процессорах Тр: Sp = —. Выпол-
Тр
ним теоретические оценки возможного ускорения каждого из параллельных алгоритмов, по возможности учитывая время, затрачиваемое при выполнении алгоритма на обмены и возможные накладные расходы на вычисления в перекрывающихся подобластях:
■Sp = 7-, Тт-----• (И)
х 1 / у “Г -L over | comm
Здесь Т\ - время вычислений на одном процессоре, Tover - время на дополнительные вычисления, связанные с декомпозицией области, Tcomm - время, требуемое для обменов.
Как следует из принятой нами схемы распределения данных, на каждой итерации метода Якоби требуется обмен граничными элементами вектора решений и порождаемые распределенностью данных дополнительные вычисления. В обоих случаях декомпозиции области для вычисления невязки на каждом разрезе количество точек, в которых необходим обмен данными, соответствует количеству точек сетки на разрезе. Однако при декомпозиции области с теневыми гранями согласно шаблону метода Якоби (см. рис. 2) требуется хранить не только точки на разрезе, но и все точки, принадлежащие треугольникам, в которые входят точки разреза. Обозначим через Nbnd количество точек сетки, в которых необходим обмен данными, далее такие точки будем называть граничными в подобласти.
Объем дополнительных вычислений, связанных с декомпозицией области зависит от количества процессоров р, на которых выполняется алгоритм, и количества дополнительных вычислений д, которые необходимы для одной граничной в подобласти точки. Источники этих вычислений при различных способах декомпозиции области отличаются. При декомпозиции без теневых граней каждый процесс, вычисляя невязку в граничной в подобласти точке, ’’обрабатывает” свою часть треугольников, не дублируя вычисления. В этом случае после обмена требуется дополнительное суммирование частей невязок, насчитанных в соседнем процессоре. При декомпозиции области с теневыми гранями соседние процессы дублируют друг друга при расчете невязки для точек перед разрезом. Таким образом, объем дополнительных вычислений оценивается следующим образом:
Vover ^ s9^bnd(P !)•
При этом время, затрачиваемое на дополнительные вычисления в Nbnd граничных точках подобласти можно оценить так:
Tover ~ SQNbndt'
Пусть также для каждой граничной точки требуется пересылка m значений. Тогда объем пересылаемых данных соседним процессам можно оценить следующим образом:
Vcomm
~ 2smNbnd.
Для оценки времени, необходимого для обменов, рассмотрим их реализацию в MPI. При использовании библиотеки MPI возможно два принципиально разных способа организации обменов - с использованием блокирующих или неблокирующих функций приема и передачи данных. Оценим ускорение в обоих случаях.
Блокирующие передачи. На рис. 6 отображено влияние блокирующих передач на производительность на примере восьми процессов [8]. На первом этапе шестой процесс передает данные седьмому процессу, в то время как процессы с нулевого по пятый простаивают в ожидании приема данных своими правыми соседями. На втором этапе пятый процесс передает данные шестому, а процессы с нулевого по четвертый простаивают. Таким образом, все обмены завершатся за семь этапов.
В общем случае обмены с использованием блокирующих передач потребуют Тсотт(р— 1) единиц времени, где Тсотт = 2зтЩп^ ~ время на обмен данными с одной границы, р -количество используемых процессоров.
Время
РО P1 Р2 РЗ Р4 Р5 Р6 Р7 Р
Рис. 6. Влияние блокирующих передач на ускорение
С учетом (11) ускорение при использовании блокирующих передач с учетом всех оценок будет определяться следующим соотношением:
1 о ^ Л/ га ’
- + -Д + 2(р-1)-Д
1
1 о „ ,.т„’
жл + кв+2(з’-1)тв
1 <Р< Nnd, Р > Nnd.
(12)
Здесь через R обозначено отношение количества граничных точек в подобласти к общему числу точек расчетной области: R = N}jnci/Nn(i.
Неблокирующие передачи. В общем случае время, затрачиваемое на обмены с использованием неблокирующих передач, не зависят от количества участвующих в обменах процессов: Tcomm = const. Оценка (11) в случае неблокирующих обменов имеет вид:
1
1 д _ 2т ’
—Ь "г-R + —
р к kNnd
1 д „ 2т + vR +
1 <р< Nnd, Р > Nnd.
(13)
Для получения численных оценок ускорения оценим количество операций к на одну расчетную точку. Каждый шаг по времени включает в себя: 1) вычисления трения, правой части, диагональных элементов локальных матриц жесткости ~ 60ЖП^ операций; 2) вычисление диагонали глобальной матрицы жесткости ~ 457Уе/ операций (Л^ - количество треугольных элементов в триангуляции расчетной области); 3) решение системы (8)
методом Якоби (9). Реализация метода Якоби со сборкой невязки по треугольникам требует 1) вычисления невязки ~ 180операций; 2) З^Уе| проверки на принадлежность точки границе; 3) ~ 60А^ дополнительных операций для вычисления интегралов в Щ граничных точках расчетной области; 4) Мпд операций для умножения на 5) ЗЛУ^ проверок для вычисления глобального максимума.
В табл. 1 приведены результаты теоретических оценок ускорения параллельного алгоритма в зависимости от способа декомпозиции области, типа обменов и количества процессов для модельной прямоугольной области с сеткой 401x401 точек. Здесь Ып(л = 160801, N^1 = 320000. Из табл. 1 видно, что наиболее эффективной является реализация, основанная на разбиении области без перекрытия подобластей. Полученные результаты показывают, что алгоритм обладает значительным объемом потенциального параллелизма и хорошей с точки зрения распараллеливания структурой, что дает ускорение близкое к линейному в зависимости от количества используемых процессоров.
Таблица 1
Значения оценки ускорения в зависимости количества процессов
Тип декомпозиции Р 1 2 4 8 12 16
Без перекрытия Без перекрытия Я? СипЫ Р 1.0000 1.0000 1.9999 1.9999 3,9986 3,9995 7,9866 7,9981 11,9526 11,9957 15,8855 15,9923
С теневыми гранями С теневыми гранями S« gunbl 1.0000 1.0000 1.9962 1.9962 3,9830 3,9849 7,9173 7,9399 11,7814 11,8653 15,5554 15,7615
3. Вычислительный эксперимент
3.1. Модельная область
Для численного исследования ускорения параллельного алгоритма рассмотрим следующую модельную задачу. Пусть (2 — «квадрат» на сфере: Q = [0,7г/10] х [7г/2,7г/2 + 7г/10]. Границы fl считаются «твердыми». В О рассмотрена задача с известным точным решением [3]. В расчетной области построены две равномерные квадратные сетки 401x401 и 801x801 точек с соответствующей согласованной триангуляцией. В вычислительных экспериментах было сделано 1000 шагов по времени.
Вычисления выполнялись на 96-процессорном кластере ИВМ СО РАН. Кластер МВС-1000/ИВМ (собственная сборка ИВМ СО РАН, [11]) содержит 24 вычислительных узла AMD Athlon64/3500+/irb (однопроцессорные, одноядерные); 12 вычислительных узлов AMD Athlon64 Х2 Dual Соге/4800+/2Гб (однопроцессорные, двухъядерные); 12 вычислительных узлов 2 X Dual-Core AMD Opteron Processor 2216/4Г6 (двухпроцессорные, двухъядерные). Управляющий узел, сервер доступа и файловый сервер Athlon64/3500+/1Gb с общей дисковой памятью 400 Гб под управлением ОС Gentoo Linux. Управляющая сеть FastEthernet (100 Мбит/сек), сеть передачи данных GigaEthernet (1000 Мбит/сек). Необходимо отметить, что поскольку кластер является гетерогенным, то временные характеристики выполнения программы осреднялись по результатам нескольких десятков расчетов.
На рис. 7 представлена зависимость ускорения вычислений от количества используемых процессов для модельных сеток. Чтобы не загромождать рисунок, не отображены графики ускорений, полученных для декомпозиции с теневыми гранями с блокирующими обменами, поскольку они практически совпадают с графиками ускорения для декомпозиции с теневыми гранями с неблокирующими обменами.
Анализ вычислительных экспериментов показывает, что при сравнительно небольшом количестве процессоров (Р < 12) значение ускорения остается одинаковым для всех рассматриваемых вариантов, при этом можно отметить незначительное преимущество декомпозиции без теневых граней.
В алгоритме пока не используется возможность совмещения вычислений и обменов, поэтому влияние блокирующих и неблокирующих вариантов пересылок должно быть минимальным. В то же время, на сетке размерностью 401x401 блокирующие обмены для декомпозиции без теневых граней не дают рост ускорения, начиная с вычислений на 12-ти процессах. В остальных случаях на сетке 401x401 рост ускорения наблюдается вплоть до использования 22-х процессов. Наличие эффективности больше единицы в районе 18-ти -22-х процессов объясняется, по-видимому, попаданием в кэш. На сетке 801x801 при общем увеличении ускорения вплоть до использования 20 процессов отмечается «провал» при 16-ти процессах. Неустойчивость при количестве процессов больше 22-х говорит о преобладании в этих случаях времени обменов над вычислениями, архитектурными особенностями кластера и реализацией МР1.
Модельная сетка 401x401 -
6 8 10 12 14 16 18 20 22 24 26 28 30 32
Количество процессов
—\-----Г"
2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32
Декомпозиция без перекрытия, неблокирующее обмены Количество процессов
■ .....« Декомпозиция без перекрытия, обмены вепйВесу
—»">*— Декомпозиция с теневыми гранями, неблокируюи^е обмены
Рис. 7. Зависимость ускорения вычислений от количества доступных процессоров
В работах [5, 6] приведены некоторые результаты, характеризующие для этой же задачи эффективность распараллеливания кода с использованием языка Гох^гап-БУМ при проведении численных экспериментов на равномерной модельной сетке и в акватории Мирового океана. Декомпозиция области в работе проводилась автоматически и использовала технологию теневых граней. Проводя сравнение с этими результатами, заметим, что при использовании 8-ми и более процессов программирование параллельных алгоритмов при помощи функций МР1 показывает несколько более высокие значения ускорения вычислений. Однако, данных для добросовестного сравнения возможностей МР1 и БУМ на данной задаче не достаточно.
3.2. Вычислительные эксперименты для акватории Охотского моря
Тестовые расчеты для акватории Охотского моря проводились на сетках, подготовлены* С.Ф.Пятаевым и И.В. Киреевым на основе открытой батиметрической базы данных ЕТ0Р02 [6, 10].
В рассмотренной сетке число узлов 7\ГП^ = 43 768, а число треугольных элементов = 78 929 (рис. 8 а).
Численный эксперимент соответствует начальным данным с локальным подъемом уровня, описываемым гауссовской функцией
£(0, А, <р) = Аехр(—(А - А0)2/(2В)2 щ)2/(20)2).
В расчетах было принято А = 10, Ао = 149,1° восточной долготы, а щ = 53,1° северной широты, В = 0,005, коэффициент трения г* = 0,0026. Для скоростей предполагались нулевые начальные данные: и(0, Л, ф) = ?;(0, А, ф) = 0.
В построенной сетке из общего числа граничных участков около 6.4% проходят по морю. На таких участках Ш2 = 1, ив краевом условии (2) полагали ^ = 1и й(£, А, ф) = 0.
Результаты численного моделирования показаны на рис. 8 б) - г), где представлена функция £(£, А, ф) для некоторых моментов времени. Линиями показаны границы, получившиеся при декомпозиции области на 8 частей. Поскольку сетка сильно сгущается в приграничных областях, то ширина полосы, которая отводится одному процессору, в прибрежной зоне значительно уже, чем в море. В табл. 2 приведены результаты численных экспериментов на различном числе процессоров. Здесь верхние индексы «Ы» и «ипЫ» соответствуют экспериментам с блокирующими и неблокирующими операциями обмена соответственно.
Таблица 2
Результаты численных экспериментов для задачи в акватории Охотского моря
Количество процессоров 1 2 4 8
Время с Ускорение Эффективность Еы 6,8729 1,0000 1,0000 4,0690 1,6891 0,8446 2,3030 2,9843 0,7461 1,6324 4,2103 0,5263
Время £“п6*, с Ускорение 5'гшЫ Эффективность ЕипЫ 7,0290 1,0000 1,.0000 4,1209 1,7057 0,8528 2,3389 3,0053 0,7513 1,5304 4,5928 0,5741
Из табл. 2 видно, что ускорение вычислений увеличивается пропорционально количеству доступных процессоров. При этом время вычислений при использовании неблокирующих передач незначительно меньше по сравнению со временем вычислений при наличии блокирующих передач.
Следует отметить, что декомпозиция области с теневыми гранями в этом случае является трудоемким процессом, поэтому не проводилась.
Сравнение с результатами, полученными при применении технологии БУМ [6], дает преимущество наших алгоритмов при большом количестве процессоров.
4. Заключение
На примере численного решения краевой задачи для уравнений мелкой воды в работе рассмотрены технологические аспекты разработки масштабируемых параллельных алгоритмов для кластерных вычислительных систем с использованием библиотеки МРІ.
Рис. 8. Численные эксперименты в акватории Охотского моря на 8 процессорах: а) общий вид сетки; б) начальное возмущение; в) через 42 мин; г) через 67 мин
Поскольку при дискретизации задачи использовался метод конечных элементов с организацией вычислений по треугольным элементам, было рассмотрено два естественных подхода к декомпозиции области - без перекрытий и с теневыми гранями.
Теоретические оценки показали, что алгоритм обладает значительным объемом потенциального параллелизма и хорошей, с точки зрения распараллеливания, структурой, что дает ускорение, в зависимости от количества используемых процессоров, теоретически близкое к линейному.
Численные эксперименты показали, что использование неблокирующего режима обменов, которое допускается алгоритмом, является безусловно более эффективным. В дальнейшем планируется повысить эффективность алгоритма, в том числе и за счет совмещения вычислений с неблокирующими операциями обмена.
Следует отметить как явное преимущество простоту организации параллельного алгоритма при декомпозиции без перекрытия подобластей на неструктурированных триангуляциях реальных акваторий.
Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант 08-01-00621-а) и Президентской программы «Ведущие научные школы РФ» (грант НШ-3431.2008.9). Статья рекомендована к печати программным комитетом международной научной конференции «Параллельные вычислительные технологии 2009>.
Литература
1. Марчук, Г.И. Динамика океанских приливов / Г.И. Марчук, Б.А. Каган. - Л.: Гидро-метиздат, 1983.
2. Agoshkov, V.I. Inverse problems of the mathematical theory of tides: boundary-function problem / V.I. Agoshkov // Russ. J. Numer. Anal. Math. Modelling. - 2005. - Vol. 20, № 1.
- P. 1 - 18.
3. Kamenshchikov, L.P. Simulation of surface waves in basins by the finite element method / L.P. Kamenshchikov, E.D. Karepova , V.V. Shaidurov // Russian J. Numer. Anal. Math. Modelling. - 2006. - Vol. 21, №4. - P. 305 - 320.
4. Karepova, E.D. Numerical Solution of the Boundary Problem foe Shallow Water Equations for Modelling Surface Waves in World Ocean by Finite Elements Methods / L.P. Kamenshchikov, E.D. Karepova , V.V. Shaidurov // Finite Difference Methods: Theory and Applications. Proceedings of Fourth International Conference FDM:T&A’06. - Bulgaria, 2007. -P. 227 - 233.
5. Каменщиков, Л.П. Моделирование поверхностных волн в водоемах методом конечных элементов на вычислительном кластере / Л.П. Каменщиков, Е.Д. Карепова, В.В. Шайдуров // Избранные матералы Четвертой школы-семинара «Распределенные и кластерные вычисления». - Красноярск, 2005. - С. 114 - 125.
6. Моделирование гравитационных волн в Мировом океане методом конечных элементов с распараллеливанием / Л.П. Каменщиков, Е.Д. Карепова, С.Ф. Пятаев, В.В. Шайдуров // Избр. матералы Шестой школы-семинара «Распределенные и кластерные вычисления». - Красноярск, 2006. - С. 52 - 64.
7. Ортега, Дж. Введение в параллельные и векторные методы решения линейных систем: пер. с англ. / Дж. Ортега. - М.: Мир, 1991.
8. Efficient Management of Parallelism in Object-Oriented Numerical Software Libraries, Modern Software Tools in Scientific Computing / S. Balay, W.D. Gropp, L.C. Mclnnes and others. - Birkhauser Press, 1997. - P. 163 - 202.
9. McBryan, O.A. An overwiev of message passing environments / O.A. McBryan // Parallel Computing. - 1994. - V 20. - P. 417 - 441.
10. National Geophysical Data Center, http://www.ngdc.noaa.gov/ngdc.html
11. Исаев, С.В. Развитие Красноярского центра параллельных вычислений / С.В. Исаев, А.В. Малышев, В.В. Шайдуров // Вычислительные технологии. - Новосибирск, 2006. -Т. 11, спецвып. - С. 28 - 33.
Отдел вычислительной математики,
Институт вычислительного моделирования СО РАН j апе<Шст. кг авп. га
Поступила в редакцию 2 марта 2009 г.
