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

Сахненко Александр Константинович

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

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

Общественные науки в целом

Дипломы

Вуз: Санкт-Петербургский государственный университет (СПбГУ)

ID: 587d36365f1be77c40d5898d

UUID: e2e3a4e3-d0ff-4f7d-8813-c759e1d8c712

Язык: Русский

Опубликовано: почти 8 лет назад

Просмотры: 16

Сахненко Александр Константинович

Источник: Санкт-Петербургский государственный университет

Комментировать 0

Рецензировать 0

Скачать - 729803 bytes

Поделиться работой

САНКТ-ПЕТЕРБУРГСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ Кафедра компьютерного моделирования и многопроцессорных систем Сахненко Александр Константинович Выпускная квалификационная работа бакалавра Оптимизация и распараллеливание численных алгоритмов для сильно нелинейных волновых процессов Направление 010300 Фундаментальная информатика и информационные технологии Научный руководитель, доктор физ.-мат. наук, профессор Богданов А.В. Санкт-Петербург 2016

Содержание Введение..................................................................................................... 3 Постановка задачи .................................................................................... 5 Глава 1. Описание предметной области ................................................. 6 Глава 2. Разностные схемы «предиктор-корректор» ............................ 8 2.1. Описание ............................................................................................ 8 2.2. Метод Мак-Кормака ......................................................................... 9 Глава 3. Обзор архитектуры Kepler....................................................... 11 Глава 4. Построение разностной схемы ............................................... 15 4.1. Конечно-разностная схема для уравнения КдВБ......................... 15 4.2. Начальные распределения.............................................................. 17 4.3. Чистый перенос ............................................................................... 18 4.4. Уравнение Бюргерса ....................................................................... 19 4.5. Уравнение Кортевега – де Вриза ................................................... 20 Глава 5. Реализация алгоритма .............................................................. 22 5.1. Работа с памятью............................................................................. 22 5.2. Функции ядра .................................................................................. 25 5.3. Запуск алгоритма............................................................................. 26 Выводы ..................................................................................................... 28 Заключение .............................................................................................. 29 Список литературы ................................................................................. 30 2

Введение Нелинейные явления появляются в природе повсюду: от волн на водной поверхности до магнитных, от оптики до прогнозов погоды. Следовательно, их описание и понимание имеет принципиальное значение, как с теоретической, так и прикладной точки зрения. Нелинейные явления, как правило, описываются дифференциальными уравнениями, решение которых часто является сложной проблемой. Тем не менее, существует специальный класс дифференциальных уравнений, которые разрешимы (в некотором смысле) - интегрируемые системы. Многие понятия современной математической физики, такие как солитоны, инстантоны и квантовые группы имеют свое происхождение в теории интегрируемых систем. Когда физическое явление описывается интегрируемой системой, ее поведение может быть понято во всем мире и его часто можно предсказать. Красота этой теории заключается в его универсальности: многие фундаментальные нелинейные уравнения оказываются как широко применимы, так и интегрируемы. Кроме того, в ряде случаев приближены, неинтегрируемые при некоторых нелинейные уравнения предположениях, могут быть нелинейными интегрируемыми уравнениями, что позволяет лучше понять явления, моделируемые ими. Однако апроксимация неинтегрируемых нелинейных уравнений не всегда положительно сказывается на получаемых результатах, и найденное решение отличается от истинного довольно сильно. Решение же самих неинтегрируемых нелинейных уравнений очень трудно даже численно и практически невозможно с помощью стандартных аналитических методов. Особенно трудно получить надежные результаты в асимптотической области. 3

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

Постановка задачи Целью данной работы является оптимизация и распараллеливание алгоритма численного решения дифференциальных уравнений, описывающих сильно нелинейные волновые процессы в гидродинамике, на примере уравнения Бюргерса. Для достижения поставленной цели предполагается решить задачу распараллеливания алгоритма, осонованного на двухшаговом методе численного интегрирования «предиктор-корректор» Мак-Кормака. Решение данной задачи требует реализации следующих этапов: 1. Анализ предметной области; 2. Обзор архитектуры Kepler графического ускорителя NVIDIA; 3. Построение разностной схемы; 4. Реализация параллельного алгоритма. 5

Глава 1. Описание предметной области В данной работе рассматривается процесс непрерывно- корпускулярного моделирования волновых процессов в жидкостях. Макроскопическое поведение любого материала определяется, в конечном счете, его микроструктурой, которая на молекулярном уровне носит принципиально межмолекулярных дискретный характер. взаимодействиях, О связанных микроструктуре с и определенными материалами, известно довольно много, и поэтому втаёт вопрос: может ли, и до какой описании степени, непрерывных направлении сущностей такая является (например, объединение корпускулярными информация моделей. Первым обработка молекул, теории использована шагом дискретных атомов и непрерывного соображениями. быть Подробнее в этом фундаментальных ионов), как частиц моделирования это в описано в и с книге И.Мёрдока [1]. Весь объём жидкости, при таком подходе, представляется набором крупных «частиц» воды, движение каждой из которых в каждый момент времени характеризуется набором внутренних состояний, инерцией и внешними силами, в том числе, со стороны других частиц (сила натяжения) [2]. Непрерывно-корпускулярный подход основан на численных схемах первого порядка с постоянным шагом в интегрировании законов движения. Каноническое представление о законах механики жидкости позволяет строго и однозначно связать численные объекты с арифметическими и логическими операциями и сложными геометрическими алгоритмами, в том числе с использованием быстрой интерполяции для нерегулярных сеточных пространств. Основной вектор математики в этой области сейчас направлен на создание прямых вычислительных экспериментов при решении практических 6

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

Глава 2. Разностные схемы «предиктор-корректор» 2.1. Описание Схемы «предиктор-корректор» - семейство методов, относящихся к классу алгоритмов, предназначенных для интегрирования обыкновенных дифференциальных уравнений. Все такие методы включают в себя два шага:  На первом шаге (предиктор) вычисляется некоторая функция от значений, посчитанных на предыдущем шаге, чтобы получить апроксимированное значение искомой функции в следующей точке.  На втором шаге (корректор) совершенствуется полученное приближение при помощи спрогнозированного значения функции и другого оператора, чтобы интерполировать значение искомой функции в той же самой точке. Эту схему можно представить следующей системой: 1 𝑢𝑛+2 = 𝑓 1 (𝑢𝑛 , … ) , { 1 𝑢𝑛+1 = 𝑓 2 (𝑢𝑛+2 , … ) 1 где 𝑢𝑛+2 – значение, полученное на шаге «предиктор», а 𝑢𝑛+1 – уточнённое искомое значение. Т.о. видно, что первый шаг реализуется с помощью явных методов, а второй шаг основан на применении формул неявных методов, в правую часть которых вместо неизвестного значения 𝑢𝑛+1 подставляется результат 1 предсказания 𝑢𝑛+2 . Методы, использующие схему предиктор-корректор:  Метод Милна – для ОДУ;  Метод Хойна (предиктор – Метод Эйлера, корректор – Метод трапеций); 8

При использовании схемы предиктор-корректор для решения ОДУ отмечают высокую точность расчета и отсутствие свойства самостартуемости (то есть для начала вычислений по схеме предиктор-корректор требуется предварительно воспользоваться другим, самостартующим методом).  Метод Адамса-Башфорта – параллельный п.-к. для решения нежестких краевых задач (используется корректор Адамса-БашфортаМултона);  Формулы Хемминга.  Метод Мак-Кормака. 2.2. Метод Мак-Кормака Метод Мак-Кормака является измененной двухшаговой схемой Лакса-Вендрофа, но он гораздо проще в применении. Рассмотрим следующее гиперболическое уравнение первого порядка: 𝜕𝑢 𝜕𝑢 +𝑎 =0 𝜕𝑡 𝜕𝑥 Предиктор: На этом шаге прогнозируемое значение 𝑢 в момент времени k + 1 (𝑢̅𝑖𝑘+1 ) оценивается следующим образом: 𝑢̅𝑖𝑘+1 = 𝑢𝑖𝑘 − 𝑎 Δ𝑡 𝑘 (𝑢𝑖+1 − 𝑢𝑖𝑘 ) Δ𝑥 Корректор: На этом шаге предсказанное значение 𝑢̅𝑖𝑘+1 корректируется в соответствии с уравнением: 1 Δ𝑡 𝑘+1 𝑘+1 𝑢𝑖𝑘+1 = [𝑢𝑖𝑘 + 𝑢̅𝑖𝑘+1 − 𝑎 (𝑢̅ − 𝑢̅𝑖−1 )] 2 Δ𝑥 𝑖 Отметим, что промежуточное значение 𝑢̅𝑖𝑘+1 никакого физического смысла не имеет. 9

Схема имеет второй порядок точности с погрешностью аппроксимации О(Δ𝑡 2 , Δ𝑥 2 ), при этом она устойчива при |С| ≤ 1, а дифференциальное приближение схемы имеет вид: Δ𝑡 2 Δ𝑥 3 2 (1 − 𝐶 )𝑢𝑥𝑥𝑥 − 𝑎 𝑢𝑡 + 𝑎𝑢𝑥 = −𝑎 𝐶(1 − 𝐶 2 )𝑢𝑥 4 + . .. , Δ𝑥 8 где С = 𝑎 Δ𝑡 Δ𝑥 – число Куранта. Для линейного волнового уравнения схема Мак-Кормака эквивалентна схеме Лакса-Вендроффа. Схема Мак-Кормака часто применяется в силу ряда своих достоинств [3]. В частности, она оперирует только величинами в основных узлах сетки и легко обобщается на многомерные задачи. Также, это схема второго порядка точности. Она обладает лишь малой диффузией, но имеет большую дисперсию и чувствительна к нелинейной неустойчивости. В областях больших градиентов могут возникать нефизические осцилляции. Эти численные колебания часто приводят к отрицательным значениям заведомо положительных величин. 10

Глава 3. Обзор архитектуры Kepler Компания NVIDIA вот уже почти целое десятилетие выпускает видеокарты для общих вычислений. Нельзя сказать, что раньше было невозможно производить вычисления на GPU, но этот процесс был довольно неудобен и требовал большого количества времени и определённых навыков. В 2007 году была предложена универсальная архитектура параллельных вычислений (CUDA), позволяющая благодаря использованию GPU повысить производительность вычислений. Это дало возможность наращивать мощность устройств путём увеличения количества их составляющих, при этом не изменяя основные особенности архитектуры. Таким образом, объём памяти и количество процессоров постоянно увеличиваются, а вот разделение памяти на глобальну, разделяемую, текстурную и регистровую остаётся однообразным, неизменным в каком-то смысле (позже появились и другие структуры памяти, такие как surface или разные типы кэша, но это только расширило возможности CUDA). Однако всё должно развиваться и архитектура и набор команд не стоят на месте и со временем модифицировались – иногда они изменялись значительно, иногда не очень. Для классификации семейств GPU с идентичной архитектурой было введено понятие «Compute сapability» (СС). Для сравнения: процессоры с CC 1.0 не поддерживали атомарных операций в принципе, с СС 1.1 – могли производить их в глобальной памяти, а с СС 1.2 – ещё и в shared. Полный список возможностей разных СС традиционно приводится в конце CUDA C Programming Guide [4]. В 2012 году NVIDIA представила архитектуру Kepler (СС 3.0). Она рассчитана на энергоэффективность, программируемость и производительность, в отличае от предыдущей архитектуры, Fermi, основная направленность которой была сконцентрирована производительности. 11 на чистой

Уменьшение энергопотребления происходит благодаря тому, что шейдерные блоки работают на одной частоте с ядром. Это позволило снизить энергопотребление даже при том, что необходимо использовать большее количество шейдерных ядер для получения такой же производительности, которая достигалась Энергоэффективность с помощью достигается не предыдущих только тем, разработок. что новая архитектурапотребляет меньше энергии, чем архитектура предыдущего поколения (двум шейдерным ядрам Kepler требуется приблизительно 90% питания, потребляемого одним ядром Fermi), но и благодаря тому, что унификация тактовой частоты приводит к снижению частоты шейдерных блоков, что в свою очередь серьёзно снижает энергопотребление[5]. Улучшенная программируемость была достигнута за счёт введения новой модели обработки текстур, которая не требует привязки к CPU. Улучшение же производительности было достигнуто за счёт внедрения абсолютно новых контроллера памяти и шины. В свою очередь это позволило поднять тактовую частоту памяти до 6 ГГц, что всё ещё ниже, чем теоретически максимальные для GDDR5 7 ГГц, но значительно больше, чем частота памяти в 4 ГГц при архитектуре предыдущего поколения [6]. Также, архитектуру Kepler отличает от предшественников новый огромный мультипроцессор. Если раньше мультипроцессоры имели 8 (CC 1.x), 32 (CC 2.0) или 48 (СС 2.1) потоковых процессоров, то в Kepler используется новый чип на 192 процессора. Другие характеристики в сравнении с архитектурой Fermi представлены ниже (Таблица.1). 12

Таблица 1. Сравнение характеристик процессоров архитектуры Kepler и Fermi FERMI FERMI GF104 GF100 KEPLER KEPLER GK104 GK110 Версия CC 2.0 2.1 3.0 3.5 Потоков в Warp'e 32 32 32 32 Число warp'ов на 48 48 64 64 1536 1536 2048 2048 8 8 16 16 32768 32768 65536 65536 63 63 63 255 Конфигурации 16К 16К 16К 16К количества shared- 48К 48К 32К 32К 48К 48К мультипроцессор Потоков на мультипроцессор Блоков на мультипроцессор 32-битные регистры на мультипроцессор Максимальное количество регистров на поток памяти Максимальный 2^16‐1 2^16‐1 2^16‐1 2^32‐1 Hyper‐Q Нет Нет Нет Есть Динамический Нет Нет Нет Есть размер grid по оси X параллелизм Чтобы извлечь наибольшую выгоду из использования потокового мультипроцессора (SMX), т.е. максимизировать количество одновременно загруженных блоков/потоков, нужно использовать гораздо меньше регистров, чем при использовании архитектуры предыдущего поколения. С другой стороны, Kepler позволяет хранить до 255 регистров на поток, и если раньше было справедливо утверждение «лучше заново вычислить, чем хранить», то сейчас оно подвергается сомнению. Таким образом, имеет место некая «развилка» на пути оптимизации – можно записывать результаты 13

вычислений в регистры, но тогда придётся пожертвовать количеством блоков, которые могут поместиться на SMX, что, вполне вероятно, приведёт к простоям на чтение/запись в RAM. Ещё несколько заметных нововведений: 1. Новый Warp Scheduler. На каждый мультипроцессор приходятся 4 планировщика, и каждый из них умеет исполнять 2 инструкции варпа за такт; 2. Dynamic Parallelism – позволяет ядру CUDA создавать и синхронизировать вложенные задачи, используя CUDA во время выполнения API для запуска других ядер, управлять памятью устройства, а также создавать и использовать потоки и события, и все это без участия центрального процессора; 3. GPUDirect – данная технология повышает скорость передачи данных между GPU и другими устройствами, обходя использование CPU, в частности – обеспечивает между графическими процессорами, находящимися в одной системе, peer-to-peer (P2P) соединение; 4. Shuffle Instruction – благодаря этим инструкциям параллельные потоки могут считывать значения регистров друг друга; 5. Read-only кэш данных – Kepler предоставляет возможность отмечать данные как read-only, чтобы компилятор помещал их в соответствующий кэш; 6. Hyper‐Q – блок, позволяющий сразу нескольким ядрам процессора работать с GPU. 7. Grid Management Unit – модуль, призванный утилизировать ресурсы GPU, что приводит к увеличению одновременно выполняемых работ. Подробнее это описано в [4]. 14

Глава 4. Построение разностной схемы 4.1. Конечно-разностная схема для уравнения КдВБ Рассмотрим уравнение Кортевега – де Вриза – Бюргерса (КдВБ), которое является универсальным для описания одномерных нелинейных волн в средах с дисперсией и диссипацией: 𝑢𝑡 + 𝑢𝑢𝑥 + 𝛼𝑢𝑥𝑥 + 𝛽𝑢𝑥𝑥𝑥 = 𝛾𝐼(𝑢), (1) где 𝛼 – мера диссипации, 𝛽 – мера дисперсии, 𝛾 – мера межфазных взаимодействий, а 𝐼(𝑢) – интегральный оператор. При интегрировании уравнения КдВБ, вместо исходного уравнения мы используем его эквивалент, написанный в дивергентной форме: 𝑢𝑡 + 0.5(𝑢2 )𝑥 + 𝛼𝑢𝑥𝑥 + 𝛽𝑢𝑥𝑥𝑥 = 𝛾𝐼(𝑢), (2) Решение данного уравнения в полуплоскости 𝑡 ≥ 0 ищется для начального 𝑢(𝑥, 0) = 𝑓(𝑥), 𝑥 ∈ (−∞; +∞). распределения Граничные условия – условия Дирихле: 𝑢(−∞, 𝑡) = а, 𝑢(+∞, 𝑡) = 𝑏. В основном, численное моделирование уравнения (2) проводилось с использованием двухступенчатой явной разностной схемы Мак-Кормака [7][8] вместе с методом коррекции потоков [9]. Также есть и другие численные подходы, в частности обобщение Залесака для гиперболических систем [10]. Вычисления выполняются в области [𝑥𝑚𝑖𝑛 , 𝑥𝑚𝑎𝑥 ] × [0, 𝑇]. Сетка конечных разностей по времени t и пространству координаты x строится следующим образом: 𝑢𝑗𝑛 = 𝑢(𝑥𝑗 , 𝑡 𝑛 ), 𝑥𝑗 = 𝑗Δ𝑥, 𝑡 𝑛 = 𝑛Δ𝑡, 𝑗 ∈ [1, 𝑀], 𝑥1 = 𝑥𝑚𝑖𝑛 , 𝑥𝑀 = 𝑥𝑚𝑎𝑥 , 𝑛 = 0,1, …, 15

где Δx является шагом пространственной координаты; Δt – шаг по времени. Граничные условия: 𝑢(𝑥𝑚𝑖𝑛 , 𝑡) = а, 𝑢(𝑥𝑚𝑎𝑥 , 𝑡) = 𝑏. Разностная схема Мак-Кормака для уравнения (2) имеетследующий вид (𝑟 = 𝛼𝛥𝑡 𝛥𝑥2 , 𝑠= 𝛽𝛥𝑡 2𝛥𝑥3 ): Шаг 1. Предиктор. 𝑢𝑗∗ = 𝑢𝑗𝑛 − Δ𝑡 2 2 𝑛 𝑛 𝑛 [(𝑢𝑗+1 − 2𝑢𝑗𝑛 + 𝑢𝑗−1 ) − (𝑢𝑗𝑛 ) ] − 𝑟(𝑢𝑗+1 )− 2Δ𝑥 𝑛 𝑛 𝑛 𝑛 −𝑠(𝑢𝑗+2 − 2𝑢𝑗+1 + 2𝑢𝑗−1 + 𝑢𝑗−2 ) + 𝛾Δ𝑡𝐼𝑗𝑛 Шаг 2. Корректор. 𝑢𝑗𝑛 + 𝑢𝑗∗ Δ𝑡 𝑟 ∗ 2 2 ∗ ∗ 𝑢= − [(𝑢𝑗∗ ) − (𝑢𝑗−1 − 2𝑢𝑗∗ + 𝑢𝑗−1 ) ] − (𝑢𝑗+1 )− 2 4Δ𝑥 2 𝑠 ∗ Δt ∗ ∗ ∗ − (𝑢𝑗+2 − 2𝑢𝑗+1 + 2𝑢𝑗−1 + 𝑢𝑗−2 ) + 𝛾 𝐼𝑗𝑛 2 2 Приближение переносного члена в шаге «предиктор» для нечетных шагов по времени выполняется разностями вперед, а для четных шагов – разностями назад. Таким образом, ошибка фазы уменьшается в расчетах, что приводит к движению волновых фронтов в правильном направлении. Известно, что схема Мак-Кормака, повышая при этом порядок аппроксимации по сравнению с монотонными схемами, порождает в окрестности сильных скачков колебаня, которые носят нефизический смысл. Метод коррекции потоков, согласно [8], включает в себя следующие этапы: 1. Вычисление диффузионных потоков: 𝑢𝑑 1 𝑗+ 2 𝑛 = 𝜈𝑗+1 (𝑢𝑗+1 − 𝑢𝑗𝑛 ) 2 2. Вычисление антидиффузионных потоков: 𝑢𝑎𝑑1 = 𝜇𝑗+1 (𝑢𝑗+1 − 𝑢𝑗 ) 𝑗+ 2 2 Коэффициенты 𝜈𝑗+1 и 𝜇𝑗+1 вычисляются следующим образом: 2 2 16

𝜈𝑗+1 2 1 1 = + 𝐶 2 1, 6 3 𝑗+2 После 𝜇𝑗+1 2 1 1 = − 𝐶2 1, 6 3 𝑗+2 корректирующего шага мы 𝐶𝑗+1 2 𝑛 𝑢𝑗𝑛 + 𝑢𝑗+1 Δ𝑡 = . 2 Δ𝑥 получаем отрегулированное значение 𝑢𝑗𝑐𝑜𝑟𝑟 . Окончательное решение на 𝑛 + 1 шаге имеет вид: 𝑐𝑜𝑟𝑟 𝑢𝑗𝑛+1 = 𝑢𝑗 + 𝑢𝑑 1 − 𝑢𝑑 1 − 𝑢𝑗𝑐𝑜𝑟𝑟 + 𝑢𝑗−1 . 𝑗+ 2 𝑗− 2 4.2. Начальные распределения При расчетах в качестве начальных профилей используются три распределения: 1. "Ударная волна" (прыжок): 𝑢(𝑥 ≥ 𝑥0 , 0) = 𝑢1 ; 𝑢(𝑥 < 𝑥0 , 0) = 𝑢2 . Скорость сдвига скачка: 𝑈 = (𝑢1 + 𝑢2 )/2 = 𝐶𝑜𝑛𝑠𝑡. 2. Солитон. Известное стационарное решение уравнения КдВ: 𝑢(𝑥, 0) = 𝑓(𝑥) = 𝜀 cosh−2 𝑥/𝛿 , 𝛿 = √12𝛽/𝜀 . 𝑤 𝑤 3. Прямоугольный профиль: 𝑢(𝑥, 0) = ℎ × 𝑤, при 𝑥 ∈ [− , ]. Выбор 2 2 между h и w продиктован законом сохранения «массы» для исходного профиля: ∞ ∞ −∞ −∞ ∫ 𝑓(𝑥)𝑑𝑥 = ∫ 𝜀 cosh−2 𝑥/𝛿 = 4√3𝜀𝛿 = ℎ × 𝑤 (3) При рассмотрении уравнения КдВ, выражение (3) имеет вид: ℎ × 𝑤 = 12. Для поддержания этого равенства в расчете, конкретное значение 𝑤 также должно быть связано с шагом по координате 𝑥. 17

4.3. Чистый перенос Рассмотрим частные случаи исходного уравнения (1): чистый перенос, уравнение Бюргерса, уравнение КдВ. Базовые свойства разностной схемы (устойчивость, порядок апроксимации и т.п.) описаны в работах [7], [8] В случае чистого переноса уравнение (2) записывается в виде: 𝑢𝑡 + 0.5(𝑢2 )𝑥 = 0. Конечно-разностная схема для (4) чистого переноса традиционно тестируется «ударной волной». Пусть скачок имеет единичную амплитуду: 𝑢1 = 1, 𝑢2 = 0, 𝑥0 = 0. Тогда скорость сдвига скачка 𝑈 = 1/2. Уравнение (4) запишется в линейной форме как: 𝑢𝑡 + 𝑈𝑢𝑥 = 0. На Рисунке 1 показаны результаты численных вычислений для 125 шагов, С=0.25 и 𝑥 = 𝑈𝑡 = 62.5. Результаты численных вычислений и аналитическое решение для «ударной волны» доволно похожи. Рисунок 1. Сравнение результатов для распределения "ударной волны" Для линейного случая (4), с использованием метода скорректированных потоков, существует известное условие устойчивости: 𝐶 ≤ 1/2, где 𝐶 = 𝑈𝛥𝑡/Δ𝑋 – это число Куранта, а Δ𝑡 и Δ𝑥 – шаги по времени и координате. Для нелинейного случая условие устойчивости преобразуется в следующее: 18

𝑢𝑗𝑛 (Δ𝑡) |max | ≤ 1/2. ∀𝑗,𝑛 Δ𝑥 Из формулы (3) хорошо видно влияние ошибки дисперсии в зависимости от числа Куранта. Без алгоритма коррекции потоков, мы не можем получить правильное численное решение. Вычисления показывают, что при шаге Δ𝑥 = 0.1 предельная величина шага по времени Δ𝑡 = 0.071 (𝐶 = 0.355). 4.4. Уравнение Бюргерса При 𝛽 = 𝛾 = 0 уравнение (2) записывается в виде уравнения Бюргерса: 𝑢𝑡 + 𝑢𝑢𝑥 + 𝛼𝑢𝑥𝑥 = 0 (5) Для уравнения (5) существует множество точных решений [11],[12]. В данном случае для тестирования схемы будет использоваться стационарное аналитическое решение: 𝑢= 𝑢2 − 𝑢1 , 𝑢 −𝑢 1 + exp (− 2 2𝛼 1 𝜉) 𝜉 = 𝑥 − 𝑈𝑡. Как и в случае чистого переноса для начального распределения, мы возьмем единичный скачок: 𝑢1 = 1, 𝑢2 = 0, 𝑥0 = 0. На Рисунке 2 аналитическое решение сравнивается с эволюцией единичного скачка и со скачком без влияния вязкости. Рисунок 2. Сравнение числовых и аналитических решений для уравнения Бюргерса 19

Доминирующий диссипативный эффект – физическая диффузия, которая подавляет числовую ошибку дисперсии. 4.5. Уравнение Кортевега – де Вриза При 𝛼 = 𝛾 = 0 уравнение (2) записывается в виде уравнения КдВ: 𝑢𝑡 + 0.5(𝑢2 )𝑥 + 𝛽𝑢𝑥𝑥𝑥 = 0 (6) Уравнение (6) не имеет диссипативных членов, и поэтому численные результаты очень чувствительны к ошибкам аппроксимации, особенно фазовым ошибкам. И не только члены с производными первого порядка имеют эффект ошибки аппроксимации, но и с производными высших порядков. Непосредственное применение описанного выше метода конечных разностей приводит к ошибочным результатам. Это связано с тем, что процедура коррекции сама имеет диффузию сетки. Результаты, полученные при применении процедуры коррекции, как и ожидалось, являются вполне удовлетворительными. Таким образом, простое уменьшение шага по времени может существенно повлиять на получаемое решение. Скорость солитона распространения U = ε / 3 для стационарных решений (ε = 12) равна 4 (𝑟 = 𝛽𝛥𝑡/(2𝛥𝑥 3 ). Скорость движения солитона в случае 2 на рис. 3 составляет около 3,97, относительная погрешность равна 0,64%, что является вполне приличным результатом. Относительная погрешность расчетной амплитуды солитона составляет 0,58%. 20

Рисунок 3.Сравнение решений для различных временных шагов: 1 − 𝑟 = 6.25 ∙ 10−3 ; 2 − 𝑟 = 6.25 ∙ 10−4 21

Глава 5. Реализация алгоритма В данной главе мы рассмотрим реализацию алгоритма решения уравнения Бюргерса: 𝑢𝑡 + 𝑢𝑢𝑥 + 𝛼𝑢𝑥𝑥 = 0 Для реализации алгоритма была выбрана технология CUDA и язык программирования С++. Также, при работе с этой технологией полезными могут оказаться такие работы, как [4], [13], [14], [15]. 5.1. Работа с памятью Константы, необходимые для нахождения решения, помещаются в конмтантную память (constant) – кэшируемую область DRAM, доступ к которой у GPU есть только на чтение, а у хоста (CPU) – на чтене и запись при помощи функции CUDA API: cudaError_t cudaMemcpuToSymbol( const char * symbol, const void * src, size_t count, size_t offset, enum cudaMemcpuKind kind); Параметр kind принимает одно из следующих значений, задающих направление копирования:  cudaMemcpuHostToHost;  cudaMemcpuHostToDevice;  cudaMemcpuDeviceToHost;  cudaMemcpuDeviceToDevice Несмотря на то, что константная память находится в DRAM, она имеет высокую скорость доступа благодаря наличию независимого кэша. В рассматриваемом алгоритме используются следующие константы (здесь dt – шаг по времени, dx – шаг по координате x):  tx = dt/(2*dx) – отношение шага по времени к шагу по координате, делённое на 2 (используется на этапе «предиктор»); 22

 tx2 = dt/(4*dx) – отношение шага по времени к шагу по координате, делённое на 4 (используется на этапе «корректор»);  atx=alf*ta/pow(dx,2) – отношение шага по времени к шагу по координате в квадрате, умноженное на параметр 𝛼 (используется на этапе «предиктор»)  atx2=alf*ta/(2*pow(dx,2)) отношение шага по времени к шагу по координате в квадрате, делённое на 2 и умноженное на параметр 𝛼 (используется на этапе «корректор»);  istep – номер шага по времени. Константная память выделяется через статическое объявление в коде с добавлением атрибута __constnat__: __constant__ __constant__ __constant__ __constant__ float float float float tx; tx2; atx; atx2; float host_tx = dt / (2*dx); float host_tx2 = dt / (4*dx); float host_atx = alf * dt / pow(dx,2); float host_atx2 = alf * dt / (2*pow(dx,2)); cudaMemcpuToSymbol(tx, host_tx, sizeof(data), 0, cudaMemcpuHostToDevice); cudaMemcpuToSymbol(tx2, host_tx2, sizeof(data), 0, cudaMemcpuHostToDevice); cudaMemcpuToSymbol(atx, host_atx, sizeof(data), 0, cudaMemcpuHostToDevice); cudaMemcpuToSymbol(atx2, host_atx2, sizeof(data), 0, cudaMemcpuHostToDevice); В ходе выполнения алгоритма для хранения полученных промежуточных и окончательных значений используются следующие массивы:  u0 – массив значений функции 𝑢 на входе итерации;  us – массив значений 𝑢∗ послевыполнения шага «предиктор»; 23

 u1 – массив промежуточных значений функции 𝑢 и на выходе из итерации;  cu_d – массив значений диффузионных коэффициентов на этапе коррекции потоков;  cu_ad – массив значений антидиффузионных коэффициентов на этапе коррекции потоков;  u_d – массив значений диффузионого потока на этапе коррекции потоков;  u_ad – массив значений антидиффузионого потока на этапе коррекции потоков;  d_u – массив первых разностей. Все данные массивы хранятся в глобальной памяти GPU, выделение памяти и копирование в неё данных происходит следующим образом (здесь и далее n – количество точек): int nb = n * sizeof(float); float* u0 = NULL; cudaError_t cuerr = cudaMalloc( (void**)&u0, nb); cuerr = cudaMemcpy(u0, host_u0, cudaMemcpuHostToDevice); При этом значение переменной cuerr будет содержать информацию, об ошибках при выполнении функции (значение cudaSuccess соответствует успешному выполнению операции). 24

5.2. Функции ядра Фунции ядра описываются с атрибутом __global__, ниже описаны некоторые из них: 1. Предиктор: __global__ void predictor_kernel (float* u0, float* us) { int j = threadIdx.x + blockIdx.x * blockDim.x + 1; int jl = j-1; int jr = j+1; us[j] = u0[j] - tx * (pow(u0[j],2) - pow(u0[jl],2)) - atx*(u0[jr]-2*u0[j]+u0[jl]); } 2. Корректор: __global__ void corrector_kernel (float* us, float* u1) { int j = threadIdx.x + blockIdx.x * blockDim.x + 1; int jl = j-1; int jr = j+1; u1[j] = 0.5 * (u0[j]+us[j]) tx2 * (pow(us[jr], 2) - pow(us[j], 2)) atx2 * (us[jr] - 2*us[j] + us[jl]); } 3. Функция вычисления коэффициентов для метода коррекции потоков: __global__ void coef_fct_kernel (float* u0, float* cu_d, float* cu_ad) { int j = threadIdx.x + blockIdx.x * blockDim.x + 1; int jr = j+1; float ur = u0[j] + u0[jr]; float txx = pow(tx*ur, 2)/2; cu_d[j] = 1.0/6.0 + 1.0/3.0 * txx; cu_ad[j] = 1.0/6.0 + 1.0/3.0 * txx; } 25

4. Вычисление диффузионного и антидиффузионного потоков: __global__ void fct_kernel(float* u0, float* u1, float* cu_ad, float* cu_d, float* u_ad, float* u_d) { int j = threadIdx.x + blockIdx.x * blockDim.x + 1; int jr = j+1; u_d[j] u_ad[j] = cu_d[j] * (u0[jr] - u0[j]); = cu_ad[j] * (u1[jr] - u1[j]); } 5. Вычисление первых разностей для метода коррекции потоков: __global__ void d_fct_kernel(float* u1, float* u_d) { int j = threadIdx.x + blockIdx.x * blockDim.x + 1; int jl = j-1; u1[j] = u1[j] + u_d[j]-u_d[jl]; } 5.3. Запуск алгоритма Далее следует задать конфигурацию запуска n нитей и вызывать функции ядра для обработки данных согласно алгоритму, описанному в Главе 4: dim3 threads = dim3(BLOCK_SIZE,1); dim3 blocks = dim3((n-2) / BLOCK_SIZE, 1); predictor_kernel<<<blocks, threads>>> (u0, us); cuerr = cudaDeviceSynchonize(); corrector_kernel<<<blocks, threads>>> (us, u1); cuerr = cudaDeviceSynchonize(); coef_fct_kernel<<<blocks, threads>>> (u0, cu_d, cu_ad); cuerr = cudaDeviceSynchonize(); coef_fct_kernel<<<blocks, threads>>> (us, u1); cuerr = cudaDeviceSynchonize(); 26

fct_kernel<<<blocks, threads>>> (u0, u1, cu_ad, cu_d, u_ad, u_d); cuerr = cudaDeviceSynchonize(); d_fct_kernel<<<blocks, threads>>> (u1, u_d); cuerr = cudaDeviceSynchonize(); du_fct_kernel<<<blocks, threads>>> (u1, d_u); cuerr = cudaDeviceSynchonize(); correct_ad_fct_kernel<<<blocks, threads>>> (u_ad, d_u); cuerr = cudaDeviceSynchonize(); ad_fct_kernel<<<blocks, threads>>> (u1, u_ad); cuerr = cudaDeviceSynchonize(); После окончания алгоритма результаты выгружаются из памяти GPU в память CPU. 27

Выводы Основной целью данной работы была реализация параллельного алгоритма для сильно нелинейных волновых процессов. Для достижения поставленной цели был изучен непрерывно- корпускулярный подход к моделированию сильно-нелинейных волновых процессов в гидродинамике, проанализирована архитектура Kepler графического ускоритля NVIDIA и построена разнстная схема решения дифференциального уравнения в частных производных. Применение различных разностных схем в задачах решения дифференциальных уравнений в частных производных требует огромного количества операций над данными одинаковой структуры, а значит, эти задачи обладают высокой степенью параллелилизма, что побуждает использовать GPGPU для их решения. В ходе проделанной работы был реализован алгоритм решения дифференциального уравнения в частных производных Бюргерса, которое моделирует нелинейные волновые процессы в средах с диффузией и диссипацией. Алгоритм основан на разностной схеме Мак-Кормака и использует технологию CUDA для высокопроизводительных вычислений на GPU с архитектурой Kepler. Реализованный алгоритм является примером использования параллельных вычислений для задач моделирования нелинейных волн и может быть использован для построения алгоритмов решения более сложных уравнений. 28

Заключение В настоящее время графические ускорители приобретают всё большую популярность, причём они используются не только для решения специализированных задач компьютерной графики или обработки видео, но и для общих вычислений. Одной из задач, где возможно использование вычислений на GPU является решение дифференциальных уравнений в частных производных. Отдельно стоит подчеркнуть, что это очень трудоёмкий процесс, а производительности CPU недостаточно для решения задач моделирования нелинейных волновых процессов в реальном времени, что делает более предпочтительным использование высокопроизводительных систем и параллельных алгоритмов. В данной работе был показан процесс построения такого алгоритма для решения дифференциального уравнения в частных производных, описывающего нелинейные волновые процессы в гидродинамике. 29

Список литературы 1. Murdoch A. I. A corpuscular approach to continuum mechanics: basic considerations //Analysis and Thermomechanics. – Springer Berlin Heidelberg, 1987. – С. 81-111. 2. Bogdanov A., Khramushin V. Tensor Arithmetic, Geometric and Mathematic Principles of Fluid Mechanics in Implementation of Direct Computational Experiments //EPJ Web of Conferences. – EDP Sciences, 2016. – Т. 108. – С. 02013. 3. Pulliam T. H., Zingg D. W. Fundamental algorithms in computational fluid dynamics. – Switzerland : Springer, 2014. – С. 84-85. 4. Cuda C. Programming guide. – 2012. 5. Круглов В. Н., Папуловская Н. В., Чирышев А. В. Преимущества совместного использования CPU И CUDA-устройства //Фундаментальные исследования. – 2014. – №. 8-2. 6. Kepler (microarchitecture) на Wikipedia https://en.wikipedia.org/wiki/Kepler_(microarchitecture) 7. Pletcher R. H., Tannehill J. C., Anderson D. Computational fluid mechanics and heat transfer. – CRC Press, 2013, 774 p. 8. Fletcher C. Computational techniques for fluid dynamics 2: Specific techniques for different flow categories. – Springer Science & Business Media, 2012, 496 p. 9. Boris J. P., Book D. L. Flux-corrected transport. I. SHASTA, A fluid transport algorithm that works //Journal of computational physics. – 1973. – Т. 11. – №. 1. – С. 38-69. 10. Zalesak S. T. Fully multidimensional flux-corrected transport algorithms for fluids //Journal of computational physics. – 1979. – Т. 31. – №. 3. – С. 335-362. 30

11. Benton E. R., Platzman G. W. A table of solutions of the onedimensional Burgers equation //Quarterly of Applied Mathematics. – 1972. – С. 195-212. 12. Bogdanov A., Stankova E., Mareev V. High performance algorithms for multiphase and multicomponent media //14th Ship stability workshop, UTMSPACE, Malaysia. – 2014. – С. 242-245. 13. Боресков А. В. и др. Параллельные вычисления на GPU //Архитектура и программная модель CUDA. М.: Изд-во Московского университета. – 2012. 14. Боресков А., Харламов А. Основы работы с технологией CUDA. – Litres, 2015. 15. Сандерс Д., Кэндрот Э. Технология CUDA в примерах: введение в программирование графических процессоров //М.: ДМК Пресс. – 2011. – Т. 232. 31

Рецензии:

Авторизуйтесь, чтобы добавить рецензию

- у работы пока нет рецензий -

Отзывы:

Авторизуйтесь, чтобы оставить отзыв