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

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

Оглавление диссертации кандидат наук Карпенко, Антон Геннадьевич

Оглавление

Введение

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.3.5 Сравнение времени выполнения на ЦПУ и ГПУ

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Заключение

Литература

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

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

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

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

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

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

Введение диссертации (часть автореферата) на тему «Численное решение задач гидроаэромеханики на графических процессорах»

Введение

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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. Структура и состав ГПУ

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

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

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

Выполнение программы, использующей ресурсы ЦПУ и ГПУ, поясняет рис. 1.4. Последовательная часть приложения (serial code) работает на ЦПУ, а часть приложения, на которую приходится наибольшая вычислительная нагрузка (kernel), выполняется на ГПУ. Скалярные процессоры в потоковом мультипроцессоре синхронно исполняют одну и ту же команду, следуя модели исполнения ОПМД. Мультипроцессоры не являются синхронизированными, в связи с чем общая схема исполнения следует модели одиночный поток команд и множественный поток нитей ОПМН (Single Instruction Multiple Threads, SIMT ).

ЦПУ

Послед код

Ядро 1

Послед код

Ядро 1

ГПУ

Сетка 1

Блок

(0,0)

Блок '

(0,У

Блок (1.0)

Блок (1.1)

Сетга

Блок (2.0)

Блок

(2.1)

\ \

I \

1 I I 11 |'г

Нить (0.0) Нить (1,0) Нить (2,0) Нить (3,0) Нить (4,0)

Нить (0,1) Нить (1.1) Нить (2Д) Нить (4,1) Нить (4,1)

Нить (0,2) Нить (1.2) Нить (2,2) Нить (3.2) Нить (4.2)

Рис. 1.4. Выполнение программы, использующей ресурсы ЦПУ и ГПУ

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

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

Иерархия выполнения программы на ГПУ, производимых компанией NVIDIA, делится на сетки блоков, блоки нитей и нити (рис. 1.5), что связывается со структурой и организацией памяти. До выполнения ядра размеры блоков нитей и сетки блоков задаются программистом в явном виде.

Рис. 1.5. Иерархия исполнения программы на ГПУ

Верхним уровнем иерархии является сетка блоков (grid), которая соответствует всем нитям, выполняющим данное ядро.

Сетка состоит из одномерного или двумерного массива блоков (block). Все блоки, образующие сетку, имеют одинаковую размерность и размер. Каждый блок в сетке имеет уникальный адрес, состоящий из одного или двух неотрицательных целых чисел (индекс блока в сетке). Во время исполнения блоки являются неделимыми и укладываются в один потоковый мультипроцессор. Блоки выполняются асинхронно и не имеется механизма, обеспечивающего синхронизацию между ними. Группировка блоков в сетки помогает при масштабировании программного кода (блоки выполняются параллельно, а при ограниченных ресурсах — последовательно).

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

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

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

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

В отличие от ЦПУ, где имеется один тип памяти с несколькими уровнями кэширования в самом процессоре, ГПУ обладает более сложной организацией памяти [1,9], связанной с их назначением (обработка графической информации). Часть памяти располагается непосредственно в каждом из потоковых мультипроцессоров (регистровая и разделяемая памяти), а часть памяти размещается в ОЗУ (локальная, глобальная, константная и текстурная памяти).

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

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

Хост

Блок(0,0) Блок(1,0)

Разделяемая память Разделяемая память

Регистры Реестры Регистры Регистры

I I 1 t

Нить(0 0) Нить(1 0) Нитцад) Нитьи о)

1 -4— —i-1—

Глобальная память

Константная память

Рис. 1.6. Структура памяти ГПУ

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

Глобальная память (global memory) обладает высокой пропускной способностью (более 100 Гб/с для некоторых типов процессоров) и используется для обменов данными между ЦПУ и ГПУ. Размер глобальной памяти составляет от 256 Мб до 1.5 Гб в зависимости от типа процессора и доходит до 6 Гб на платформе NVIDIA Tesla.

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

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

ной памяти в один (coalescing), что позволяет получить почти 16-кратное ускорение (обращения к глобальной памяти происходят посредством 32 или 128-битовых слов, а адрес, по которому происходит доступ, выравнивается по размеру слова).

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

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

Разделяемая память (shared memory) размещается в самом мультипроцессоре и выделяется на уровне блоков (каждый блок получает в свое распоряжение одно и тоже количество разделяемой памяти), являясь доступной на чтение и запись всем нитям блока. Разделяемая память объемом от 16 Кб до 48 Кб поровну делится между всеми блоками сетки, исполняемыми на мультипроцессоре.

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

При надлежащем использовании разделяемой памяти дополнительное ускорение счета составляет от 18 до 25% (в зависимости от технологии и методов оптимизации), а при использовании константной и текстурной памяти ускорение счета достигает 35-40% [3].

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

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

Наибольшая производительность достигается при вычислениях с одинарной точностью (single precision). При использовании вычислений с двойной точностью (double precision) производительность снижается теоретически на порядок (по сравнению с пиковой производительностью, заявленной производителем) и примерно в два раза на практике [3]. Снижение быстродействия вычислений на ГПУ существенно больше, чем на ЦПУ, что сказывается в расчетах, требующих повышенной точности.

Расчеты, проведенные в работе [39], показывают, что вычисления с двойной точностью оказываются примерно на 46-66% (в зависимости от размера сетки) медленнее вычислений с одинарной точностью. При этом пиковая производительность вычислений с двойной точностью (78 ГФлопс для видеокарты NVIDIA GT200) в 12 раз меньше, чем производительность вычислений с одинарной точностью (936 ГФлопс для видеокарты NVIDIA GT200).

Потенциал использования вычислений с различной точностью (mixed precision) обсуждается в работах [23,56]. Одинарная точность используется для расчета левой части дискретного уравнения, связанной с дискретизацией по времени, а двойная точность — для расчета правой части, представляющей собой невязку, обусловленную дискретизацией невязких и вязких потоков [56] (в последовательной части кода, запускаемой на ЦПУ, исполь-

зуются вычисления с двойной точностью). Такой подход позволяет получить выигрыш в производительности вычислений без потери их точности.

1.6. Схема решения задачи и технология программирования СИБА

Схему решения задачи на ГПУ поясняет рис. 1.7 (шаг 6). Стрелки —> соответствуют командам, а стрелки => — командам и обменам данными. Подготовка исходных данных (чтение с жесткого диска, формирование и инициализация массивов, сортировка) осуществляется на ЦПУ. Массивы данных копируются на ГПУ (пропускная способность шины, через которую производится копирование данных, составляет около 5 Гб/с). Интенсивные вычисления производятся на ГПУ. На ГПУ над массивами данных производится заданная последовательность действий. Каждая нить обрабатывает один узел сетки. Результирующие массивы копируются с ГПУ на ЦПУ, после чего производится обработка и визуализация данных, а также их сохранение на жесткий диск.

ЦПУ

Считать

параметры *

Выделить память на ЦПУ

Отправить данные на ГПУ

Получить данные с ГПУ - | | Отправить данные на ЦПУ

~ 1 1

* 1 Записать 1 | > данные в файл | ■ 1

Освободить | 1 Освободить память на ГПУ

1амять на ЦПУ | .

Выделить память на ГПУ

Получить данные с ЦПУ

Рис. 1.7. Схема решения задачи при использовании ресурсов ГПУ

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

Для поддержки неграфических приложений на ГПУ компании NVIDIA используется унифицированная архитектура компьютерных вычислений CUDA, основанная на расширении языка С и позволяющая получить доступ к набору инструкций ГПУ для управления его памятью [1]. При этом ГПУ рассматривается как вычислительное устройство (device), способное поддерживать параллельное исполнение большого числа нитей или потоковых программ (ГПУ является сопроцессором ЦПУ хост-компьютера, host). Реализуется оптимизированный обмен данными между ЦПУ и ГПУ. Поддерживаются как 32-битовые, так и 64-битовые операционные системы. В язык программирования добавляются спецификаторы типа, встроенные переменные и типы, директивы запуска ядра, а также другие функции.

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

Список литературы диссертационного исследования кандидат наук Карпенко, Антон Геннадьевич, 2013 год

Литература

1. Боресков A.B., Харламов А. Основы работы с технологией CUDA. М.: Изд-во ДМК-Пресс, 2010. 232 с.

2. Боярченков A.C., Поташников С.И. Использование графических процессоров и технологий CUDA для задач молекулярной динамики // Вычислительные методы и программирование. 2009. Т. 10. № 1. С. 9-23.

3. Галимов М.Р., Биряльцев Е.В. Некоторые технологические аспекты применения высокопроизводительных вычислений на графических процессорах в прикладных программных системах // Вычислительные методы и программирование. 2010. Т. 11. № 1. С. 77-93.

4. Геллер О.В., Васильев М.О., Холодов Я.А. Построение высокопроизводительного вычислительного комплекса для моделирования задач газовой динамики // Компьютерные исследования и моделирование. 2010. Т. 2. № 3. С. 309-317.

5. Годунов С.К., Забродин A.B., Иванов М.Я., Крайко А.Н., Прокопов Г.П. Численное решение многомерных задач газовой динамики. М.: "Наука 1976. 400 с.

6. Горобец A.B., Суков С.А., Железняков А.О., Богданов П.Б., Четве-рушкин Б.Н. Применение GPU в рамках гибридного двухуровневого распараллеливания MPI+OpenMP на гетерогенных вычислительных системах // Параллельные вычислительные технологии. Челябинск: Издательский центр ЮУрГУ, 2011. С. 452-460.

7. Губайдуллин Д.А., Никифоров А.И., Садовников Р.В. Библиотека шаблонов итерационных методов подпространств Крылова для численного решения задач механики сплошных сред на гибридной вычислительной системе // Вычислительные методы и программирование. 2010. Т. 11. № 1. С. 351-359.

8. Евстигнеев Н.М. Интегрирование уравнения Пуассона с использованием графического процессора технологии CUDA // Вычислительные методы и программирование. 2009. Т. 10. № 1. С. 268-274.

9. Казенное A.M. Основы технологии CUDA // Компьютерные исследования и моделирование. 2010. Т. 2. № 3. С. 295-308.

10. Клосс Ю.Ю., Черемисин Ф.Г., Шувалов П.В. Решение уравнения Больцмана на графических процессорах; // Вычислительные методы и программирование. 2010. Т. 11. № 1. С. 144-152.

11. Колган В.П. Применение принципа минимальных значений производных к построению конечноразностных схем разрывных решений газовой динамики // Уч. зап. ЦАГИ. 1972. Т. 3. № 6. С. 68-77.

12. Лебедев A.C., Чёрный С.Г. Практикум по численному решению уравнений в частных производных. // Учебное пособие. - Новосибирск: Новосиб. гос. ун-т., 2000.

13. Лойцянский Л.Г. Механика жидкости и газа // Государственное издательство технико-теоретически литературы, 1950. 678 с.

14. Моисеев В.И., Денисов Г.И. Подавление естественной конвекции, как способ сохранения текучести горячих вязких жидкостей при перевозках. // Восьмой сборник академии: Актуальные проблемы обеспечения войсковых Сил в транспортном отношении. С.Пб.: «ВАТТ», 2009., С.124-133.

15. Рыбакин Б.П., Шидер Н.И. Построение параллельных алгоритмов для решения задач гравитационной газовой динамики // Вычислительные методы и программирование. 2010. Т. 11. № 1. С. 388-394.

16. Соркин Р. Е. Газотермодинамика ракетных двигателей на твёрдом топливе.// Наука; Глав. ред. физико-математической лит-ры, 1967, 368 с.

17. Спенс А., Скотт Н., Гиллан Ч. Новые технологии высокопроизводительных вычислений: численное интегрирование на графических процессорах // Вычислительные методы и программирование. 2009. Т. 10. № 1. С. 385-388.

18. Чадов С.Н. Реализация алгоритма решения несимметричных систем линейных уравнений на графических процессорах // Вычислительные методы и программирование. 2009. Т. 10. № 2. С. 135-140.

19. Шлихтинг Г. Теория пограничного слоя // М.: Наука, 1974. 712 с.

20. Antoniou A.S., Karantasis K.I., Polychronopoulos E.D., Ekaterinaris J. A. Acceleration of a finite-difference. WEN О scheme for large-scale simulations on many-core architectures // AIAA Paper. 2010. No. 20100525.

21. Appleyard J., Drikakis D. Higher-order CFD and interface tracking methods on highly-parallel MPI and GPU systems // Computers and Fluids. 2011. Vol. 46. No. 1. P. 101-105.

22. Asuncion M., Mantas J.M., Castro M.J., Fernandez-Nieto E.D. An MPI-CUDA implementation of an improved Roe method for two-layer shallow water systems // Journal of Parallel and Distributed Computing. 2012.

23. Baboulin M., Buttari A., Dongarra J., Kurzak J., Langou J., Langou J., Luszczek P., Tomov S. Accelerating scientific computations with mixed precision algorithm // Computer Physics Communications. 2009. Vol. 180. No. 12. P. 2526-2533.

24. Bailey P., Myre J., Walsh S.D.C., Lilja D.J., Saar M.O. Accelerating Lattice Boltzmann fluid flow simulations using graphics processors // Proceedings of the 38th International Conference on Parallel Processing, 22-25 September 2009, Vienna, Austria. 2009. P. 550-557.

25. Bell N., Garland M. Efficient sparse matrix vector multiplication on CUDA // NVIDIA Technical Report. 2008. No. NVR-2008-004.

26. Ben-Dor G. Shoek Wave Reflection Phenomena // Springer-Verlag, 1992.

27. Bernardon F.F., Callahan S.P., Comba J.L.D., Silva C.T. An adaptive framework for visualizing unstructured grids with time-varying scalar fields // Parallel Computing. 2007. Vol. 33. No. 6. P. 391-405.

28. Brandvik T., Pullan G. Acceleration of a 3D Euler solver using commodity graphics hardware // AIAA Paper. 2008. No. 2008-607.

29. Brandvik T., Pullan G. An accelerated 3D Navier-Stokes solver for flows in turbomachines // AS ME Paper. 2009. No. GT2009-60052.

30. Buatois L., Caumon G., Levy B. Concurrent number cruncher: an efficient sparse linear solver on the GPU // High Performance Computing and Communications. 2007. Vol. 4782. P. 358-371.

31. Buluc A., Gilbert J.R., Budak C. Solving path problems on the GPU // Parallel Computing. 2010. Vol. 36. No. 5-6. P. 241-253.

32. Castro M.J., Ortega S., de la Asuncion M., Mantas J.M., Gallardo J.M. GPU computing for shallow water flow simulation based on finite volume schemes // Comptes Rendus Mecanique. 2011. Vol. 339. No. 2-3. P. 165184.

33. Chambon G., Bouvarel R., Laigle D., Naaim M. Numerical simulations of granular free-surface flows using smoothed particle hydrodynamics // Journal of Non-Newtonian Fluid Mechanics. 2011. Vol. 166. No. 12-13. P. 698-712.

34. Chaudhary R., Vanka S.P., Thomas B.G. Direct numerical simulations of magnetic field effects on turbulent flow in a square duct // Physics of Fluids. 2010. Vol. 22. No. 7. 075102 (15 pages).

35. Che S., Boyer M., Meng J., Tarjan D., Sheaffer J.W., Skadron K. A performance study of general-purpose applications on graphics processors using CUDA // Journal of Parallel and Distributed Computing. 2008. Vol. 68. No. 10. P. 1370-1380.

36. D. Choi, C. L. Merkle Application of time-iterative schemes to incompressible flow // AIAA Journal, 1985, Vol. 23, No. 10, pp. 15181524.

37. Y.H. Choi, C.L. Merkle The Application of Preconditioning in ViscousFlows // Journal of Computational Physics, 1993,Vol. 105, Issue 2, P. 207-223

38. Chorin A.J. A numerical method for solving incompressible viscous flow problems // Journal of Computational Physics. 1967. Vol. 2. No. 1. P. 12-26.

39. Cohen J.M., Molemaker M.J. A fast double precision CFD code using CUDA // Proceedings of Parallel CFD Conference, 18-22 May 2009, Moffet Field, California, USA. 2009. 17 p.

40. Corrigan A., Camelli F., Lohner R., Wallin J. Running unstructured grid-based CFD solvers on modern graphics hardware // AIAA Paper. 2009. No. 2009-4001.

41. Corrigan A., Camelli F., Lohner R., Mut F. Semi-automatic porting of a large-scale Fortran CFD code to GPUs // International Journal for Numerical Methods in Fluids. 2011.

42. Dick C., Georgii J., Westermann R. A real-time multigrid finite hexahedra method for elasticity simulation using CUDA // Simulation Modelling Practice and Theory. 2011. Vol. 19. No. 2. P. 801-816.

43. Elsen E., LeGresley P., Darve E. Large calculation of the flow over a hypersonic vehicle using a GPU // Journal of Computational Physics. 2008. Vol. 227. No. 24. P. 10148-10161.

44. Emery A.E. An evaluation of several differencing methods for inviscid fluid flow problems// Journal of Computational Physics, Vol. 2 (1968). 306-331.

45. Fu W.-S., Wang W.-H., Huang Y., Li C.-G. An investigation of compressible forced convection in a three dimensional tapered chimney by

CUDA computation platform // International Journal of Heat and Mass Transfer. 2011. Vol. 54. No. 15-16. P. 3420-3430.

46. Galiano V., Migallon H., Migall'on V., Penad'es J. GPU-based parallel algorithms for sparse nonlinear systems // Journal of Parallel and Distributed Computing. 2011.

47. Goddeke D., Strzodka R., Mohd-Yusof J., McCormick P., Wobker H., Becker C., Turek S. Using GPUs to improve multigrid solver performance on a cluster // International Journal of Computational Science and Engineering. 2008. Vol. 4. No. 1. P. 36-55.

48. Goodnight N., Woolley C., Lewin G., Luebke D., Humphreys G. A multigrid solver for boundary value problems using programmable graphics hardware // Proceedings of the Graphics Hardware, 26-27 July 2033, San Diego, USA. 2003. 11 p.

49. Griebel M., Zaspel P. A multi-GPU accelerated solver for the three-dimensional two-phase incompressible Navier-Stokes equations // Computer Science — Research and Development. 2010. Vol. 25. No. 1. P. 65-73.

50. N. Hakimi Preconditioning methods for time dependent navier-stokes equations application to environmental and low speed flows / / Ph.D.Thesis, 1997

51. Harten A. High resolution schemes for hyperbolic conservation laws. Journal of Computational Physics. 1983. Vol. 49. No. 3. P. 357-393.

52. Hirsch C. Numerical Computation of Internal and External Flows // Wiley, 1994

53. Hofler M. Real-time visualization of unstructured volumetric CFD data sets on GPUs // Proceedings of the 10th Central European Seminar on Computer Graphics, 23-26 April 2006, Vienna, Austria. 2006. 8 p.

54. Jacobsen D., Thibault J., Senocak I. An MPI-CUDA implementation for massively parallel incompressible flow computations on multi-GPU clusters // AIAA Paper. 2010. No. 2010-522.

55. Jespersen D.C. Acceleration of a CFD code with a GPU // NAS Technical Reeport. 2009. No. NAS-09-003. 10 p.

56. Kampolis I.C., Trompoukis X.S., Asouti V.G., Giannakoglou K.C. CFD-based analysis and two-level aerodynamic optimization on graphics processing units // Computer Methods in Applied Mechanics and Engineering. 2010. Vol. 199. No. 9-12. P. 712-722.

57. Klôckner A., Warburton T., Bridge J., Hesthaven J.S. Nodal discontinuous Galerkin methods on graphics processors // Journal of Computational Physics. 2009. Vol. 228. No. 21. P. 7863-7882.

58. Komatsu K., Soga T., Egawa R., Takizawa H., Kobayashi H., Takahashi S., Sasaki D., Nakahashi K. Parallel processing of the building-cube method on a GPU platform // Computers and Fluids. 2011. Vol. 45. No. 1. P. 122-128.

59. Kruger J., Westermann R. Linear algebra operators for GPU implementation of numerical algorithms // ACM Transactions on Graphics. 2003. Vol. 22. No. 3. P. 908-916.

60. Kuo F.-A., Smith M.R., Hsieh C.-W., Chou C.-Y., Wu J.-S. GPU acceleration for general conservation equations and its application to several engineering problems // Computers and Fluids. 2011. Vol. 45. No. 1. P. 147-154.

61. Kuznik F., Obrecht C., Rusaouen G., Roux J.-J. LBM based flow simulation using GPU computing processor // Computers and Mathematics with Applications. 2010. Vol. 59. No. 7. P. 2380-2392.

62. Lastra M., Mantas J.M., Urena C., Castro M.J., Garcia-Rodriguez J.A. Simulation of shallow-water systems using graphics processing units //

Mathematics and Computers in Simulation. 2009. Vol. 80. No. 3. P. 598618.

63. Majda. A., Osher S. Numerical viscosity and entropy condition // Communication on Pure and Applied Mathematics. 1979. 32. 797-838.

64. Mavriplis D.J. Unstructured Mesh Discretizations and Solvers for Computational Aerodynamics // 18th AIAA Computational Fluid Dynamics Conference, 2007, Miami

65. Mehra R., Raghuvanshi N., Savioja L., Lin M.C., Manocha D. An efficient GPU-based time domain solver for the acoustic wave equation // Applied Acoustics. 2012. Vol. 73. No. 1. P. 83-94.

66. Michalakes J., Vachharagani M. GPU acceleration of numerical weather prediction // Parallel Processing Letters. 2008. Vol. 18. No. 4. P. 531-548.

67. Micikevicius P. 3d finite difference computation on GPUs using CUDA // Proceedings of 2nd Workshop on General Purpose Processing on Graphics Processing Units (GPGPU-2), 8 March 2009, Washington, DC, USA. New York: ACM, 2009. P.. 79-84.

68. Obrecht C., Kuznik F., Tourancheau B., Roux J.-J. Multi-GPU implementation of the lattice Boltzmann method // Computers and Mathematics with Applications. 2011. Vol. 62.

69. Peng L., Nomura K., Oyakawa T., Kalia R., Nakano A., Vashishta P. Parallel Lattice Boltzmann flow simulation on emerging multi-core platforms // Lecture Notes in Computer Science. 2008. Vol. 5168. P. 763777.

70. Phillips E.H., Zhang Y., Davis R.L., Owens J.D. Rapid aerodynamic performance prediction on a cluster of graphics processing units // AIAA Paper. 2009. No. 2009-565.

71. Richter C., Abdel-Hay J., Panek L., Schonwald N., Busse S., Thiele F. Time domain impedance modelling and applications // Procedia IUTAM. 2010. Vol. 1. P. 133-142.

72. Roe, P. L. Approximate Riemann solvers, parameter vectors and difference schemes // Journal of Computational Physics. 1981. Vol. 43. No. 2. P. 357372

73. Rossinelli D., Bergdorf M., Cottet G.-H., Koumoutsakos P. GPU accelerated simulations of bluff body flows using vortex particle methods // Journal of Computational Physics. 2010. Vol. 229. No. 9. P. 3316-3333.

74. Scheidegger C.E., Comba J.L.D., da Cunha R.D. Practical CFD simulations on programmable graphics hardware using SMAC // Computer Graphics Forum. 2005. Vol. 24. No. 4. P. 715-728.

75. Schenk O., Christen M., Burkhart H. Algorithmic performance studies on graphics processing units // Journal of Parallel and Distributed Computing. 2008. Vol. 68. No. 10. P. 1360-1369.

76. Schive H.-Y., Tsai Y.-C., Chiueh T. Gamer: a graphic processing unit accelerated adaptive-mesh-refinement code for astrophysics // The Astrophysical Journal Supplement Series. 2010. Vol. 186. No. 2. P. 457484.

77. Senocak I., Thibault J., Caylor M. Rapid-response urban CFD simulations using a GPU computing paradigm on desktop supercomputers // Proceedings of the 8th Symposium on the Urban Environment, 10-15 January 2009, Phoenix, Arizona, USA. 2009. No. J19.2. 7 p.

78. Shams R., Sadeghi P. On optimization of finite-difference time-domain (FDTD) computation on heterogeneous and GPU clusters // Journal of Parallel and Distributed Computing. 2011. Vol. 71. No. 4. P. 584-593.

79. Shi Y., Green W.H., Wong H.-W., Oluwole O.O. Redesigning combustion modeling algorithms for the graphics processing unit (GPU): chemical kinetic rate evaluation and ordinary differential equation integration // Combustion and Flame. 2011. Vol. 158. No. 5. P. 836-847.

80. Shinn A.F., Vanka S.P. Implementation of a semi-implicit pressure-based multigrid fluid flow algorithm on a graphics processing unit //

Proceedings of the ASME International Mechanical Engineering Congress and Exposition (IMECE-2009), 13-19 November 2009, Lake Buena Vista, Florida, USA. No. IMECE2009-11587.

81. Shinn A.F., Vanka S.P., Hwu W.W. Direct numerical simulation of turbulent flow in a square duct using a graphics processing unit (GPU) // AIAA Paper. 2010. No. 2010-5029.

82. Stock M.J., Gharakhani A. A GPU-accelerated boundary element method and vortex particle method // AIAA Paper. 2010. No. 2010-5099.

83. Strzodka R., Doggett M., Kolb A. Scientific computation for simulations on programmable graphics hardware // Simulation Modelling Practice and Theory. 2005. Vol. 13. No. 8. P. 667-680.

84. Tanno I., Morinishi K., Satofuka N., Watanabe Y. Calculation by artificial compressibility method and virtual flux method on GPU // Computers and Fluids. 2011. Vol. 45. No. 1. P. 162-167.

85. Thibault J.C., Senocak I. CUDA implementation of a Navier-Stokes solver on multi-GPU desktop platforms for incompressible flows // AIAA Paper. 2009. No. 2009-758.

86. Tolke J., Krafczyk M. TeraFLOP computing on a desktop PC with GPUs for 3D CFD // International Journal of Computational Fluid Dynamics. 2008. Vol. 22. No. 7. P. 443-456.

87. Toro Eleuterio F. Rieman Solvers and Numerical Methods for Fluid Dynamics: A Prectial Introduction - 2.ed.Springer, 1999. 606 p.

88. E. Turkel Preconditioned Methods for Solving the Incompressi and Low Speed Compressible Equations // Journal of Computational Physics, 1987, Vol. 72, p. 277-298

89. E. Turkel, A. Fiterman, B. Van Leer, Preconditioning and the Limit to the Incompressible Flow equations // NASA contractor report ; NASA CR-191500, 1993

90. Vanka S.P., Shinn A.F., Sahu K.C. Computational fluid dynamics using graphics processing units: challenges and opportunities // Proceedings of the AS ME International Mechanical Engineering Congress and Exposition (IMECE-2011), 11-17 November 2011, Denver, Colorado, USA. 2011. No. IMECE2011-65260.

91. Weiss, J.M., and Smith, W.A. Preconditioning Applied to Variable and Constant Density Flows //AIAA J.,1995, Vol. 33, No. 11., pp. 2050-2057.

92. Wesseling P. Principles of computational fluid dynamics // Springer, 2000. 664 p.

93. Woodward P., Colella P. The numerical simulation of two-dimensional fluid flow with strong shocks // Journal of Computational Physics. 1984. V.54, PP.115-173.

94. Wong H.-C., Wong U.-H., Feng X., Tang Z. Efficient magnetohydrodynamic simulations on graphics processing units with CUDA // Computer Physics Communications. 2011. Vol. 182. No. 10. P. 2132-2160.

95. Xian W., Takayuki A. Multi-GPU performance of incompressible flow computation by lattice Boltzmann method on GPU cluster // Parallel Computing. 2011. Vol. 37. No. 9. P. 521-535.

96. Yokota R., Narumi T., Sakamaki R., Kameoka S., Obi S., Yasuoka K. Fast multipole methods on a cluster of GPUs for the meshless simulation of turbulence // Computer Physics Communications. 2009. Vol. 180. No. 11. P. 2066-2078.

97. Yokota R., Barba L.A. Comparing the treecode with FMM on GPUs for vortex particle simulations of a leapfrogging vortex ring // Computers and Fluids. 2011. Vol. 45. No. 1. P. 155-161.

98. Zuo W., Chen Q. Fast and informative flow simulations in a building by using fast fluid dynamics model on graphics processing unit // Building and Environment. 2010. Vol. 45. No. 3. P. 747-757.

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

Годунов в своей работе [5] показал, что давление на контактом разрыве р* можно найти, решив относительного него нелинейное уравнение:

/(р, Vl, Vr) = fL(р, VL) + fR(p, Vñ) + Au = 0, Au = uR - uL, (A.l)

где функция fk имеет следующий вид (здесь и далее соотношения взяты из работы Topo [87]):

fkip, У к) = <

(р ~ Рк)

2 Ск (7-1)

-Л

27

, если р > Рк (ударная волна),

, если р ^ Рк (волна разрежения),

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

2

Ак =

в,=.

(7 + 1W (7+1)'

Скорость контактного разрыва п* находится из соотношения:

U* = + uR) + ^[/я(р*) - /l(p*)] •

Таким образом, решив нелинейное алгебраическое уравнение (А.1) относительно давления на контактом разрыве р* , мы можем вычислить и скорость контактного разрыва и*.

Вычислим первую и вторую производную функции /(р) = ¡ь{р) + /д(р) + Au по переменной р:

fk —

\Вк+р,

p—pk 2 (Вк+р)

-(7+i) 27

, если р > рк (ударная волна),

27 ' если^ ^ Рк (волна разрежения).

/Г =

_1 f^iLj

4 \ Bfc+W

(7+l)cfc / v2,n2 \Pfc у

4Bfc+3p+pfc

-(37+1)

27 2PÍ

, если р > рк (ударная волна), , еслир ^ рк (волна разрежения). 163

р

Рис. А.1. График функции /(р)

Видно что при всех р > О первая производная /' = /ь + > 0, а вторая производная / = + /д < 0. Поэтому /(р) является монотонно возрастающей выпуклой функцией аргумента р.

Значение функции /(р) в точках ртЫ = тт(рь,ря) и ртах = тах(р1:рц) позволяет еще до решения уравнения (А.1) определить какая конфигурация возникает при распаде разрыва. Для заданных р/, и рд скачок скорости А и определяет положения корня уравнения р* см. рис. А.1. Так как Аи входит линейно в уравнение (А.1), то для различных Аи график функции /(р) будет смещаться либо вверх, либо вниз. Поэтому существует четыре варианта решения :

1. р* лежит в интервале Ь = (0,рт{п), если }{ртт) > 0 и /(ртах) > 0. Этот случай характеризуется существенно большими значениями Ащ . Так как р* < ртгП < Ртах в обе стороны расходятся волны разрежения.

2. р* лежит в интервале /2 = \ртт,ртах], если /(рт{п) ^ 0 и /(Ртах) ^ о. В этом случае , так как ргтп ^ р* ^ ртах в одну сторону распространяется волна разрежения, а в другую ударная волна. Если Рь = Ртт, то ударная волна будет распространятся влево, а при Ря = Ртт наоборот.

3. р* лежит в интервале /3 = (;ртах, оо), если /(рто»„) < 0 и ¡(ртах) < 0. В этом случае, так как ртгП < ртах < р* в обе стороны расходятся ударные волны.

4. В случае когда значение Аи превысит критическое значение (Аи)^ решение уравнения (А.1) не имеет вещественного корня. В этом случае образуется зона вакуума.

Критическое значение перепада скорости на разрыве можно определить, положив в уравнении (А.1) р = 0:

/д ч 2аь 2ак

{Аи)крит =-- Н---.

7—1 7—1

Поэтому условие, при котором существует действительное решение р* и не образуется зона вакуума имеет следующий вид:

2аь 2ад

Н--7 > иЯ - иь-

7—1 7—1

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

/(р»-0

Рг = Рг-1 - 777-

/ \Рг-1)

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

Ьг ~Рг-11

\\Рг "Ь Рг—1

165

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

Для нахождения первого приближения Topo в работе [87] использовал подход основанный на аналитических и приближенных решениях.

Запишем систему уравнений (2.9) в физических переменных в квазилинейной форме:

dV dV

Сложность решения этого уравнения в том, что матрица якобиана a(V) зависит от решения V. Если заменить матрицу якобиана на мало отличающуюся матрицу зависящую от осредненных параметров а, то систему уравнений (А.2) можно решить аналитически, так как она преобразуется в линейную гиперболическую систему дифференциальных уравнений с постоянными коэффициентами. Давление на контактном разрыве в этом случае записывается следующим образом:

Р* = ^(Pl + Pr) + - uR)(pc),

где рис осредненые значения плотности и скорости звука:

1 1

Р = -j(PL + Pr), с = ~{cl + cR).

Используя это приближенное значение = р* можно определить вариант картины течения и использовать то или иное приближение.

Если Ppyrs < Pmin то образуются две волны разрежения. В таком случае имеется точное аналитическое выражение для давления на контактом разрыве:

CL + CR- ^(UR ~ Ub)

Р* =

£L. _i_ £К pí^ pZR

где комплекс z имеет следующий вид:

7-1

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

Если Рругз > Ртах то образуются две ударные волны и аналитического выражения для давления на контактом разрыве не существует. В этом случае уравнение (А.1) можно переписать так:

1{р) = {р - рь)дь(р) + (р - рк)дя{р) + Агг = О, Аи = ия- ии

где

i

9к(р) =

Ai

2

_р + вк

Это уравнение остается нелинейным. В работе Годунова [5], а затем в работе Topo [87] предложено заменить значение функции дк{р) на ее приближенное значение ди{ро), где давление ро некоторое приближенное давление. Поэтому первое приближение давления для решение задачи о распаде разрыва в этом случае запишется так:

= 9l(po)pl + 9r(po)pr ~ (uR - uL) 9ь(Ро) + gR(po)

За приближение ро в этом случае следует брать решение Ppws линеаризованной системы уравнений :

Ро = max(0,ppws).

Если же pmin < Ppvrs < Ртах то в одну сторону будет распространятся волна разрежения, а в другую ударная волна. Уравнение (А.1) не имеет аналитического решения и за первое приближение для решение задачи о распаде разрыва принимают РруГ8.

Решив уравнение (А.1) можно определить скорость и давление на контактном разрыве . Значение же плотности слева и справа от контактного разрыва зависит от типа прилегающей волны.

Если давление на контактном разрыве больше давления слева р* > Рь, то влево распространяется ударная волна. Тогда плотность слева от контактного разрыва записывается так:

Р*ь = рь

2*.

Рь 7+1

£• + 1

. 7+1 рь

Скорость же распространения левой ударной волны выражается следующим соотношением:

= сь

7 + 1Р* 7

27 рь 27

Если же р* < рь, то влево распространяется веер волн разрежения, и плотность слева от контактного разрыва определяется исходя из того, что течение является изоэнтропным:

Р*ь = Рь

р*

:рь.

а скорость звука за волной разрежения:

Р*

с*ь = сь — \РЬ,

2=1 27

При этом скорость первой или "головной" волны Бнь равна:

Янь = иь ~ сь,

А последней или "замыкающей" Эть '■

Бть = и* — с*ь-

Значение параметров внутри веера волн разрежения Уь/С (рь/ап,иь/ап:Рь/ап)Т определяется на основе инвариантов Римана:

Р = Рь

, (7-1) ( _ х\ 7+1 (7+1)с/, ^ г)

сь + ^иь + |'

Уь/ап = < и =

7+1 = РЬ

2

7-1

_2_ (7-1) ( _ х\ 7+1 ^ (7+1 )сь ^иЬ Ь )

7-1

(А.З)

где наклон характеристик у лежит в пределах от головной волны до замыкающей:

х

Знь < - < Зть-

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

- Ря

Скорость распространения которой:

= ия + Сд

2* +1=1' Ря 7+1

1=1 Е*. + 1

[_ 7+1 рн ^

7 + 1Р* 7-1

27 рд 27

Если р* < ря, то вправо распространяется волна разрежения и ее характеристики имеют следующий вид:

р* 4 7

_'* \ I У*

Р*Я = РЯ i — , с*д = сд -

\Ря) \Ря

Скорости головной и замыкающей волны:

рЛ

Зяя = ия + сд, Бтя = и* - с*д.

Значение параметров внутри волны разрежения:

Р = Ря

Л___(т-1) л,.р _ ап

7+1 (7+1 )сс\иЯ г)

-СЦ + ^ия + |

Уя/ап = и =

7+1 Р== РД

_2___(7-1) (и _ х)

7+1 (7+1)ся^Л ^

2 7-1

2т 7-1

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

х

¿>тя < ~ < Яня-

Описанный выше алгоритм позволяет находить решение задачи о распаде разрыва в любой точке х, £. На практике необходимо знать значение

параметров после разрыва t > 0 в точке х — 0. В зависимости от скорости контактного разрыва и* > 0 возможны два варианта.

Если скорости контактного разрыва больше нуля и* > 0, то значение параметров на прямой (0, t) определяется по характеристикам левой воны, а характеристики правой волны нет необходимости рассчитывать см. рис. 2.1. Для случая, когда влево распространяется ударная волна р* > pl значение параметров в месте где произошел разрыв определяется следующим образом:

( у shock

I Vl, если Sl ^ О,

где Sl - скорость распространения левой ударной волны , V^ock - параметры потока за левой ударной волной. Для случая когда влево распространяется волна разрежения возможны три варианта:

Vl> если Shl ^ О, V(0,t) = VLfan, если Shl < 0 ^ STL, V/lLr\ если Stl ^ 0 ^ u*,

где Shl и Stl скорости первой и последней волны разрежения, V/£n -параметры за волной разрежения (слева от контактного разрыва), VLfan - параметры внутри левой волны разрежения, которые определяются по формуле (А.З) если принять f = 0.

Если же скорость контактного разрыва меньше нуля и* < 0, то все расчеты необходимо производить по формулам для правой волны в зависимости от ее конфигурации. Для правой ударной волны :

1/(0 t) = { К*ОСк' eC™U* ^ ° ^ Vr, если SR ^ 0,

а для правой волны разрежения:

V(0,t) = l

если щ < Str ^ о, vrfan, если STr ^ 0 ^ SHR, VR, если SHr < 0.

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

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

Уравнение для собственных чисел (3.8) можно переписать в следующем виде:

аг\ - 0.

<9(2

Откуда найдем:

ррух - А© рру\ + 1 - Х&ух Рр^х уу - А

РрУхУг - Авиг

ртух - Арт Рт^х - ЬртУх

РТ^хУу - \pTVy РТ'Чх^г ~ \pTVz

= 0.

р оо

2 рух - А р 0 0

рЬу рух — А р 0

риг 0 рьх — А р

ррНух — ХОН + А р(Н + г>1) - Хрух рухуу - Хруу рухуг - Хру2 (ртН + рСр)(ух - А)

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

ррух - А© р 0 0 Рт(их - А)

1 р{ух - А) 0 0 0

det о о р{ух - А) о о =0.

0 0 0 р(ух - А) 0

А рух(ух ~ А) руу(ух - А) руг{ух - А) рСр{ух - А)

Раскрывая этот определитель получим уравнение относительно собственных чисел:

(^х — А)3 (рСр—у\(рт+рРрСр) + \ух(2рт+рррСр+ОрСр) — \2(рт+ОрСр)) = 0.

Видно что три корня этого уравнения Ах^з = ух , а для А45 уравнение можно переписать если ввести параметр 0:

А2 - ух( 1 + (317?)А - и2г{1 - Ру2х) = 0, (Б.1)

где параметр /3 равен:

Р = Рр +

рт рСр

Корни уравнения (Б.1) равны:

А4)5 = ухТС, где у'х и с выражаются следующим образом:

у'х = ^С1 - а)

с — \/ о? у2х + и?-Параметр а имеет следующий вид:

1 - ри?

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

Найдем матрицы составленные из правых Мт и левых Мр 1 собственных векторов для Аг. Соотношение для правых собственные векторов гг по определению можно записать так:

Если подставить выражение для А-р, то получим:

'дК

х

дЯ

АгГ г< = 0.

г

г

Подставив выражение для матриц ^ и Г получим:

оЯ

( ррУх - Аг0 р 0 0 ртьх - \грт

ррух + 1 - \гОьх 2рух - \гр 0 0 рту\ - АхртУх

ррухъу - Аг8уу руу рух - \р 0 Ртухуу - КртУу

РрУХУг - рУг 0 рух - Агр РтУхУг ~ КРТ^г

\ РрНух - А,9Я + А,. р{Н + у1) - \грух рухуу - хгруу рухуг - \гру2 {ртН + рСр){ух - Аг) )

Вычтим из второй строки первую умноженную на ух, из третей строки первую домноженную на уу, из четвертой строки первую умноженную на

уг, а из пятой строки первую умноженную на Н получим: ( . ... .. . .. \

РРУХ - Л9 р 0 0 РТ(ух - А)

1 р{ух - А) 0 0 0

0 0 р(ух - А) 0 0

0 0 0 р(ух - А) 0

л рух(ух - А) руу(ух - А) ру2(ух - А) рСр{ух -А)

г 2

Г4 \Г5/

В случае если положить Л = ух для правых собственных векторов г1, г2, г3 тогда уравнение примет вид:

/ . Л Л Л\

V

РрЩ - ух9 р О О О

1 0 0 0 0

0 0 0 0 0

О 0 0 0 0

гл,. 0 0 0 0

/

Г2 Гз

\гч

= 0.

Откуда очевидно г\ = 0,г2 — 0, так как только при таких значениях выполнены соотношения. Компоненты вектора г3, г4, Г5 любые, но они должны быть такие, чтобы вектора гг были линейно независимые. Выбираем г3,г4,г5 следующим образом:

гх =

М м

0 0 0

0 , Г2 = 0 , г3 = 1

0 1 0

к1)

Для правых собственных векторов г4, г5 с соответствующими собственными числами Л4 = у'х — с , Л5 = у'х + с из второго уравнения можно получить что:

г 2 =

п

/о(Л - Ух)

а из третьего и четвертого уравнения получим гз = 0,74 = 0. Из пятого уравнения значение 75 равно:

П =

г\

рСр

Так если Л = и' — с , то получаем А — и = — (аи + с) и выбираем значение для г\ = аи + с . Тогда правый собственный вектор соответствующий Л4 равен:

г4 =

(

аи + с

_! р

\

О О

(■аи+с ) \ РСр /

Для Л5 = и + с , Л — и = — (аи — с) и выбираем г\ — аи — с .

г5 =

( ' ^

аи — с

_ 1 р

О О

(аи—с )

\ —¿сГ /

Тогда матрицу правых собственных векторов можно записать так:

^ 0 0 0 (аи + с) (ски - с') ^

Мг =

1 р

ООО-1

р

0 0 10 о 0 10 0 о

10 0 (аи+с) (ац~с) рСр рСр

Матрица левых собственных векторов является обратной матрице Мр и имеет следующий вид:

0

рСр

О О О о

Мг1

_1_ 2с'

V

_1_ '2с'

__

2с

О 0 1 \

0 1 о

1 о о

ООО