Numerical solution of Navier-Stokes equations on computerswith parallel architecture.pdf Создание новых эффективных параллельных вычислительных систем являетсястратегическим направлением развития компьютерной техники. Это обстоятель-ство вызвано как ограничением максимально возможного быстродействия обыч-ных последовательных ЭВМ, так и существованием вычислительных задач, длярешения которых ресурсов существующих средств вычислительной техники невсегда оказывается достаточным. Такие современные проблемы науки и техникикак моделирования климата, генная инженерия, проектирование интегральныхсхем, анализ загрязнения окружающей среды, создание лекарственных препара-тов, требуют для анализа ЭВМ с производительностью более 1000 миллиардовопераций с плавающей точкой в секунду (1 TFLOPS).К настоящему времени создано большое количество вычислительных алго-ритмов, ориентированных на последовательную модель программирования. Од-нако далеко не всегда удается разработать эффективную параллельную реализа-цию для многих из них. Поэтому актуальной в настоящее время становится про-блема создания новых параллельных алгоритмов решения задач науки и техникина суперкомпьютерах, в частности, для имеющих большие теоретическое и при-кладное значение задач гидродинамики.Основными уравнениями гидродинамики являются уравнения Навье - Стокса.С помощью уравнений Навье - Стокса описываются такие процессы, как движе-ние волн на поверхности жидкости, движение крови по кровеносным сосудам че-ловека и животных, течения в мантии Земли. Задача о движение жидкости в ка-верне также описывается уравнениями Навье - Стокса, и в ней отражаются сле-дующие характерные черты гидродинамических процессов: конвективная нели-нейность, различные соотношения между инерционными и вязкими силами, од-новременное существование областей малых и больших градиентов и т.п., поэто-му задача о каверне широко распространена в качестве «тестовой» при выбореэффективного численного алгоритма для моделирования гидродинамическихпроцессов [1]. Зачастую для получения приближенного решения с высокой точно-стью используют сетку с большой плотностью узлов, при этом существенно воз-растает время расчета и тем самым становится актуальным применение много-процессорных ЭВМ для уменьшения времени расчетов.Физическая постановка задачиРассматривается стационарное, рециркуляционное, ламинарное течение не-сжимаемой, вязкой жидкости в прямоугольной полости (каверне) с верхней стен-кой, движущейся в своей плоскости со скоростью U (рис. 1).Рис. 1. Движение жидкости в кавернеЖидкость, целиком заполняющая каверну, вовлекается в движение силами вяз-кости. Действие массовых сил (сил тяжести, инерционных, центробежных или ко-риолисовых сил) незначительное. Имеет место случай изотермического движениявязкой несжимаемой жидкости с постоянными значениями плотности и коэффи-циентов вязкости.Кроме того, требуется выполнение условия прилипания частиц жидкости ктвердой стенке. Это означает отсутствие скольжения жидкости по поверхности.Таким образом, выполняется граничное условие равенства нулю скорости жидко-сти на поверхности неподвижных стенок и совпадение скорости жидкости со ско-ростью движущейся стенки, с которыми частицы жидкости соприкасаются (ра-венство касательной компоненты скорости жидкости и скорости стенки).Математическая постановка задачиРассмотренный физический процесс можно описать с помощью следующихуравнений:- уравнения движения в проекции на ось Ox:2 22 21x x y x x xv v v v p v vx y x x y∂ ∂ ∂ ⎛∂ ∂ ⎞+ =− +μ⎜ + ⎟∂ ∂ ρ∂ ⎝∂ ∂ ⎠;- уравнения движения в проекции на ось Oу:2 22 21x y y y y y v v v v p v vx y y x y∂ ∂ ∂ ⎛∂ ∂ ⎞∂ + ∂ =−ρ∂+μ⎝⎜⎜∂ +∂ ⎠⎟⎟;- уравнения неразрывности:0 x yv vx y∂ ∂+ =∂ ∂.Здесь ( , ) x yv = v v - скорость движения жидкости, p - давление, ρ - плотностьжидкости, μ - коэффициент кинематической вязкости.Граничные условия в соответствии с физической постановкой будут следую-щими:- из условия непротекания:0 : 0; : 0;0 : 0; : 0;x x xy y yx v x L vy v y L v= = = == = = =- из условия прилипания:0 : 0; : 0;0 : 0; : ,y x yx y xx v x L vy v y L v U= = = == = = =где xL и yL - соответственно длина и ширина каверны, U - скорость движениякрышки.Получение разностной схемыПриближенное решение ищется на равномерной шахматной сетке hΩ (рис. 2)с помощью метода конечных объемов [2].hx(vx)ij hy(vy)ij(vx)i+1 j(vy)i j+1pijРис. 2. Фрагмент используемой шахматной сетки(штриховка нанесена на конечный объем)При получении разностной схемы используется следующее:1. Значения сеточных неизвестных ( , x yv v ) рассматриваются на соответст-вующих гранях конечных объемов, а давление ( ) pij - в их центрах.2. Приближенное интегрирование членов дифференциальных уравнений про-водится с помощью квадратурной формулы средних прямоугольников.3. Для устранения возникающей нелинейности в разностной схеме значениеодного из множителей следующих произведений: x xv v , x yv v , y yv v - полагаетсяравным значению, полученному на предыдущей итерации.4. Конвективные члены аппроксимируются с помощью схемы «против пото-ка», диффузионные и производные от давления - центрально-разностной схемой.5. Значения сеточных функций с дробными индексами ( 1/2),( 1/i± j± 2), из-вестные с предыдущей итерации, в точках, не совпадающих с узлами сетки, рас-считываются из значений в ближайших соседних узлах с помощью линейной ин-терполяции.В итоге получается, что разностная схема, аппроксимирующая уравнение не-разрывности, имеет вид1j 1 [( ) ( ) ] [( ) ( ) ] 0, 0, ; 0, ; x i x ij y y i j y ij x v v h v v h i Nx j Ny + +− + − = = =разностная схема для проекции уравнения движения по оси Ох1 j 1 j 1 2 1 1 1 +11 j 1 1 1 1( ) ( ) ( ) ( )( ) ( ), 0, 1, 0, ;i x i i j x i j i j x ij i j x i ji x i j i j ij i jaP v aE v aW v aN vaS v d p p i Nx j Ny+ + + + + + ++ + − + += + + ++ + − = − =разностная схема для проекции уравнения движения по оси Оyj+1 j+1 1 1 1 1 1 1 1 21 1 1( ) ( ) ( ) ( )( ) ( ), 0, , 0, 1 .i yi ij y i j ij yi j ij y iji j y ij i j ij i jbP v bE v bW v bN vbS v l p p i Nx j Ny+ + + + − + + ++ + += + + ++ + − = = −Здесь x h и yh - шаги сетки по горизонтальной и вертикальной оси соответствен-но, 1/i j xl h+= ρ и 1/i j yd h+ = ρ, x N и yN - размер сетки по горизонтали и вертика-ли соответственно. Коэффициенты ij, ij, ij , ij , ij , ij , ij , ij , ij , ijaP aE aW aN aS bP bE bW bN bS- известные сеточные функции, зависящие от поля скорости, полученного на пре-дыдущей итерации. У этих коэффициентов имеются следующие свойства [2]:- неотрицательность всех коэффициентов;- ij ij ij ij ijaP ≥aE +aW + aN + aS и ij ij ij ij ijbP ≥bE +bW + bN + bS , причем строгоенеравенство выполняется на границе.Предполагая, что функции ( , ) x x v =v xy , 蕔( , )y y v v xy = и p p(x,y) = , представ-ляющие точное решение, имеют достаточное число непрерывно дифференцируе-мых частных производных, и, используя формулу Тейлора, можно доказать, чтопогрешность аппроксимации разностной схемы для уравнения неразрывностиимеет второй порядок относительно каждого координатного направления, а дляуравнений движения - первый.Алгоритм численного метода решенияВ силу нелинейности уравнений Навье - Стокса становится невозможным пря-мое решение системы и рассматриваемый процесс приобретает итерационный ха-рактер. Для совместного решения общей системы из трех систем разностныхуравнений был применен алгоритм SIMPLE [2], который можно записать в сле-дующем виде:1. Задаем поле давления { } *pij ;2. Решаем системы уравнений11* * * *1 1 1* * * *1 1 1( ) ( ) ( ),( ) ( ) ( )i ji ji j x i j nb x nb i j ijnbi j y i j nb y nb i j ijnbaP v a v d p pbP v b v l p p+++ + ++ + +⎧ = + −⎪⎨= + − ⎪⎩ΣΣ(1)92и находим{( ) } *xijv и {( ) } *yijv ,где*( ) nb x nbnbΣa v - означает суммирование по соседним узлам сеточного шаблона;3. Решаем систему уравнений* * * * 1 j1j 1 11+11 1 11[( ) ( ) ] [( ) ( ) ] [ ( )( )] [ ( ) ( )] 0ix i x ij y y i j y ij x ij i ji jij i j iji j ij y ij ij ij ij xij i j ijdv v h v v h p paPd l lp p h p p p p haP bP bP++ + ++− + −+− + − + ′ − ′ −′ − ′ + ′ − ′ − ′ − ′ =(2)и определяем поправку давления{ } pij′ ;4. По соотношениям11*1 1 1 1*1 1 1 1( ) ( ) ( )/ ,( ) ( ) ( )/i ji jx i j x i j i j ij i jy i j y i j i j ij i jv v d p p aPv v l p p bP+++ + + ++ + + +⎧ ЃЊ ЃЊ = + − ⎪⎨ЃЊ ЃЊ = + − ⎪⎩вычисляем сеточные функции {( ) } x ijv и {( ) } y ijv , удовлетворяющие разностномууравнению неразрывности;5. Корректируем давление по формуле*p p p ′ = +ω (0< ω Σa =b + c + d + e имеет место на границе рассматриваемойобласти, что является достаточным для сходимости итерационного процесса[3, с. 89].При 0 < ω < 1 получается метод нижней релаксации. Использование методанижней релаксации замедляет изменение зависимой переменной. С помощьюданного метода решались уравнения для скоростей.При 1 < ω < 2 имеет место метод верхней релаксации. Он увеличивает ско-рость сходимости итерационного процесса. Уравнение для поправки давлениярешалось данным методом.Из рекуррентной формулы (3) видно, что производится строго последователь-ный пересчет значений неизвестной сеточной функции. Иными словами, чтобырассчитать k 1ij+ Φ , нужно знать значения 11ki j+−Φ и 11kij+−Φ , которые, в свою очередь,вычисляются через значения сеточной функции в левом и нижнем узлах сеточно-го шаблона разностной схемы. А это означает, что вычисления ведутся последо-вательно в заранее предопределенном порядке. Распараллеливание вычисленийпо рекуррентным формулам представляет определенную проблему для компью-тера с параллельной архитектурой, поэтому рассмотрим другой, отличный от ис-пользуемого в (3), способ обхода узлов сетки, который основывается на их «крас-но-черном» упорядочивании [3]. В соответствии с этим приемом (рис.3) все узлысетки делятся на два подмножества R Bh h hΩ =Ω ∪Ω (сеточные узлы «раскрашива-ются» в черные и красные цвета в шахматном порядке), а поскольку при получе-нии разностных схем использовался шаблон «крест», то, как видно из рис. 3, длярасчета значения сеточной функции k 1ij+ Φ в узле черного цвета нужны значенияknbΦ в узлах красного цвета и наоборот. В таком случае, используя порядок«Черный» узел «Красный» узелРис. 3. «Красно-черное» упорядочивание узлов сеткивычислений значений сеточной функции сначала во всех узлах черного цвета, азатем в оставшихся красного, метод релаксации (3) преобразуется в двухэтапныйвычислительный процесс с использованием явных (не рекуррентных) формул:( )( ) ( ) ( ) ( ( ) ( ) )( )( )(( ) ( ) ( ) ( ) ( ))( )1 1 1 111 1 1 11 1 1 11(1 ) ; ( , ) ;(1 ) ; ( , ) .k k k kk ij i jR ij i jR ij i j R ij i j R ij Bij Bijk Bij B hk k k kk ij i jB ij i jB ij i j B ij i j B ij Rij Rijk Rij R hb c d e fai jb c d e fai j− + − +++ + + +− + − ++⎛Φ + Φ + Φ + Φ +⎞Φ =ω⋅⎜ ⎟+⎜ ⎟⎝ ⎠+ −ω Φ ∈Ω⎛Φ + Φ + Φ + Φ +⎞Φ =ω⋅⎜ ⎟+⎜ ⎟⎝ ⎠+ −ω Φ ∈Ω(4)Здесь символ «R» означает обход узлов красного цвета, а «B» - черного. При вы-бранном способе обхода узлов одна итерация в методе релаксации сводится кдвум этапам использования явных расчетных формул (типа формул метода Яко-би), применяемого сначала к расчету сеточных значений в ‹‹черных›› узлах, а за-тем в ‹‹красных››, при этом объем вычислений на каждой итерации последова-тельного алгоритма остается прежним.Другой важной особенностью метода ‹‹красно-черного›› упорядочивания яв-ляется то, что количество итераций последовательного и параллельного алгорит-мов при проведении вычислений с заданной точностью для задачи одного и тогоже размера совпадают. Это позволяет утверждать, что при увеличении количествапроцессоров для решения рассматриваемой задачи, вычислительная нагрузка накаждый из них будет прямо пропорционально уменьшаться.Время работы последовательной программыВ таблице показано время работы программы, построенной по последователь-ному алгоритму с «красно-черным» упорядочиванием, в зависимости от размеровиспользуемой сетки. Расчеты проводились на одном вычислительном узле кла-стера ТГУ СКИФ Cyberia [4]. Использовался компилятор Intel 12.Время работы программы при использовании различных сеток(точность - ε = 0,001)Сетка Итераций Время, с128.128 1004 21256.256 3176 294512.512 12160 50371024.1024 53713 107655Из таблицы видно, что при уменьшении каждого шага равномерной сетки в2 раза число итераций увеличивается нелинейно, возрастает объем вычислений ивремя работы программы при каждом изменении размеров сетки увеличиваетсяболее чем на порядок. Тем самым становится актуальным сокращение времениработы программы за счет использования ЭВМ с параллельной архитектурой.Для проверки правильности работы построенного алгоритма и программы былипроведены сравнения расчетных данных с результатами, представленными в [5].Слева на рис. 4 для квадратной каверны показан профиль продольной скоростиv в среднем поперечном сечении xx=L , получаемый в результате вычислений,справа построен график поперечной скорости y v в среднем продольном сеченииy y = L ( Re / 1000x=UL μ= ).yLyvUyv U x /x L / x0,20,60,8-0,4 0 0,4 0,8-0,4-0,6-0,20,200,2 0,4 0,8Рис. 4. График продольной скорости (при x=Lx/2) и поперечной (при y = Ly/2)соответственно, полученный на сетке 256.256, значки - расчеты [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 в силу меньшего энергопотребления и стоимости.Рис. 6. Ускорение параллельных программ на различных сеткахНа сетке 1024.1024 время расчета последовательной программы составляетоколо 30 часов, согласно рис. 6, время работы программы с использованием сер-вера с видеокартой NVIDIA GTX260 сократилось до ~ 37 мин, и до ~ 24 мин сиспользованием 64 вычислительных ядер кластера СКИФ Cyberia.ЗаключениеДля решения задач вычислительной гидродинамики разработаны версии алго-ритма SIMPLE, ориентированные на использование высокопроизводительной вы-числительной техники с параллельной архитектурой, а именно многопроцессор-ного linux-кластера с распределенной памятью или сервера с графическими уско-рителями. В основе построенных параллельных алгоритмов лежит применениепринципа неодномерной геометрической декомпозиции, «красно-черного» упоря-дочивания при обходе узлов сетки и метода релаксации для решения сеточныхуравнений. В результате полученные параллельные алгоритмы численного реше-ния уравнений Навье - Стокса обладают замечательным свойством прямо про-порционального уменьшения числа арифметических операций, выполняемых од-ним процессором/ядром, при увеличении общего количества используемых вы-числительных процессоров/ядер. При решении задачи о течении вязкой несжи-маемой жидкости в каверне с движущейся верхней крышкой полученное свойствопараллельных алгоритмов позволило обеспечить ускорение в вычислениях насетках с более чем 106 узлов в несколько десятков раз. Кроме того, анализ резуль-татов вычислений показал высокую эффективность построенного параллельногоалгоритма на графических процессорах, что существенно расширяет возможностипри чиссленном исследовании задач механики жидкости и газа.

Скачать электронную версию публикации