CC BY
68
Вычислительные технологии
Том 14, № 6, 2009
Параллельная реализация метода конечных элементов для
kj «_> >к
начально-краевой задачи мелкой воды
Е. Д. КАРЕПОВА, В. В. Шлйдуров Учреждение Российской академии наук Институт вычислительного моделирования СО РАН, Красноярск, Россия e-mail: [email protected], [email protected]
Проведено исследование эффективности нескольких параллельных реализаций алгоритма численного решения начально-краевой задачи для уравнений мелкой воды, выполненных с помощью библиотеки MPI для языка Си. Рассмотрены два подхода к декомпозиции вычислительной области и две схемы реализации двухточечных обменов в алгоритме. Приведены сравнительные результаты ускорения вычислений в зависимости от количества процессов, способа реализации коммуникаций, способа декомпозиции вычислительной области, архитектуры кластерной системы.
Ключевые слова: модели мелкой воды, декомпозиция согласованной триангуляции, метод конечных элементов, блокирующие/неблокирующие передачи, ускорение вычислений.
1. Моделирование поверхностных волн в водоемах методом конечных элементов
Модели мелкой воды хорошо описывают большой круг природных явлений, таких как крупномасштабные поверхностные волны, возникающие в морях и океанах, цунами, приливные течения, поверхностный и русловой сток, гравитационные колебания поверхности океанов [1, 2]. В работах [2-4] рассмотрено численное моделирование поверхностных волн в больших акваториях с учетом сферичности Земли и ускорения Кориолиса на основе уравнений мелкой воды. В [2] для дифференциальной постановки задачи выведены априорные оценки, обеспечивающие устойчивость решения и однозначную разрешимость задачи. В [3, 4] для этой же задачи построен метод конечных элементов (МКЭ), для которого получены необходимые априорные оценки. Там же приведены результаты численных экспериментов на модельных сетках для акваторий Охотского моря и Мирового океана.
Поскольку для расчетов на реальных акваториях объем вычислений велик, то оправданно применение высокопроизводительных ЭВМ. В [7] приведены некоторые данные по ускорению параллельного алгоритма для данной задачи, требующие некоторого уточнения. Цель данной работы — тщательное исследование эффективности реализации параллелизма по данным для рассматриваемой задачи и особенностей поведения кластерных систем различной архитектуры при расчетах.
*Работа выполнена при поддержке РФФИ (грант № 08-01-00621).
© ИВТ СО РАН, 2009.
Рассмотрим задачу в следующей постановке. Для стандартной сферической системы координат (г, Л, в) с началом в центре земного шара будем использовать вместо угла в географическую широту ^ = п — в, так что 0 < ^ < п. Через Л обозначена географическая долгота 0 < Л < 2п. Полагаем всюду г = Яе, где Яе — радиус Земли, который считается постоянным. Пусть П — заданная область на сфере с границей Г=Гі и Г2, где Г ! — часть границы, проходящая вдоль берега, а Г2 = Г \ Г! — часть границы, проходящая по морю. Обозначим через т! и т2 характеристические функции соответствующих участков границы. Для простоты считается, что полюса ^ = 0 и ^ = п не входят в П. Относительно неизвестных функций и = и(і, Л, <^), V = г>(£, Л, <^) и £ = £(£, Л, <^) запишем в П х (0,Т) уравнения баланса импульсов и уравнение неразрывности [2]:
ди
— = Ь + т^ дЛ — Я/и + дv
— = —1и + — Я/V + ^
где u, v — компоненты вектора скорости U по осям Л и ^ соответственно; £ — отклонение свободной поверхности от невозмущенного уровня; H(Л, <^) > 0 — глубина водоема в точке (Л, <^); функция Rf = r*|U|/H учитывает силу трения о дно, г* — коэффициент трения; l = —2w cos ^ — параметр Кориолиса; m = 1/(RE sin <^); n = 1/Re; g — ускорение силы тяжести; / = /i(t, Л, <^), /2 = /2(t, Л, <^) и /3 = f3(t, Л, <^) — заданные функции внешних воздействий.
Граничные условия рассмотрим в следующем виде [2]:
HUn + em^\/gH£ = m2\/gHd на Г х (0,T), (2)
где Un = U • n, n = (ni,—n2) — вектор внешней нормали к границе; в £ [0,1] —
m
заданный параметр, d = d(t, Л, <^) — заданная граничная функция, определенная на границе Г2.
Зададим также начальные условия
u(0, Л, = u0(Л, ^ v(0, Л, ^ = v0(Л ^ £(0, Л, = £0(Л ^)- (3)
Рассмотрим некоторую согласованную триангуляцию T = (wi }|^=£1 области П, состоящую из невырожденных треугольников с прямолинейными сторонами в координатах Л и ^ и содержащую область П. Сетка в общем случае может быть неструктурированной, но обязательно является согласованной, т. е. каждое ребро является либо стороной двух треугольников, либо стороной одного треугольника, находящейся на границе области (рис. 1). Пусть Пh — множество узлов (т. е. вершин треугольных элементов) общим числом Nnd, П — множество внутренних узлов. Для каждого узла Zj £ Пh введем базисную функцию Фj• (Л, <^), которая равна единице в Zj, равна нулю во всех остальных узлах Пh и линейна в каждом треугольнике. Обозначим линейную оболочку этих функций через Hh(nh) = span{ Ф,}^ •
Для действительных вектор-функций Фh = (uh,vh,£h), Фh = (uh, vh, £h) £ Hh(nh) = (Hh(nh))3 рассмотрим [2-4] дискретный аналог скалярного произведения
Nei i 2
(ФЬ., Ф h)h ^ ] 3 Si ^ ] RE sin(^ij ) (^Hij (uij uij + vij vij ) + g£ij£j • (4)
i=1 j=0
Рис. 1. Фрагмент триангуляции
Здесь через Бг обозначена площадь ¿-го треугольного элемента, вершины которого занумерованы через 0, 1, 2; следовательно /у = / (Ау, ) — значение функции в ]-й
вершине ¿-го элемента.
Для дискретизации по времени системы (1)-(2) разобьем временной отрезок [0,Т] на К интервалов: 0 = ¿0 < ¿1 < ••• < ¿к = Т с шагом т = Т/К. Аппроксимируем производные по времени левыми разностями. В результате на каждом временном интервале (¿^, ¿£+1) получим полудискретный аналог системы (1)-(2), для которого в терминах скалярного произведения (4) сформулируем метод Бубнова—Галеркина [3, 4]. Для аппроксимации интегралов при получении билинейной и линейных форм используем формулу трапеций и ее двумерный аналог.
Занумеровав узлы П^ от 1 до ^„¿, запишем задачу в векторно-матричной форме: для фиксированного момента времени ¿£+1 найдем вектор Vfc+1 = (и1,..., и^пЛ, ^1,... , ,
£1,... , ), удовлетворяющий системе линейных алгебраических уравнений
В работе [3] показан второй порядок сходимости решений в норме, порождаемой скалярным произведением (4) на равномерной сетке.
Для решения системы (5) использовался итерационный метод Якоби, который обладает хорошим параллелизмом, а диагональное преобладание для его сходимости легко обеспечивается выбором шага по времени т. В векторно-матричной форме метод Якоби запишется в следующем виде:
Здесь V — номер итерации, индексы (к +1) для шага по времени опущены, однако компоненты глобальной матрицы жесткости А и вектора Е правой части выражения зависят от времени и должны пересчитываться в начале каждого временного шага.
Отметим некоторые особенности реализации, диктуемые методом конечных элементов. Глобальная матрица жесткости А зависит от времени и должна пересчитываться на каждом временном шаге. Однако для реализации метода Якоби на конечных элементах не требуется явного хранения глобальной матрицы жесткости. В программе насчитываются только элементы локальных матриц жесткости (причем лишь их диагональные элементы зависят от времени и перевычисляются на каждом временном шаге). Вычис-
А£+^£+1 = е^+1
(5)
(6)
ление невязки АфМ - F в (6) производится по треугольникам с использованием элементов локальных матриц. Сходимость метода определялась малостью разности Ф(^ для двух соседних итераций в дискретном аналоге равномерной нормы:
max |Ф,(^+1) - Ф,М|< е. (7)
2. Параллельный алгоритм
Используя явный параллелизм по данным, исходную расчетную область можно разбить на несколько частично перекрывающихся подобластей. Расчеты в каждой подобласти выполняются независимо друг от друга в рамках итерации Якоби. После каждой итерации Якоби необходимо проводить согласование данных в перекрытиях. Имеют место по крайней мере два следующих варианта разбиения.
1. Декомпозиция с теневыми гранями. Исходная область включает взаимно перекрывающиеся подобласти, определяемые шаблоном. При этом невязка в граничных точках подобласти i-го процесса насчитывается в подобластях соседних процессов. С учетом семиточечного шаблона и согласованности триангуляции достаточно перекрытия областей в два слоя расчетных точек. Поскольку алгоритм вычисления невязки требует хранения в каждой граничной точке подобласти семи коэффициентов матрицы жесткости, трех значений вектора решения текущей и предыдущей итераций и одного значения правой части, то избыточность хранения информации для р-процессов составляет 56(р — 1)Nbnd* SizeOfDouble байт. Здесь через Nbnd обозначено количество точек на разрезе, SizeOfDouble — количество байт, занимаемых переменной, хранимой с двойной точностью.
2. Декомпозиция без перекрытий. Исходная область разрезается на подобласти, пересекающиеся только по границам разреза. Для каждой граничной точки подобласти невязка насчитывается частично только по тем треугольникам, которые лежат в подобласти. При обмене данными после каждой итерации Якоби требуется дополнительное суммирование для значений невязки в граничных для подобласти точках. Этот способ декомпозиции более экономичен по памяти, прост в программировании, однако предполагает дополнительные арифметические операции на каждой итерации Якоби. Очевидно его достоинство для неструктурированных сеток, когда границы подобластей не являются последовательным множеством точек.
В рамках выбранной схемы распределения данных все процессы осуществляют одни и те же вычисления, но над разными подобластями. Структура обменов, за исключением первого и последнего процессов, также является однородной. Обмены необходимы на каждой итерации метода Якоби.
Реализация параллельной программы осуществлялась на языке программирования Си с применением функций библиотеки передачи сообщений MPI.
Выполним теоретические оценки возможного ускорения каждого из параллельных алгоритмов, следуя [8].
Обозначим время выполнения одной арифметической операции t^, а время выполнения пересылки одного значения — tcomm. Пусть Nnd — общее количество точек сетки расчетной области, s — количество операций, выполняемых в одной расчетной точке на каждой итерации Якоби, k — количество шагов по времени, v — среднее количество итераций Якоби на каждом временном шаге. Тогда общий объем вычислений Vcaic в ал-
горитме определяется соотношением Уса1с = а время выполнения алгоритма на
одном процессоре соответственно можно оценить следующим образом: Т к^Ы„^ор.
Потенциальное ускорение алгоритма оценивается как отношение времени вычисле-
Т\
ния на одном процессоре Т к времени вычислений на р-процессах Тр: Бр = —. Вы-
Тр
полним теоретические оценки ускорения для каждого рассматриваемого случая декомпозиции области, учитывая время, затрачиваемое на выполнение обменов, а также на возможные накладные расходы при вычислениях в перекрывающихся подобластях.
Предположим, что при декомпозиции области нам удается равномерно распределить весь объем вычислений Уса1с по процессорам. Тогда для ускорения параллельного алгоритма можно записать следующее:
Т±
Т1 /р + Таует + Тс<
Здесь Т1 — время вычислений на одном процессоре, Тоуег — время на дополнительные вычисления, связанные с декомпозицией области, Тсотт — время, требуемое для обменов.
Как следует из принятой нами схемы распределения данных, на каждой итерации метода Якоби требуются:
1) глобальная операция приведения для вычисления равномерной нормы разности между двумя последовательными приближениями итераций Якоби;
2) обмен граничными элементами вектора решений на разрезах;
3) дополнительные вычисления, порождаемые распределенностью данных.
Пусть д — количество дополнительных вычислений, связанных с декомпозицией области, которые необходимы для одной граничной в подобласти точки. Декомпозиция с теневыми гранями не требует прямых дополнительных вычислений, т. е. д = 0. При декомпозиции без теневых граней каждый процесс, вычисляя невязку в граничной в подобласти точке, “обрабатывает” свою часть треугольников. В этом случае после обмена требуется дополнительное суммирование частей невязок, насчитанных в соседнем процессоре. Поскольку в каждой точке вычисляется невязка для трех компонент вектора (и, г>,£), то д = 3. В результате объем дополнительных вычислений Уоуег ~ 2&^дЛгьга^(р — 1). Поскольку суммирование процессы выполняют асинхронно и независимо, то время, затрачиваемое на дополнительные вычисления в Ыьпй граничных точках подобласти, можно оценить следующим образом:
op•
Для вычисления критерия остановки итерационного процесса на каждой итерации требуется вычислить глобальный максимум по всем процессам. Предположим, что реализация глобальных операций приведения в MPI выполняется по оптимальному алгоритму сдваивания, что дает время выполнения одной операции приведения Tc1omm = (top + tcamm) log2 p. Пусть также для каждой граничной точки требуется пересылка m значений (в нашем случае m =3). Объем пересылаемых данных каждым процессом соседу можно оценить как
Vcomm ~ kvmNbnd.
При использовании библиотеки MPI возможно два принципиально разных способа организации обменов — с использованием блокирующих или неблокирующих функций передачи данных. Оценим время, необходимое для обменов в обоих случаях.
Блокирующие передачи. На рис. 2 показана схема реализации обменов по цепочке процессов с помощью функций совмещенных приема-передачи МР!_БепНгесу(...), отражающая влияние блокирующих передач на производительность на примере восьми процессов. На схеме операции отправки данных правому и левому соседу обозначены через БпН_Р и БпН_Ь соответственно, а операции приема от правого и левого соседа — через Рсу_Р и Рсу_Ь. Сначала все процессы, кроме последнего, отправляют данные своим правым соседям и от них же ожидают поступления данных. Затем все процессы, кроме первого, посылают данные своим левым соседям и ожидают поступления данных от них же. Поскольку на первом этапе последний процесс не передает данные вправо, то он успешно выполняет операцию приема от шестого (левого) процесса и передает ему свои данные, в то время как процессы с нулевого по пятый простаивают в ожидании приема данных своими правыми соседями. На втором этапе шестой процесс получает данные от пятого и осуществляет передачу ему своих данных, а процессы с нулевого по четвертый простаивают. Таким образом происходит разрешение всех обменов по цепочке. Все обмены завершатся за р + 1 этап.
В общем случае обмены с использованием блокирующих передач потребуют 2 (р — 1) пересылок Усотт единиц данных, где р — количество используемых процессов. Тогда Тсотт = 2(р — 1)китЩ^¿сотт — общие затраты на обмены в блокирующем режиме.
С учетом (8) и всех проведенных оценок ускорение при использовании блокирующих передач будет определяться следующим соотношением:
SbJ =
1
1 д
+ 2- Я + р 8 вМш1
1082 Р(1 + к) + 2(р - 1) ™Я*
8
(9)
Здесь Я — отношение количества граничных точек в подобласти к общему числу точек расчетной области: Я = Щпл/Мпл, к — отношение времени пересылки одного значения к времени выполнения одной арифметической операции: к = Ьсотт/Ьор.
Неблокирующие передачи. В общем случае время, затрачиваемое на обмены с использованием неблокирующих передач, не зависит от количества участвующих
Рис. 2. Схема реализации блокирующих обменов
в обменах процессов: ТС2отт = 2 Уотт^сотт = 2кут^ЬпЛЬССутт. Оценка (8) в случае неблокирующих обменов имеет вид
5Г“ = 1----------------------------------------------------------1-1-■ (10)
' 1 + 2 »Я +1°|2Р (1 + к)+ 2 тЯк
р 5 в^па в
Из (9), (10) следует, что для достаточно мелких сеток потенциальное ускорение близко к линейному на достаточно большом диапазоне количества процессов. Величина ускорения определяется двумя параметрами. Первый из них Я = N^/^4 характеризует декомпозицию расчетной области. Приемлемое ускорение обеспечивается малостью Я, поэтому при построении декомпозиций сложных вычислительных областей наряду с требованием равенства вычислительной нагрузки на процесс необходимо обеспечивать минимальную протяженность границы подобласти для каждого процесса. Второй параметр к = 1сотт/1ор характеризует коммуникационную среду. Этот параметр описывает вычислительную сеть очень условно, но он показывает, что для приемлемого ускорения следует выбирать архитектуру кластера с небольшими значениями к. Расчеты показали, что неоднородность архитектуры высокопроизводительной системы приводит к непредсказуемому поведению к в зависимости от количества процессов, конкретных вычислительных узлов, на которых будет выполняться задача, сложившейся общей ситуации на кластере (т.е. задач других пользователей). Отметим также, что величина Тоюег на несколько порядков меньше, чем другие слагаемые знаменателя в (8) и не зависит от количества процессов. С учетом экономии памяти и легкости реализации на неструктурированных сетках это дает преимущество декомпозиции без перекрытий. Кроме того, неблокирующие передачи, даже в ситуации, когда нет совмещения вычислений и обменов, несколько предпочтительней блокирующих пересылок.
3. Вычислительный эксперимент
3.1. Модельная область
Для численного исследования ускорения параллельного алгоритма рассмотрим следующую модельную задачу. Пусть П — “квадрат” на сфере: П = [0,п/10] х [п/2, п/2 + п/10]. Границы П считаются “твердыми”. В П рассмотрена задача с известным точным решением [3]. В расчетной области построены две равномерные квадратные сетки 401x401 и 801x801 точек с соответствующими согласованными триангуляциями. В вычислительных экспериментах было сделано 1000 шагов по времени.
Вычислительный эксперимент проведен на трех высокопроизводительных комплексах.
Во-первых, вычисления выполнялись на 99-ядерном кластере Института вычислительного моделирования СО РАН. Кластер МВС-1000/ИВМ (собственная сборка ИВМ СО РАН [12]) содержит 27 вычислительных узлов AMD Athlon64/3500+/1K (однопроцессорные, одноядерные); 12 вычислительных узлов AMD Athlon64 X2 Dual Соге/4800+/2Гб (однопроцессорные, двухъядерные); 12 вычислительных узлов 2XDual-Core AMD Opteron Processor 2216/4Гб (двухпроцессорные, двухъядерные). Управляющий узел, сервер доступа и файловый сервер — Athlon64/3500+/1K с общей дисковой памятью 400 Гб под управлением ОС Gentoo Linux. Управляющая сеть — FastEthernet, сеть передачи данных — GigaEthernet.
Вторым кластером, на котором проведены серии расчетов, является кластер Томского государственного университета SKIF Cyberia, содержащий 283 двухпроцессорных двухъядерных узла (1132 ядра) IntelXeon5150 2.66 Ггц с суммарным объемом дисковой памяти 1136 Гб и объемом дискового пространства 22.56 Тб. Внешняя дисковая система хранения данных имеет объем 10 Тб. Все узлы объединены высокопроизводительной сетью передачи данных InfiniBand.
Наконец, эксперименты были частично повторены на высокопроизводительном вычислительном комплексе Информационного вычислительного центра Новосибирского государственного университета, на базе кластера из 64 восьмиядерных вычислительных узлов, основанных на блэйд-серверах Hewlett-Packard BL460c. Все вычислительные узлы объединены высокопроизводительной сетью InfiniBand DDR 4x. В кластере используется параллельная файловая система хранения SFS емкостью 24 Тб, основанная на технологии Lustre, скорость чтения СХД — более 1 Гб/с, записи — 800 Мб/с.
Отличительной чертой первого кластера является его гетерогенная архитектура. Кластер SKIF Cyberia и HP-кластер ИВЦ НГУ отличаются в интересующем нас плане количеством ядер на узел.
Исследования проводились на грубой и мелкой сетке — 401 х 401 и 801 х 801 точек соответственно. Ясно, что потенциальное ускорение параллельной программы зависит от размерности сетки. Исследования зависимости ускорения проводились вплоть до 32 процессов, где на рассматриваемых сетках алгоритм хорошо масштабируем. Временные характеристики замерялись средствами MPI каждым процессом в отдельности, за измеряемый показатель принималось максимальное значение. Все временные характеристики осреднялись по результатам нескольких десятков расчетов, исключая экстремальные значения. Дисперсия измерений, кроме некоторых расчетов на кластере ИВМ СО РАН с гетерогенной архитектурой, как правило, была незначительна. Численные эксперименты показали несомненное преимущество однородной архитектуры высокопроизводительной вычислительной системы.
Результаты серии расчетов на всех кластерных системах на грубой модельной сетке 401 х 401 точек начиная с девяти процессов показали на неблокирующих режимах обмена для обоих вариантов декомпозиции ускорение выше линейного, что можно объяснить эффектом “попадания в кэш”. Однако для кластера ИВМ СО РАН при блокирующих обменах эффект попадания в кэш также наблюдается, но начиная с 12 процессов рост ускорения прекращается. Существенных отличий ускорений для разных видов декомпозиции не наблюдается.
На рис. 3 приведены зависимости ускорения вычислений от количества используемых процессов для мелкой сетки, полученные на кластерах ИВМ СО РАН и SKIF Cyberia. Для сравнения представлен график потенциального ускорения согласно оценке (10). Чтобы не загромождать рисунок, не показаны графики ускорений, полученных для декомпозиции с теневыми гранями, поскольку они практически совпадают с представленными. В экспериментах на кластере ИВМ СО РАН общий тренд ускорения совпадает с теоретическими оценками, однако его рост носит сильно неустойчивый характер. Расчеты, проведенные на кластере SRIF Cyberia, показывают классическую картину ускорения, подтверждающую линейный характер его роста с увеличением количества процессов с эффективностью около единицы (эффективность расчета для 32 узлов « 0.85). Эксперименты на кластере НГУ подтверждают эти результаты. Таким образом, неустойчивый характер ускорения в расчетах, проведенных на кластере ИВМ СО РАН, объясняется скорее всего его неоднородной архитектурой.
-■о--
351 —■
30 - —X—
25 - --°--
20 -
15 -
10 -
5 -
П -
без перекрытий, неблокирующие обмены (SKIF Cyberia)
без перекрытий, обмены SendRecv (SKIF Cyberia)
без перекрытий, неблокирующие обмены (кластер ИВМ СО РАН)
без перекрытий, обмены SendRecv (кластер ИВМ СО РАН)
- теория
14 16 18 20
Количество процессов
Рис. 3. Зависимость ускорения вычислений (ось ординат) от количества доступных процессов; сетка 801 х 801
Отметим, что тип декомпозиции существенно не влияет на ускорение параллельной программы, а неблокирующие операции следует признать более надежными и эффективными.
Поскольку время выполнения программы складывается из времени вычислений, времени коммуникаций и дополнительного времени для операций, связанных с распределенностью данных, то были проведены серии расчетов, в которых названные составляющие были измерены отдельно. На рис. 4 приведены результаты этих экспериментов для случая декомпозиции без перекрытий и неблокирующего режима двухточечных обменов. Для того чтобы была возможность сравнить результаты расчетов с теоретическими оценками, рассмотрены зависимости безразмерных величин — отношения времени выполнения вычислений и времени выполнения обменов для р-процессов к времени выполнения всего алгоритма для одного процесса: Рса1с = Тра1с/Т и всотт = Тротт/7\. Время на дополнительные операции не зависит от количества участвующих в вычислениях процессов и на пять-шесть порядков меньше суммарного времени вычислений и коммуникаций. Время выполнения вычислений (см. рис. 4 вверху) для гомогенных архитектур совпадает с теоретическими оценками. Ввиду неоднородности узлов кластера ИВМ СО РАН время вычислений в этом случае уменьшается немонотонно.
Схема двухточечных обменов алгоритма, диктуемая декомпозицией, начиная с трех процессов неизменна: каждый процесс, за исключением первого и последнего, на каждой итерации Якоби пересылает для компонент скорости и возвышения свободной поверхности невязку в граничных точках своим левым и правым соседям (шесть операций посылок и шесть операций получения). В результате время выполнения обменов в неблокирующем режиме не зависит от количества процессов, участвующих в расчетах. Для кластера с неоднородной архитектурой наблюдается ожидаемое немонотонное поведение исследуемой величины. Однако немонотонное поведение наблюдается и в случае расчетов для гомогенных архитектур. Исследование времени выполнения обменов показало следующее: 1) время обменов минимально и не зависит от количества процессов, участвующих в обменах, если загружены все ядра на узле по одному процессу на ядро (для БК^ СуЬепа количество процессов кратно четырем, для НР-кластера НГУ — восьми); 2) при существовании в расчетах узлов, не полностью загруженных,
Рис. 4. Исследование времени счета и времени коммуникаций; неблокирующие обмены, декомпозиция без перекрытий
время обменов тем больше, чем больше количество простаивающих ядер; 3) время, затраченное на обмены, уменьшается с ростом количества задействованных процессов (т. е. е уменьшением времени вычислений).
3.2. Эксперименты для акватории Охотского моря
Тестовые расчеты для акватории Охотского моря проводились на сетках, подготовленных С.Ф. Пятаевым и И.В. Киреевым на основе открытой батиметрической базы данных ЕТ0Р02 [6, 11].
Численный эксперимент соответствует начальным данным с локальным подъемом уровня, описываемым гауссовской функцией
£(0,А,р) = Аехр(-(А - Ао)7(2С)2 - (р - Ро)7(2В)2).
В расчетах было принято А = 10, А0 = 149.1° восточной долготы, р0 = 53.1° северной широты, П = 0.005, коэффициент трения г* = 0.0026. Для скоростей предполагались нулевые начальные данные: и(0, А, р) = г>(0, А, р) = 0.
Рис. 5. Численные эксперименты в акватории Охотского моря на восьми процессах: а — общий вид сетки; б — начальное возмущение; уровень свободной поверхности через 42 мин (в), через 67 мин (г)
В построенной сетке из общего числа граничных участков около 6.4% проходят по морю. На таких участках m2 = 1, и в краевом условии (2) полагали в = 1, d(t, Л, = 0.
Результаты численного моделирования показаны на рис. 5, б — г, где представлена функция £(t, Л, ^) для некоторых моментов времени. Линиями показаны границы, получившиеся при декомпозиции области на восемь частей. Поскольку сетка сильно сгущается в приграничных областях, ширина полосы, которая отводится одному процессу, в прибрежной зоне значительно уже, чем в море. Следует отметить, что декомпозиция области с теневыми гранями в этом случае является трудоемким процессом, поэтому не проводилась.
Заключение
На примере численного решения краевой задачи для уравнений мелкой воды в работе рассмотрены технологические аспекты разработки масштабируемых параллельных алгоритмов для кластерных вычислительных систем с использованием библиотеки MPI.
Поскольку при дискретизации задачи использовался метод конечных элементов с организацией вычислений по треугольным элементам, были рассмотрены два естественных подхода к декомпозиции области — без перекрытий и с теневыми гранями.
Теоретические оценки показали, что алгоритм обладает значительным объемом потенциального параллелизма и хорошей с точки зрения распараллеливания структурой, что в зависимости от количества используемых процессов дает ускорение, близкое к линейному.
Численные эксперименты показали, что использование неблокирующего режима обменов, которое допускается алгоритмом, является безусловно более эффективным. Как явное преимущество следует отметить простоту организации параллельного алгоритма при декомпозиции без перекрытия подобластей на неструктурированных триангуляциях реальных акваторий.
Поскольку расчеты были проведены на трех высокопроизводительных ВС двух различных архитектур, то показано преимущество однородного устройства кластера (SRIF Cyberia) над гетерогенным. Кроме того, продемонстрирована неустойчивость ускорения при неоднородном составе кластера, которая не присуща ни алгоритму, ни реализации.
Авторы благодарят коллективы МВЦ ТГУ и ИВЦ НГУ, а также проф. А.В. Стар-ченко и Д.Л. Чубарова за предоставленную возможность проведения серии вычислительных экспериментов на кластере SRIF Cyberia и кластере ИВЦ НГУ.
Список литературы
[1] Марчук Г.И., Каган Б.А. Динамика океанских приливов. Л.: Гидрометеоиздат, 1983.
[2] ÂGOSHKOV V.I. Inverse problems of the mathematical theory of tides: boundary-function problem // Rus. J. Numer. Anal. Math. Modelling. 2005. Vol. 20, N 1. P. 1-18.
[3] Kamenshchikov L.P., Karepova E.D., Shaidürov V.V. Simulation of surface waves in basins by the finite element method // Ibid. 2006. Vol. 21, N 4. P. 305-320.
[4] Kamenshchikov L.P., Karepova E.D., Shaidurov V.V. Numerical solution of the boundary problem foe shallow water equations for modelling surface waves in World ocean by finite elements methods // Finite Difference Methods: Theory and Appl. Proc. of Fourth Intern. Conf. FDM:T&A’06. Bulgaria: Rousse, 2007. P. 227-233.
[5] Каменщиков Л.П., Карепова Е.Д., Шайдуров В.В. Моделирование поверхностных волн в водоемах методом конечных элементов на вычислительном кластере // Избр. материалы Четвертой школы-семинара “Распределенные и кластерные вычисления”. Красноярск: ИВМ СО РАН, 2005. С. 114-125.
[6] Каменщиков Л.П., Карепова Е.Д., Пятаев С.Ф., Шайдуров В.В. Моделирование гравитационных волн в Мировом океане методом конечных элементов с распараллеливанием // Избр. материалы Шестой школы-семинара “Распределенные и кластерные вычисления”. Красноярск: ИВМ СО РАН, 2006. С. 52-64.
[7] Карепова Е.Д., Шайдуров В.В., Вдовенко М.С. Параллельные реализации метода конечных элементов для краевой задачи для уравнений мелкой воды // Вестник ЮУрГУ. Серия “Математическое моделирование и программирование”. 2009. Т. 17(150), вып. 3. С. 73-85.
[8] Ортега Дж. Введение в параллельные и векторные методы решения линейных систем: Пер. с англ. М.: Мир, 1991.
[9] BALAY S. Efficient management of parallelism in object-oriented numerical software libraries, modern software tools in scientific computing / Eds. S. Balay, W.D. Gropp, L.C. McInnes and others // Modern Software Tools for Scientific Computing. Birkhauser Press, 1997. P. 163-202.
[10] McBRYAN O.A. An overwiev of message passing environments // Parallel Computing. 1994. Vol. 20. P. 417-441.
[11] National Geophysical Data Center. http://www.ngdc.noaa.gov/ngdc.html
[12] Исаев С.В., Малышев А.В., Шлйдуров В.В. Развитие Красноярского центра параллельных вычислений // Вычисл. технологии. 2006. Т. 11, спецвыпуск: Избранные докл. X Российской конф. “Распределенные информационно-вычислительные ресурсы”. С. 27-33.
Поступила в редакцию 27 октября 2009 г.
