CC BY
45
УДК 519.688:523.035.1
0 развитии двух параллельных алгоритмов численного моделирования взаимодействия излучения с веществом
1 2 3
Л. И. Рубина , О. Н. Ульянов , М. А. Чащин
1 [email protected], 2 [email protected], 3 [email protected]
1’2’3 ФГБУН Институт математики и механики им. Н. Н. Красовского УрО РАН (ИММ УрО РАН)
Поступило в редакцию 22.12.2013
Аннотация. Представлена методика численного моделирования взаимодействия ионизирующего излучения с веществом для ряда физических и математических постановок задачи радиационного переноса в смесях веществ. Методика предполагает расчет вклада каждой учитываемой спектральной линии; вычисление населенностей уровней в соответствии с уравнениями кинетики населенностей; нахождение интенсивности излучения в соответствии со спектральным уравнением переноса с коэффициентами, вычисляемыми по населенностям уровней и электронной температуре, которая удовлетворяет уравнению энергобаланса. Решение задачи в достаточно полных постановках сопряжено с большими объемами вычислений и требует применения параллельных технологий. Авторы разрабатывают методику, основанную на двух различных методах решения задачи и на развитии двух параллельных алгоритмов. Предложенные алгоритмы и разработанные программы для суперкомпьютеров позволили проводить моделирование переноса излучения в плоских слоях вещества с учетом до тысячи спектральных линий за приемлемое время с приемлемой точностью. Приводятся результаты численных расчетов.
Ключевые слова. Перенос излучения; численное моделирование; параллельные алгоритмы
ВВЕДЕНИЕ
Перенос излучения и процессы, связанные с переносом излучения в однородных и неоднородных средах, вызывают большой интерес исследователей. Численное моделирование и исследование этих процессов используются как при изучении многих явлений природы, так и при создании новой техники [ 1-4]. Детальное изучение процесса переноса излучения и связанных с ним процессов - непростая задача, требующая учета многих тонких эффектов и приводящая к проведению большого объема вычислений, к обработке больших многомерных массивов информации.
Статья рекомендована к публикации программным комитетом Международной научной конференции «Параллельные вычислительные технологии 2012». Работа выполнена в рамках программы межрегиональных и межведомственных фундаментальных исследований УрО РАН (проект 12-С-1-1001) и программы фундаментальных исследований Президиума РАН «Информационные, управляющие и интеллектуальные технологии и системы» (проект 12-П-1-1023) при поддержке УрО РАН
Разработкой методов, алгоритмов и программ для численного моделирования радиационного переноса занимаются многие исследователи (например, [5]).
Авторы на протяжении ряда лет проводили работы в указанном направлении, основанные на реализации двух развиваемых ими методов (МАПИ и МПЛЧ) численного решения задачи. Метод МАПИ (метод аналитического представления излучения) использует явное представление интенсивности излучения, метод МПЛЧ (метод полиномов Лагранжа-Чебышева) основан на использовании интерполяционных полиномов Лагранжа по узлам Чебышева. Методика разрабатывается для случая взаимодействия ионизирующего излучения с веществом, состоящим из водородоподобных и гелиоподобных ионов, а также ионов, у которых на внешней оболочке отсутствуют электроны.
На первом этапе изучалась такая физическая модель явления, которую адекватно представляет случай, когда рассматриваемые вещества имеют доплеровские профили излучения и поглощения спектральных линий. В качестве математической модели была взята система стационарных уравнений, включающая систему уравнений кинетики населенностей уровней
и систему интегро-дифференциальных уравнений, состоящую из уравнения переноса излучения в «непрерывном спектре» и уравнений переноса излучения в каждой из учитываемых спектральных линий с соответствующими граничными условиями [6]. Пересечением спектральных линий пренебрегалось, энергобаланс не учитывался. При наличии доплеровских контуров было возможно не учитывать пересечения резонансных линий. Первое предположение — отсутствие пересечения спектральных линий вполне естественно, что и показали дальнейшие исследования. Со вторым — не учетом энергобаланса — дело обстоит, очевидно, не так. Но, как показали дальнейшие исследования, в некоторых случаях можно изучать явление, пренебрегая энергобалансом [7, 8], хотя, конечно, более полные модели должны содержать уравнение энергобаланса для электронной температуры. Была разработана методика для численного моделирования процессов переноса излучения в неоднородном плоском слое смеси веществ, содержащем несколько подслоев, имеющих разные плотности и температуры, при учете до ста резонансных линий за приемлемое время с приемлемой точностью [9—11]. Развитие методики состояло как в развитии оригинальных, так и в подборе известных, но отвечающих особенностям задачи математических и вычислительных методов решения задачи, совершенствовании и модификации численных алгоритмов и реализующих их программ.
Благодаря использованию параллельных технологий и все более мощных суперкомпьютеров удалось решать задачи в достаточно полных постановках. В настоящее время возможно рассчитывать интенсивность излучения, учитывая не несколько десятков, а сотни линий. Слой, через который проходит излучение, может содержать до десяти подслоев, отличающихся как составом вещества, так и плотностью. Электронная температура в слое меняется согласно уравнению энергобаланса или может быть задана.
В данной статье представлены результаты развития методики при переходе к физическим моделям с фойгтовскими профилями излучения и поглощения в спектральных линиях. Этот класс задач с физической точки зрения является существенно более содержательным, чем ранее решавшиеся задачи. С математической и вычислительной — значительно более сложным. Переход к таким моделям потребовал пересмотреть математическую модель. Ранее в случае допле-
ровских контуров было возможно не учитывать пересечение резонансных линий. Рассмотрение моделей с фойгтовскими профилями сделало этот учет необходимым, что привело к отказу от раздельного вычисления интенсивности излучения в линиях и в непрерывном спектре и, следовательно, к разработке нового подхода к решению единого интегро-дифференциаль-ного
уравнения переноса излучения.
Кроме того, на этом этапе развития методики и применения новых параллельных технологий для современных суперкомпьютеров была поставлена задача оценки влияния учета энергобаланса на решение в различных постановках. Для более полного описания физического процесса в математической модели температуру требуется определять, используя информацию
0 спектре излучения в среде. Таким образом, при решении задачи необходимо совместно с решением уравнения переноса излучения (УП) и системы уравнений кинетики (УК), описывающих долю ионов веществ смеси, пребывающих в определенном состоянии, которое называется населенностью уровня, решать уравнение энергобаланса (УБ) для электронной температуры вещества. Результаты развития методики представлены в данной работе.
ПОСТАНОВКА ЗАДАЧИ
Рассматривается стационарный случай задачи переноса излучения в плоском слое смеси веществ. Слой конечной толщины Ь может состоять из отдельных подслоев От (1 < т < М, М < 10). В каждом подслое смесь может обладать своими характеристиками (состав, плотность). Смесь состоит из Г компонент, у — номер компоненты. Атомы каждой компоненты могут быть в нескольких состояниях с разной степенью возбуждения и ионизации. Каждое состояние определяется парой индексов (/7), где
1 — степень ионизации, / — номер (дискретного) возбужденного состояния для данной степени ионизации.
Относительная доля ионов у-го вещества, находящегося в состоянии (/I), обозначается \х) и называется населенностью этого состояния. Все населенности для каждой конкретной компоненты вещества можно считать некой векторной величиной С(У) (х) = (с1у) (х),
с2у)(х),...,С;у){х}) . Они удовлетворяют системе уравнений кинетики:
Ж(у) с(у) = Ь(у). (1)
Элементы кинетических матриц Ж(г)(х) и компоненты вектора Ь(:1\х) зависят от средних по углу и энергии интенсивностей излучения в каждой точке сетки по пространственной переменной. Кроме того, имеется условие нормировки ^ с(Ц) = 1 .
и
Вещества в смеси связаны через поле излучения. Интенсивность излучения зависит от векторов населенностей для каждого вещества в смеси. Если І(х, ц, е, Те) — интенсивность излучения вдоль некоторого луча в точке среды х с энергией е и электронной температурой Те, то в среде выполняется уравнение переноса излучения (ц — косинус угла между направлением луча и осью х, ось х направлена перпендикулярно поверхности слоя):
т — +КІ = 5. (2)
йх
Здесь К = к(ц, е, Те, С((р ) — коэффициент поглощения, 5 = 5(ц, е, Те, с\р) — функция источника, цє [-1,1], е є (-да, да), хє [0,Х]; если ц > 0, І|х=о = Іо(ц, е); если ц < 0, І\х=ь = Іі(ц, е).
Будем рассматривать случай, когда Те = Те(х) удовлетворяет уравнению энергобаланса:
- йх (Кт: ——Т-) = 2р|_+1йт 0К/ - 5 ]йе. (3)
Здесь к0 — коэффициент, являющийся функцией эффективного заряда атомного ядра смеси, п — показатель молярной теплоемкости.
Ранее для доплеровского контура излучения рассматривалось уравнение энергобаланса, когда к0 = 0. В настоящей работе приведены результаты решения системы уравнений (1—3) для случая фойгтовского профиля излучения и произвольного коэффициента к0 Ф 0.
МЕТОД РЕШЕНИЯ УРАВНЕНИЯ ЭНЕРГОБАЛАНСА
В уравнении энергобаланса (3) сделаем замену, полагая Q = Те(п+1). Для определения функции Q получаем уравнение:
й 2Q +1 +м
—— = -[2я(п +1)/к0]| йц | (КІ - 5)йе. (4)
—х і
— 1 — ¥
Сравнивая (2) и (4), получаем уравнение:
— = [2р(п + 1)/ к0]Г+1тф [+ і—е, (5)
йх •м •'-¥
к0 Ф 0 .
Решение уравнения (5) может быть сведено к последовательному вычислению трех интегралов
+1 +¥
Зт (х, е) = 17|1ф, 1Ш(х) = | Jj.de,
-1
х/+1
Q(хм) = [2р(п +1)/к0] 11м(x)dx,
х/
(/ = 0, 1, ... , /),
если известно Те(х0), то Те (х/+1) = Q1/{”+1)(х/+1),
{х0, х;, ., хг} — сетка по пространственной переменной.
В методике МАПИ для вычисления инте-
+1
грала JT (х, е) = 11у^т отрезок [-1,+1] разбива-
-1
ется на пять частей: [-1, ц1], [ць ц2], [ц2, ц3], [ц3, ц4], [ц4, 1]. Точки ц и ц совпадают, как правило, с точками разрыва заданной начальной интенсивности излучения. Отрезок [ц2, ц3] содержит точку ц = 0. Здесь интенсивность излучения имеет свои особенности, связанные с тем, что точка ц = 0 — особая точка уравнения переноса. Чтобы учесть особенности поведения подынтегральной функции на этом участке, требуется достаточно частая сетка. Для вычисления интегралов на каждой части отрезка используется полином Лежандра своей специальным образом выбранной степени.
Для вычисления несобственного интеграла
7и/(х) = | JTde используется метод Лобатто,
позволяющий варьировать пределы интегрирования и достигать в расчете заданной высокой гарантированной точности на каждом интервале, предварительно полученном после выделения точек разрыва подынтегральной функции.
Для вычисления значений Q(x/+1) в узлах сетки по пространственной переменной также используются полиномы Лежандра.
В методике МПЛЧ все интегралы вычисляются построением полиномов Лагранжа специально выбранной степени с узлами Чебышева. Совпадение результатов расчетов по двум значительно отличающимся алгоритмам расчета принимается за критерий истины.
О РАЗВИТИИ МЕТОДИКИ МПЛЧ
Метод МПЛЧ (метод полиномов Лагранжа-Чебышева) существует в нескольких модификациях.
Программа £^1;_псот39_8 осуществляет расчеты вариантов с фойгтовскими профилями для смесей, содержащих до четырех компонент
в неоднородных слоях (10 слоев). Слои могут отличаться плотностью, температурой (постоянной) и составу вещества, максимальное число линий — 39.
Программа Ьа1480в1 рассчитывает варианты с фойгтовскими профилями, смеси могут содержать до 4 компонент, неоднородные по плотности и составу вещества слои (10 слоев), решается уравнение энергобаланса. Размерности массивов по пространственной и по энергетической переменным могут изменяться от варианта к варианту. Выбор размерности зависит от того, какую задачу надо решать (либо много слоев, либо много линий, максимально 39 линий и 10 слоев). Ограничение на количество точек в сетке по пространственной переменной (480 точек на все подслои и 120 точек на одном слое).
Для задач с учетом нескольких десятков линий предложен и реализован в программных кодах новый высокоэффективный алгоритм, на порядок уменьшено время расчета одного варианта (сравните времена счета варианта на одинаковом числе процессоров в табл. 1 и 2). В табл. 1 и 2 представлены результаты расчета одинакового варианта: 2 линии, 2 компоненты вещества в слое, фойгтовский профиль. Результат, полученный по старой методике (табл. 1, 2008 г.). Результат, полученный по усовершенствованной программе метода МПЛЧ (табл. 2, 2010 г.).
Т аблица 1
N T(N) (с) H(N) %
2 14,7 1,00 100
4 7,9 1,86 93
8 5,3 2,77 69
16 4,0 3,68 45
32 3,1 4,74 29
Т аблица 2
N T(N) (с.) H(N) %
1 6,12 1 100
16 0,39 15,7 98
32 0,26 23,5 73
64 0,20 30,6 47
12 0,16 38,3 30
15 0,15 40,8 27
Примечание. N — число процессоров, T(N) — время счета на данном числе процессоров (с.), H(N) — ускорение при использовании данного числа процессоров, Q(N =
= 100%T(1)/{N*T(N)} — эффективность параллельной реализации.
Разработанная параллельная программа 1ш§£_псот_ос3 позволяет проводить численное
моделирование радиационного переноса в смесях веществ с учетом нескольких сотен спектральных линий (до 1000), (см. рис. 1).
3.00E+01 -| 1
2.50E+01 -
2.00E+01 -
1.50E+01 -
1.00E+01 -
5.00E+00 -
0.00E+00 -
0 5 10 15
Рис. 1. Интенсивность излучения при
р = 0,48, х = 0,4, плотность 0,03, температура Te = 1,5, смесь из 4-х веществ, 988 резонансных линий
Учет при численном моделировании значительного количества спектральных линий является весьма актуальным, поскольку моделирование с учетом большого числа спектральных линий наиболее адекватно отражает физическую сущность явления. Решение таких задач сопряжено с очень большими объемами вычислений, которые под силу только мощным современным компьютерам. Более того, считается общепринятым, что для большинства задач и в обозримом будущем решение полного уравнения переноса слишком дорогостояще как по времени, так и по памяти. Современные базы данных содержат сотни тысяч линий. Поэтому в других подходах при численном моделировании, как правило, используются различные физические и математические приближения, в частности многогрупповой метод и т. п.
Но развитие параллельных технологий, использование современных суперкомпьютеров, разработка новых математических методов и вычислительных подходов позволяют надеяться на продвижение в решении таких задач в достаточно полных постановках (метод МПЛЧ).
ПАРАЛЛЕЛЬНАЯ ПРОГРАММА МПЛЧ
Параллельная программа по методике МПЛЧ разработана в системе DVM (Distributed Virtual Memory) или (Distributed Virtual Machine). Распараллеливание программы ведет-
ся по угловой переменной. Все циклы по угловой переменной (а также данные (массивы)) автоматически распределяются на заданное число процессоров. Программа легко масштабируется. Возможно задание любого числа процессоров от 1 до N .
В связи с разработкой программы с фойг-товским профилем излучения и поглощения и использованием новой математической модели (единое уравнение переноса), параллельная программа по методике МПЛЧ существенно изменилась. Ранее в распараллеливании участвовало три подпрограммы. Это подпрограмма решения уравнения кинетики (DCWR), подпрограмма решения уравнения переноса в непрерывном спектре (DKDUC) и подпрограмма решения уравнения переноса в дискретном спектре (DKDUK). При этом подпрограмма решения уравнения переноса в дискретном спектре выполнялась для каждой линии отдельно. Цикл по линиям был распараллелен на произвольное (1 <p < L, L — число линий) число процессоров. При этом был реализован так называемый «параллелизм задач». Внутри подпрограммы DKDUK было добавлено распараллеливание циклов по угловой переменной. В новой программе по методике МПЛЧ отсутствует «параллелизм задач», т. е. подпрограмма DKDUK решает уравнение переноса для всех линий, а также все необходимые интегралы. При этом у нас решается единое уравнение переноса в непрерывном и дискретном спектре. Также подпрограмма DKDUK решает уравнение энергобаланса. Таким образом, распараллеливается только подпрограмма DKDUK по угловой переменной. Подпрограмма DKDUC исключена. Модификация программы по методике МПЛЧ привела к десятикратному ускорению программы.
Для распараллеливания данных в подпрограмме DKDUK система DVM использует модель параллелизма по данным (МПД) [15].
При работе с удаленными данными модель параллелизма по данным преобразуется в модель передачи сообщений.
В модели передачи сообщений каждый процесс имеет собственное локальное адресное пространство. Обработка общих данных и синхронизация осуществляется посредством передачи сообщений. В нашем случае с процессора на процессор (каждый каждому) передаются недостающие части интегралов по угловой переменной. Передача с процессора на процессор интеграла намного экономичнее, чем передача интенсивности излучения.
О РАЗВИТИИ МЕТОДИКИ МАПИ
Методика МАПИ (метод аналитического представления излучения) существует в нескольких модификациях: МАПИ_ЯКОБИ, МА-ПИ_Л, МАПИ_2Л, МАПИ_3Л.
В модификации МАПИ_ЯКОБИ для вычисления интегралов используются разные полиномы (Эрмита-Чебышева, Якоби, Лагерра, Лежандра). Выбор того или иного полинома зависит от асимптотики подынтегральной функции. Естественно, что выбор степени каждого полинома при изменении задачи требует отдельной работы.
В модификации МАПИ_Л все интегралы по энергетической переменной е в линиях вычисляются с использованием метода Лобатто, который позволяет задавать точность вычисления интеграла, варьировать размер области интегрирования при расчете интегралов с бесконечными пределами. Выбор сетки интегрирования в данном методе зависит только от значений подынтегральной функции на сетке и итоговой точности расчета интеграла, заданной для метода. В нашей задаче подынтегральная функция зависит от параметра ц (угловой переменной). Сетка для каждого значения параметра ц строится своя, число точек в сетке для разных значений параметра может значительно отличаться, если поведение подынтегральной функции разное при разных значениях параметра. В нашей задаче поведение интенсивности излучения при малых значениях угловой переменной более сложное, чем при других значениях этого параметра. Соответственно сетка при малых значениях ц содержит больше точек, чем при других значениях этого параметра.
Следует заметить, что подынтегральные функции при решении нашей задачи негладкие. Они имеют разрывы и «пики», поэтому предварительно область интегрирования разбивается на участки, где функция достаточно гладкая, что позволяет вычислять интегралы с высокой точностью. Если заданная высокая точность не может быть достигнута, то это, как правило, является указанием на то, что какая-то особенность функции не выделена, что иногда обнаруживается при счете новых задач.
В модификации МАПИ_2Л метод Лобатто используется для вычисления интегралов по энергии не только в линиях, но и в непрерывном спектре. Разработка этого комплекса программ связана с тем, что при переходе к фойгтовскому профилю излучения и проведению расчетов для широкого диапазона постоянных электронных
температур оказалось, что результаты расчета населенностей уровней по двум методикам не совпадают. Изучение причины расхождений результатов показало, что результат вычисления одного из несобственных интегралов в непрерывном спектре был разным при расчете по методу МПЛЧ с использованием полиномов Лагранжа и по методу МАПИ с использованием полиномов Лагерра. Расчет «спорного» интеграла по методу Лобатто с гарантированной точностью позволил разрешить «спорный» вопрос (пункт 8). Были установлены необходимые для требуемой точности вычислений границы интегрирования, что позволило получить совпадение результатов расчета населенностей уровней по обеим методикам.
В модификации МАПИ_3Л метод Лобатто используется для вычисления интегралов по энергетической переменной во всех случаях, в том числе при вычислении интегралов, необходимых при решении уравнения энергобаланса. Применение этого комплекса программ оправдано при решении новых задач, использующих уравнение энергобаланса, когда поведение функций значительно меняется. Следует заметить, что у нас нет физических тестовых вариантов расчетов. Критерий истины один — совпадение результатов при расчете по разным методикам. Если результаты не совпадают, то большую роль в установлении истины играет возможность расчетов на сетках с любой заданной точностью, не связанных с какими-то определенными полиномами. Такие возможности заложены в комплексе МАПИ_3Л.
Программа МПЛЧ ориентирована на сокращение времени счета варианта, но этот метод при переходе на новые задачи требует дополнительной настройки. Для получения быстрого результата нужна гарантия, что полученный результат правильный. Такую гарантию дает «настройка» с использованием программы МА-ПИ_3Л. Существуют понятия: «методический счет» и «серийный счет». Комплекс МПЛЧ способен проводить серийный счет, а комплексы МАПИ_Л, МАПИ_2Л, МАПИ_3Л — методические расчеты. Они требуют на несколько порядков большего времени для расчета варианта, чем МПЛЧ или достаточно большого числа процессоров, но зато дают картину процесса с заданной гарантированной точностью.
В комплексе МАПИ нет ограничений на выбор сетки по пространственной переменной. Она может задаваться любой, в том числе по узлам полиномов Чебышева, как в методе
МПЛЧ. Это создает дополнительные удобства при сравнении результатов, полученных по разным методикам.
ПАРАЛЛЕЛЬНЫЙ КОМПЛЕКС МАПИ_3Л
В комплексе МАПИ_3Л для распараллеливания используется МР1-технологии параллельного программирования. Как было отмечено ранее [7], комплекс МАПИ разбит на задачи, каждая из которых распараллеливается отдельно по тем параметрам, которые используются в данной задаче. Внутреннее наполнение таких задач зависит от постановки глобальной задачи (ГЗ), но параметры, по которым производится распараллеливание, сохраняются. Иногда меняется диапазон их изменения, если этого требует (ГЗ). Так, задача ЛуШ вычисляет средние по энергии от интенсивности излучения в дискретном спектре. Она распараллелена по точкам (/) в сетке по пространственной переменной (/ <241), по компонентам (/) в смеси веществ (/ <5), по числу линий к в дискретном спектре (к<40), по точкам (/) в сетке по угловой переменной (/<96). В реальных расчетах заложенное в этой задаче полное распараллеливание удается осуществить только тогда, когда диапазон изменения параметров в расчетном варианте достаточно мал. В сложных вариантах на одном из предоставленных р процессоров вычисляется (/х/'хкх/)/р точек задачи.
Другая задача С/, например, вычисляет интегральные характеристики, осуществляющие связь кинетической матрицы с интенсивностью излучения вне линий (в непрерывном спектре), населенности уровней и температуру на сетке по пространственной переменной. Она распараллелена по точкам (/) в сетке по пространственной переменной и по точкам (/) компонент в смеси веществ. Следует заметить, что большая эффективность распараллеливания достигается благодаря тому, что вычисление многих величин повторяется на каждом из используемых процессоров, а передаются с процессора на процессор только сравнительно небольшие массивы переменных. Число используемых процессоров р может быть задано любым. Естественно, что в сложных вариантах, чем больше процессоров используется в счете, тем меньше времени требуется на расчет варианта.
ИЗУЧЕНИЕ ВЛИЯНИЯ ТОЧНОСТИ ВЫЧИСЛЕНИЯ ИНТЕГРАЛА
| е2 ехр(-е / Те )а(е / е л )йе
етт
НА ПОВЕДЕНИЕ НАСЕЛЕННОСТЕЙ УРОВНЕЙ
Приведенный в заголовке раздела 7 интеграл используется при вычислении компонент матрицы в уравнении кинетики. Для расчета данного интеграла в программе МАПИ_ЯКОБИ использовались полиномы Лагерра вычисления несобственных интегралов, а в программе МПЛЧ — полиномы Лагранжа. В методе МПЛЧ пределы интегрирования задаются конечными (верхний предел интегрирования выбирался Я=30). Когда рассчитывались варианты с невысокой постоянной электронной температурой, итоговые результаты расчетов населенностей и интенсивностей излучения по обеим методикам совпадали (что принималось за критерий истины). При расчете варианта, в котором была задана электронная температура Те = 10, обнаружились значительные расхождения в результатах расчетов населенностей по двум разным методикам (рис. 2, кривая 1 получена по методу МПЛЧ, кривая 2 получена по программе МА-ПИ_ЯКОБИ). Анализ расчетов показал, что результаты расчета интенсивности излучения одинаковы. В кинетической матрице компоненты, зависящие от интенсивности излучения, совпадали в результате расчетов по обеим методикам. Отличия обнаружились при вычислении выписанного в заголовке раздела несобственного интеграла. Причина могла быть в неточности счета как по методике МАПИ_ЯКОБИ, так и по методике МПЛЧ. Возможно, неправильно выбрана степень используемого полинома. Потребовался такой метод вычисления, который не связан со степенью полинома. Так возник комплекс МАПИ_2Л, который позволил провести расчеты несобственного интеграла с разной заданной точностью и с разными пределами интегрирования и установить, когда увеличение верхней границы интегрирования Я не влияет на величину интеграла, что привело к уточнению счета данного интеграла по методу МПЛЧ. В результате расчеты по двум методикам совпали (рис. 2, Я = 1000, линия 2). На рис. 2 изображена населенность уровня с для двух значений верхней границы вышеуказанного интеграла в методе МПЛЧ((1) — Я = 30 и (2) — Я = 1000), вещество — цинк. При Я = 1000 графики насе-
ленности, сосчитанные по разным методикам (МАПИ_2Л, МАПИ_ЯКОБИ, МПЛЧ) совпали.
Рис. 2. Населенность уровня с1
РЕЗУЛЬТАТЫ ЧИСЛЕННОГО МОДЕЛИРОВАНИЯ
Разработана методика на основе подхода МАПИ (комплекс программ МАПИ_3Л) позволила решать задачу с заданной (гарантированной) точностью, что дает возможность определить необходимые параметры и сетки, оценить точность вычислений по «быстрой» методике, основанной на методе МПЛЧ. Сравнительные расчеты по методам МПЛЧ и МАПИ приведены на рис. 2, 3. На рис. 3 представлен расчет варианта для слоя, состоящего из двух подслоев, содержащих вещество кальций (п = 25,9; к0 = 0,5). В одном подслое (0 < х < 0,4) плотность вещества равна 0,03; в другом подслое (0,4 < х <0,8) плотность вещества равна 0,1. Учитывались две резонансные линии. На рис. 3, а — населенность сь сосчитанная по методике МАПИ (сплошная линия) и по методике МПЛЧ (маркеры). На рис. 3, б — электронная температура, полученная по двум методикам для данного варианта.
Исследована зависимость электронной температуры и населенностей уровня от плотности вещества в слое (рис. 4, а, б). Наблюдалось незначительное уменьшение электронной температуры в слое при увеличении плотности вещества.
Оценка влияния учета уравнения энергобаланса (рис. 5). Изучение влияния показало, что качественно картина не меняется в случае переменной электронной температуры.
Рис. 3. Сравнение результатов, полученных по двум методикам
б
а
б
а
Рис. 4. Сравнение результатов, полученных для веществ с различной плотностью: а - электронная температура; б - населенности уровней с1 в случае Те= Те(х), п = 25,9, к0 = 0,5, 1 — плотность 0,03, 2 — плотность 0,25, 3 — плотность 0,35
Рис. 5. Интенсивность излучения при:
^ = 0,48, х = 0,4, плотность — 0,03, смесь из 4 компонент вещества, 156 линий: 1 — интенсивность излучения для температуры Те = 1,5; 2 — интенсивность излучения;
Те = Те(х), п = 25,9, к0 = 0,5
ЗАКЛЮЧЕНИЕ
В статье представлены результаты развития методики численного моделирования взаимодействия ионизирующего излучения с веществом. Рассматриваемая модель предполагает, что изучается случай вещества, имеющего фойгтов-ские профили излучения и поглощения, энергобаланс учитывается. Методика предполагает расчет вклада каждой учитываемой спектральной линии; вычисление населенностей уровней в соответствии с уравнениями кинетики населенностей; нахождение интенсивности излучения в соответствии со спектральным уравнением переноса с коэффициентами, вычисляемыми по населенностям уровней и электронной температуре, которая удовлетворяет уравнению энергобаланса. Методика основана на двух методах решения задачи и, соответственно, на разработке двух реализующих эти методы параллельных алгоритмов.
Основными результатами являются:
• разработаны параллельные алгоритмы и программы, реализующие методы МАПИ и МПЛЧ для новой математической постановки задачи, отражающей новую (и более полную) физическую модель;
• разработан алгоритм решения уравнения энергобаланса (пункт 2);
• создан комплекс программ для решения задачи переноса излучения в плоском слое смеси веществ с учетом нескольких сотен линий (рис. 1);
• разработаны новые алгоритмы и программы МПЛЧ, позволившие на порядок сократить время расчетов по сравнению с ранее применявшимся параллельным алгоритмом МПЛЧ;
• достигнуто приемлемое время расчетов достаточно сложных и физически содержательных задач;
• получена возможность проведения серийных инженерных и научных расчетов;
• проведены серии расчетов на суперкомпьютерах URAN (ИММ УрО РАН, г. Екатеринбург) и МВС-100К (МСЦ РАН, г. Москва), в том числе для определения электронной температуры в однородных и неоднородных средах;
• проведено сравнение результатов, полученных по методикам МАПИ и МПЛЧ, установлены допустимые пределы изменения параметров задачи.
СПИСОК ЛИТЕРАТУРЫ
1. Михалас Д. Звездные атмосферы. Т. 1, 2. М.: Мир, 1982.
2. Никифоров А. Ф., Новиков В. Г., Уваров В. Б.
Квантово-статистические модели высокотемпературной плазмы и методы расчета росселандовых пробегов и уравнений состояния. М.: Физ.-мат. лит-ра, 2000.
3. Сушкевич Т. А Математические модели переноса излучения. М. БИНОМ. Лаборатория знаний, 2006.
4. Gonzalez M., Audit E. Numerical treatment of radiative transfer // Astrophysics and Space Science 298: 357-362, 2005.
5. Manual of the FLYCHK code (National Institute of Standards and Technology, USA) // URL: http://nlte.nist.gov/FLY/Doc/Manual_FLYCHK_Nov08.pdf.
6. Параллельные вычислительные технологии в задаче о переносе радиационного излучения / Е. Ф. Леликова [и др.]. // Вопросы атомной науки и техники. Сер. Математическое моделирование физических процессов. М. 2002. Вып. 3. С. 3-13.
7. Параллельные вычисления в задачах, возникающих при математическом моделировании переноса излучени / Е. Ф. Леликова [и др.]. // Автоматика и телемеханика. 2007. № 5. C. 126-140.
8. Параллельные вычисления в задачах, возникающих при математическом моделировании переноса излучения / Е. Ф. Леликова [и др.]. // В кн.: Энциклопедия низкотемпературной плазмы, Серия "Б", т. VII-1 «Математическое моделирование в низкотемпературной плазме». Часть 3. [Под ред. чл.-корр. РАН Ю.П. Попова. Гл. ред. энцикл. cер. акад. В.Е. Фортов]. М.: ЯНУС-К, 2008. 698 с. С. 514-522.
9. Чащин М. А., Рубина Л. И., Ульянов О. Н. Развитие алгоритмов для задач моделирования радиационного переноса // Тезисы докладов Всероссийской школы молодых исследователей и V Всероссийской конференции «Актуальные проблемы прикладной математики и механики»,
посвященной памяти академика А.Ф. Сидорова (Абрау-Дюрсо, 13-18 сентября 2010 г.). Екатеринбург: УрО РАН, 2010. С. 86.
10. Численное моделирование переноса излучения в смесях веществ / О. Н. Ульянов [и др.] // Тезисы международной конференции «X Забабахинские научные чтения», РФЯЦ-ВНИИТФ, 2010 г. Снежинск. С.272-273.
11. Параллельные вычислительные технологии в задаче о переносе излучения / М. А. Чащин [и др.] // Параллельные вычислительные технологии (ПаВТ’2010): Труды международной научной конференции (Уфа, 29 марта - 2 апреля 2010 г.) [Электронный ресурс] - Челябинск: Издательский центр ЮУрГУ, 2010. С. 692. ISBN 978-5-69603987-9 [Электронное издание].
12. Bowen C. Lee R. W. Ralchenko Yu. Comparing plasma population kinetics codes: Review of the NLTE-3 Kinetics Workshop // Journal of Quantitative Spectroscopy and Radiative Transfer, vol. 99, no1-3, pp. 102-119 (2006).
13. The HITRAN Database. URL: http:// www.hitran.com .
14. Rothman L. S. The HITRAN 2008 molecular spectroscopic database // Journal of Quantitative Spectroscopy and Radiative Transfer, vol. 110, pp. 533-572 (2009).
15. URL: http://www.keldysh.ru/dvm.
ОБ АВТОРАХ
Рубина Людмила Ильинична, ст. науч. сотр. отдела прикладных задач Ин-та математики и механики им. Н.Н. Кра-совского УрО РАН, г. Екатеринбург. Дипл. математик по вычислительн. математике (Уральск. гос. ун-т,
г. Екатеринбург, 1961). Канд. физ.-мат. наук (1973). Иссл. в обл. методов получения решений нелинейных уравнений и систем в частных производных.
Ульянов Олег Николаевич, ст. науч. сотр. того же отдела, доц. каф. вычислительн. методов и уравнений матем. физики Уральск. федерального ун-та, г. Екатеринбург. Дипл. математик (Уральск. гос. ун-т, г. Екатеринбург, 1979). Канд. физ.-мат. наук (ИММ УрО РАН, 1992). Иссл. в обл. численных методов и параллельных алгоритмов решения прикладных задач.
Чащин Михаил Александрович, гл. программист того же отдела. Дипл. математик (Уральск. гос. ун-т, г. Екатеринбург, 1987). Иссл. в обл. численных методов решения прикладных задач и параллельных вычислений на многопроцессорных системах.
METADATA
Title: Development of two parallel algorithms for numerical simulation of interaction of radiation with substance.
Authors: L. I. Rubina, O. N. Ulyanov, M. A. Chashchin Affiliation:
N.N. Krasovskii Institute of Mathematics and Mechanics of the Ural Branch of the Russian Academy of Sciences (IMM UB RAS), Ekaterinburg, Russia.
Email: [email protected]; [email protected]; cma@imm. uran.ru .
Language: Russian.
Source: Vestnik UGATU (Scientific journal of Ufa State Aviation Technical University), 2013, Vol. 17, No. 2 (55), pp. 6474. ISSN 2225-2789 (Online), ISSN 1992-6502 (Print). Abstract: We present the technique of modeling the interaction of ionizing radiation with a substance, which aimed to solve radiation transfer problems in case of mixture of sub-
stances for a number of physical and mathematical problem statements. The technique involves the calculation of the contribution of each accounted spectral line, computation of populations levels in accordance with the equations of the kinetics; the finding of the radiation intensity in accordance with the spectral transfer equation with coefficients calculated on populations levels and electronic temperature, which satisfies the equation of the energy balance. Solving the problem in a comprehensive statements is connected with large amounts of computing and requires the application of parallel technologies. The authors design the technique which is based on implementing two different methods of the numerical solution and developing two parallel algorithms. The proposed algorithms and designed codes are allowed us to carry out simulation of radiation transfer for case of plane layer of substance taking into consideration a large number of spectral lines (up to 1000) in acceptable time with reasonable accuracy. The results of numerical experiments are presented.
Key words: radiative transfer, numerical simulation, parallel algorithms.
References (English Transliteration):
1. Mihalas D., Stellar Atmospheres, San Francisco: Freeman,
1978. Translated under the title Zvezdnye atmosfery, Moscow: Mir, 1982, vols. 1, 2.
2. Nikiforov A.F., Novikov V.G., and Uvarov V.B. Kvantovo-
statisticheskie modeli vysokotemperaturnoi plazmy i meto-dy rascheta rosselandovykh probegov i uravnenii sostoya-niya (Quantum-Statistical Models of High-Temperature Plasma and Computing Methods for Rosseland Length and State Equations), Moscow: Fizmatlit, 2000. (In Russian).
3. Sushkevich T.A. Matematicheskie modeli perenosa izluche-
niya (Mathematical models of Radiation Transfer), Moscow: BINOM, 2006. (In Russian).
4. Gonzalez M., Audit E. Numerical treatment of radiative trans-
fer // Astrophysics and Space Science, 298: 357-362, 2005.
5. Manual of the FLYCHK code (National Institute of Stan-
dards and Technology, USA) // URL: http ://nlte.nist. gov/FLY/Doc/Manual_FLY CHK_Nov08. pdf.
6. Lelikova E.F., Rubina L.I., Ul’yanov O.N., and Chashchin
M.A. «Parallel Computer Technologies in the Problem of Radiation Transfer» in Vopr. Atom. Nauki Tekh., Ser. : Mat. Model. Fiz. Prots., 2002, no. 3, pp. 3-13. (In Russian).
7. Lelikova E.F., Rubina L.I., Ul’yanov O.N., and Chashchin
M.A. «Parallel computing in problems that occur during mathematical simulation of radiation transfer» in Automation and Remote Control. May 2007, Volume 68, Issue 5, pp 860-873. (Original Russian Text © E.F. Lelikova, L.I. Rubina, O.N. Ul’yanov, M.A. Chashchin, 2007, published in Avtomatika i Telemekhanika, 2007, No. 5, pp. 126140.) (In Russian).
8. Lelikova E.F., Rubina L.I., Ul’yanov O.N., and Chashchin
M.A. «Parallel computing in problems that occur during mathematical simulation of radiation transfer» in Encyclopedia of Low-Temperature Plasma [Editor-in-Chief V. E. Fortov]. Series B, Volume VII-1, Part 3 Mathematical modeling of low-temperature plasmas [Editor-in-Chief Yu. P. Popov], pp. 514-522. (In Russian).
9. Chashchin M.A., Rubina L.I., Ulyanov O.N. «Development
of algorithms for radiation transfer simulation problems» in Talk abstracts of young scientists School-Conference and V All-Russian conference «Modern Problems of applied
mathematics and mechanics» dedicated to the memory of academician of A.F. Sidorov (September 13-18, 2010, Ab-rau-Dyurso), p. 86. (In Russian).
10. Ulyanov O.N., Lelikova E.F., Rubina L.I., and Chashchin M.A. «Numerical modeling of radiation transfer in mixtures of substances» in Abstracts of Zababakhin Scientific Talks (International conference. March 15-19, 2010), pp.272-273.
11. Chashchin M.A., Lelikova E.F., Rubina L.I., Ulyanov O.N. «Parallel computational technologies for radiation transfer problem» in Proceedings of International conference «Parallel computational technologies (PCT) 2010» (March 29-Aprile 02, 2010, Ufa). P. 692. [Electronic publication]. ISBN 978-5-696-03987-9. (In Russian). 12. Bowen C. Lee R.W. Ralchenko Yu. «Comparing plasma population kinetics codes: Review of the NLTE-3 Kinetics Workshop» in Journal of Quantitative Spectroscopy and Radiative Transfer, 2006, vol. 99, no1-3, pp. 102-119.
13. The HITRAN Database URL: http: // www.hitran.com .
14. Rothman L.S. «The HITRAN 2008 molecular spectroscopic database» in Journal of Quantitative Spectroscopy and Radiative Transfer, vol. 110, pp. 533-572 (2009).
15. URL: http://www.keldysh.ru/dvm.
About authors:
1. Rubina, Ludmila Il’inichna, Senior Research Assistant, Dept. of Applied Problems (IMM UB RAS), Dipl. mathematician (Ural State Univ., 1961), Cand. of Sci. (Ph.D.) in physics and mathematics (Ural State Univ., 1973), Senior Researcher (1985).
2. Ulyanov, Oleg Nikolaevich, Senior Research Assistant, Dept. of Applied Problems (IMM UB RAS), Dipl. mathematician (Ural State Univ., 1979), Cand. of Sci. (Ph.D.) in physics and mathematics (IMM UB RAS, 1992), Assoc. Prof., Dept. of Numerical Methods and Equations of Mathematical Physics (UrFU, 2001), Senior Researcher (1999).
3. Chashchin, Mikhail Aleksandrovich, chief programmer, Dept. of Applied Problems (IMM UB RAS), Dipl. mathematician (Ural State Univ., 1987).
