Численное решение уравнений Навье Стокса на компьютерах с параллельной архитектурой Текст научной статьи по специальности «Математика»

ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА

2012 Математика и механика № 2(18)

УДК 519.6

Д.В. Деги, А.В. Старченко ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЙ НАВЬЕ - СТОКСА НА КОМПЬЮТЕРАХ С ПАРАЛЛЕЛЬНОЙ АРХИТЕКТУРОЙ

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

Ключевые слова: уравнения Навье - Стокса, численное решение, библиотека МР1, технология СиВЛ.

Создание новых эффективных параллельных вычислительных систем является стратегическим направлением развития компьютерной техники. Это обстоятельство вызвано как ограничением максимально возможного быстродействия обычных последовательных ЭВМ, так и существованием вычислительных задач, для решения которых ресурсов существующих средств вычислительной техники не всегда оказывается достаточным. Такие современные проблемы науки и техники как моделирования климата, генная инженерия, проектирование интегральных схем, анализ загрязнения окружающей среды, создание лекарственных препаратов, требуют для анализа ЭВМ с производительностью более 1000 миллиардов операций с плавающей точкой в секунду (1 TFLOPS).

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

Основными уравнениями гидродинамики являются уравнения Навье - Стокса. С помощью уравнений Навье - Стокса описываются такие процессы, как движение волн на поверхности жидкости, движение крови по кровеносным сосудам человека и животных, течения в мантии Земли. Задача о движение жидкости в каверне также описывается уравнениями Навье - Стокса, и в ней отражаются следующие характерные черты гидродинамических процессов: конвективная нелинейность, различные соотношения между инерционными и вязкими силами, одновременное существование областей малых и больших градиентов и т.п., поэтому задача о каверне широко распространена в качестве «тестовой» при выборе эффективного численного алгоритма для моделирования гидродинамических процессов [1]. Зачастую для получения приближенного решения с высокой точностью используют сетку с большой плотностью узлов, при этом существенно возрастает время расчета и тем самым становится актуальным применение многопроцессорных ЭВМ для уменьшения времени расчетов.

Физическая постановка задачи

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

Рис. 1. Движение жидкости в каверне

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

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

Математическая постановка задачи

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

- уравнения движения в проекции на ось Ох:

дух Ух дх

дууУх 1 др

■—-— =--------—+ ц

ду р дх

(я1

д

1

д V-

дх1 ду1

- уравнения движения в проекции на ось Оу:

д^ д^ 1 др

---- + —=----—+ ц

дх ду р ду

( д V

дх2 ду

- уравнения неразрывности:

д^х дуу п

—х + —у = 0.

дх ду

Здесь V = (ух , Уу) - скорость движения жидкости, р - давление, р - плотность жидкости, р - коэффициент кинематической вязкости.

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

- из условия непротекания:

х =0 :ух = 0х = 1х:х =0;

У = О-.Уу = 0;у = Ьу -.уу = 0 ;

- из условия прилипания:

х = 0: у = 0; х = Ьх : у = 0; у = 0:ух = °у = :ух = и,

где Ьх и Ьу - соответственно длина и ширина каверны, и - скорость движения крышки.

Получение разностной схемы

Приближенное решение ищется на равномерной шахматной сетке (рис. 2) с помощью метода конечных объемов [2].

Рис. 2. Фрагмент используемой шахматной сетки (штриховка нанесена на конечный объем)

п

При получении разностной схемы используется следующее:

1. Значения сеточных неизвестных ( ух , уу) рассматриваются на соответствующих гранях конечных объемов, а давление (р у) - в их центрах.

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

3. Для устранения возникающей нелинейности в разностной схеме значение одного из множителей следующих произведений: ухух , ухуу , уууу - полагается

равным значению, полученному на предыдущей итерации.

4. Конвективные члены аппроксимируются с помощью схемы «против потока», диффузионные и производные от давления - центрально-разностной схемой.

5. Значения сеточных функций с дробными индексами (i ± 1/2),( j ± 1/2), известные с предыдущей итерации, в точках, не совпадающих с узлами сетки, рассчитываются из значений в ближайших соседних узлах с помощью линейной интерполяции.

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

[(vx)i+ij- (vx)у ]hy+[(vy\j+i- (vy)y i = °>Nx ; J = °> Ny;

разностная схема для проекции уравнения движения по оси Ох

aPi+1j (vx )i+l j = aEi+lj (vx )/'+2 j + aWi+1j (vx \j + aNi+lj (vx ) i+1 j+1 +

+aSi+1 j (vx )i+1 j-1 + di+1 j(Pij - Pi+1 j X 1 = 0, Nx -1, j = 0, Ny;

разностная схема для проекции уравнения движения по оси Оу

bPi j+l (vy )i j+l = bEi j+l(vy )i+l j+l + bWij+l(vy X-lj+l + bNij+l (vy )i j+2 +

+ bSi j+l (vy )ij + li j+l (Pij - Pi j+l X i = 0,, j = 0, -1 .

Здесь hx и hy - шаги сетки по горизонтальной и вертикальной оси соответственно, /г- j+1 = hx / р и di+1 j = / р, Nx и - размер сетки по горизонтали и вертика-

ли соответственно. Коэффициенты aPj, aE y ,aWy, aNy, aSy, bPy, bE у, bWy, bNy, bSy

- известные сеточные функции, зависящие от поля скорости, полученного на предыдущей итерации. У этих коэффициентов имеются следующие свойства [2]:

- неотрицательность всех коэффициентов;

- aPy > aE у + aW у + aNy + aSy и bPy > bE у + bW у + bNу + bSy , причем строгое

неравенство выполняется на границе.

Предполагая, что функции vx = vx (x, y), vy = vy (x, y) и p = px, y), представляющие точное решение, имеют достаточное число непрерывно дифференцируемых частных производных, и, используя формулу Тейлора, можно доказать, что погрешность аппроксимации разностной схемы для уравнения неразрывности имеет второй порядок относительно каждого координатного направления, а для уравнений движения - первый.

Алгоритм численного метода решения

В силу нелинейности уравнений Навье - Стокса становится невозможным прямое решение системы и рассматриваемый процесс приобретает итерационный характер. Для совместного решения общей системы из трех систем разностных уравнений был применен алгоритм SIMPLE [2], который можно записать в следующем виде:

1. Задаем поле давления {p};

2. Решаем системы уравнений

aPi+1 j(v* Х-+1 j = Z anb(v* )nb + di+1 j (P* - p+1 y).

* nb * * * (1)

bPi j+1(v* )i j+1 = Z bnb (v* )nb + li j+1 (p - p],+1 )

nb

и находим

(2)

где ^ апЬ (у* )пь - означает суммирование по соседним узлам сеточного шаблона;

пЬ

3. Решаем систему уравнений

>к >к >к >к I 1 \ .

[О* );+1] - О*) 9 ]Ау + [(У у\]+1 - 0^ ) д ]А* + [-р-(Рд - Р+1 д ) -

ар+17

(й'-19 - Р )]^ + [;ь^^ (Рд - Р'ц+1) - Ьр- (р'и-1- Р'д Ж = О

и определяем поправку давления {р};

4. По соотношениям

[о*)г-+1 / = О*)г-+1 / + 4+1 у(Р/ -Р\+1,Vар+1 у>

1(^ )Ц+1 = (у*)г/+1 + 1г/+1(Р/ - Р' /+1 )/Щ/+1

ij+1

вычисляем сеточные

функции {(vx)„} и |(vy )„|, удовлетворяющие разностному

уравнению неразрывности;

5. Корректируем давление по формуле р = р* + юр' (0 <ю< 1);

6. Проверяем условия завершения итерационного процесса, если они не вы-

* о

полняются, то принимаем р = р и переходим к пункту 2.

В качестве критерия завершения итерационного процесса выбрано выполнение следующих условий:

Жх

i) \\р'н \\= max Ё\(р)•• \ <є;

J=l,.,Ny^-l V >4

Nx

2) max ^ [(vxХ+ij - (vx)у]hy + [(vy)ij+i - (vy)у ]hx 1 < є.

J=1> Ny i=i

Основные вычислительные трудности данного алгоритма - это решение уравнений для компонентов скорости (І), а также для поправки давления (2) на каждой глобальной итерации.

Один из перспективных методов решения систем линейных алгебраических уравнений вида

а Ф' - Ьу Фі-1 y - Су Фі+1 у - dy Фу_1 - ву Фу+1 = ау Фу - Z апьф nb =fij > * = І> nJ = І>m;

nb

на многопроцессорных вычислительных системах является метод релаксации [З]. Метод релаксации имеет вид

Фк-+1=®-

(ь фк+1 + с фк, + d Фк+1 + е Фк , + f ^

bj І-1 j Cj І + 1 j Ij I j-1 ij І j + 1 f IJ

aij

+ (1 -ю) Фк-,

i = 1,n, j = 1,m, k = 0,..,K. (3)

Здесь к - номер итерации. Для каждой из полученных разностных схем для определения |(ух).)..} и {р у} выполняется равенство

ау ^ Ъу + Су

+ еу = Е апЪ ,

пЬ

причем

ау > Е а

пЬ

пЬ Ь/у

■ву имеет место на границе рассматриваемой

области, что является достаточным для сходимости итерационного процесса [3, с. 89].

При 0 < ю < 1 получается метод нижней релаксации. Использование метода нижней релаксации замедляет изменение зависимой переменной. С помощью данного метода решались уравнения для скоростей.

При 1 < ю < 2 имеет место метод верхней релаксации. Он увеличивает скорость сходимости итерационного процесса. Уравнение для поправки давления решалось данным методом.

Из рекуррентной формулы (3) видно, что производится строго последовательный пересчет значений неизвестной сеточной функции. Иными словами, чтобы

рассчитать Фк+1, нужно знать значения Фк+1

и Ф

к+1

которые, в свою очередь,

у , -1 j и '*'у_1 ,

вычисляются через значения сеточной функции в левом и нижнем узлах сеточного шаблона разностной схемы. А это означает, что вычисления ведутся последовательно в заранее предопределенном порядке. Распараллеливание вычислений по рекуррентным формулам представляет определенную проблему для компьютера с параллельной архитектурой, поэтому рассмотрим другой, отличный от используемого в (3), способ обхода узлов сетки, который основывается на их «красно-черном» упорядочивании [3]. В соответствии с этим приемом (рис.3) все узлы

сетки делятся на два подмножества Ой = и (сеточные узлы «раскрашива-

ются» в черные и красные цвета в шахматном порядке), а поскольку при получении разностных схем использовался шаблон «крест», то, как видно из рис. 3, для расчета значения сеточной функции фк+1 в узле черного цвета нужны значения

ФкпЬ в узлах красного цвета и наоборот. В таком случае, используя порядок

«Черный» узел

«Красный» узел

Рис. 3. «Красно-черное» упорядочивание узлов сетки

вычислений значений сеточной функции сначала во всех узлах черного цвета, а затем в оставшихся красного, метод релаксации (3) преобразуется в двухэтапный вычислительный процесс с использованием явных (не рекуррентных) формул:

к+1 \

и ) = Ю'

и >Б

(ЬИ (<1.,), + С, (<,), + ), + е. (фк,+, ), + , )

(ф:

(*е(

)б + С,

+ (1 -ю) (Ф (Ф

к

у /Б

к+1 ) V

;(и Л) епБ;

<и(Ф

к+1

iI--}б

«1 )б + (Л'),

(4)

+(1 -ю) (фк.),;(/, Л) еП;

Здесь символ «Л» означает обход узлов красного цвета, а «В» - черного. При выбранном способе обхода узлов одна итерация в методе релаксации сводится к двум этапам использования явных расчетных формул (типа формул метода Якоби), применяемого сначала к расчету сеточных значений в <<черных>> узлах, а затем в <<красных>>, при этом объем вычислений на каждой итерации последовательного алгоритма остается прежним.

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

Время работы последовательной программы

В таблице показано время работы программы, построенной по последовательному алгоритму с «красно-черным» упорядочиванием, в зависимости от размеров используемой сетки. Расчеты проводились на одном вычислительном узле кластера ТГУ СКИФ СуЬепа [4]. Использовался компилятор Ше! 12.

Время работы программы при использовании различных сеток (точность - е = 0,001)

Сетка Итераций Время, с

128x128 1004 21

256x256 3176 294

512x512 12160 5037

1024x1024 53713 107655

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

2 раза число итераций увеличивается нелинейно, возрастает объем вычислений и время работы программы при каждом изменении размеров сетки увеличивается более чем на порядок. Тем самым становится актуальным сокращение времени работы программы за счет использования ЭВМ с параллельной архитектурой.

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

ух в среднем поперечном сечении х = Ьх , получаемый в результате вычислений, справа построен график поперечной скорости V в среднем продольном сечении у = Ь (Яе = иъх / ц = 1000 ).

X. п Ly 0,8 -

0,6 -

Гл

/ 0,2Н

V

/

/

-0,4

0

0,4

0,8 Vx /U

Vy .

U 0,2 ■

0 І -0,2 --0,4 --0,6 -

x/Lx

1 і ч

0,2 0,4

0,8 j X t

YJ

♦

Рис. 4. График продольной скорости (при х=£х/2) и поперечной (при у = £/2) соответственно, полученный на сетке 256x256, значки - расчеты [5]

Видно, что даже при использовании для конвективных членов противопотоко-вой аппроксимации первого порядка на подробной сетке имеет место хорошее согласование с расчетами, полученными на основе использования разностной схемы более высокого порядка на относительно более грубой сетке [5].

Ускорение программы на многопроцессорных ЭВМ

Построение параллельной версии алгоритма SIMPLE осуществлялось на основе принципа геометрической декомпозиции [6] сеточной области, когда вся область исследования делится на равные по площади (или по количеству сеточных узлов) части, вычисления в которых следует проводить одновременно и независимо. Как указывалось выше, основная вычислительная сложность алгоритма SIMPLE - это решение на каждой глобальной итерации систем сеточных уравнений для компонент скорости и поправки давления. Для решения таких систем в данной работе предлагается использовать метод релаксации, поскольку другие более быстро сходящиеся при последовательных вычислениях методы (GMRES, BiCGStab, CG и др.) при параллельной реализации на основе метода геометрической декомпозиции показывают увеличение общего объема вычислительной работы (числа итераций для обеспечения сходимости с заданной точностью) при возрастании количества применяемых в вычислениях процессоров по сравнению с последовательной версией [7].

В данной работе использовался более эффективная по сравнению с одномерной двумерная декомпозиция сеточной области [8], когда вся расчетная область делится на двумерные блоки одинаковых размеров (рис. 5), в каждом из которых значения сеточных функций вычисляются одновременно и независимо по явным формулам (4). Однако для обеспечения таких идеально параллельных вычислений необходимо обеспечить каждую подобласть дополнительными значениями сеточ-

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

Рис. 5. Геометрическая декомпозиция данных

Разработанный параллельный алгоритм был реализован на параллельных ЭВМ с различной архитектурой: вычислительный кластер ТГУ СКИФ Cyberia [4] и сервер с видеокартой NVIDIA GTX260, которая включает 192 процессора с небольшой тактовой частотой. При создании параллельной программы для кластера использовался стандарт Message Passing Interface [8], а для сервера с графическими процессорами - технология Compute Unified Device Architecture [9]. Параллельная версия программы для каждого случая получается на основе последовательной, при этом применение библиотеки MPI изменяет код последовательной программы лишь частично, в то время как для технологии CUDA требуется полное перестроение последовательной программы.

На рис. 6 показано ускорение (отношение времени работы последовательной программы к времени работы параллельной) параллельных программ, написанных с использованием CUDA или MPI-интерфейса, на различных сетках. MPI-программа была запущенна на 32 и 64 вычислительных ядрах кластера СКИФ Cyberia. Из рисунка видно, что ускорение параллельных программ растет с увеличением размеров сетки, что связано с увеличением объема вычислительной работы по сравнению с комммуникационными затратами. Для решения рассматриваемой задачи на сетках с различной плотностью узлов производительность сервера с видеокартой оказалась лучше 32 вычислительных ядер кластера, но хуже чем 64, однако с экономической точки точки зрения расчеты на видеокарте GTX 260 предпочтительнее использования 64 вычислительных ядер кластера СКИФ Cyberia в силу меньшего энергопотребления и стоимости.

CUDA

MPI-32 PROC MPI-64 PROC

80 70 60 „ 50 I 40 §

& 30 20 10 0

128x128 256x256 512x512 1024x1024

Размер сетки

Рис. 6. Ускорение параллельных программ на различных сетках

На сетке 1024x1024 время расчета последовательной программы составляет около 30 часов, согласно рис. 6, время работы программы с использованием сервера с видеокартой NVIDIA GTX260 сократилось до ~ 37 мин, и до ~ 24 мин с использованием 64 вычислительных ядер кластера СКИФ СуЬепа.

Заключение

Для решения задач вычислительной гидродинамики разработаны версии алгоритма SIMPLE, ориентированные на использование высокопроизводительной вычислительной техники с параллельной архитектурой, а именно многопроцессорного linux-кластера с распределенной памятью или сервера с графическими ускорителями. В основе построенных параллельных алгоритмов лежит применение принципа неодномерной геометрической декомпозиции, «красно-черного» упорядочивания при обходе узлов сетки и метода релаксации для решения сеточных уравнений. В результате полученные параллельные алгоритмы численного решения уравнений Навье - Стокса обладают замечательным свойством прямо пропорционального уменьшения числа арифметических операций, выполняемых одним процессором/ядром, при увеличении общего количества используемых вычислительных процессоров/ядер. При решении задачи о течении вязкой несжимаемой жидкости в каверне с движущейся верхней крышкой полученное свойство параллельных алгоритмов позволило обеспечить ускорение в вычислениях на сетках с более чем 10б узлов в несколько десятков раз. Кроме того, анализ результатов вычислений показал высокую эффективность построенного параллельного алгоритма на графических процессорах, что существенно расширяет возможности при чиссленном исследовании задач механики жидкости и газа.

ЛИТЕРАТУРА

1. Лойцянский Л.Г. Механика жидкости и газа: учебник для вузов. 6-е изд., перераб. и доп.

М.: Наука, 1987. 840 с.

2. Патанкар С. Численные методы решения задач теплообмена и динамики жидкости. М.:

Энергоатомиздат, 1984. 124 с.

3. Ортега Дж. Введение в параллельные и векторные методы решения линейных систем / под. ред. Х.Д. Икрамова. М.: Мир, 1991. 367 с.

4. http://skif.tsu.ru

5. U., K.N. Ghia and C.T. Shin. High-Resolutions for incompressible flow using the Navier-Stokes equations and a multigrid method // J. Comp. Phys. 1982. V. 48. P. 387-411.

6. Богословский Н.Н., Есаулов А.О., Старченко А.В. Параллельная реализация алгоритма вычислительной гидродинамики SIMPLE // Сибирская школа семинар по параллельным вычислениям. Томск: Изд-во Том. ун-та, 2002. С. 118-124.

7. Старченко А.В., Данилкин Е.А. Параллельная реализация численного метода решения системы уравнений Навье - Стокса при моделировании крупных вихрей турбулентных течений // Вестник Новосибирского государственного университета. Серия: Информационные технологии. Новосибирск: НГУ, 2009. Т. 7. № 2. С. 49-61.

8. Беликов Д.А. Говязов И.В. и др. Высокопроизводительные вычисления на кластерах: учеб. пособие / под ред. А.В. Старченко. Томск: Изд-во Том. ун-та, 2008. 198 с.

9. Боресков А.В., Харламов А.А. Основы работы с технологией CUDA. М.: ДМК Пресс, 2010. 232 с.

Статья поступила 13.02.2012 г.

Degi D.V., Starchenko A.V. NUMERICAL SOLUTION OF NAVIER-STOKES EQUATIONS ON COMPUTERS WITH PARALLEL ARCHITECTURE. In this paper, we numerically solve the Navier-Stokes equations using high-performance computing technology, such as systems with distributed memory and graphic accelerators. As a test problem, we consider the classical problem of computational fluid dynamics, namely, the motion of a fluid in a rectangular cavity.

Keywords: Navier-Stokes equations, OpenMP system, MPI library, CUDA technology.

DEGI Dmitrii Vladimirovich (Tomsk State University)

E-mail: [email protected]

STARCHENKO Alexandr Vasil’evich (Tomsk State University)

E-mail: [email protected]