автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.18, диссертация на тему:Разработка численных методов для решения уравнения переноса излучения и их реализация с использованием графических ускорителей

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

Цыбулип Иван Владимирович

РАЗРАБОТКА ЧИСЛЕННЫХ МЕТОДОВ ДЛЯ РЕШЕНИЯ УРАВНЕНИЯ ПЕРЕНОСА ИЗЛУЧЕНИЯ И ИХ РЕАЛИЗАЦИЯ С ИСПОЛЬЗОВАНИЕМ ГРАФИЧЕСКИХ УСКОРИТЕЛЕЙ

Специальность 05.13.18 —Математическое моделирование, численные методы и комплексы программ

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

11 ноя 2015

005564372

Москва —2015

005564372

Работа выполнена на кафедре информатики и вычислительной математики Московского физико-технический института (государственного университета)

Научный руководитель: Скалько Юрий Иванович,

кандидат физико-математических наук

Официальные оппоненты: Боговалов Сергей Владимирович,

доктор физико-математических наук, Национальный исследовательский ядерный университет «МИФИ», кафедра молекулярной физики, профессор

Лукин Владимир Владимирович,

кандидат физико-математических наук, Институт прикладной математики им. М.В. Келдыша Российской академии наук, отдел вычислительных методов и математического моделирования, научный сотрудник

Ведущая организация: Институт автоматизации проектирования Российской академии наук

Защита состоится «03» декабря 2015 г. в 11— часов на заседании диссертационного совета Д 212.156.05 на базе Московского физико-технического института (государственного университета) по адресу: 141700, Московская обл., г. Долгопрудный, Институтский пер., д. 9, ауд. 903 КПМ.

С диссертацией можно ознакомиться в библиотеке МФТИ и на сайте университета https://www.mipt.ru.

Автореферат разослан У*» £-5/7-9 2015 г.

Ученый секретарь

диссертационного совета Д 212.156.05 ^ /

— Федько Ольга Сергеевна

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

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

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

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

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

Целями данной работы являются:

1. Построение и исследование численных методов решения уравнения переноса излучения.

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

3. Моделирование линейчатого спектра излучения звезды типа Т Тельца при наличии конического ветра.

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

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

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

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

4. Создана параллельная, программная реализация метода, позволяющая проводить расчеты на графическом ускорителе.

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

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

7. Показана корректность данного алгоритма в случае использования тетраэдральной сетки, удовлетворяющей условию Делоне.

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

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

10. Выведено трассировочное соотношение и доказана лемма об устойчивой трассировке.

11. Метод реализован в многопроцессорном МР1-варианте, а также в варианте, использующем кластер из графических ускорителей.

12. Изучены вопросы ускорения и эффективности параллельной реализации.

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

14. По результатам гидродинамического моделирования звезды типа Т Тельца проведен расчет спектра излучения в линии Н-а при различных ориентациях плоскости акреционного диска.

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

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

2. Впервые предложен маршевый алгоритм для решения уравнения переноса на неструктурированной тетраэдральной сетке.

3. Предложена оригинальная распределенная многопроцессорная реализация метода длинных характеристик.

Теоретическая и практическая значимость работы: Предложены и исследованы новые методы решения уравнения переноса излучения. Представ-

лены алгоритмы распараллеливания предложенных методов. Оценена скорость сходимости методов по пространственным и угловым переменным.

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

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

Работа выполнялась при поддержке проекта Министерства образования и науки РФ №3.522.2014/К «Исследование процессов, происходящих в веществе, и изменения его свойств при импульсном вводе энергии высокой плотности».

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

Положения, выносимые на защиту соответствуют основным результатам диссертации, приведенным в заключении.

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

Основные результаты диссертации изложены в 8 публикациях [1—8], 2 из которых изданы в журналах, рекомендованных ВАК [1,2].

Основные результаты работы докладывались и получили одобрение специалистов на следующих научных конференциях:

1. 53-й научной конференции МФТИ «Современные проблемы фундаментальных и прикладных наук», Долгопрудный, 2010;

2. 54-й научной конференции МФТИ «Проблемы фундаментальных и прикладных естественных и технических наук в современном информационном обществе», Долгопрудный, 2011;

3. 55-й научной конференции МФТИ «Проблемы фундаментальных и прикладных естественных и технических наук в современном информационном обществе», Долгопрудный, 2012;

4. 56-й научной конференции МФТИ «Актуальные проблемы фундаментальных и прикладных наук в современном информационном обществе», Долгопрудный, 2013;

5. 57-й научной конференции МФТИ «Актуальные проблемы фундаментальных и прикладных наук в области физики», Долгопрудный,

2014.

Основные результаты работы докладывались и получили одобрение специалистов на следующих научных семинарах:

1. семинар лаборатории математического моделирования нелинейных процессов в газовых средах, МФТИ, Москва, 2012;

2. научная сессия VII школы по высокопроизводительным вычислениям, Университет Иннополис, Казань, 2015;

3. семинар лаборатории флюидодинамики и сейсмоакустики, МФТИ, Москва, 2015;

4. семинар института автоматизации проектирования РАН, Москва,

2015.

Основное содержание работы

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

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

Первая глава

В первом параграфе описывается математическая модель процесса переноса излучения. В работе рассматривается стационарное уравнение переноса для функции интенсивности излучения /(г, ÎÎ),

(nV)/ + *(r)/(r,i2) = *(r)/p(r),

которое справедливо в предположении отсутствия рассеяния и переизлучения на других частотах. Уравнение выражает баланс лучистой энергии, заключенной в элементарном объёме окружающем луч г — го = fis. Приводятся несколько типов граничных условий для уравнения переноса излучения. Показана роль излучения в системе уравнений радиационной газовой динамики.

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

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

оо

J х„dv

— оо

Здесь >cv — коэффициент поглощения на частоте v, av — сечение рассеяния, х)п — доля атомов, находящихся в состоянии п, N — концентрация -атомов, $пп> — сила осциллятора для перехода п —> п'. Для вычисления доли 0п использовалось ионизационное уравнение Саха

Jo^f «+ Шр / ^ч 1-Q+ uo рА3 V кг)'

ОО

1 j OvL

7Г6

= №п I crvdu = Ndn-fnn,.

тпес

о

где а+ — доля ионов водорода, и+, щ — электронные статистические суммы для иона и нейтрального атома, Л — тепловая дебройлевская длина волны электрона, с — энергия ионизации атома.

Вторая глава

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

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

по полусфере направлений, содержащихся в области. Описывается подход к построению такой квадратурной формулы.

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

ЩП) = ехр(-е<||П - Шг\\2) + ехр(-е4||П + с^||2).

(Пп)<0

£<р = — сИУ

d^ag(cl,c2,c3)

§гас! ср + >яр,

где среди чисел С\,С2, Сз не более двух различных. Аналогичное предобуслав-ливание производится и для метода с радиальными базисными функциями.

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

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

Третья глава

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

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

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

где 1\, ¡2 — значения интенсивности в концах ребра, 1тц — значение интенсивности в середине ребра, полученное до монотонизации, /тм — значение после монотонизации, с1атр(а:, [а, Ь]) = тах(а, тт(Ь, ж)).

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

clamp Imii -

h+h Г \h + h\ \h + h

2 ' 4 ' 4

/ / \

1 I

I1 1

Рис. 1 — Направление распространения излучения в численной схеме (контурная стрелка) и истинное направление П в случае нарушения условия размещения узлов (точки) на сторонах прямоугольной сетки

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

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

Четвертая глава

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

7(г0 + Пз, П) = а(0, з)1{г0, П) + /3(0, я),

справедливо трассировочное соотношение

а(во,з) = »(во^^а^ьв),

0Ы,= /3(в0, яОа^!, в) + в),

где точка го + Г^х расположена на луче между г0 + ГЬ0 и г0 + $15, см. рисунок 2. Данное соотношение крайне удобно, так как позволяет вычислить а(0, в) и /3(0, в), последовательно проходя по сеточным элементам от конца луча к началу.

Во втором параграфе вводится понятие входной, выходной и касательной к направлению излучения грани тетраэдра. Если п — нормаль к грани, а (I — направление излучения, то в случае

- если (пГ2) < —6, то грань называется входной (луч входит в тетраэдр);

- если (пГ2) > е, то грань называется выходной (луч выходит из тетраэдра);

- иначе, грань называется касательной.

В этом определении е 1 — малое число, позволяющее различить входные и выходные грани устойчиво к ошибкам округлений. Доказывается лемма об устойчивой трассировке, в которой утверждается, что луч, проходящий через точку в многограннике Т в направлении —Л, гарантированно выйдет через входную (в смысле определения выше) грань при условии, что точка <3 находится на расстоянии более 8 = е - <1(Т) от каждой грани. Здесь с1(Т) обозначает диаметр многогранника Т. Иллюстрация данного условия приведена на рисунке 3.

Чтобы гарантировать условие леммы, при трассировке начальная точка (5 смещается внутрь многогранника. Это также решает проблему возможного зацикливания луча при прохождении через вершину, когда из-за ошибок округления точка <3 при трассировке через несколько шагов возвращается в тетраэдр, который был уже пройден. При использовании смещения точки <3, она гарантированно сдвигается хотя бы на е2й(Т) в направлении — Л. Для работы алгоритма

втв = -(П,п) < е п эт а = е

/

Рис. 3 — Многогранник Т, запрещенная область (серая) и гарантированно допустимая область из леммы (пунктир). В трехмерном случае запрещенная область имеет более сложную форму

требуется, чтобы пунктирная область на рисунке 3 была непуста, что накладывает определенные ограничения на качество расчетной сетки.

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

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

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

После решения системы линейных уравнений и определения интенсивности на границах подобластей можно провести трассировку внутренних то-

чек каждой подобласти в отдельности. Это действие производится автономно в каждой подобласти.

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

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

CPU

GPU _

I-1,-, I-1 Трассировка

□Подготовка Нс^раниРц

Рис. 4 — Временная диаграмма работы распределенного алгоритма метода длинных характеристик

■Решение ■СЛАУ

щ Пересылки s данных

Трассировка подобластей

Пятая глава

В главе содержатся результаты исследования построенных численных методов. В первом параграфе на примере задачи об излучающем оптически плотном шаре изучается сходимость метода, основанного на вариационном принципе. Использовалось шесть первых приближений метода сферических гармоник (с числом угловых функций К = 1,6,15,26,45,66 соответственно). Показано, что в областях гладкости решения пространственная сходимость имеет второй порядок, характерный для используемого метода Ритца. В окрестности границы шара, где коэффициент поглощения скачкообразно меняется и

градиент плотности энергии излучения имеет логарифмическую особенность, сходимость по пространственным переменным имеет порядок 0,4 (см. графики, представленные на рисунке 5). Порядок сходимости метода в норме ¿2 по количеству угловых функций составляет 0,29 (е ~ Л'-0-29), см. рисунок 6.

Рис. 5 — Ошибка сходимости е в центре шара (слева) и на краю шара (справа) в зависимости от мелкости сетки ./V-1/3

к

Рис. 6 — Ошибка £ в норме Lz в зависимости от используемого числа угловых

функций К

Метод сферических гармоник не обладает «эффектом луча», поскольку линейная оболочка набора сферических функций инвариантна относительно вращений пространства и не содержит выделенных направлений. Отмечается, что решения, полученные методом сферических гармоник, могут быть нефизичны. Например, тензор направленности излучения (тензор квазидиффузии) Daß в задаче об излучающем шаре содержит отрицательные компоненты в областях, удаленных от источника излучения (см. график на рисунке 7, в области, где Drr > 1 остальные главные компоненты отрицательны, поскольку tr Daß = 1).

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

задаче с излучением шара. Пунктиром приведено аналитическое значение

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

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

Рис. 8 — Вид излучающего тела в задаче для метода коротких характеристик

Рис. 9 — Решение в различные моменты марша вдоль одного направления. Показаны лишь те элементы, в которых решение уже известно

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

Рис. 10 — Решения вдоль направления к наблюдателю, полученные методами первого порядка, второго порядка и второго порядка с ограничителем

Сравнение плотности излучения показало, что методы первого и второго порядка с ограничителем дают решения, отличающиеся не более чем на 2%. В методе второго порядка более выражен «эффект луча», в то время как в методе первого порядка он сглажен численной диффузией лучей.

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

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

Рис. 11 — Решение вдоль одного направления с одной подобластью (слева) и с 8 подобластями (в центре). Справа представлены отдельные подобласти

В таблице 1 приведено сравнение времени работы алгоритма при различных разбиениях расчетной области. Данные результаты получены для сетки с 65 ■ 103 узлов (363 ■ 103 тетраэдров) с использованием 230 направлений из квадратурной формулы Лебедева (25-й порядок) на высокопроизводительном стенде кафедры информатики и вычислительной математики МФТИ. Незначительное ускорение обусловлено хаотическим доступом к памяти при трассировке лучей, плохой балансировкой нагрузки на графическом ускорителе, а также значительным количеством условных операторов в коде трассировки.

р Устройства ¿GPU, с ¿CPU-С ¿CPU /£ GPU

1 1 х Tesla С2075 117 439 3,75x

2 2 х Tesla С2075 52 196 3,77x

3 3 х Tesla С2075 32,7 125 3,82x

4 3 х Tesla C2075 + GTX 680 28,1 108 3,84x

8 3 x Tesla C2075 + GTX 680 25,5 69 2,70x

Таблица 1 — Время работы алгоритма в зависимости от числа подобластей Р для MPI реализации (CPU) и MPI+CUDA реализации (GPU)

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

порождает конический ветер, который можно обнаружить по синему смещению поглощения в линии Н-а. В этом расчете вместо квадратурной формулы Лебедева использовалась квазиравномерная угловая сетка из 1091 направления. Это позволило получить усредненные по периоду обращения звезды спектры при различных углах наклонения г (угол между плоскостью диска и лучом зрения). Использовалось 64 частотные группы в интервале, соответствующем доплеровскому смещению линии Н-а от —600 км/с до 600 км/с. Зависимость усредненной за период обращения звезды относительной интенсивности /0 от частоты излучения при разных углах наклонения г представлена на рисунке 12. Для удобства смещение частоты Дг/ переведено в относительную скорость Дг> = с^. Можно заключить, что при углах наклонения г < 45° поглощение в диапазоне синего смещения Ау ~ 150-^-200 км/с существенно. Нормализованный профиль линии излучения имеет существенный провал в области красного смещения на скорости Дъ> и 100 км/, который объясняется значительным поглощением излучения веществом диска, аккрецирующего на звезду.

Рис. 12 — Профиль линии Н-а в зависимости от угла наклонения i

Для расчета задачи о спектре излучения использовался кластер лаборатории флюидодинамики и сейсмоакустики МФТИ, состоящий из двух узлов, на каждом из которых имеется 8 ускорителей Tesla K40m. На сетке из 5 • 10э узлов (около 3 ■ 10(l тетраэдров) при использовании 110 направлений и 64 частотных групп удалось получить ускорение в 4,5 раза при использовании 16 графических ускорителей по отношению к MPI реализации с 16 процессами.

В заключении приведены основные результаты работы.

Основные результаты диссертации

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

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

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

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

5. Разработанные вычислительные алгоритмы реализованы в виде программного комплекса. В рамках модели локального термодинамического равновесия вычислен коэффициент поглощения частично ионизованной плазмы. Для задачи моделирования спектра излучения звезды типа Т Тельца построен спектральный профиль линии Н-а в зависимости от ориентации плоскости аккреционного диска.

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

1. Маршевый алгоритм решения задачи переноса излучения методом коротких характеристик / Ю. И. Скалько, И. В. Цыбулин, Р. Н. Карасев и др. // Компьютерные исследования и моделирование. — 2014. — Т. 6, № 2. — С. 203-215.

2. Цыбулин И. В., Скалько Ю. И., Павлова Е. С. Распределенный матод длинных характеристик для решения уравнения переноса излучения // Труды Московского физико-технического института. — 2015. — Т. 7, № 2. — С. 51—59.

3. Цыбулин И. В., Скалько Ю. И. Конечно-элементный метод для решения уравнения переноса излучения // Математическое моделирование информационных систем: сб. науч. тр. — М.: МФТИ, 2015. — С. 38-44.

4. Цыбулин И. В. Решение систем обыкновенных дифференциальных уравнений большой размерности с использованием графических ускорителей // Современные проблемы фундаментальных и прикладных наук. Часть VII. Управление и прикладная математика: Труды 53-й научной конференции МФТИ. — Т. 3. — М.: МФТИ, 2010. — С. 33-34.

5. Цыбулин И. В., Сколько Ю. И. Численное построение квадратурной формулы для полусферы методом продолжения по параметру // Проблемы фундаментальных и прикладных естественных и технических наук в современном информационном обществе. Часть VII. Управление и прикладная математика: Труды 54-й научной конференции МФТИ. — Т. 2. — М.: МФТИ, 2011.

- — С. 27-28.

6. Цыбулин И. В., Сколько Ю. И. Вариационный метод решения задачи переноса излучения в трехмерных областях // Проблемы фундаментальных и прикладных естественных и технических наук в современном информационном обществе. Часть VII. Управление и прикладная математика: Труды 55-й научной конференции МФТИ. — Т. 2. — М.: МФТИ, 2012. — С. 96-97.

7. Цыбулин И. В., Сколько Ю. И. Использование метода длинных характеристик для решения задачи переноса излучения на системах с распределенной памятью // Актуальные проблемы фундаментальных и прикладных наук в современном информационном обществе. Часть VII. Управление и прикладная математика: Труды 56-й научной конференции МФТИ. — Т. 2. — М.:

МФТИ, 2013,—С. 89-90.

8. Цыбулин И. В., Сколько Ю. И. Реализация распределенного метода длинных характеристик с использованием GPU ускорителей // Актуальные проблемы фундаментальных и прикладных наук в области физики. Часть VII. Управление и прикладная математика: Труды 57-й научной конференции МФТИ. — Т. 2. — М.: МФТИ, 2014. — С. 103-105.

Личный вклад автора в публикации с соавторами.

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

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

В работах [3,5-8] автором построены численные методы для решения уравнения переноса, выполнена алгоритмическая и программная реализация на графических ускорителях методов численного решения систем линейных алгебраических и систем обыкновенных дифференциальных уравнений. Сформулирована задача построения квадратурной формулы для полусферы, которая решена предложенным численным методом продолжения по параметру.

Цыбулин Иван Владимирович

РАЗРАБОТКА ЧИСЛЕННЫХ МЕТОДОВ ДЛЯ РЕШЕНИЯ УРАВНЕНИЯ ПЕРЕНОСА ИЗЛУЧЕНИЯ И ИХ РЕАЛИЗАЦИЯ С ИСПОЛЬЗОВАНИЕМ ГРАФИЧЕСКИХ УСКОРИТЕЛЕЙ

Автореферат

Подписано в печать 02.10.2015. Формат 60 х 84 У16. Усл. печ. л. 1,3. Тираж 100 экз. Заказ №411. Федеральное государственное автономное образовательное учреждение высшего профессионального образования «Московский физико-технический институт (государственный университет)» Отдел оперативной полиграфии «Физтех-полиграф» 141700, Московская обл., г. Долгопрудный, Институтский пер., 9