CC BY
31
УДК 519.63
И. В. Цыбулин, Ю. И. Сколько, Е. С. Павлова
Московский физико-технический институт (государственный университет)
Распределенный метод длинных характеристик для решения уравнения переноса излучения
В работе предлагается модификация метода длинных характеристик, используемого для решения задачи переноса излучения. Предлагаемый метод позволяет решать уравнение переноса на декомпозированных неструктурированных тетраэдральных сетках. Коллективное взаимодействие процессов требуется только для решения задачи на границах области. Благодаря тому, что в основе метода лежат длинные характеристики, вспомогательная задача допускает простое представление и сводится в решению небольшой разреженной системы линейных уравнений.
Ключевые слова: перенос излучения, метод длинных характеристик, неструктурированные сетки, распределенная реализация, графические ускорители.
I. V. Tsybulin, Y. I. Skalko, E. S. Pavlova
Moscow Institute of Physics and Technology (State University)
A distributed version of the long characteristics method for the equation of radiative transfer
A modified version of the long characteristics method used for the radiative transfer problem is proposed. The method is suitable for solving the equation on a decomposed unstructured tetrahedral grid. Interprocess communication is required only for solving an auxiliary problem on the domain interfaces. Utilization of the long characteristics admits a simple representation of the auxiliary problem which reduces to solving a small sized sparse linear system of equations.
Key words: radiative transfer, long characteristics method, unstructured grid, distributed implementation, GPU.
Введение
Уравнение переноса излучения является составляющей частью системы уравнения динамики излучающего газа. Излучение может играть определяющую роль в высокотемпературных газодинамических процессах. Без учета излучения невозможно моделирование сильноточных разрядов, например, рентгеновских источников на основе z-пинчей, проектирование теплозащиты спускаемых аппаратов. Излучение играет ключевую роль в астрофизике при решении задач динамики звездных фотосфер, а также спектр излучения является важнейшей характеристикой, с помощью которой возможно сравнивать астрофизические наблюдения с теоретическими моделями.
Одним из существенных отличий уравнения переноса излучения от остальных уравнений системы радиационной газовой динамики является его многомерность. Даже в простейшем стационарном одночастотном случае уравнение переноса содержит в дополнение к трем пространственным измерениям еще два угловых измерения. Эта особенность свойственна всем кинетическим уравнениям, описывающим функцию распределения, в данном случае такой функцией является интенсивность излучения.
Довольно часто при решении уравнения переноса излучения совместно с вычислительными комплексами для моделирования динамики излучающего газа тип расчетной сетки и способ ее декомпозиции на подобласти уже определены. Обычно такая декомпозиция производится геометрическим образом и отражает локальность уравнений балансов основных
физических величин. В уравнении переноса излучения данная локальность отсутствует, в связи с чем даже весьма удаленные участки области могут напрямую взаимодействовать, обмениваясь излучением вдоль соединяющих их прямой. Эта проблема осложняет численное решение уравнения переноса, так как взаимодействующие участи могут находиться в различных подобластях, а соединяющий их луч может проходить через несколько других подобластей.
1. Математическая постановка
Поле излучения можно описать функцией интенсивности излучения I, которая с точностью до умножения на энергию фотона совпадает с функцией распределения фотонов по координатам, направлениям и частотам. Для случая нерелятивистских течений газа можно пренебречь нестационарностью процесса переноса излучения, то есть считать, что излучение распространяется от источника мгновенно. Пренебрежем также рассеянием излучения и изменением частоты фотонов. Таким образом, процесс переноса можно рассматривать независимо в каждом направлении и на каждой частоте излучения.
Перенос излучения в стационарном одночастотном приближении описывается уравнением [1]
(«V)/(х, О) + к(х)1(х, О) = к(х)1р(х),
где О — единичный вектор направления излучения, к(х) — коэффициент поглощения излучения в точке х, учитывающий вынужденное испускание, 1р(х) — интенсивность равновесного излучения в точке х.
Рассмотрим излучение в выпуклой области С с граничным условием отсутствия источников излучения во внешности области С. Граничное условие для интенсивности принимает вид
I(х, О) = 0, х е дС, (п, О) < 0,
то есть интенсивность излучения, входящего в область через границу, равна нулю.
Данное уравнение переноса является гиперболическим и имеет простое семейство характеристик прямых:
х — хо = О(в — во).
Вдоль каждой характеристики уравнение превращается в обыкновенное дифференциальное уравнение:
^+ К(8)1 (8) = к(8)1р (8) и может быть легко проинтегрировано:
I(8) = I(8о)е-+ [а к(01р(0е-т(1)
J «о
В последнем выражении т(а,Ь) = ^ к(в) йв — оптическая толщина характеристики от точки 8 = а до точки в = Ь.
Выражение (1) лежит в основе метода длинных характеристик. Чтобы вычислить интенсивность в точке х в направлении О, из точки х выпускается луч в направлении —О до пересечения с границей области в точке хо (см. рис. 1). Так как в точке хо интенсивность входящего в область излучения известна из граничных условий, вычисление интенсивности в точке х в направлении О сводится к вычислению интегралов вдоль отрезка.
Заметим, что интенсивность 1(8) в формуле (1) является линейной функцией интенсивности на границе 1(8о). Введем
а(8о, 8) = е-т(зо'з), Р(8о, 8)= 13 к(01р(0а((, 8)<%.
' «о
Т(50,5)
При 5 > 5о для величин , $), Р(«о, 5) справедливо
0 < а(во,в) < 1,
0 ^ р(зо, з) ^ тах 1р(£) • Г к(£)е-= тах 1р(£) • [ е-и¿и = (2)
Js0 С ^[«0,з] у о
= (1 - а(во,в)) тах 1р(£).
С ^0,5]
В этих обозначениях выражение для интенсивности (1) принимает особенно простой вид:
I(з) = а(во, (во) + Р(¿о, 5). (3)
С учетом (2) для I(в) верен принцип максимума:
I(в) < а(во, в)1 (во) + (1 - а(^о, «)) тах 1р(£) < тах| I(«о); тах 1р(£Н .
Лемма 1. Заметим, что функции а и /3 удовлетворяют следующему трассировочному соотношению при 5о ^ 51 ^ 5:
а(во, в) = а(во, 51 )а(^1, в), (4)
/3(во, в) = Р(5о, 51 )а(^1, 5) + р(51, 5).
Данное соотношение позволяет найти значения коэффициентов а, Р для целого отрезка, если известны коэффициенты для двух его частей.
Доказательство. Из аддитивности оптической толщины т(5о,5) = т(5о,в1) + т(в1 ,в) следует первое соотношение
а(5о ,5) = е-т (30 ^ = е-т (30,31 )-т ^ = е-т (з0,3!) е-т ^ = а(5о ,51 )а(*1,5).
Разобьем интеграл в Р(5о, 5) на сумму двух:
Гв 1 гв
р (зо ,з)= к а )1р (£ )а(£,з№ + / к (Щ )1р (£ )а(£,8)<%.
«/ 50 «/в!
Второй интеграл в точности является Р(в1,5). В первом интеграле воспользуемся доказанным свойством: а(£, в) = а(£,в1 )а(в1 ,в). Его можно применить, так как £ ^ 51. Величина а(^1,5) является константой в первом интеграле и может быть вынесена за его пределы:
Рв!
Р(5о, 5) = а(^1, в) к(£)1р(£)а(£, 51)^ + Р(51, 5) = Р(во, 51 )а(в1, 5) + Р(51, 5). □ ■3 50
2. Предлагаемый метод
Реализация метода длинных характеристик осложняется, если вычислительная область разбита на подобласти. Обычно каждая подобласть принадлежит своему вычислительному процессу и трассировка характеристики становится сложной коллективной операцией.
Заметим, что формула (3) позволяет не только вычислить неизвестное значение 1(в) через известное значение 1( в о), но и просто связать линейным соотношением два значения интенсивности в двух точках на одном луче. Если интенсивность излучения каким-то образом оказалась известной на границе вычислительных подобластей, то трассировку лучей достаточно провести в каждой подобласти, независимо от других. В этом случае точка хо — это точка, в которой характеристика пересекает границу подобласти, а не всей вычислительной области (см. рис. 2).
Рис. 2. Трассировка луча до границы подобласти
Предположим, что границы подобластей проходят по граням тэтраэдральной сетки, построенной в области. Зафиксируем направление О. Будем говорить, что узел сетки (вершина тетраэдра) с координатой х принадлежит подобласти, если точка х — Ог лежит в подобласти при £ ^ +0. Если таких подобластей несколько, например, если вектор О лежит в плоскости грани, разделяющей подобласти, выберем подобласть с минимальным номером. Заметим, что определенное таким образом отношение принадлежности зависит не только от геометрических параметров сетки, но и от направления О также. Если узел не принадлежит ни одной подобласти, будем говорить, что он принадлежит внешней подобласти, то есть М3 \ С. Заметим, что интенсивность в граничных узлах, принадлежащих внешней подобласти, задана граничными условиями.
Для каждого граничного узла с координатой х, принадлежащего подобласти С С, выпустим характеристику до пересечения с границей подобласти . В общем случае характеристика выйдет из области через некоторую грань / в точке хо = х — О 5. Тогда, пользуясь соотношением (3), заключаем, что
1(х) =а(0, з)1(хо )+Р (0, 8).
Однако точка х является узлом расчетной сетки, а точка хо в общем случае — нет. Вычислим 1(хо) приближенно интерполяцией интенсивности по вершинам грани /:
к
1(хо) I(xj),
3=1
где К — число вершин грани / (для тетраэдральных сеток К = 3), а х3- — координаты вершин этой грани, которые, в отличие от хо, являются узлами вычислительной сетки. Величины 73 являются обобщенными (для К > 3) барицентрическими координатами точки хо в грани /. При этом £К=1 73 = 1.
Таким образом, получается система линейных уравнений, в которые входят значения интенсивности только на границах подобластей и не входят значения внутри подобластей. Уравнения записаны для каждого граничного узла, принадлежащего некоторой вычислительной подобласти. В узлах, принадлежащих внешней подобласти, интенсивность уже задана граничными условиями.
Отметим, что данная система уравнений имеет вид
и обладает диагональным преобладанием, поскольку на диагонали стоит коэффициент 1, а сумма модулей внедиагональных коэффициентов не превосходит 1 — а(0, з) < 1 (было учтено, что коэффициенты имеют один знак, но некоторые из неизвестных I(х) могли быть исключены, если в х задано граничное условие).
После коллективного решения системы в узлах на границах подобластей становится известными значения интенсивности. Далее, выпуская характеристики из внутренних точек подобластей, восстанавливается решение в каждом узле каждой подобласти всей расчетной области. Эта операция не является коллективной и проводится автономно в каждой подобласти. Также она является самой трудоемкой.
Необходимо подчеркнуть, что решение, полученное данным методом, отличается от решения, полученного методом длинных характеристик. В частности, решение зависит от способа декомпозиции на подобласти. Предельный случай декомпозиции области на одну подобласть соответствует методу длинных характеристик, а гипотетический случай декомпозиции на отдельные тетраэдры приводит к методу, где характеристика выпускается от одной грани тетраэдра до другой его грани, то есть интерполяция производится на каждой грани, а не только на границе подобласти, как в предлагаемом методе. Такой метод называется методом коротких характеристик [2]. Увеличение количества областей приводит к более частой интерполяции, которая, в свою очередь, вызывает численную диффузию луча. С другой стороны, при этом средняя длина трассируемых лучей уменьшается, что существенно снижает время работы алгоритма.
3. Вычислительный алгоритм
Уравнение переноса решается в некотором наборе направлений. В качестве такого набора были выбраны узлы квадратурной формулы для сферы Лебедева [3], что позволяет получать значения интегральных характеристик (плотности и потока излучения), суммируя найденные значения интенсивности с весами Wk той же квадратурной формулы:
Для каждого направления О граничные узлы каждой вычислительной подобласти разбиваются на те, которые ей принадлежат (в смысле определения, данного выше) и которые не принадлежат. Из тех, что принадлежат, производится трассировка лучей вдоль направления —О до пересечения с границей подобласти.
Рассмотрим алгоритм трассировки подробнее. Трассировка начинается из некоторой точки х, которая является вершиной тетраэдра Т. Выбирается любой тетраэдр, содержащий в себе точку х — Ое при е ^ +0. Далее определяются барицентрические координаты точки пересечения луча х — О,в, в > 0, со всеми плоскостями, содержащими грани тетраэдра. Определяется грань с минимальным значением суммы модулей барицентрических координат. В точной арифметике сумма модулей координат равна единице в случае, если луч пересекает грань, и больше единицы, если проходит вне грани. Далее определяется тетраэдр, смежный с данным по найденной грани, и трассировка повторяется уже в нем из
(5)
з=1
точки, в которой луч пересекает грань. Если смежный тетраэдр отсутствует, трассировка завершена — достигнута граница подобласти. Барицентричекие координаты запоминаются и используются в качестве весов интерполяции ^^ в уравнениях (5).
Во избежание зацикливания алгоритма, например, когда луч проходит в малой окрестности вершины, после попадания луча в тетраэдр точка входа смещается от границы в направлении центра масс тетраэдра на расстояние порядка 6 = 10-6 диаметра тетраэдра.
£
Рис. 3. Трассировка луча по триангуляции в подобласти
Одновременно с трассировкой луча при помощи соотношений (4) производится вычисление величин а(0, 8), /(0, 8). В начале трассировки полагается а := 1, / := 0. При прохождении очередного тетраэдра для него вычисляются значения а = а(80, 81), 3 = /3(8о, 81), где ^ о, — координаты точки выхода и точки входа вдоль луча (см. рис. 3). В предположении однородности параметров к, 1р в тетраэдре можно использовать явные выражения
а = е-к£, 3 = (1 - а) 1р,
где I — длина отрезка луча, заключенного в тетраэдре. После этого значения а, /3 модифицируются по правилу:
/ := / + /За, а := а а.
К моменту завершения трассировки луча значения а и / равны соответственно а(0, 8) и /(0, 8).
Таким образом, на каждый процесс формирует в своей подобласти часть общей системы уравнений, в которые входят неизвестные значения интенсивности на границе области. Далее эта система собирается воедино на одном из процессов, где решается прямым методом для разреженных матриц. Выбор такого метода вместо итерационного обусловлен двумя факторами. Во-первых, размер системы существенно меньше общего количества узлов сетки, так как неизвестные находятся лишь на границах подобластей. Во-вторых, матрицу системы можно превратить в близкую к треугольной, если упорядочить неизвестные по значению их проекции на луч П.
Чтобы эффективно использовать остальные процессы, пока на одном из них происходит решение системы, направления обрабатываются не по одному, а по несколько за раз. Это позволяет загрузить все процессы решением систем линейных уравнений, причем каждый процесс будет решать систему для своего направления.
Алгоритм решения состоит из следующих этапов:
1) Выбирается очередной пакет направлений по одному на процесс.
2) Каждый процесс определяет принадлежащие ему точки и производит трассировку от границ в каждом направлении из пакета.
3) На каждом процессе собирается система уравнений для его направления из пакета.
4) Система решается на каждом процессе и решение рассылается по процессам.
5) В каждой подобласти производится вычисление решения во всех внутренних узлах автономно от остальных процессов.
6) Если остались направления, в которых решение еще не строилось, алгоритм повторяется с первого пункта.
Для решения разреженной системы линейных уравнений применялась библиотека иМРРАСК.
Так как пятый этап алгоритма - самый трудоемкий, он был реализован с использованием графических ускорителей. Это позволило не только сократить время расчета задачи, но и произвести перекрытие подготовительных этапов (2-4), производимых на центральном процессоре с этапом 5, производимым на графическом ускорителе (см. рис. 4).
Рис. 4. Временная развертка алгоритма
4. Результаты
Данный алгоритм был реализован и тестировался на модельной задаче с излучающим оптически плотным шаром в оптически тонкой области. Использовалась сетка с 66 тысячами узлов (363 тысячи тетраэдров) и 230 направлений из квадратурной формулы Лебедева 25-го порядка.
В случае нескольких подобластей метод имеет незначительную диффузию луча, проходящего через поверхность раздела подобластей. На рис. 5 приведены решения для случая одной подобласти и восьми подобластей. Разбиение на подобласти (см. рис. 6) производилось с помощью библиотеки METIS.
Рис. 5. Срез решения в одном направлении. Слева Рис. 6. Используемая декомпозиция об-
одна подобласть, справа — 8 ласти на подобласти
В экспериментах наблюдалось сверхлинейное ускорение алгоритма (см. таб. 1), которое можно обяснить различной сложностью задачи в зависимости от числа подобластей. То есть данный метод не является распараллеленной версией некоторого последовательного алгоритма, а параллельным алгоритмом, результат работы которого напрямую зависит от числа используемых процессов. Реализация с использованием графических ускорителей (Nvidia Tesla C2075) показала незначительное ускорение, что, видимо, объясняется хаотичностью доступа к памяти в алгоритме трассировки и недостаточной его вычислительной
сложностью. В эксперименте с 8 подобластями использовались только 4 графических ускорителя.
Т а б л и ц а 1
Время работы алгоритма
Число процессов tори С ¿СРИ С tcPlJ/¿ори
1 117 439 3.75х
2 52 196 3.77х
3 32.7 125 3.82х
4 28.1 108 3.84х
8(4) 25.5 69 2.70х
Так как рассматриваемая задача сферически симметрична, хорошим тестом вычислительного метода является тест на сферичность изоповерхностей решения (см. рис. 7). Видно, что искажения в направлениях координатных осей довольно сильны. Так как используется неструктурированная пространственная сетка, эти искажения могут быть связаны только с угловой сеткой, в которой эти направления выделены. В данном методе, как и в методе длинных характеристик, сильно выражен «эффект луча», или численная анизотропия переноса излучения в различных направлениях. Для этой задачи малая численная диффузия луча является скорее недостатком метода. Для сглаживания «эффекта луча» требуется брать больше направлений.
Рис. 7. Изоповерхности плотности излучения
5. Заключение
В работе предложена модификация метода длинных характеристик для решения уравнения переноса излучения. Метод допускает эффективную реализацию на вычислительных системах с распределенной памятью и может использоваться совместно с имеющимися программами для численного газодинамического моделирования.
Самая трудоемкая часть алгоритма была распараллелена с использованием графических ускорителей. Однако полученное ускорение свидетельствует о плохой пригодности графических ускорителей для этой задачи. Вероятно, причиной низкого ускорения является хаотический доступ к памяти при трассировке лучей и недостаточная вычислительная сложность трассировки.
К основным недостаткам метода можно отнести существенно больший объем вычислений, по сравнению, например, с методом коротких характеристик. Данный недостаток
можно преодолеть, если использовать внутри области приближение ближайшего луча, выпущенного из точки на границе.
Работа выполнена при финансовой поддержке проекта Министерства образования и науки РФ №3.522.2014/К в лаборатории флюидодинамики и сейсмоакустики МФТИ.
Литература
1. Зельдович Я. Б., Райзер Ю.П. Физика ударных волн и высокотемпературных гидродинамических явлений. М.: Физматгиз, 1963. 632 с.
2. Четверушкин Б. Н. Математическое моделирование задач динамики излучающего газа. М.: Наука, 1985. 304 с.
3. Лебедев В. И., Лайков Д. Н. Квадратурная формула для сферы 131-го алгебраического порядка точности // Докл. РАН. 1999. Т. 366, № 6. С. 741-745.
References
1. Zeldovich Ya. B., Raiser Yu. P. Physics of shock waves and high-temperature hydrodynamic phenomena. V. 1. N.-Y.: Academic Press, 1966. 464 p.
2. Chetverushkin B.N. Mathematical modelling in radiative gas dynamics problems. M.: Nauka, 1985. 304 p. (in Russian).
3. Lebedev V. I., Laikov D. N. A quadrature formula for the sphere of the 131st algebraic order of accuracy. Doklady Mathematics. 1999. V. 53, N 3. P. 477-481.
Поступила в 'редакцию 26.06.2015.
