Параллельные алгоритмы численного моделирования динамики турбулентного следа в линейной стратифицированной жидкости Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 16, № 4, 2011

Параллельные алгоритмы численного моделирования динамики турбулентного следа в линейной стратифицированной жидкости*

Н.П. Мошкин

Институт вычислительных технологий СО РАН, Новосибирск, Россия Institute of Science, Suranaree University of Technology, Nakhon Ratchasima, Thailand e-mail: nikolay.moshkin@gmail. com

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

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

Введение

Исследование турбулентных следов и еледоподобных образований в стратифицированных средах имеет важное значение для метеорологии (течения над островами, горами и т.д.). Следоподобные образования играют важную роль в формировании тонкой микроструктуры гидрофизических полей в океане [1]. Турбулентные следы за буксируемыми и самодвижущимися телами в стратифицированных жидкостях рассмотрены в целом ряде теоретических, экспериментальных и численных исследований [2-13]. Анализ математических моделей и их численных реализаций показывает, что использование современных полуэмпирических подходов для моделирования динамики турбулентного следа в стратифицированной среде требует существенных затрат времени вычислений на однопроцессорной вычислительной системе. Цель настоящей работы — разработка и анализ параллельных алгоритмов, уменьшающих время вычислений с использованием модели турбулентности второго порядка. Рассмотрены два подхода к распараллеливанию задачи. В одном из них учитывается возможность независимого вычисления отдельных компонент модели (функциональный параллелизм). Во втором случае вычисления проводятся параллельно в различных подобластях (геометрический параллелизм). Теоретические оценки ускорения (speed up) подтверждаются численными экспериментами. Представленный материал является продолжением предыдущих исследований автора [14].

* Работа выполнена при финансовой поддержке РФФИ (грант № 10-01-00435а) и СО РАН (Интеграционный проект 26).

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

Данная работа посвящена реализации примененных ранее численных моделей турбулентного следа в стратифицированной среде для широкого класса вычислительных многопроцессорных систем, условно называемых кластерами. Чтобы обеспечить переносимость программ на ЭВМ с различной вычислительной архитектурой, используется коммуникационная библиотека MPI (Message Passing Interface), Вычисления проводились на двух кластерах со следующими параметрами: кластер 1 управлялся операционной системой Solaris и состоял из PC-сервера с шестью процессорами (один двухъ-ядерный CPU, Opteron (2 cores), 1,6 ГГц, 8 GB RAM, один четырехъядерный CPU, Intel (4 cores), 2,0 ГГц, 8 GB RAM) и двух машин на базе двухъядерного CPU, AMD (2 cores), 2.4 ГГц, 2 GB of RAM.

Кластер 2 под управлением операционной системы Linux состоял из семи PC с двумя четырехъядерными CPU, Quad cores Xeon, 2.33 ГГц, 8 GB of RAM и семи PC с двумя двухъядерными CPU Dual cores Xeon, 3.0 ГГц, 4 GB of RAM. На обоих кластерах использовался компилятор Sun FORTRAN на базе Sun Studio 11с библиотекой MPI 2.0. Латентность и пропускная способность канала передачи данных оценивались экспериментально с помощью специально написанной программы. Для кластера 1 латентность составила (5.2)-5 с, пропускная способность капала 118 МБ/с, для кластера 2 — соответственно (5.9)-5 с и 260 МБ/с.

В качестве исходной рассмотрим модель, примененную в работах [4, 13]. Для описания течения в дальнем турбулентном следе за осеснмметрпчным телом в стратифицированной среде используется трехмерная параболизованная система осредненных уравнений Навье — Стокса в приближении Обербека — Буссипеска:

TT 9Ud , vdUd , u/dUd 9 / ' /\ , 9 / ' /\ m

U°-te+V~dï + W-dr = dï{uv) + Tz{uw)> (1)

и — v— W— - 1 9^ A (v'2\ A (v'w'\ (2)

0 dx dy dz po dy dy^ ' dz

TTdW dW dW 1 d (pi) d d , m (pi) /0.

UQir + V— + W— =---f-L - — (v'w') - — (w'2) - gZ-L, 3)

dx dy dz p0 dz dy dz p0

TT ^(Pi) , T^d(pi) , , iXrdps d , , d , ,

U°— + V^ + W— + W^ = -d-y{vp)-d-z{wp)> (4)

dV dW dUd

d

ду дг дх ^

где и4 = и0 — и — дефект осредненной продольной компоненты скорости; и,У,Ш — компоненты скорости осредненного движения в направлении осей х,у,г соответственно; (р1) — отклонение давления от гидростатического, обусловленного стратификацией р«(г); ио — скорость набегающего невозмущенного потока; д — ускорение силы тяжести; (Р1) — осредненный дефект плотности: р1 = р — р3, р3 = р3(г) — плотность невозмущенной жидкости, которая полагается линейной: р3(г) = р0(1 — аг), а = eonst > 0;

()

жидкости считается линейной функцией температуры, стратификация предполагается слабой. В уравнениях (1) - (4) в предположении малости отброшены члены с молекуляр-

х

Система уравнений (1)-(5) незамкнута. Ниже рассмотрена математическая модель, которая вместе с (1)-(5) образует замкнутую модель течения.

Величины компонент тензора рейнольдеовых напряжений (и^Ц) (кроме (и'2и3) = (у'ш'), (и']2) = (и'2), (и'22) = (у'2), (и'32) = (ш'2)), турбулентных потоков (и[р') и дисперсии флуктуаций плотности (р'2) аппроксимируются алгебраическими соотношениями [15, 16] (по повторяющимся индексам предполагается суммирование):

Ки;)

-(и'гр')

5и +

1~С2 ГРг С1 Vе

3

е

С'

2 „ С

е

С'т е

/ / /\5(р) . (л ч Л / 9г , /2\

2 е. ' д(р)

<р') = —Ли'кр')

Ст е

дхк

. ' ' ди^ ' ' дЩ

РгЗ = ~ 1 +

дхь

(6)

(7)

(8) (9)

1

= — «ЧР + (и^ }дг), г,], к =1,2, 3, ро ^

g =(0, 0, -д), 2Р = Ргг, 2С = Сгг.

По аналогии с [4, 5] упростим выражения (6)-(9) с учетом физических особенностей рассматриваемого течения — спутного струйного турбулентного течения в поле силы тяжести на больших расстояниях от тела; при этом интересующие нас слагаемые порождения (9) заменяются приближенными соотношениями

р о Л ' >\д17(1 , / ' >\д17(1 Рп = 2^{иь)— + {ит) —

Р22 = Р3

33

0, Р

12

(у'2)

дЩ

ду

Р1

13

Выражения (6)-(8) упрощаются следующим образом:

(Ш2)

{иь} =---— = К„

С1

дУ

дУ

дЩ

дг

(10)

(и'ш ')

(1 - С2)е(ш'2) -

(1 - Сз)(1 - С2Т) д е2

С1т

ро е

(ш'р')

С1е 1 -

(1 ~ сз) 9 е2 д{р) \ с\С\т ро£2 дг )

/ /2ч 2 е/ / !\д{р)

дил

дг

К.

дЩ

дг

/ / л 1 е

"ИР = —т

С1т е

-(у'р')

, ' Л д(р) ... ,.ди

{ит}-^- + (1 - с2Т){и) р ) —

1 ел./2\д(р)

С1Г е <9у

К

ад

(П)

(12)

(13)

(14)

е

е

е

-{w'p')

CiT £

dz po

e{w'2) d(p) d(p)

с1т£ ( 1 _ о1 ~ °2Т 9 е2д(рА дх дх

1Т С1ТСТ р0£2 дх )

Для определения значений касательного рейнольдеова напряжения (у'ы') привлекается дифференциальное уравнение переноса [15, 16]

d(v' w') т d(v'w') rrrd(v'w') d dlv'w')

» я + V \ + W \ = —Кеу \ dx dy dz dy dy

+ (1 _ C2)p23 + (1 _ Сз)С28 _ Cl£(„v>, (16)

dW t ._dV\ ^ g

Коэффициенты турбулентной вязкости в уравнении (16) Key = Ky, Kez = Kz, величины Ci = 2.2, C2 = 0.55, C3 = 0.55, cit = 3.2, cei = 1.44, C£2 = 1.92, ct = 1.25, C2T = 0.5, a = 1.3 являются общепринятыми эмпирическими конетантами. Величины [u'2) (i = 1, 2, 3)

£

соответствующих дифференциальных уравнений переноса с упрощенными представлениями коэффициентов турбулентной вязкости:

d(u'2) d(u'2) d(u'2) d ~ d(u'2) d ~ d(u'2)

Uo^.--b Vх—--b w -- = —key—--b --b n 1 + Сгц —

dx dy dz dy dy dz dz

-\s - ^ Uu>2) - \e) - c2 (pn - ¡p) - c2 (gll - ¡g) , (17)

d(v>2} d(v>2} d(v>2} d 2 d(v>2} d 2 d(v>2}

)-^--^ V --^ И/—- = —Key—--b ^-iiez^,-

dx dy dz dy dy dz dz

~ cA (V2> - \e\ - c2 (p22 - IA - c2 (g22 - |g ) , (18)

d(w'2) d(w'2) d(w'2) d ~ d(w'2) d ~ d(w'2)

Uo—--b V—--b W—-- = -¿-Key—---b —---b P33 + G33-

dx dy dz dy dy dz dz

-Cl£- Uw'2) - |e) - c2 (P33 - - C2 (G33 ~ ^g) , (19)

Здесь

TT d£ T d£ TTd£ d^d£ d^d£ £.n £2

+ + W7T = 7ГК^7Г + + Cei~(P + G - ce2 —,

dx dy dz dy dy dz dz e e

e = «u'2) + [v'2) + [w'2))/2.

ee Key = cs-{v12), Kez = cs-(w'2), Key = Key/(J, Kez = Kez/(7,

(20) (21)

Р22 = о, с22 = О, Рзз = 0, С33 = -2^-(^У) = 2С,

ро

р / I !\ди<1 пдиа (диЛ2 (диЛ2 д , , д др

В результате выполненных построений модель дальнего турбулентного следа представляет собой систему дифференциальных уравнений (1)-(5), (16)—(20).

В дополнение к уравнениям модели используются уравнения переноса осредненной концентрации пассивного скаляра в и дисперсии флуктуаций концентрации пассивного скаляра (в'2) [10, 13] :

дв дв дв д дв д дв и01г + у-гг- + iv— = т-к&у— + 7— , 22

дх ду дг ду ду дг дг

тт 4. 4. д к д^ к д^ и0^;--Ь к^г--Ь УУ —— = —Ккэу—--Ь тг^-Юг

I I 1 I I л/ I 1 Л I 1 Л I л/

о I о I ' ' о о 1~у о 1 о о

дх ду дг ду ду дг дг

дв дв -2т-(23)

где

^ _Юе „ _К>е _ с,{у'2)е С>*)е

—-) —-, —-, —-,

С1те С1те е е

дв дв в 2 е

ду дг е

Величины С1Т = 3.2, СТ = 1.25, С^ = 0.13 являются известными эмпирическими константами [15, 16],

Применимость рассматриваемой и других подобных моделей к расчету безымпульсных следов и следов за буксируемыми телами детально анализировалась в [4, 5, 10],

х

(16)-(20), (22), (23) играет роль времени. На расстоянии х = х0 от тела задаются следующие начальные условия:

е(хо,у,г) = $1(0, е(xо, y, г) = ф2(r), Ща(хо,у,г) = фз(r), в(хо,у,г) = ф4(r),

(в'2)(х0,у,г) = Ф5(г), (у'ш') = (р1) = V = Ш = 0, г2 = у2 + г2, 0 <г< то.

Здесь Ф1(г), Ф2(г), Ф3(г) — функции, согласующиеся с экспериментальными данными об эволюции турбулентного следа в однородной жидкости; Ф4(г), Ф5(г) — некоторые финитные колоколообразные функции [3, 4, 10], Дополнительно задавались начальные

/ 2

условия для нормальных рейнольдеовых напряжений (и2) = —е (г = 1, 2, 3). При г2 = у2 + г2 ^ то ставились условия невозмущенного потока

и = V = Ш = (р1) = е = е = (у'ш') = в = (в'2) = 0, х > хо.

При численном решении задачи нулевые краевые условия, соответствующие г — то, переносятся на границы достаточно большого прямоугольника. Из соображений симметрии решение отыскивается лишь в первом квадранте плоскости (у, г), На осях симметрии принимаются граничные условия

д(рг) ду

V

дЦ де де

= <р1) = W

дУ

ду_

дг

ду ду ду

дЦ*

дг

де дг

де дг

50 _ д{9'2) ду ду

д<0/2)

0,

у

0,

г > 0,

дв

дг

дг

0, у > 0.

Дополнительно задавались условия симметрии па осях координат для нормальных рей-нольдеовых напряжений (пг 2) , г =1, 2, 3,

Переменные задачи могут быть обезразмерены с применением масштаба длины Д — диаметра тела и масштаба Ц _ скорости невозмущенного потока. Безразмерные переменные вводятся следующим образом:

г

ъ'

Ш

х у

!>' у =ъ,

/ /3\

/3\* {ш )

/

иг

иг С/о' ы* <Р1) / ' '\* = ТТ2П > те) = ио Ро {Ки'з Л2

еБ <р)* = аДро \Роу \ (1ра 1

©*

©

$4(0);

<Г)* =

{о'2) ф42(0):

т

Аор

©оо^о'

При этом в обезразмереппых уравнениях вместо $ появится вели чина 4п2/Г,, где плотностное число Фруда, определяемое равенством

и0т в '

т

2тт

л/од'

Для удобства интерпретации результатов расчетов ниже используется время ¿, связанное с расстоянием от тела соотношением

i = — i*

и0'

т

г

т

хД

х

ЦоДТ Г,'

*

*

г

*

е

2. Реализация алгоритма решения задачи на многопроцессорной системе

Несколько численных моделей динамики турбулентного следа в стратифицированной жидкости успешно использовались для численного моделирования [4, 5], При этом в [4] модель основана на введении переменных завихренность — функция тока, в работе [5] применялся явный метод расщепления по физическим процессам для нахождения компонент вектора скорости V, W,

Используемая в настоящей работе модель состоит из 12 дифференциальных уравнений и одного алгебраического уравнения для энергии турбулентности.

2.1. Последовательный алгоритм

Последовательный алгоритм решения задачи основан на неявном методе расщепления по пространственным переменным [17] для уравнений (1), (4), (16)-(20), (22), (23) и на явном методе расщепления по физическим процессам [18] для системы уравнений (2), (3), и (5), Предположим, что все неизвестные функции известны при х = хп. Процедура решения задачи состоит из шагов, позволяющих определить все неизвестные функции в следующем сечении х = хп+1. Детальное описание последовательного алгоритма приведено в [14],

Алгоритм решения задачи сводился к последовательному интегрированию системы уравнений (1) — (5), (16)-(23), записанных в новой системе координат. Условимся, что верхний индекс п соответствует сечению х = хп = хп-1 + Л-П, а нижние индексы %,] отвечают узлу сетки тц) в сечении х = ^^^ (соответственно (г + 1/2, j + 1/2) — узлу сдвинутой сетки),

хп

сечении х = хп+1 находится следующим образом:

1) из уравнения (1) отыскиваются значения дефекта осевой компоненты скорости 1- При этом используется неявная схема расщепления по пространственным переменным [17];

2) для определения компонент вектора скорости V, Ш и дефекта давления (р1) используются уравнения (2), (3) и уравнение несжимаемости (5), которые интегрируются с применением явного метода расщепления по физическим процессам [18] (неявный алгоритм использовался в [19]);

3) с применением неявного метода расщепления по пространственным переменным из уравнений (4), (16)-(20), (22), (23) последовательно определяются сеточные функ-

ч™ (Р1)п+1> з (уы )п+1, (и%+1> (у%+1> (ы%+1> з (©)пг> (нп; ^

4) проверка необходимости следующего шага, В случае необходимости перейти к выполнению пункта 1,

Вычисленные на очередном этапе алгоритма переменные участвуют в определении неизвестных на последующих этапах. Реализованная таким образом идея блочного аналога метода Зейделя существенно упрощает процедуру расчетов,

2.2. Функциональная модель параллельного алгоритма

В основе принципа распараллеливания использована возможность решать некоторые уравнения независимо друг от друга. Например, если не применять блочный аналог метода Зейделя, то уравнения для компонент тензора рейнольдеовых напряжений можно решать независимо. Так как уравнения имеют одинаковую структуру, то проблем с балансировкой загрузки процессоров не возникает.

Поясним идею на примере абстрактного алгоритма. Предположим, что имеется М процессов, которые выполняются последовательно, как показано на рис, 1, о.

Обозначим ^ время, необходимое для выполнения процесса 7 для г = 1, 2,..., М.

м

Время последовательного выполнения М-процеесов Т3 = ^ Предположим также,

г=1

что последовательный процесс может быть реализован в виде последовательного выполнения к групп несвязных процессов (рис, 1,6), На рис, 1, б использованы обозначе-к

ния: тз = М, ^ — время выполнения процесса Ьв? — время передачи данных 3 = 1

а

б

Рис. 1. Схема абетрактншх) поеледовательншх) а.;п'оритма (а) и блок-схема абетрактншх) па-раллельншх) а;пх)ритма (б)

процесса, г = 1, 2,..., т3-, j = 1, 2,..., к, Пусть вычислительная система имеет Р процессоров, число которых достаточно дня независимого выполнения заданий в каждой

группе Р > тах{т3-}3=1)2)...,к. Можно оценить время параллельного выполнения прок

цессов Тр = ^ (Т3 + Т/птгп), где Т3 и Т/отт — соответственно время выполнения и

3=1

передачи данных, необходимое для завершения всех работ в группе ^ j = 1, 2,... ,к. Очевидно, что

т

тах 3 и Т/отт « V (р - 1) .

=1

Здесь [•] означает функцию взятия целой части. Время передачи данных процессо-

. 3 МеввадеБгге/

ром можно представить как = Н--;-——. где — латептпость — дли-

3 БапаШгат

тельность подготовки сообщения для передачи, МеввадеБгге/ — размер сообщения, передаваемого процессом 7/, ^ = 1, 2,... ,к, г = 1, 2,... ,т/, БапвШгаИН — пропускная способность канала передачи данных. Поэтому ускорение параллельного алгоритма, основанного па функциональной независимости отдельных заданий, можно оцепить как

Т3

Р

м

5

=1

к к

3=1

щ Р

т

тах / + > (Р - 1) и +

=1

МеввадеБгге/

Бапашгагн

(25)

к

цессов т (т3- = т, j = 1,..., к). Каждый процесс требует Ь секунд для выполнения

задания (¿г = ¿, г = 1,..., М) и передает одинаковый размер сообщения МеззадеБгге7 МеззадеБгге, Тогда ускорение

5

т х £

тп ч i МеззаоеБгге\

Р

Если время коммуникаций между процессорами нулевое, то ускорение равно Р и соответствует идеальному. Ускорение зависит от латентности, пропускной способности канала передачи данных, размера передаваемого сообщения и времени выполнения отдельного процесса. Можно установить соотношения между параметрами, при которых ускорение теоретически возможно, Б > 1:

5

к

Е

.7 = 1

М г=1

щ

Р

шах1<г<т. + Е(Р - 1) ¿г +

г=1

МеээадеЗгге^

<

<

к

Е # т7

г=1

к т

к т

Е & + Е ЕОР - 1)^ + Е ЕОР - 1)

7=1

7=1г=1

7=1г=1

МеззадеБгге7

(26)

Здесь ¿7 — максимальное время выполнения работ и МеззадеБгге7 — минимальный размер передаваемых данных в группе Если правая часть неравенства меньше 1, то, очевидно, ускорения не произойдет. Следовательно, необходимо рассматривать только случай, когда правая часть (26) больше 1, Поэтому имеем

к к

7=1

щ Р

к т

к т

7

7=1 г=1 7 = 1 г=1 4 7 г=1

ИЛИ

Е [^1 * - Е + Е - + - и* < о.

7=1

7=1

7=1

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

к

Случай 1. Е 1

7=1

Очевидно, что

-щ Р

■гщ р

V — Vт7\ < 0.

< т. для всех ] = 1,..., к, тривиальный случай.

Если ¿7 > больше 1,

Меззадевгге^ 1

(Р — 1)МеззадеБгге7

Случай 2, ^ < — ¿7т7- + т7-(Р — 1)

<0

для всех ^ = 1,..., к, то можно получить ускорение

Случай 3. Е #т7 > М(Р — 1)^.

7=1

Если все ¿7 = {, то Ь > (Р — 1)£г и можно получить ускорение больше 1,

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

1, Инициализация для создания групп коммуникаций в MPI,

2, Задание начальных данных,

3, Распределение работ между процессорами,

4, Одновременное вычисление Ud, V* и W* (m1 = 3),

5, Обновление данных для Ud, V* и W* (передача данных перед выполнением следующего шага),

6, Вычисление (pi)n+1 (m2 = 1). На этом этапе возможно использование всех свободных процессоров,

7, Вычисление Vn+1 и Wn+1 на двух процессорах одновременно (m2 = 2),

8, Обновление данных для Vn+1 и Wn+1 (передача данных перед выполнением следующего шага),

9, Вычисление (p1)n+S £n+1, (v'w')n+1, (u/2)n+1, (v/2)n+1, (w/2)n+1, en+1, (0)n+1 и (0/2)n+1 (m4 = 8), Число независимых процессов при использовании других моделей может варьироваться,

10, Обновление данных для (p1)n+S £n+1, (v' w')n+1, (u/2)n+1, (v/2)n+1; (0)n+1 и (w/2)n+1 (передача данных перед выполнением следующего шага),

11, Проверка необходимости выполнения следующего шага, В случае необходимости переход к выполнению пункта 4,

2.3. Геометрическое распараллеливание

Второй алгоритм распараллеливания основан на принципе геометрического параллелизма, Программа разбивает расчетную область на непересекающиеся подобласти, содержащие (возможно, приблизительно) одинаковое количество точек. Вычисления в каждой подобласти выполняются на отдельном процессоре.

Рассматриваемая модель турбулентности требует решения нескольких однотипных параболических дифференциальных уравнений. Эти уравнения решаются с применением метода расщепления по пространственным переменным [17], Приведем пример типичного уравнения (это может быть уравнение переноса компоненты тензора рей-

£

ции):

dF dVF dWF d / dF\ d / dF\ птто , . 0 . _

F(y,z, 0) = U0(y, z, 0), (y,z) G Q, F (y, z, t) = g(y, z, t), (y, z) G dQ x [0, T]. (28)

Для простоты рассмотрим однородную прямоугольную сетку в области Q

Qh = {(yi, zj )|yi = iAy, zj = j Az, i = 0,...,M + 1,

j = 0,..., N +1, Ay = 1/M, Az = 1/N}.

Схема расщепления записывается в виде двух полушагов во времени. На первом полушаге используются следующие разностные уравнения:

Аг 2 Ау

(Ку)» 1/2„-(*£#а - - - КУЛ

Ау2

На втором полушаге рассматриваются уравнения

Ц?1 - (Ш^ - ^Р)^ _

+ яяя

г7

(29)

2Д г

(К, )г,7+1/2 (7 — Р^1) — (К, )г,7-1/2 — 7)

(30)

Выражение (29) представляет собой систему линейных алгебраических уравнений с трех-диагональной матрицей

<1/2 = д (0,*,.),

_а..рп+1/2 + г Р _ ь.г = ь . г =-1 2 м 7 = 1 2 N

аг7 ^г—1,7 + Сг7 Рг,7 Ьг7 Рг+1,7 = Ьг,7, г = 1, 2,-.,М, 7 = 1, 2,...,^,

тП+1/2

тП+1/2 _

^М++У7 = д (ум ,г7).

(31)

Для каждого фиксированного 7 имеется система из М линейных уравнений с трехдиа-гональной матрицей. Поэтому легко организовать параллельный процесс для решения (31), приписывая каждому процессору ^Р систем. Процессор К (более точно (К — 1))

(К — ^ будет решать системы уравнений для ] = -—--Ь 1,..., . При такой стратегии параллелизма матрица коэффициентов разбивается по строкам, как показано на рис. 2, а. На втором дробном шаге необходимо решить систему

РП+1 = д (уг, 0),

-п. . гп+1 + с • Р^ - ь • Рп+1 = Ь • 7 = 12

аг,7Рг,7 — 1 + сг,7Рг,7 ьг,7Рг,7+1 = Ьг,7 , 7 = 1, 2, . .

., N г = 1, 2,..., М,

РП++1 = д (уг )■

(32)

СРЩМ)

СРи(1) СРи(О)

• • •

СРИ(О) СРи(1) СРи(Р-1)

Рис. 2. Схема разбиения матрицы по строкам вдоль оси У (а) и по столбцам вдоль оси X (б)

Для каждого фиксированного г имеем систему из N линейных уравнений с трехдиаго-нальной матрицей. Процессор (К — 1) будет решать системы уравнений, соответствую-

(К — 1)М КМ ^

щие г = -—--Ь 1,..., , Рисунок 2, б показывает направления вычислений по

столбцам на втором дробном шаге схемы расщепления.

Алгоритм, основанный на геометрическом параллелизме, подобен последовательному алгоритму. Однако неизвестные (Ц^)п+1, (р1)га+1, ёп+\ (ь'т')п+\ (и'2)п+1, (^'2)га+1, (и>'2)п+1 и др. вычисляются с применением геометрического подхода. Остальные неизвестные функции находятся параллельно на различных процессорах.

Рассмотрим сетку, состоящую из М х N узлов, и предположим что вычисления в каждом узле требуют £ секунд, а во всех узлах (М х N) х £ секунд. Пусть имеется Р процессоров и каждый из них работает с (М х N) /Р узлами сетки. Каждому процессору требуется (М х N) х Ь/Р секунд и необходимо Р х (Р — 1) раз обновить данные блоков, В этом случае можно ожидать ускорения

5 = _„_(33)

N х М х Ь „ / МеввадеБъгеЛ к '

+ Рх(Р- 1)х Н у *

P к ; V 1 Bandwidth J

Оценим ситуации, в которых ускорение теоретически возможно. Пусть T — время выполнения последовательной программы T = N х M х t. Ускорение возможно в случае

T

T ^ , . ( MessaqeSize\ р + рх(р-1)х(*'+ BandLdth)

> 1, (34)

или

'^Ннйвг)- »

л г т г,2 ( Ме88адеБ%хе\

Ускорение невозможно в случае 1 < -г ¿Н----г-т- . Другими словами, ис-

\ BandWгdth )

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

В уравнении (35) МеззадеБ1ге = -^242—— Мбайт и в случае ¿г = 0 имеем

2 2 ( N х М х 8 \

" х м х '> Р и + Р {Ш^Р^в^Шт) '

или

8

* >

10242 х BandWidth'

Следовательно, ускорение может быть больше единицы, если время выполнения зада-

8

ння на одну точку сетки / >---,

10242 х BandWidth

3. Результаты расчетов

Параллельные алгоритмы детально тестировались путем сопоставления с данными численных экспериментов [8, 11, 14].

Таблица 1. Безразмерные осевые значения энергии турбулентности ео = ео(х) = е(х, 0, 0) на различных расстояниях от тела; след за самодвижущимся и буксируемым телами; ^ = 280

Последовательный Алгоритм Алгоритм

х/Б алгоритм с функциональным параллелизмом с геометрическим параллелизмом

След за самодвижущимся телом

12.00 0.42813Е-01 0.42890Е-01 0.42813Е-01

19.00 0.33442Е-01 0.33436Е-01 0.33442Е-01

63.65 0.13778Е-01 0.13744Е-01 0.13778Е-01

119.22 0.80534Е-02 0.80220Е-02 0.80534Е-02

252.74 0.41380Е-02 0.41117Е-02 0.41380Е-02

947.74 0.13891Е-02 0.13771Е-02 0.13891Е-02

1502.74 0.96091Е-03 0.95273Е-03 0.96091Е-03

След за буксируемым телом

12.00 0.46145Е-01 0.46144Е-01 0.46145Е-01

19.00 0.33749Е-01 0.33737Е-01 0.33749Е-01

63.35 0.21323Е-01 0.21258Е-01 0.21323Е-01

118.78 0.14147Е-01 0.14093Е-01 0.14147Е-01

252.24 0.79067Е-02 0.78716Е-02 0.79067Е-02

947.24 0.29912Е-02 0.29963Е-02 0.29912Е-02

1502.24 0.21874Е-02 0.21867Е-02 0.21874Е-02

Таблица 2. Безразмерные осевые значения дефекта осевой скорости и а на различных расстояниях от тела; след за самодвижущимся и буксируемым телами; ^ = 280

Последовательный Алгоритм Алгоритм

х/Б алгоритм с функциональным параллелизмом с геометрическим параллелизмом

След за самодвижущимся телом

12.00 0.84457Е-01 0.84848Е-01 0.84457Е-01

19.00 0.50144Е-01 0.50165Е-01 0.50144Е-01

63.65 0.11243Е-01 0.11336Е-01 0.11243Е-01

119.22 0.54013Е-02 0.54826Е-02 0.54013Е-02

252.74 0.27806Е-02 0.28392Е-02 0.27806Е-02

947.74 0.14482Е-02 0.14851Е-02 0.14482Е-02

1502.74 0.12296Е-02 0.12611Е-02 0.12296Е-02

След за буксируемым телом

12.00 0.21287Е-00 0.21288Е-00 0.21287Е-00

19.00 0.17274Е-00 0.17276Е-00 0.17274Е-00

63.35 0.73696Е-01 0.73956Е-01 0.73696Е-01

118.78 0.44832Е-01 0.45083Е-01 0.44832Е-01

252.24 0.27355Е-01 0.27544Е-01 0.27355Е-01

947.24 0.16515Е-01 0.16626Е-01 0.16515Е-01

1502.24 0.14563Е-01 0.14622Е-01 0.14563Е-01

В табл. 1, 2 сравниваются результаты вычислительных экспериментов по трем алгоритмам — последовательному, с функциональным и геометрическим параллелизмом (для плотноетного числа Фруда Fd = 280), В табл. 1 приведены безразмерные осевые значения энергии турбулентности e0 = e0(x) = e(x, 0, 0) на различных расстояниях от тела в случаях безымпульсного следа и следа за буксируемым телом, в табл. 2 — безразмерные осевые значения дефекта осевой скорости Ud на различных расстояниях от тела за самодвижущимся и буксируемым телами.

Некоторые несоответствия, видимые из таблиц, можно объяснить различающимися численными реализациями. Например, последовательный алгоритм использует идею блочного аналога метода Зейделя, Приведенные сравнения указывают на достаточную надежность используемых алгоритмов,

3.1. Экспериментальная оценка эффективности параллельных алгоритмов

Время выполнения последовательного алгоритма измерялось временем реализации программы на одном процессоре после этапа генерации начальных данных. Для оценки времени выполнения использовалось "Wall-clock time", содержащее "idle time" и время обмена данными между процессорами.

Зависимость ускорения от числа процессоров и используемых алгоритмов для кластеров 1, 2 приведена в табл. 3, 4, Максимальное число процессоров равно восьми. Оценка ускорения параллельных алгоритмов проведена на трех сетках размером 400 х 750, 600 х 1150 и 800 х 1550 узлов то направлениям z и у, Видно, что в случае алгоритма с использованием функционального параллелизма ускорение возрастает линейно с увеличением числа процессоров до четырех, а далее существенно замедляется до уровня насыщения с величиной ускорения чуть больше двух, В случае геометрического параллелизма ускорение с увеличением числа используемых процессоров до восьми возрастает почти линейно. Видно также, что при увеличении числа расчетных узлов ускорение возрастает несущественно.

На рис, 3 приведено сравнение ускорений параллельных алгоритмов, полученных в результате численных экспериментов на кластере 2, с теоретическими оценками (25) (о)

Таблица 3. Зависимость ускорения от числа процессоров и алгоритмов с функциональным и геометрическим параллелизмом для кластера 1

Размер сетки

NCPU 400 х 750 600 х 1150 800 х 1550

Время Ускорение Время Ускорение Время Ускорение

Алгоритм с функциональным параллелизмом

1 6824.52 _ 14058.94 _ 27420.92 _

2 5593.86 1.22 11247.15 1.25 19727.28 1.39

4 3554.43 1.92 6994.49 2.01 13183.13 2.08

7 3174.19 2.15 6361.51 2.21 11718.34 2.34

Алгоритм с геометрическим параллелизмом

1 6824.52 _ 14058.94 _ 27420.92 _

2 4382.75 1.56 8824.99 1.59 16876.53 1.62

4 2293.28 2.98 4608.72 3.05 8624.59 3.18

7 1683.57 4.05 3426.01 4.10 6368.64 4.31

8 1564.32 4.36 3171.43 4.43 6039.14 4.54

T а б .л и ц а 4. Зависимость ускорения от числа процессоров и алгоритмов е функциональным и геометрическим параллелизмом для кластера 2

Размер сетки

NCPU 400 х 750 600 х 1150 800 х 1550

Время Ускорение Время Ускорение Время Ускорение

Алгоритм с функциональны,м параллелизмом

1 4913.07 11952.72 25420.44

2 3584.59 1.37 8417.48 1.42 17614.25 1.44

4 2516.09 1.95 5683.31 2.10 12058.83 2.11

7 2134.25 2.30 4927.48 2.43 10291.68 2.47

Алгоритм с геометрическим параллелизмом

1 4913.07 11952.72 25420.44

2 3125.54 1.57 7431.28 1.61 14869.83 1.71

4 1628.22 3.02 3850.87 3.10 7916.05 3.21

7 1186.73 4.14 2851.70 4.19 5846.71 4.35

8 1120.38 4.39 2686.21 4.45 5487.65 4.63

3.0

2.5

4>

К

я и 2.0

Р. Q

Ьй £ 1.5

1.0

0.5

0.0

"2

2 3 4 5 Число процессоров

7 6 5

и Я

S 4

I

£ з 2 1

0

1

' 2

2 3 4 5 6 Число процессоров

Рис. 3. Сравнение ускорений, полученных теоретически (1) и в результате численных экспериментов на кластере 2 (2); размер сетки 800 х 1550; а — функциональный, б— геометрический параллелизм

и (33) (б). Из графиков следует, что теоретическая оценка но сравнению с численной завышена. Это объясняется не идеально сбалансированными заданиями. Тем не менее наблюдаемое согласование результатов можно признать удовлетворительным.

Автор признателен проф. Г.Г. Черных за полезные обсуждения в ходе выполнения работы и Dr. Sa-at Muangehan за помощь при выполнении расчетов.

Список литературы

[1] МОПИП А.С., ОЗМИДОВ Р.В. Океанская турбулентность. Л.: Гидрометеоиздат, 1981. 319 с.

[21 Lix J.T., РАО Y.H. Wakes in stratified fluids /7 Annu. Rev. Fluid Meeh. 1979. Vol. 11. P. 317 336.

[3] Hassid S. Collapse of turbulent wakes in stratified media //J- Hvdronautics. 1980. Vol. 14. P. 25-32.

[4] Chernykh G.G., Voropayeva O.F. Numerical modelling of momentumless turbulent wake dynamics in a linear stratified medium // Comput. Fluids. 1999. Vol. 28. P. 281-306.

[5] Воропаева О.Ф., Мошкин Н.П., Черных Г.Г. Внутренние волны, генерируемые турбулентными следами за буксируемым и самодвижущимся телами в линейно стратифицированной среде // Матем. моделирование. 2000. Т. 12, № 1. С. 77-94.

[6] Gourlay M.J., Л нем л S.C., Fritts D.C., Werne J. Numerical modelling of initially turbulent wakes with net momentum // Phvs. Fluids. 2001. Vol. 13, No. 12. P. 3783-3802.

[7] Dommermuth D.G., Rottman J.W., Innis G.E., Novikov E.A. Numerical simulation of the wake of a towed sphere in a weakly stratified fluid //J. Fluid Mech. 2002. Vol. 473. P. 83-101.

[8] Воропаева О.Ф., Мошкин Н.П., Черных Г.Г. Внутренние волны, генерируемые турбулентными следами в устойчиво стратифицированной среде // Докл. АН. 2003. Т. 392, № 2. С. 190-194.

[91 Мег мен P., Spedding G.R. Stratified propelled wakes // J. Fluid Mech. 2006. Vol. 552. P. 229-256.

[10] Мошкин Н.П., Фомина А.В., Черных Г.Г. Численное моделирование динамики турбулентного следа за буксируемым телом в линейно стратифицированной среде // Матем. моделирование. 2007. Т. 19, № 1. С. 29-56.

[11] Chernykh G.G., Demenkov A.G., Fomina A.V. et al. Numerical modeling of some free turbulent flows // Notes on Numerical Fluid Mechanics and Multidisciplinarv Design. 2008. Vol. 101. P. 82-101.

[12] Воропаева О.Ф. Численная модель анизотропного вырождения турбулентности в дальнем безымпульсном следе в стратифицированной среде // Матем. моделирование. 2008. Т. 20, № 11).С. 23-38.

[13] Chernykh G.G., Moshkin N.P., Fomina A.V. Numerical models of turbulent wake dynamics behind towed body in linearly stratified fluid //J. Eng. Thermophvs. 2009. Vol. 18, No. 4. P. 279-305.

[14] Мошкин Н.П. Метод расщепления в задаче численного моделирования турбулентного следа за буксируемым телом в стратифицированной жидкости // Вычисл. технологии. 2009. Т. 14, № 4. С. 81-92.

[15] Роди В. Модели турбулентности окружающей среды // Методы расчета турбулентных течений: Пер. с англ. М.: Мир, 1984. С. 227-322.

[16] Rodi W. Examples of calculation method for flow and mixing in stratified fluids //J. Geophvs. Res. 1987. Vol. 92, No. C5. P. 5305-5328.

[17] Яненко H.H. Метод дробных шагов решения многомерных задач математической физики. Новосибирск: Наука, 1967.

[18] Белоцерковский О.М. Численное моделирование в механнике сплошных сред. М.: Наука, 1977.

[19] Даниленко А.Ю., Костин В.И., Толстых А.И. О неявном алгоритме расчета течений однородной и неоднородной жидкости. \!.. 1985 (Препр. ВЦ АН СССР).

Поступила в редакцию 21 мая 2010 г., с доработки — 15 августа 2010 г.