Ускорение расчетов на графических процессорах при исследовании течения Стокса методом граничных элементов Текст научной статьи по специальности «Физика»

УДК 532.5:519.6

Ускорение расчетов на графических процессорах при исследовании течения Стокса методом граничных элементов

1 3 1 3 1 2

О. А. Солнышкина , Ю. А. Иткулова , Н. А. Гумеров

[email protected], [email protected], [email protected] 1 «Центр микро- и наномасштабная динамика дисперсных систем Башкирского государственного университета»

(ЦМНДДС БашГУ)

2 Institute for Advanced Computer Studies, University of Maryland, USA 3 Институт механики им. Р. Р. Мавлютова Уфимского научного центра РАН

Поступило в редакцию 22.04.2012

Аннотация. В работе исследуется динамика двух вязких несмешивающихся жидкостей в неограниченной области и течение вязкой жидкости в канале сложной геометрии в трехмерной постановке. Численный подход основывается на методе граничных элементов, используются треугольные сетки. При увеличении масштаба задач возникает нехватка памяти вычислительной системы. Эта проблема решалась заменой прямого метода решения итерационным GMRES. Для ускорения расчетов разработан модуль матрично-векторного произведения «на лету» и распараллелен на графических процессорах ^Ри). Проведено сравнение полученных численных результатов с аналитическими решениями. Представлены результаты по эффективности использования GPU.

Ключевые слова. Уравнения Стокса; метод граничных элементов; деформируемые капли; течение эмульсии в канале; GMRES; графические процессоры

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

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

Работа выполнена при поддержке Министерства образования и науки Российской Федерации грант №11.G34.31.0040, соглашение 8861, грант РФФИ №12-01-31173, Fantalgo, LLC (Maryland, USA)

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

В работе исследуются трехмерные течения Стокса двух вязких несмешивающихся жидкостей в неограниченной области [2] и несжимаемой вязкой жидкости в канале переменного сечения. Рассматриваемые задачи будут лежать в основе численного исследования эффекта динамического запирания эмульсий в микроканалах. Метод граничных элементов для Стоксовых течений изложен в [3] и успешно применялся для расчета течений в каналах [4] и взаимодействия капель и частиц в дисперсных течениях [5]. В упомянутых работах также можно найти обзор литературы.

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

лировать процессы в геометрически сложных объектах, границы которых дискретизируются треугольными сетками. Программная реализация исследуемых задач предусматривает выбор оптимальных алгоритмов в зависимости от количества узлов сетки. Для ускорения расчетов разработан модуль матрично-векторного произведения без хранения матрицы в памяти вычислительной системы, который используется в итерационном решателе СМЯЕ8. и распараллелен как на обычном многоядерном процессоре (СР11), так и на графических процессорах (СРи) с использованием технологии С1ГОА.

ТРЕХМЕРНОЕ ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ДИНАМИКИ КАПЛИ В НЕОГРАНИЧЕННОМ ПОТОКЕ ЖИДКОСТИ

Математическая модель

Постановка задачи об обтекании капли ньютоновской жидкости плотности р2 и вязкости |12 потоком другой ньютоновской жидкости плотности р| и вязкости Ц| приведена, например, в [2] и заключается в следующем. Рассматриваемые жидкости занимают объем и У\ соответственно, Л' - граница раздела фаз. На бесконечности скорость потока обтекающей каплю жидкости 1 равна иоа(х). При числах Рейнольдса Явь Яе2 значительно меньше 1, поведение потока жидкости 7, обтекающего каплю, и поведение жидкости 2 внутри капли описывается уравнениями Стокса

V ст; = -Ург + ^У2и; = О,

V • и; = 0, г = 1,2, где ст - тензор напряжений; р - давление, которое включает гидростатический компонент.

Условия на границе двух жидкостей и, = и2 = ч, { = СТ • П = /п,

/ = у(У • п) + (р2 - рх)^ • х), х £ 5,

где I - разность векторов нормального напряжения в жидкости 1 и 2; п - нормаль к поверхности, направленная в жидкость 1; у - коэффициент поверхностного натяжения; g - ускорение свободного падения; х - радиус-вектор рассматриваемой точки.

Условие на бесконечности

11! (х) -> ию(х),

X -> СО.

Кинематическое условие, описывающее динамику капли с1х

(1)

Численное моделирование

Задача решалась методом граничных элементов. Суть метода состоит в преобразовании дифференциального уравнения в частных производных, описывающего поведение неизвестной функции внутри и на границе области, в интегральное уравнение, связывающее только граничные значения, и поиске численного решения этого уравнения. Если требуется найти значения искомой функции во внутренних точках области, то их можно вычислить, используя известные решения на границе. Таким образом, размерность задачи уменьшается на единицу. Для данной поверхности 5” скорость и в произвольной точке у может быть вычислена через граничные интегралы [3]

уеУг, и (у) — 2ию(у) у £ У2, Аи(у)

1 + А

У —— и(у) -ию(у);

= J (--Є(у,х) • ^х) -

-(1 - Л)[Т(у,х) • п(х)] • и(х))гі5(х),

где ц = ц,: X = ц2/ці - отношение вязкостей внутренней и внешней жидкостей; Є и Т - тензоры второго и третьего ранга, компоненты которых в прямоугольной декартовой системе координат имеют вид:

1 (5ц П г/

(3)

Тц( У,х) = -

3 Г і Г: Гк

(4)

1] Ап г5

П=У1~ х» и], к = 1,2,3, где 5у - символ Кронекера. Поверхность Л' разбивается на N узлов хг, по которым строятся квадратурные формулы граничных интегралов.

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

N

“т“и;+2/1 - 'и*=

г=і

N

— = (и • п)п, х Є 5.

(2)

= = 1..........................(5)

1[С)(У) = [ С(х,у) • п(х)с*5(х),

1Г(У) = [ т (Х'У) ‘ п(х)^5(х).

Вычисление сингулярных элементов производится на основе известных интегральных тождеств [3].

Рис. 1. Динамика капли

Решая систему линейных алгебраических уравнений (5), получаем компоненты скорости на границе.

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

Наиболее трудоемкой процедурой является вычисление средней кривизны поверхности в каждом узле сетки. Средняя кривизна вычислялась двумя методами: контурных интегралов и установления параболоида [5]. Ошибка вычисления кривизны находилась из сравнения с аналитическим выражением для эллипсоида. Метод установления параболоида дает более точное значение средней кривизны с погрешностью не более 3 %, но для его использования валентность всех узлов сетки должна быть не менее пяти.

Результаты тестирования для одиночной капли

Из рисунков видно, что изначально эллипсоидальная капля (рис. 1 (слева)) под действием поверхностного натяжения приобретает сферическую форму (рис. 1 (справа)). Число узлов сетки N варьировалось в пределах 600-3000.

Было проведено сравнение численного решения с аналитическим [6] в 3 этапа:

1. Компоненты скорости на границе капли, полученные аналитически, сравнивались с численными значениями, полученными методом граничных элементов. Погрешность составила 0,5 % для N = 2562.

2. Компоненты скорости внутри и вне капли, полученные соответственно из первого и второго уравнений (3), используя для вычисления правой части точные значения Г(х), и(х), сравниваются с аналитическими решениями.

3. Компоненты скорости внутри и вне капли, полученные соответственно из первого и второго уравнений (3), используя для вычисления правой части численные значения Г(х), и(х), сравниваются с аналитическими решениями.

На рис. 2 показаны результаты, полученные на 2-м и 3-м этапах сравнения с аналитическим решением для X = 2,5. Результаты сравнения показывают хорошее приближение численных расчетов к аналитическим с некоторой погрешностью, которая зависит от соотношения вязкостей рассматриваемых жидкостей X и количества узлов дискретизации. Трехмерное численное моделирование течения вязкой жидкости в канале переменного сечения

Математические модели периодического и непериодического потока

Рассматривается течение вязкой несжимаемой жидкости с вязкостью ц в ограниченном канале переменного сечения с гладкой поверхностью 5” при малых числах Рейнольдса (рис. 3).

Рис. 3. Канал переменного сечения

Движение жидкости описывается уравнениями Стокса (1).

Аналитическое решение на границе. МГЭ

Ух Линии тока Уз

т

0.1

*0.1

-4-2 0 2 4

-4-2 0 2 4

ЯП

Численно е решение на грашще. МГЭ

0.1

-4-2 0 2 4

т

-4-2 0 2 4 2

О

0.1

0.1

-0.1

о

-5

-4-2 0 2 4

Аналитическое решение

Рис. 2. Компоненты скоростей внутри и вне капли

Поверхность канала разбивалась на два подмножества: 5 = 5^ и 5^, где 5^г) - жесткая граница и - мягкая граница. Задача решается со следующими граничными условиями:

• для непериодического потока

и(х) = и°(х), х е 5,

Г(х„) = Г°(х0), х0 е 5(5).

• для периодического потока

и(х) = 0, х е 5(г),

Г(х0) = Г°(х0), х0 е 5(5), и(0, у, г) = и(1, у, г),

Г (0 ,у,г) = -ЦЬ,у,г) + 1хАр.

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

Методы воплощены для произвольной трехмерной геометрии, в тестовой задаче аналитическая форма канала осесимметрична.

На рис. 4 слева представлена модель канала для непериодического потока, справа — модель экспериментального канала для периодического потока.

Метод граничных элементов в нескольких модификациях

Уравнения Стокса (1) переписываются в граничных интегралах [3]. Сначала определяется вектор нормального напряжения Г на границе

- [ в (у, х) ■ Г(х)с15(х) =

№ (7)

= £и°(у) - 4 К (у,х) • и°(х)с*5(х), у е 5,

где К(у,х) = Т(у,х) • п(х) , С, Т находятся по

формуле (4).

Затем значения скорости и в любой точке области у вычисляются через граничные интегралы

и (у) = - [ й (у,х) • {(х)(15(х) +

(8)

+ I К (у, х) • и°(х)<±!>(х), у £ V.

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

Г = Г° + Г, и = и° + и'. (9)

Граничные условия (6) для возмущений и' и Г выглядят следующим образом

и' = 0, х е 5^, и' = —и0, х е Б(г\

Г(х0) • п(х0) = 0, х0 е 5^.

Из интегральных уравнений (7)-(8), которым удовлетворяют возмущения и' и 1^, находятся их значения. С учетом (9) определяется вектор нормального напряжения Г на границе и компоненты скорости и в любой точке области.

Рис. 4. Триангуляция канала переменного сечения

Для периодического потока уравнения Стокса (1) переписываются в граничных интегралах

^и(у) - [ К (у, х) • и(х)с?5(х) =

\ (11)

- [ С (у, х) • f(x)dS(x), у £ 5.

Интегральные уравнения (11) переписываются в матричном виде

Аи = ВГ, А = -1 — К, В = -С.

2 д

Граничные условия имеют вид

иг = 0, и2 = и3 = и5, ^ = {г,

f2 = *8, *3 = + К, fp = .

где 1 — боковая поверхность; 2 — входное сечение; 3 — выходное сечение.

Из (11) получаем систему линейных алгебраических уравнений относительно неизвестных (и5,{г~)

/•^12 ( А22 \а32

/в13^

= В2^р \B33fp/

+ А13 В13 " “ ^12 Вц\ /Щ

+ •^23 ^23 " “ ^22 -В21 Г5

+ А33 В3З " “ ^32 -В31/

Ух

о

-10 1

У:

0.2

0.1

0

Аналитическое решение

Линии тока

6 4

2

Уг

х 0.5

и 0.5

0 0.5 1

О

-10 1

0.2

0.1

О

Непериодический поток. МГЭ

6(

4

-3

; 10

И

-10 1

I

о

-10 1

Непериодический поток. Модифицированный МГЭ

6г

х 10

-18

0.2

0.1

0

□

-10 1

У

Периодический поток. МГЭ

0.2

0.1

0

У

X 10

2

1

0

•1

■2

Рис. 5. Сравнение модификаций МГЭ с аналитическим решением

Сингулярные интегралы рассчитываются, используя известные тестовые решения: постоянное и линейное течения.

Результаты тестирования для канала

Для тестирования программы использовался канал c входным и выходным сечением радиуса Я = 1, радиусом сужения г = 0,5 и длиной Ь = 5.

Проведено сравнение компонент вектора нормального напряжения Я на границе методом граничных элементов в нескольких модификациях с аналитическим решения для непериодического течения Пуазейля. Погрешность вычисления неизвестных значений на границе составляет 0,7 %. Для периодического потока погрешность компонент вектора нормального напряжения Я на границе составляет 0,24 %, для компонент вектора скорости и около 0,9 %.

На рис. 5 представлены поля скорости и линии тока, вычисленные в осевом сечении цилиндра плоскостью Охг, для течения Пуазейля. Приведены случаи для канала с N = 2598 для метода коллокаций в узлах сетки и N = 2704 для метода колокаций с центром в грани элемента. Из рисунка видно, что численные расчеты хорошо приближаются к аналитическому решению. Модифицированный метод граничных элементов, особенно для компонент скорости вдоль оси Ог, дает более точные результаты по сравнению с обычным методом граничных элементов.

Непериодический поток

Ух Метод граничных эелементои '/г

5.

► -<

Модифицированный метод граничных элементов 5.

Периодический поток

Рис. 6. Канал переменного сечения

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

В результате, относительная погрешность составила 0,1-1,7 % в зависимости от формы канала. По мере уменьшения радиуса сужения канала погрешность увеличивается и принимает максимальное значение 1,7 % на цилиндрическом канале.

На рис. 6 представлены поля скорости для каналов разной геометрии и радиусом сужения г = 0,5. Сглаженный канал состоит из 2598 вершин и 5192 граней, экспериментальный — из 1354 вершин и 2704 граней. Таким образом, количество расчетных узлов N не превышает 3000 узлов. На границе сингулярные интегралы вычисляются неточно, что ведет к большой погрешности вычисления поля скорости в точках, лежащих близко к стенке канала. Наблюдаемый эффект на входе для периодического потока качественно совпадает с экспериментальным [7] и приближенным аналитическим решением [8].

Ускорение расчетов на графических процессорах

При трехмерном численном моделировании физических процессов для областей со сложной геометрией, например течение эмульсии в микроканалах переменного сечения, необходимо построение сеток с большим количеством узлов. Решение подобных многомасштабных задач требует разработки и применения эффективных численных методов. Для сеток маленького размера при решении СЛАУ применялись прямые методы, но при увеличении масштаба задачи их использование затрудняется. Это связано с тем, что размер необходимой памяти пропорционален квадрату числа узлов сетки, также при их увеличении возрастает время вычислений. При использовании прямых методов, начиная с некоторого количества узлов, возникает нехватка памяти вычислительной системы.

Эту проблему можно решить, используя итерационные методы решения, которые существенно снижают затраты памяти и времени. Наиболее эффективными и устойчивыми среди итерационных методов решения таких систем уравнений являются так называемые проекционные методы, и особенно тот их класс, который связан с проектированием на подпространства Крылова (например, GMR.ES). Эти методы обладают целым рядом достоинств: они устойчивы, допускают эффективное распараллелива-

ние, работу с различными строчными (столбцовыми) форматами и предобуславливателями разных типов [9].

Для эффективной программной реализации итерационного метода необходимо решить две проблемы:

• разработать подпрограмму, быстро умножающую матрицу на вектор;

• ускорить сходимость метода с помощью предобуславливателя.

В рамках данной работы для решения проблем, связанных с использованием памяти, в среде Matlab разработан программный модуль умножения матрицы на вектор без хранения матрицы в памяти системы («MV product on the fly», который используется в GMRES при решении СЛАУ. Каждый элемент матриц G и К вычисляется по следующим формулам:

^тп хп)

о • 8 pi

(чхт xrl) (хт XR)

X-vvj X-V)

|X|7| X™

к

4 _

xn)

(12)

3 {Xin-Xn){xL-Xn) 4 pi 11 |xm-xn|5

m,n = 1, ... ,N, i,j, к = 1,2,3, где N - количество узлов дискретизации области. Таким образом, матрицы G и К имеют размер 3 N х 3 N и 3 N х N.

Применение модуля матрично-векторного произведения позволяет решить проблему ограничения по памяти вычислительной системы, но матричные вычисления являются вычислительно-трудоемкими, поэтому представляют собой классическую область применения параллельных вычислений. Для ускорения расчетов модуль распараллелен на центральном многоядерном процессоре (CPU) средствами Matlab Parallel Computing Toolbox и с помощью про-

граммно-аппаратной технологии CUDA на графических процессорах (GPU) [10].

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

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

В рамках данной работы распараллеливание модуля «MV product on the fly» на GPU основывается на разбиении матриц горизонтальными полосами на т частей так, что N = т / L. где N— размерность матрицы, L - число строк матрицы в блоке, т - количество потоков на GPU.

На каждой итерации каждый из т потоков вычисляет свою часть вектора решения. То есть при умножении матрицы G размера 3N х 3N на вектор b размера 3N / 1 на одной итерации поток вычисляет очередной свой блок матрицы размером 3 х 3 по формулам (12) и умножает полученный блок на соответствующую часть вектора b размером 3 х 1 и прибавляет результат к полученному на предыдущей итерации вектору (рис. 7). Каждый поток хранит часть результирующего вектора, которая по окончании вычислений копируется в глобальную память для получения полного вектора решения.

о о о; о о о:о о о о о о

о о О; о о о:о о о о о о

о о О: о о 0;0 о о о о о

о о о о о 0:0 о о о о о

о о о о о 0:0 о о о о о

о о ° о о о:о о о о о о

о о о: о о 0^0 о о о о о

о о о: о о о':<=> о о 0 о о

о о о: о о ОІО о о о о о

о о о о о о-о о о о о о

о о О о о 0*0 о о о о о

о о О: о о о-о о о о о о

о

о

о

о

о

о

"o'

о

о

о

о

о

о

о

о

о

о

о

о

о

о

о

о

о

Рис. 7. Разделение данных и организация вычислений при выполнении параллельного модуля матрично-векторного произведения на ОРИ

— □ -MVP onlhelty

—A— parallel CPU MVP Ihe liy

—МаїйЬМИ3

•—*— parallel GPU M\4= rtlhe Ity

I ю’ І

Ю1 Matrix si» 10

Рис. 8. Время вычисления матричновекторного произведения

На рис. S показано сравнение времени выполнения модуля матрично-векторного произведения на одном ядре CPU (MVP on the fly), на S ядрах CPU (parallel MVP on the fly), на GPU (parallel GPU MVP on the fly) и встроенной функции умножения матрицы на вектор Matlab (Matlab MVP).

Все расчеты проводились на вычислительной системе с CPU Intel Xeon 566О, 2.8GHz, GPU NVIDIA Tesla C2050. При проведении вычислительных экспериментов установлено, что для различных размеров задачи наилучшие результаты по времени достигаются при размере блока равном 256 потокам.

Рис. 9. Ускорение расчетов в зависимости от размера матрицы

Рис. lO. Производительность в зависимости от размера матрицы

На рис. 9 приведено ускорение расчета на GPU по сравнению с CPU для операций с двойной и одинарной точностью. При расчетах на графической карте NVIDIA Tesla C2050 для количества узлов сетки N от 1 000 до 100 000 для уравнений Стокса получено ускорение до 520 раз для операций с двойной точностью и до 700 раз для операций с одинарной точностью.

При проведении тестовых расчетов достигнута следующая производительность: до 220 Gflops для операций с двойной точностью и 300 Gflops для операций с одинарной точностью (рис. 10). Учитывая, что пиковая производительность графической карты Tesla C2050 для чисел с плавающей точкой двойной точности составляет 515 Gflops, одинарной точности -1,03 Tflops, получены хорошие результаты.

Результаты тестов на графической карте NVIDIA Tesla C2050 показали возможность решения граничных задач для уравнений Стокса размером до 100 000 элементов на одной рабочей станции.

ЗАКЛЮЧЕНИЕ

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

Для ускорения расчетов разработан модуль матрично-векторного произведения «на лету», то есть без хранения матрицы в памяти вычислительной системы, который используется в GMRES, и распараллелен на графических процессорах (GPU). Представлены результаты по эффективности использования GPU.

Решеные задачи являются начальным шагом для дальнейшего моделирования поведения водонефтяных эмульсий в микроканалах.

СПИСОК ЛИТЕРАТУРЫ

1. Ахметов А. Т., Саметов С. П. Особенности течения дисперсии из микрокапель воды в микроканалах // Письма в ЖТФ. Т. 36, вып. 22, 2010. С. 21-28.

2. Rallison J.M., Acrivos A. A numerical study of the deformation and burst of a viscous drop in an extensional flow // J. Fluid Mech. 1978. Vol. 89, part 1. P. 191-200.

3. Pozrikidis C. Boundary Integral and Singularity Methods for Linearized Viscous Flow. 1992 (Cambridge University Press, Cambridge, MA).

4. Pozrikidis C. Creeping flow in two-dimensional channels// J. Fluid Mech. 1987. Vol. 180. P. 495-514.

5. Zinchenko A. Z. and Davis R.H. An efficient algorithm for hydrodynamical interaction of many deformable drops // J. Comp. Phys. 2000. Vol. 157. P. 539-587.

6. Хаппель Дж., Бреннер Г. Гидродинамика при малых числах Рейнольдса. М.: Мир, 1976. 623 c.

7. Boger D. V. Viscoelastic flows through contractions //J. Fluid Mech. 1987. Vol. 19. Р. 157-182

8. An approximate solution to flow through a contraction for high Trouton ratio fluids / A. S. Lubansky [et al.] // Non-Newtonian Fluid Mech. 2007. Vol. 144. P. 87-97

9. Saad Y. Iterative Methods for Sparse Linear System. 2000, SIAM.

10. NVIDIA Corporation. NVIDIA CUDA Compute Unified Device Architecture Programming Guide. Version

3.2.2010.

ОБ АВТОРАХ

Солнышкина Ольга Александровна, асп. каф. прикладной информатики и численных методов БашГУ, стажер-исследователь Центра микро- и наномасштабной динамики дисперсных систем БашГУ, мл. науч. сотр. ИМех УНЦ РАН. Магистр прикладной математики и информатики (2009).

Иткулова Юлия Айратовна, асп. той же каф., стажер-исследователь того же Центра, м.н.с. ИМех УНЦ РАН. Магистр прикладной математики и информатики (2011).

Гумеров Наиль Асгатович, д-р физ.-мат. наук, руководитель вычислительной лаборатории того же Центра, проф. Ин-та передовых компьютерных исследований университета Мэриленда.

METADATA

Title: Acceleration of the calculations on the graphics processors for the study of Stokes flow by boundary element method. Authors: O. A. Solnyshkina1,3, Yu. A. Itkulova1,3, N. A. Gu-merov1,2 Affiliation:

1 Center for Micro- and Nanoscale Dynamics of Dispersed Systems, Bashkir State University, Ufa, Russia.

2 Insitute for Advanced Computer Studies University of Maryland, USA.

3 Institute of Mechanics Russian Academy of Sciences Ufa Branch, Ufa, Russia.

Email: '[email protected].

Language: Russian.

Source: Vestnik UGATU (Scientific journal of Ufa State Aviation Technical University), 2013, Vol. 16, No. 6 (51), pp. 92100. ISSN 2225-2789 (Online), ISSN 1992-6502 (Print).

Abstract: The dynamics of two immiscible liquids in an unbounded domain and viscous liquid flow in a channel with complex geometry in three dimensional case are studied in the present work. The numerical approach is based on the boundary element method and triangular meshes are used. As the problem size increase, the deficiency of the computer memory arises. This problem is solved by replacing an direct method of solution by a iterative method GMRES. To accelerate the calculations the module for a matrix-vector product "on the fly" is developed and parallelized on graphics processors (GPU). The obtained results are compared with the analytical solutions. The results of the efficiency GPU are performed.

Key words: Stokes equations, boundary element method, deformable droplets, emulsion channel flow, GMRES, graphics processors.

References (English Transliteration):

1. Аhmetov А. T, Sametov S. P. Features of flow dispersion of water microdroplets in microchannels. // Letters in JTP. 2010. Vol. 36, No. 22. P. 21-28.

2. Rallison J.M., Acrivos A. A numerical study of the deformation and burst of a viscous drop in an extensional flow // J. Fluid Mech. Vol. 89. part 1, 1978. P. 191-200.

3. Pozrikidis C. Boundary Integral and Singularity Methods for Linearized Viscous Flow. 1992(Cambridge University Press, Cambridge, MA).

4. Pozrikidis C. Creeping flow in two-dimensional channels// J. Fluid Mech. Vol. 180. 1987. P. 495-514.

5. Zinchenko A. Z. and Davis R.H. An efficient algorithm for hydrodynamical interaction of many deformable drops // J. Comp. Phys. Vol. 157. 2000. P. 539-587.

6. Happel J., Brenner H. Low Reynolds number hydrodynamics. 1976. P. 623.

7. Boger D.V. Viscoelastic flows through contractions //J. Fluid Mech. 1987. Vol. 19. Р. 157-182.

8. Lubansky A.S., Boger D.V., Servais C. Burbidge A.S. Cooper-White J.J. An approximate solution to flow through a contraction for high Trouton ratio fluids// Non-Newtonian Fluid Mech. 2007. Vol. 144. P. 87-97.

9. Saad Y. Iterative Methods for Sparse Linear System. 2000, SIAM.

10. NVIDIA Corporation. NVIDIA CUDA Compute Unified Device Architecture Programming Guide. Version

3.2.2010.

About authors:

1. Solnyshkina, Olga Aleksandrovna, Postgrad. (PhD) Student, Dept. of Applied Informatics and Numerical Methods. Junior Researcher in the Center for Micro- and Nanoscale Dynamics of Dispersed Systems (CMNDDS), BSU, Junior Researcher in the Institute of Mechanics Russian Academy of Sciences Ufa Branch, Master of Applied Mathematics and Informatics (BSU, 2009).

2. Itkulova, Yuliya Airatovna, Postgrad. (PhD) Student, Dept. of Applied Informatics and Numerical Methods. Junior Researcher in the CMNDDS, BSU, Junior Researcher in the Institute of Mechanics Russian Academy of Sciences Ufa Branch, Master of Applied Mathematics and Informatics (BSU, 2011).

2. Gumerov, Nail Asgatovich, Dr. of Phys. and Math Sci, Header of the Computational Lab in the CMNDDS, BSU, Research Prof. of the Institute for Advanced Computer Studies, University of Maryland, USA.