Численное решение задач гидроаэромеханики на графических процессорах

Карпенко, Антон Геннадьевич

Численное решение задач гидроаэромеханики на графических процессорах тема автореферата и диссертации по механике, 01.02.05 ВАК РФ

Карпенко, Антон Геннадьевич АВТОР

кандидата физико-математических наук УЧЕНАЯ СТЕПЕНЬ

Санкт-Петербург МЕСТО ЗАЩИТЫ

2013 ГОД ЗАЩИТЫ

01.02.05 КОД ВАК РФ

Диссертация по механике на тему «Численное решение задач гидроаэромеханики на графических процессорах»

Автореферат диссертации на тему "Численное решение задач гидроаэромеханики на графических процессорах"

САНКТ-ПЕТЕРБУРГСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ

На правах рукописи

00553414/

Карпенко Антон Геннадьевич

ЧИСЛЕННОЕ РЕШЕНИЕ ЗАДАЧ ГИДРОАЭРОМЕХАНИКИ НА ГРАФИЧЕСКИХ ПРОЦЕССОРАХ

Специальность 01.02.05 — Механика жидкости, газа и плазмы

3 ОКТ 2013

Автореферат диссертации на соискание учёной степени кандидата физико-математических наук

Санкт-Петербург — 2013

005534147

Работа выполнена на кафедре гидроаэромеханики математико-механического факультета Санкт-Петербургского государственного университета.

Научный руководитель: доктор физико-математических наук,

профессор Матвеев Сергей Константинович

Официальные оппоненты: доктор технических наук, профессор

Гурьев Юрий Владимирович (ВУНЦ ВМФ «Военно-морская академия им. Н.Г. Кузнецова»), заведующий кафедрой «Механика и гидромеханика»

кандидат физико-математических наук, доцент Бестужева Алла Николаевна, доцент кафедры прикладной математики (Петербургский государственный университет путей сообщения)

Ведущая организация: Федеральное государственное унитарное предприятие

Российский федеральный ядерный центр — Всероссийский научно-исследовательский институт экспериментальной физики ФГУП «РФЯЦ - ВНИИЭФ>

Защита состоится /¿¡¿-Т 2013 г. в #часов на заседании диссертационного совета Д 211.232.30 по защите диссертаций на соискание ученой степени кандидата наук, на соискание ученой степени доктора наук при Санкт-Петербургском государственном университете по адресу: 198504, Санкт-Петербург, Петродво-рец, Университетский пр., 28, математико-механический факультет, ауд. 405.

С диссертацией можно ознакомиться в Научной библиотеке им. М.Горького Санкт-Петербургского государственного университета по адресу: 199034, Сантк-Петербург, Университетская наб., 7/9.

Автореферат разослан "_"_2013 года.

Ученый секретарь диссертационного совета Д 211.232.30, д.ф.-м.н., проф.

Кустова Е.В.

Общая характеристика работы

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

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

Одним из способов ускорения вычислений является распараллеливание. В традиционном подходе производится декомпозиция расчетной области и каждая область рассчитывается на своем центральном процессорном устройстве (ЦПУ). Но для обеспечения согласованного решения необходима пересылка данных между процессорами. На пересылку данных тратится значительная часть времени. В последнее время появились графические процессорные устройства (ГПУ) ориентированные на научные вычисления. Они представляют собой векторные процессоры с большим количеством арифметико-логических устройств (АЛУ), в отличие от ЦПУ у которых всего несколько АЛУ. Для ГПУ необходимо разрабатывать новые схемы расчета учитывающие особенности архитектуры. Необходимо разрабатывать подходы к описанию течения жидкости и газа имеющие параллелизм на уровне данных. Это позволит использовать весь потенциал ГПУ и существенно ускорить вычисления. Вновь разработанные схемы становятся конкурентоспособными со схемами разработанными для ЦПУ.

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

Для достижения поставленной цели необходимо было решить следующие задачи:

1. Разработать, реализовать и протестировать схемы расчета нестационарного одномерного течения газа на графических процессорах;

2. Разработать, реализовать и протестировать схемы расчета нестационарного вязкого сжимаемого течения газа на неструктурированных сетках в многомерном случае, используя ГПУ;

3. Разработать и реализовать универсальный подход расчета течения газа при малых числах Маха и течения несжимаемой жидкости;

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

Основные результаты:

1. Разработаны, реализованы и протестированы схемы расчета нестационарного одмерного течения газа на графических процессорах. Показана возможность и преимущество использования ГПУ. Разработан, реализован и протестирован программный комплекс расчета на ГПУ нестационарного вязкого сжимаемого течения газа на неструктурированных сетках в многомерном случае.

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

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

4. Разработанный комплекс использован при решении практических задач теплогидродинамической конвекции (ТГК) при транспортировки застывающих наливных грузов (ЗНГ).

Научная новизна:

1. Разработаны модификации схем расчета нестационарных одномерных течений газа на графических процессорных устройствах ГПУ.

2. Разработаны схемы расчета нестационарного вязкого сжимаемого газа на

неструктурированных сетках в многомерном случай с использованием ГПУ.

3. Разработан универсальный подход расчета течений газа при малых числах Маха и несжимаемых течений, а также течений со свободной конвекцией. Подход основан на предобусловливании системы уравнений сжимаемого газа.

4. Впервые аналитически выведены матрицы диагонализации якобиана потока предобусловленной системы в многомерном случае.

Практическая значимость Разработка универсального программного комплекса расчета задач гидроаэромеханики на ГПУ позволит повысить быстродействие и обеспечит оперативность вычислений. Результаты работы и разработанные схемы расчета были внедрены в пакет программ трехмерного имитационного моделирования на супер-ЭВМ ЛОГОС, являющегося отечественным аналогом Пакет программ инженерного анализа ЛОГОС предназначен

для моделирования процессов тепломассопереноса (аэродинамика, газодинамика, гидродинамика, теплопроводность) и разрабатывается Российским федеральным ядерным центром — Всероссийский научно-исследовательский институт экспериментальной физики ФГУП «РФЯЦ-ВНИИЭФ» г. Саров.

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

Апробация работы. Основные результаты работы докладывались на:

• XIII Международный семинар "Супервычисления и математическое моделирование" (Саров, ФГУП «РФЯЦ-ВНИИЭФ», 3-7 октября, 2011);

• XVIII Школа-семинар молодых ученых и специалистов под руководством академика РАН А.И. Леонтьева "Проблемы газодинамики и тепломассообмена в новых энергетических технологиях" (Звенигород, 23-27 мая, 2011);

• XXII Всероссийский семинар с международным участием по струйным, отрывным и нестационарным течениям (Санкт-Петрбург, 22-25 июня, 2010);

• XVII Школа-семинар молодых ученых и специалистов под руководством академика РАН А.И.Леонтьева "Проблемы газодинамики и теплообмена в аэрокосмических технологиях" (Жуковский, 25-29 мая, 2009);

• Международная конференция по механике и баллистике «VI Окуневские чтения» (Санкт-Петербург, 23-27 июня, 2008);

• XVI Школа-семинар молодых ученых и специалистов под руководством академика РАН А.И.Леонтьева "Проблемы газодинамики и теплообмена в энергетических установках" (Санкт-Петербург, 21-25 мая, 2007);

• III Всероссийская научная конференция "Проектирование научных и инженерных приложений в среде MATLAB" (Санкт-Петербург, 23-26 октября, 2007).

Публикации. Основные результаты по теме диссертации изложены в 12 печатных изданиях, 3 из которых изданы в журналах, входящих в перечень рецензируемых научных журналов и изданий, 7 — в трудах российских, международных конференций и семинаров. В работе [1] автору принадлежит разработка и реализация метода постановки дозвуковых граничных условий для одномерной системы уравнений газовой динамики, а также расчет тестовых задач. Емельянов В.Н. в этой работе выступал в роли научного консультанта и разработчика идеи постановки характеристических граничных условий. В работе [2] автору принадлежит обзор архитектуры графических процессоров, использования их для расчетов и технологии программирования. Соавторам принадлежит реализация схемы решения уравнения Лапласа, расчет течения в каверне. В работе [3] автору принадлежит разработка и реализация метода контрольного объема на графических процессорах, решение одномерной задачи Сода и задачи об ударной трубе в трехмерной постановке. Смирнов П.Г. и Тетерина И.В. выполнили расчет течения в плоском канале. Волкову К.Н. принадлежит подробное описание алгоритма расчетов и общей архитектуры ГПУ. В работах [4], [5] автором представлен обзор возможности и перспективы использования ГПУ для расчета течений жидкости и газа. В работе [6] Емельянову В.Н. принадлежит общая постановка задач. В работах [10], [12] Азарову М.А. принадлежит разработка математической модели течения в лазерном резонаторе, а автору расчет протекающих газодинамических процессов.

Объем и структура работы. Диссертация состоит из введения, четырех глав, заключения и приложения. Полный объем диссертации 178 страниц текста с 84 рисунками и 3 таблицами. Список литературы содержит 98 наименований.

Содержание работы

Во введении обосновывается актуальность выбранной темы исследования, формулируется цель работы, ее научная новизна и практическая значимость, сформулированы результаты выносимые на защиту.

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

для вычислений. Дан обзор инструментов и технологий программирования ГПУ. Произведен обзор литературы по использованию ГПУ для расчетов течений жидкости и газа.

Вторая глава посвящена использованию графических процессоров для решения задач гидроаэромеханики. На первом этапе производится численное решение системы уравнений нестационарного одномерного течения газа на ГПУ. Рассматривается численное решение краевой задачи для системы уравнений одномерной нестационарной газовой динамики:

дИ ) п тт ( 9 ^ „(тп

+ = 0, где и = \ ри , =

Щх,0) = и<Р\х), (1)

им = щг), и{Ь,1) = иг{ь).

Рассматривается эволюция начального распределения функции консервативных переменных и^\х) в области х € [О, Ц во времени. Вначале и в конце рассматриваемой области задаются значения функций от времени. Такой постановкой можно описывать течение в ударной трубе после разрыва мембраны. Поэтому часто эту задачу называют задачей об ударной трубе.

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

Ограничение на шаг по времени вытекает из условия, что волны, образовавшиеся после распада разрыва на грани контрольного объема, не достигли его центра, а при более слабом ограничении - другой грани. Схема Эйлера первого порядка обеспечивает положительность решения если Ьсрь = Дг— < х> где Дг = т;Д:гг - расстояние между центром ячейки и ближайшей гранью. Это условие должно быть выполнено для всех ячеек, поэтому за шаг интегрирования берется наименьший.

Различные версии метода контрольного объема различаются способом

вычисления потоков на грани 1'\, х. В работе Годунова предложено вычислять

поток путем решения задачи о распаде произвольного разрыва на грани между

соседними ячейками, так как решение является автомодельным. Для расчета по-

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

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

А+;

(ил — и£) можно рассматривать как некую диссипативную добавку в схему с центральными разностями.

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

На рис. 1 представлено распределение плотности для тестовых задач. Сплошной линией показано точное аналитическое решение, красной линей отмечено решение с первым порядком аппроксимации по пространству, синей со вторым.

Таблица 1: Тестовые задачи.

Тест PL «L PL PR UR PR tfin

1 Задача Сода 1.0 0.75 1.0 0.125 0.0 0.1 0.2

2 Две волны разрежения 1.0 -2.0 0.4 1.0 2.0 0.4 0.12

3 Сильная ударная волна 1.0 0.0 1000.0 1.0 0.0 0.01 0.012

4 Движущийся контактный разрыв 1.4 0.1 1.0 1.0 0.1 1.0 2.0

Рис. 1: Распределение плотности.

Расчеты производились на одном из модулей Tesla SI070 с частотой ядер 1.44 ГГц (количество ядер 256) и на одном ядре ЦП AMD Phenom 2 с частотой ядра 3 ГГц. Для сравнения скорости вычисления расчеты проводились на сетках с различным числом ячеек. При переходе от сетки к сетке число ячеек увеличивается на порядок (1024 ячеек для сетки 1, 30720 ячеек для сетки 2 и 307200 ячеек для сетки 3). Сетка наилучшей разрешающей способности содержит около 3 миллионов ячеек (сетка 4). Время, необходимое для расчета одного шага (производилось осреднение по 10 шагам), а также ускорение S решения задачи Сода приводятся в табл. 2 (время приводится в миллисекундах). Вариант 1 соответствует расчету по схеме Годунова, вариант 2 — расчету по схеме Рое.

Density. t=0.1200? dh=0.01

Таблица 2: Задача Сода. Время и ускорение вычислений.

№ Сетка 1 Сетка 2 Сетка 3 Сетка 4

ЦПУ ГПУ 5 ЦПУ ГПУ 5 ЦПУ ГПУ 5 ЦПУ ГПУ 5

1 2 1.63 0.14 0.13 0.07 12.43 1.87 47.70 5.51 0.20 0.17 245.25 33.17 460.64 43.58 0.92 0.57 502.50 76.00 4627.61 436.09 8.06 5.22 574.39 83.48

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

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

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

Взаимодействие прямой ударной волны с клином. В рассматриваемой задаче прямая ударная волна с числом Маха ударной волны М$ = 1.7 и интенсивностью I = 3.2 распространяется вдоль оси ОХ рис. 2 слева. Ударная волна набегает на клин, расположенный под углом ф = 25 В общем случае, в зависимости от и ф существует четыре возможных типа взаимодействия. При заданных параметрах происходит нерегулярное отражение с образованием прямого скачка уплотнения.

І__ х

Рис. 2: Слева схема расчетной области. Справа поле градиента плотности вблизи передней

кромки клина.

Для визуализации ударно волновых структур в вычислительной газодинамики часто вместо полей параметров используют поле модуля градиента плотности рис. 2 справа.

Расчет трансзвукового обтекания профиля №АСА0012 под углом атаки. Выполняется расчет обтекания плоского крыла трансзвуковым потоком воздуха. Схема расчетной области показана на рис. 3 слева внизу. На внешней границе расчетной области задаются граничные условия типа Римана, число Маха потока Моо = 0.75, поток натекает на профиль под углом атаки а = 4

Рис. 3: Слева вверху прямоугольная расчетная сетка. Справа вверху треугольная расчетная сетка. Слева внизу схема расчетной области. Справа внизу поле чисел Маха.

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

данном течении локально образуется сверхзвуковая зона над профилем. На рис. 3 справа внизу представлено поле чисел Маха.

Обтекание бесконечно тонкой пластинки вязким потоком. Рассчитывалось обтекание бесконечно тонкой пластинки длиной Ь = 1и вязким потоком воздуха. Скорость набегающего потока равна V = 10 м/с. Для разрешения пограничного слоя к пластинке выполнено сгущение. Для того чтобы сетка у краев пластинки не была сильно вытянута производилось сгущение к краям.

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

Рис. 4: Слева поле модуля скорости. Справа распределение безразмерной скорости в сечениях х = 0.01 м и х = 0.1 м.

Произведено сравнение времени расчета на ГПУ и центральном процессоре. Расчеты производились на сетках различного размера от 10000 ячеек до 10 млн., с учетом и без учета вязкости, с различными порядками аппроксимации по времени и по пространству. Расчеты производились на одном из модулей Tesla S1070 с частотой ядер 1.44 ГГЦ (количество ядер 256) и на одном ядре ЦП AMD Phenom 2 с частотой ядра 3 ГГц. На данном оборудовании удалось получить ускорение в 42 раза, но для реальных вязких задач необходимо производить расчет со вторым порядком аппроксимации по пространству. В этом случае удалось достигнуть ускорения лишь в 22 раза. Графические процессоры постоянно

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

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

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

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

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

Стандартные методы расчета сжимаемых течений, описанные в предыдущих главах, оказываются неэффективными при расчете течений с малыми числами Маха. Это происходит из за жесткости системы уравнений для сжимаемого

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

Первой работой описывающей подход видоизменения основной системы уравнений или предобусловливание была работа СЬопп. Применение сжимаемых алгоритмов для расчета системы уравнений несжимаемаемой жидкости выполнено путем добавления в уравнение неразрывности производной от дав-

с др г,

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

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

3<3 дЬ дх

д и

где - якобиан перехода от консервативных

дУ

(2)

переменных II =

эи

(р,ри,ру,ри>,рЕ) к физическим С? = {р,и, и, ш,Т). Матрица Г = Р— называется матрицей предобусловливания и выбирается таким образом, чтобы устранить жесткость системы при малых числах Маха. На первый взгляд изменение системы уравнения приведет к решению отличному от решения исходной системы, но в случае расчета стационарных течений первый член стремится к нулю. В случае расчета нестационарных течений используется метод двойных шагов. Все существующие методы отличаются выбором матрицы предобусловливания. В данной работе матрица предобусловливания Г имеет вид:

Г =

Параметр 0 водится следующим образом:

\и? РСР.

С его помощью происходит модификация исходной системы уравнений. Параметр 11г имеет смысл скорости распространения возмущений давления. Дополнительный член -4?- введен для упрощения аналитических выкладок.

е 0 0 0 Рт

9 их Р 0 0 РТУх

6 Уу 0 р 0 РТ'"у

Qvz 0 0 Р

ен-1 рух ртН + рСр

Для совершенного газа регулирующий параметр (7Г можно ввести следующим образом:

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

Метод контрольного объема в многомерном случае с использованием матрицы предобусловливания можно записать так:

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

где в диссипативной добавке производится диагонализация якобиана потока пре-

др

добусловлештой системы упавнений Аг = Г-1-. В паботе впервые выведены

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

Часто бывает интересно стационарное решение задачи течения в канале. Решение такой задачи может быть получено методом установления. Суть метода состоит в том, что решение стационарной задачи отыскивается как предел к которому стремится решение нестационарных уравнений при £ —> оо. Интегрирование системы уравнений по времени ведется до тех пор, пока не будет

ЄС, если |у| < ЕС 17г — |у|, если ЕС < |у| < с с, если |у| > с

= 2 № + - зГМгІЛгІМ^ДО,

достигнут стационарный режим. О достижении стационарного режимов можно судить по тому, как отличается решение на слое п от решения на слое п + 1. Это можно делать, оценивая Ь-2 норму приращения сеточного решения. Решение можно считать сошедшимся, если Ь-2 < где е некоторая заданная малая величина.

Будем рассматривать течение в канале круглого сопла, профиль которого

/ \ 1/2

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

При перепаде давления Ар = 1600 Па, число Маха в критическом сечении приближенно будет равно М = 0.3. При таких скоростях потока, подходы с использованием предобусловливания и без предобусловливания дают одинаковую скорость сходимости. Это говорит о том, что для чисел Маха М > 0.3 нет необходимости использовать предобусловливание.

При уменьшении перепада давления до Ар =175 Па, число Маха в критическом сечении приближенно будет равно М = 0.1. Решение полученное по схеме Роэ при таких числах Маха начинает отличатся от решения полученного с использованием предобусловливания. Это связано с тем что в схеме Роэ используются абсолютные значения давления и жесткостью системы уравнений при малых числах Маха. На рис. 5 представлены графики сходимости решения при использовании 20 ячеек. Линия 1 соответствует невязке скорости, а линия 2 невязке давления. Видно, что решение без использования предобусловливания не сходится, тогда как с использованием сходится. Стоит отметить, что распределения параметров в данных случаях отличаются. Это говорит о плохой сходимости алгоритма без предобусловливания при числах Маха М < 0.1.

Рис. 5: Сходимость решения. Слева без использования предобусловливания, справа с

прсдобусловливанисм

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

Для демонстрации возможности подхода был произведен расчет свободно конвективного течения около вертикально поставленной нагретой пластины. При описании свободно-конвективных течений используется приближение Бу-синесска. Основная идея которого состоит в том, что зависимость плотности от температуры учитывается только при вычислении массовых сил. Расчетная область является двумерной и замкнутой, течение в области является нестационарным, но в окрестности пластины, после некоторого промежутка времени, параметры слабо меняются и можно считать течение установившимся. В начальный момент времени воздух покоится имеет равномерную температуру Т = 273 К, затем пластина приобретает температуру Т = 300 К и из за разности плотности нагретого и холодного газа возникает конвективное течение. На пластинке выставлялись условия прилипания и температура стенки.

распределение безразмерной температуры У=0.03

Рис. 6: Профиль безразмерной температуры в сечении х = 0.03 м.

•Эта задача может быть решена аналитически в приближении пограничного слоя. В этом подходе вводят безразмерные температуру скорость £ и координату г]. Стоит отметить, что данная постановка задачи отличается от постановки в приближении пограничного слоя. В данном случае решаются полные уравнения Навье-Стокса. Поэтому профили скорости у кромки пластинки будут отличаться от автомодельного решения. На рис. 6 представлено сравнение профилей безразмерной температуры в сечении у = 0.03 м. полученного численно

с автомодельным решением. Из графиков видно, что численное решение хорошо совпадает с автомодельным решением.

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

В четвертой главе произведен расчет практической задачи о конвекции мазута в железнодорожной цистерне.

Жидкие нефтепродукты, в том числе и мазут перевозятся в цистернах но железной дороге в холодное время года. При охлаждении у мазутов и масел происходит резкий рост вязкости. Такие вязкие нефтепродукты при перевозках охлаждаются и загустевают настолько, что их выгрузка без разогрева становится невозможной. При температурах, близких к температурам застывания, большинство застывающих наливных грузов (ЗНГ) становятся структурированными реологическими средами. Вязкость таких сред зависит от напряжения сдвига. При высоких же температурах, когда вязкость меньше 0.3 • итах можно применять модель Ньютоновской жидкости.

На станции загрузки жидкие топлива заливаются горячими с температурой около Тж = 70 °С. Температура окружающей среды в зимнее время в некоторых регионах может опускаться до Тв = —40 °С. При этом, время ожидания цистерн к отправке достигает несколько десятков часов. При таких условиях в цистерне из за разности температур развивается естественная или теплогид-родинамическая конвекция (ТГК). При этом имеет место существенный отток тепла из цистерны. В работе Моисеева В.И. предлагается подавлять ТГК около стенки цистерны встречным восходящим потоком жидкости. Этот поток жидкости предлагается организовывать с помощью тепловых аккумуляторов (ТА). Тепловые аккумуляторы представляют собой тела с большой теплоемкостью.

Тепловые аккумуляторы предлагается установить по контуру в цистерну. В начале загрузки температура нефтепродуктов равна Тж = 70 "С и ТА нагреваются до той же температуры. В процессе охлаждения цистерны температура нефтепродуктов уменьшается, но температура ТА остается почти такой же и снижается намного медленнее. Тем самым, около ТА образуется ТКГ направленная вверх, против течения около стенки цистерны рис. 7. Этим предлагается снижать интенсивность ТГК около стенки цистерны и уменьшать теплоотдачу.

Рис. 7: Схема расчетной области для цистерн без тепловых аккумуляторов (слева) и с ТА

(справа).

Будем рассматривать течения мазута марки М40. Это один из самых распространенных видов топлива. Зависимость вязкости нефтепродуктов от температуры хорошо описываются экспоненциальной зависимостью:

где для мазута марки М40 параметр А = 4.2127 • 10~4, а параметр В = 0.0325. Так как цистерна достаточно длинная, будем рассматривать двумерную задачу. Расчетная сетка является неструктурированной и к стенке цистерны и ТА выполнено сгущение ячеек для разрешения пограничного слоя. Очевидно, задача имеет плоскость симметрии, поэтому рассматривается только половина области, а на границе выставляются условия симметрии. Так как время охлаждение достаточно длительное, будем рассматривать задачу в квазистационарной постановке. Будем рассчитывать ТГК в некоторый момент времени, когда жидкость уже остыла до температуры Тж = 50 °С, а тепловые аккумуляторы до температуры Тта = 55 °С. Поэтому, в начале расчета в обоих случаях жидкость покоится и имеет температуру Тж = 50 °С. На стенке цистерны выставляется температура равная Тст = 0°С, а на стенке тепловых аккумуляторов Тта = 55 °С. Производится расчет установившегося течения при таких параметрах. Это течение соответствует некоторому моменту времени охлаждения. В результате расчета выполняется сравнение интенсивности течения около холодной стенки цистерны в обоих случаях.

I/(Т) = Ае вт,

На рис. 8 представлено поле скорости для варианта без установки тепловых аккумуляторов и с установкой. Видно, что восходящий поток от ТА направлен в противоположную сторону от течения около стенки. Встречный поток замедляет течение около стенки, тем самым снижая теплоотдачу.

Рис. 8: Векторное поле скорости для цистерн без тепловых аккумуляторов (слева) и с ТА

(справа)

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

1. Разработаны, реализованы и протестированы схемы расчета нестационарного одномерного течения газа на графических процессорах. Показана возможность и преимущество использования ГПУ. Разработан, реализован и протестирован программный комплекс расчета на ГПУ нестационарного вязкого сжимаемого течения газа на неструктурированных сетках в многомерном случае.

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

4. Разработанный комплекс использован при решении практических задач теплогидродинамической конвекции (ТГК) при транспортировке застывающих наливных грузов (ЗНГ).

В Приложении А подробно описан метод решения задачи о распаде произвольного разрыва методом последовательных приближений.

В Приложепии Б приведены аналитические выводы собственных чисел и матриц диагонализации якобиана потока предобусловлеиной системы уравнений. Также приведен вид матрицы обратной к матрице предобусловливания в многомерном случае.

Публикации автора по теме диссертации

1. Емельянов В.Н., Карпенко А.Г. Метод постановки дозвуковых граничных условий для системы уравнений газовой динамики // Вестник тихоокеанского государственного университета. 2013. N 2 (29). С. 29-38. (Из перечня ВАК)

2. Волков КН., Емельянов В.Н., Карпенко А.Г., Курова И.В., Серов А.Е., Смирнов П.Г. Численное решение задач гидродинамики на графических процессорах общего назначения // Вычислительные методы и программирование. 2013. Т. 14. N 1. С. 82-90. (Из перечпя ВАК)

3. Волков КН., Емельянов В.Н., Карпенко А.Г., Смирнов П.Г., Тетерина И.В. Реализация метода конечных объемов и расчет течений вязкого сжимаемого газа на графических процессорах // Вычислительные методы и программирование. 2013. Т. 14. N 1. С. 183-194. (Из перечня ВАК)

4. Emelyanov V.N., Karpenko A.G., Volkov K.N. Development of advanced CFD tools and their application to simulation of internal turbulent flows // Proceedings of the 5th European Conference for Aeronautics and Space Sciences (EUCASS 2013), Munich, Germany. P. 15.

5. Волков K.IL, Емельянов В.H., Курова КВ., Карпенко А.Г. Современные вычислительные технологии в разработке и оптимизации технических устройств различного назначения // Труды XIII международного семинара «Супервычисления и математическое моделирование». Издательско-полиграфический комплекс ФГУП «РФЯЦ - ВНИИЭФ». г. Саров. 2011. С. 161-171.

6. Емельянов В.Н., Карпенко А.Г. Реализация алгоритма решения трехмерных уравнений нестационарной газовой динамики на графических процессо-

pax // Тезисы XVIII Школы-семинара молодых ученых и специалистов под руководством академика РАН А.И. Леонтьева: Проблемы газодинамики и тепломассообмена в новых энергетических технологиях. М: Издательский дом МЭИ. 2011. С. 56-58.

7. Жигалко Е.Ф., Карпенко А.Г., Цветков А. И. Деструкция вихря // Труды семинара «Компьютерные методы в механике сплошной среды». Издательство С.Петерб. ун-та. 2011. С. 3-6.

8. Карпенко А.Г. Решение задачи о распаде произвольного разрыва используя векторные алгоритмы на графических процессорах // Тезисы XXII Всероссийский семинар с международным участием по струйным, отрывным и нестационарным течениям. Балт. гос. техн. ун-т. СПб. 2010. С. 84-86.

9. Карпенко А.Г. Разработка алгоритма расчёта течений сжимаемого газа на неструктурированных сетках // Труды XVII Школы-Семинара молодых ученых и специалистов под руководством академика РАН А.И.Леонтьева: Проблемы газодинамики и теплообмена в аэрокосмических технологиях. М: Издательский дом МЭИ. 2009. Т. 1. С. 353-356.

10. Азаров М.А., Карпенко А.Г. Гидродинамика и теплообмен в импульсно-периодических лазерных системах // Сборник трудов Международной конференции по механике и баллистике «VI Окуневские чтения». Балт. гос. техн. ун-т. 2008. Т. 3. С. 79-84.

11. Карпенко А.Г. Математическое моделирование нестационарных газодинамических процессов и теплопереноса в рабочих полостях импульсно-периодических систем // Труды XVI Школы-семинара молодых ученых и специалистов под руководством академика РАН А.И.Леонтьева: Проблемы газодинамики и теплообмена в энергетических установках. М: Издательский дом МЭИ. 2007. Т.2. С. 29-32.

12. Азаров М.А., Емельянов В.Н„ Карпенко А.Г. Математическое моделирование нестационарных газодинамических процессов и теплопереноса в рабочих полостях импульсно-периодических систем // Труды III научной конференции "Проектирование Инженерных и научных приложений в среде Matlab". 2007. С. 5-24.

Текст научной работы диссертации и автореферата по механике, кандидата физико-математических наук, Карпенко, Антон Геннадьевич, Санкт-Петербург

САНКТ-ПЕТЕРБУРГСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ

КАРПЕНКО Антон Геннадьевич

ЧИСЛЕННОЕ РЕШЕНИЕ ЗАДАЧ ГИДРОАЭРОМЕХАНИКИ НА ГРАФИЧЕСКИХ

ПРОЦЕССОРАХ

Специальность 01.02.05 — Механика жидкости, газа и плазмы

Диссертация на соискание ученой степени кандидата физико-математических наук

Научный руководитель: доктор физико-математических наук, профессор

Матвеев С.К.

Санкт-Петербург - 2013

Оглавление

Введение......................................................................5

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

1.1 Основные отличия графических процессоров....................9

1.2 Устройство графических процессоров............................11

1.3 Модель программирования........................................14

1.4 Структура памяти..................................................16

1.5 Вычисления с различной точностью..............................19

1.6 Схема решения задачи и технология программирования С1ША 20

1.7 Области применения ..............................................22

1.8 Выводы и заключение по главе ..................................36

Глава 2. Использование графических процессоров для решения задач

гидроаэромеханики......................................................37

2.1 Система уравнений газовой динамики............................37

2.2 Численное решение задачи о распаде произвольного разрыва для системы уравнений одномерной нестационарной газовой динамики на графических процессорах..........................40

2.2.1 Постановка задачи и основные уравнения................40

2.2.2 Численная реализация метода контрольного объема . . 43

2.2.3 Реализация векторного алгоритма........................52

2.2.4 Расчет тестовых задач ....................................54

2.3 Расчет трехмерных сжимаемых течений газа на неструктурированных сетках с использованием графических процессоров 63

2.3.1 Численная схема и метод контрольного объема .... 63

2.3.2 Неструктурированная сетка и расчет ее геометрических параметров..................................................73

2.3.3 Особенности реализации векторного алгоритма .... 74

2.3.4 Решение тестовых задач ................. 77

2.3.5 Сравнение времени выполнения на ЦПУ и ГПУ .... 94 2.4 Выводы и заключение по главе ................. 96

Глава 3. Расчет течений несжимаемой жидкости и течений газа при

малых числах Маха........................... 98

3.1 Проблема расчета течений газа при малых числах Маха ... 98

3.2 Численная реализация метода предобусловливания......102

3.2.1 Вывод матрицы предобусловливания..........102

3.2.2 Метод контрольного объема и расчет потоков на грани 108

3.2.3 Схема для расчета стационарных течений методом установления ..........................110

3.2.4 Схема для расчета для не стационарных течений с использованием метода двойных шагов по времени . . . 111

3.3 Расчет стационарного квазиодномерного течения в канале при малых числах Маха........................116

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

3.3.2 Предобусловливание системы уравнений и численная схема............................117

3.3.3 Численная реализация граничных условий.......120

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

3.4 Расчет многомерных тестовых задач..............135

3.4.1 Расчет стационарного двумерного течения в канале при малых числах Маха....................135

3.4.2 Расчет свободно-конвективного течения около вертикально поставленной нагретой пластины........139

3.5 Выводы и заключение по главе .................143

Глава 4. Расчет свободной конвекции в цистерне с мазутом.......145

4.1 Обзор проблемы..........................145

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

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

4.4 Заключение и выводы по главе .................150

Заключение..................................151

Литература..................................152

Приложение А. Решение задачи о распаде произвольного разрыва методом последовательных приближений................163

Приложение Б. Предобусловливание системы уравнений газовой динамики ..................................172

Б.1 Вывод формул для собственных чисел предобусловленной системы уравнений..........................172

Б.2 Вывод матриц диагонализации якобиана потока предобуслов-

ленной системы уравнений....................173

Б.З Вид обратных матриц к матрицам предобусловливания ... 177

Введение

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

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

арифметико-логических устройств АЛУ, в отличии от ЦПУ у которых всего несколько АЛУ. Для ГПУ необходимо разрабатывать новые схемы расчета учитывающие особенности архитектуры. Необходимо разрабатывать подходы к описанию течения жидкости и газа имеющие параллелизм на уровне данных. Это позволит использовать весь потенциал ГПУ и существенно ускорить вычисления. Вновь разработанные схемы становятся конкурентоспособными со схемами разработанными для ЦПУ.

Для достижения поставленной цели необходимо было решить следующие задачи:

2. Разработать, реализовать и протестировать схемы расчета нестационарного вязкого сжимаемого течения газа на неструктурированных сетках в многомерном случае используя ГПУ;

4. Продемонстрировать возможности разработанного подхода на примере решения практической задачи теплогидродинамической конвекции ТГК при транспортировке застывающих наливных грузов ЗНГ. Основные результаты:

нестационарного вязкого сжимаемого течения газа на неструктурированных сетках в многомерном случае;

4. Разработанный комплекс использован при решении практических задач теплогидродинамической конвекции ТГК при транспортировки застывающих наливных грузов ЗНГ.

Научная новизна:

2. Разработаны схемы расчета нестационарного вязкого сжимаемого газа на неструктурированных сетках в многомерном случай с использованием ГПУ.

3. Разработан универсальный подход расчета течений газа при малых числах Маха и несжимаемых течений, а так же течений со свободной конвекцией. Подход основан на предобусловливании системы уравнений сжимаемого газа.

4. Впервые аналитически выведены матрицы диагонализации якобиана потока предобусловленной системы в многомерном случае. Разработка программного универсального программного комплекса

расчета задач гидроаэромеханики на ГПУ позволит повысить быстродей-

ствие и обеспечит оперативность вычислений. Результаты работы и разработанные схемы расчета были внедрены в пакет программ трехмерного имитационного моделирования на супер-ЭВМ ЛОГОС, являющегося отечественным аналогом АМБУЗ. Пакет программ инженерного анализа ЛОГОС предназначен для моделирования процессов тепломассопереноса (аэродинамика, газодинамика, гидродинамика, теплопроводность) и разрабатывается Российским федеральным ядерным центром — Всероссийский научно-исследовательский институт экспериментальной физики ФГУП «РФЯЦ-ВНИИЭФ» г. Саров.

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

1.1. Основные отличия графических процессоров

Графические процессоры представляют собой программируемые вычислительные устройства, предназначенные для обработки графической информации. Графический процессор занимается расчетами выводимого изображения, освобождая от этой обязанности ЦПУ, а также производит расчеты для обработки команд трехмерной графики (геометрическое преобразование объектов, моделирование освещения). Ряд процедур обработки изображений и видео с переносом на GPU стали работать в режиме реального времени. Графические процессоры используют концепцию одиночный поток команд и множественный поток данных ОПМД (Single Instruction Multiple Data, SIMD).

Графические процессоры общего назначения ГПУ (General Purpose GPU, GPGPU) обладают собственной динамической памятью ОЗУ (Dynamic Random Access Memory, DRAM) и содержат множество мультипроцессоров, которые управляют высокоскоростной памятью, что делает их использование эффективным как для графических, так и для неграфических вычислений (например, платформы NVIDIA Tesla и Fermi специально разрабатываются для вычислительных приложений).

Несмотря на то, что тактовые частоты ГПУ ниже, чем у ЦПУ, высокие показатели производительности ГПУ обеспечиваются благодаря большому количеству потоковых процессоров. По сравнению с традиционными высокопроизводительными архитектурами вычислительных систем (параллельные и векторные компьютеры), ГПУ и кластеры ГПУ обладают лучшими показателями энергоэффективности и производительности в расчете на сто-

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

Структурные отличия ЦПУ и ГПУ приводят к тому, что теоретическая производительность ГПУ, измеряемая количеством арифметических и логических операций в единицу времени, значительно превосходит теоретическую производительность ЦПУ (рис. 1.1). Значки • соответствуют вычислениям с двойной точностью на ЦПУ компании Intel, а значки о и □ — вычислениям с одинарной и двойной точностью на ГПУ компании NVIDIA.

Несмотря на то, что потенциальные возможности ГПУ при решении широкого круга прикладных задач не подлежат сомнению, технологические вопросы реализации вычислительных задач на ГПУ общего назначения требуют дальнейшего развития (адаптация существующих программ к использованию вычислительных ресурсов ГПУ, исследование производительности на определенных классах задач, разработка средств и систем программирования, оптимизация кода). Основной проблемой, препятствующей массовому внедрению ГПУ в вычислительную практику, является отсутствие высокоуровневых средств программирования. Появление нескольких технологий доступа к ГПУ и интерфейсов к ним (например, CUDA, OpenCL) ставит вопрос о сравнительной эффективности и особенностях реализации вычислительных задач с использованием различных сред и технологий. Распространение кластеров ГПУ приводит к необходимости выбора технологии организации кластера [47,96]. Немаловажным представляется также вопрос, связанный с оптимизаций программного кода, предназначенного для выполнения на ГПУ [1].

Рис. 1.1. Пиковая производительность ЦПУ и ГПУ

В данной работе обсуждаются применение технологии С1ЮА для решения задач, связанных с моделированием течений жидкости и газа. Рассматриваются особенности реализации программного кода на ГПУ и ряд вопросов, связанных с его оптимизацией. Сравнивается ускорение решения задачи на ГПУ по сравнению с расчетами на ЦПУ при использовании сеток различной разрешающей способности.

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

Упрощенную организацию, а также устройство ЦПУ и ГПУ поясняют схемы, приведенные на рис. 1.2.

□I II

□l I I I I I

□ii i ii__

MINN-

0_____

01 I I I I I-

0 I I I I

1_II

0_____

01 I | I I I-

□i i i i i__

ОЗУ

6) ГПУ

Рис. 1.2. Условные схемы ЦПУ (а) и ГПУ (б)

У современных ЦПУ (рис. 1.2 а) имеется небольшое количество арифметико-логических устройств АЛУ (Arithmetic Logical Unit, ALU), контроллер управления, отвечающий за передачу следующей машинной инструкции АЛУ и ряд других функций, контроллер доступа к внешней памяти (не показан на рисунке) и внешняя ОЗУ-память. Каждое АЛУ содержит набор регистров и свой сопроцессор для расчета сложных функций. В кэш-память, обычно располагающейся на одном кристалле с АЛУ, загружаются те данные, к которым в исполняемой программе происходит обращение несколько раз. При первом обращении к данным время считывания является таким же, как и при обычном обращении к внешней памяти. Данные считываются в кэш-память, а потом поступают в АЛУ. При последующих обращениях данные уже считываются из кэш-памяти, что происходит намного быстрее, чем обращение к внешней памяти.

Структура и организация памяти ГПУ является более сложной, чем ЦПУ В ГПУ находится большое количество АЛУ, несколько контроллеров управления и внешняя память (рис. 1.2б).

Каждый ГПУ состоит из потокового мультипроцессора ПМ (Streaming Multiprocessor, SM), содержащего несколько скалярных процессоров СП (Scalar Processor, SP). Число скалярных процессоров зависит от типа ви-

Модуль управления АЛУ АЛУ

АЛУ АЛУ

КЭШ память

ОЗУ

а) ЦПУ

деокарты и изменяется в широких пределах. Потоковый мультипроцессор имеет ряд регистров и общую память, используемую для коммуникаций между скалярными процессорами (рис. 1.3). Система команд скалярных процессоров включает арифметические команды для вещественных и целочисленных вычислений с 32-разрядной точностью, команды управления (ветвления и циклы), а также команды обращения к памяти. Команды доступа к оперативной памяти выполняются асинхронно. Для уменьшения задержек в очереди выполнения ГПУ переключение осуществляется за один такт. За переключение потоков отвечает диспетчер потоков, который не является программируемым.

Рис. 1.3. Структура и состав ГПУ

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