CC BY
14
УДК 51.73
DOI: 10.14529/cmse210106
ПРОГРАММНЫЙ комплекс радуга-т ДЛЯ МОДЕЛИРОВАНИЯ ПОЛЕЙ НЕЙТРОНОВ В ЯДЕРНО-ЭНЕРГЕТИЧЕСКИХ УСТАНОВКАХ
© 2021 О.В. Николаева, С.А. Гайфулин, Л.П. Басс
Институт прикладной математики им. М.В. Келдыша РАН (123047 Москва, Миусская пл., д. 4) E-mail: nika@kiam.ru, sagarp@yandex.ru, bass@kiam.ru Поступила в редакцию: 14.07.2020
При проектировании и сопровождении эксплуатации ядерно-энергетических установок (ЯЭУ) необходимо выполнять моделирование в этих установках потоков нейтронов. При задании геометрии ЯЭУ необходимо учитывать границы разномасштабных конструктивных элементов, состоящих из материалов с существенно различными свойствами. Из-за больших размеров ЯЭУ для расчетов желательно использовать параллельные компьютеры. Для выполнения такого моделирования развиваются алгоритмы и программы численного решения краевой задачи для интегро-дифференциального уравнения переноса нейтронов на неструктурированных сетках. В статье приводится описание реализованных в программном комплексе РАДУГА-Т алгоритмов решения такой задачи. Представлены сетки, сеточные схемы, итерационные методы решения систем сеточных уравнений. Рассмотрены методы распараллеливания вычислений на гибридных компьютерах (используются технологии MPI и OpenMP). Представлены методы работы с пространственными сетками (построение, улучшение качества, декомпозиция, визуализация). Описаны особенности программной реализации. Проведено сравнение используемых в программном комплексе РАДУГА-Т алгоритмов с алгоритмами в других аналогичных программных комплексах. Приведены результаты исследования эффективности распараллеливания вычислений в задаче расчета коэффициента размножения нейтронов в модели легководного реактора. Исследования выполнены на многопроцессорном компьютере МВС-10П (МСЦ РАН). Приведены значения ускорения вычислений каждого из используемых в расчете алгоритмов и суммарного ускорения всего расчета.
Ключевые слова: уравнение переноса, неструктурированные сетки, сеточные схемы, итерационные методы, параллельные вычисления.
ОБРАЗЕЦ ЦИТИРОВАНИЯ
Николаева О.В., Гайфулин С.А., Басс Л.П. Программный комплекс РАДУГА-Т для моделирования полей нейтронов в ядерно-энергетических установках // Вестник ЮУрГУ. Серия: Вычислительная математика и информатика. 2021. Т. 10, № 1. С. 75-89. DOI: 10.14529/cmse210106.
Введение
При проектировании и сопровождении эксплуатации ядерно-энергетических установок (ЯЭУ) необходимо выполнять моделирование в этих установках потоков нейтронов. Нейтроны в работающих ЯЭУ возникают в результате ядерных реакций; в дальнейшем они вызывают такие же реакции и тем самым порождают энерговыделение в активной зоне (АЗ) ЯЭУ. Кроме того, нейтроны повреждают конструктивные элементы как в АЗ, так и в окружающих ее слоях радиационной защиты. Нейтроны, покинувшие ЯЭУ, влияют на радиационный фон в ее окрестности.
Для определения потоков нейтронов в АЗ часто используется уравнение диффузии. Такое приближение применимо в предположении, что потоки слабо зависят от направления скорости нейтрона. Однако такое предположение часто не выполнено даже в АЗ из-за наличия щелей и зазоров между конструктивными элементами [10]. Тем более оно не выполнено в слоях защиты и в окрестности ЯЭУ. В этом случае для определения потоков нейтронов используется интегро-дифференциальное уравнение переноса.
Для задания распределения материалов в 3Б области расчета с декартовой геометрией используются либо структурированные сетки с ячейками-параллелепипедами, либо неструктурированные сетки с ячейками-тетраэдрами, гексаэдрами, призмами. Структурированные сетки легче в построении. Кроме того, сеточные схемы на ячейках-параллелепипедах обладают меньшей погрешностью аппроксимации, чем на ячейках-тетраэдрах и призмах [9]. С другой стороны, грани ячеек-параллелепипедов не могут точно передать границы между материалами. Только с помощью неструктурированных сеток могут быть корректно описаны локальные неоднородности в области расчета [7]. Однако при переходе от структурированных сеток к неструктурированным возникают дополнительные трудности, связанные с построением сеточных схем, итерационных алгоритмов, распараллеливанием вычислений и обработкой результатов расчетов.
В настоящей работе описаны возможности и алгоритмы программого комплекса РАДУГА-Т, предназначенного для решения интегро-дифференциального уравнения переноса в 3Б областях на неструктурированных сетках на гибридных многопроцессорных системах. Статья имеет следующую структуру. В разделе 1 представлена постановка краевой задачи для уравнения переноса нейтронов. В разделе 2 обсуждаются методы ее решения; описаны методы, используемые в комплексе РАДУГА-Т и методы, применяемые в других комплексах. В разделе 3 перечислены возможности и алгоритмы программного комплекса РАДУГА-Т. В разделе 4 приведены результаты исследования эффективности распараллеливания вычислений в программном комлексе РАДУГА-Т. В заключении приводится краткая сводка результатов, полученных в работе, и указаны направления дальнейших исследований.
1. Краевая задача для уравнения переноса нейтронов
Уравнение переноса нейтронов имеет следующий вид
64/60. + а(г,Е)4(г,П,Е) = ^4(г,а',Е') &',Е,Е')йП' йЕ'+ (1)
+ 7^ 1(ЮЦПгМ',Е') ог(т,Е'Ш' <1Е' + Р(г,а,Е\ Е=п <Е< Етах, г ЕС, ПеП.
Здесь решение 1¥(г,0.,Е) — плотность потока нейтронов кинетической энергии Е из интервала (£=;?, Етах) в пространственной точке г выпуклой области С, движущихся в направлении единичного вектора О, О — единичная сфера.
Величина о(г,Е) — полное сечение взаимодействия со средой в точке г нейтрона энергии Е; левая часть уравнения (1) отвечает уменьшению числа нейтронов в единичном объеме фазового пространства (г, О., Е) за счет пролета нейтронов без столкновения с атомами среды и за счет их поглощения и рассеяния атомами.
Первый интеграл в правой части уравнения (1) отвечает за увеличение числа нейтронов в единичном объеме фазового пространства при изменении энергии Е' и направления движения О' нейтрона за счет рассеяния на атомах среды. Коэффициент ст5(г, П • 0.',Е,Е') есть сечение рассеяния. Отметим, что сечение рассеяния зависит от скалярного произведения П • О' векторов О и О' , то есть от косинуса угла рассеяния. Второй интеграл в правой части уравнения (1) отвечает увеличению числа нейтронов за счет деления атомов. Коэффициент 0;(г,Е') есть сечение деления, функция х(Е) задает спектр образующихся при делении нейтронов, ке;; — коэффициент размножения. Функция Е(г,0.,Е) есть плотность источника нейтронов.
Краевые условия для уравения (1) задаются на границе 6С области С
4(г*,П,Е) = /(г*,Ф(г*,а,Е)) при Г*едв, п(г*) • П < 0. (2)
Здесь п(г*) — внешняя нормаль к границе Зб в точке г*, / — кусочно-линейная функция Таким образом, краевые условия в каждой точке г* границы Зб задаются для тех направлений скорости О, которые имеют тупой угол с внешней нормалью, то есть отвечают входящим в область С нейтронам.
Отметим, что существует два варианта задачи (1), (2).
- В предположении А:е// = 1 имеем задачу с источником F(r, П, £").
- В предположении F(r, П, £") = 0 имеем задачу определения коэффициента А:е//.
Для решения задачи (1), (2) вводится сетка (как правило, неравномерная) по энергии Е, состоящая из Ц интервалов (групп)
^тах = > ¿2 > > > £9+1 = -^тт. (3)
Задача (1), (2) принимает следующий вид
^К] = Е^АЛ4?] + ЭД^Ае// + ^(г,«), (4)
= зчузп + ^(г)%(г,п), ^[Чр] = /^р(г,П') ^,р,„(г,П • П')^п', (5)
■?/№>]= ¿г ^Ч2и=1^/,р(г)/Чр(г,П') ЙП', (6)
%(г*,П) = /(г*,Чч(г*,П)) при г*еЗС, п(г*) • П < 0, д = 1.....(7)
Здесь Чч(г, П) — плотность группового потока нейтронов с энергией из интервала < Е < функция ^(г, П) — плотность источника нейтронов с энергией из д-го интер-
вала. Величины стч(г), П • П'), о/,р(г) суть групповые сечения взаимодействия, рас-
сеяния и деления, ^ — групповой спектр образующихся при делении нейтронов. Решение задачи (4)—(7) может быть найдено методом Зейделя
^к+ч=+^?/[ч3р]/^К//+^(Г,п), = ^ ^?/[ч3р+1]/^?/[ч3р], (8)
Ч^+4(г*,П) = /(г*,Ч^+Чг*,П)) при г*еЗС, п(г*) • П < 0, д = 1.....(9)
где г — индекс итерации. Эти итерации первого уровня. Они также называются внешними итерациями. Для решения задачи (8), (9) тоже используется метод Зейделя
ад+Ч = Ер^^^+Ч+Е^+ЛА7] + ■^/[Чр]/^е// + ^(г,п), (10)
Ч^+ЧЛП) = /(г*,Ч'7+1(г*,П)) при г*еЗС, п(г*) • П < 0, д = 1.....(11)
где к — индекс итерации. Это итерации второго уровня.
Поскольку только медленные нейтроны при взаимодействии со средой могут увеличить свою кинетическую энергию, а группы нумерованы по убыванию энергии, см. (3), то при р > ц большинство сечений рассеяния = 0. Если при р > ц все сечения = 0,
то в методе (10), (11) выпоняется одна итерация. В любом случае в методе (10), (11) на каждой итерции последовательно для всех энергетических групп решаются задачи
=$,,,[%]+ ^0-,"), (12) %(г*,П) =/(г*,%(г*,П)) при г*еЗС, п(г*) • П < 0, д = 1.....(13)
Здесь правая часть ^(г, П) содержит как плотность источника ^(г, П), так и интегралы от плотностей потока в других группах. Базовой при моделировании потоков нейтронов является краевая задача (12), (13).
2. Методы решения базовой задачи
2.1. Сетки и сеточные алгоритмы
При решении сеточным методом задачи (12), (13) на единичной сфере W вводится набор узлов (Пт), m=1,...,M, каждая функция Ч^Сг,П) заменяется функциями 4qm(r) = Wq(r,Q.m), интеграл в правой части (12) заменяется квадратурной формулой [1].
Далее в пространственной области G вводится неструктурированная сетка из ячеек-выпуклых многогранников. Используются ячейки-тетраэдры [2, 4, 12, 14, 20], произвольные гексаэдры, треугольные и четырехугольные призмы [4]. Тетраэдры универсальны, они позволяют аппроксимировать любые границы между материалами. Треугольные призмы хорошо приспособлены для задания тепловыделяющих элементов (ТВЭЛ) в АЗ и удобны для моделирования движения стержней в ТВЭЛах. Использование произвольных гексаэдров позволяет уменьшить число ячеек в однородных подобластях.
В каждой ячейке функция 4qm(r) представляется своими значениями в каких-либо выбранных точках, например, в вершинах. Аналогично эта функция представляется на каждой грани каждой ячейки. Причем лучше использовать так называемые «разрывные» сеточные аппроксимции, когда значения решения в вершинах грани, общей для двух ячеек, предполагаются различными в каждой из этих ячеек.
Сеточная схема для уравнения переноса (12) строится отдельно в каждой ячейке. При этом используется та особенность уравнения переноса (12), что для каждого вектора направления скорости Q.m грани каждой ячейки разделяются на «входные», для которых Лт • п < 0, и «выходные», для которых Q.m • п > 0, где n — внешняя нормаль к грани. Для уравнения (12) в ячейке краевые условия вида (13) задаются на «входных» гранях.
При построении сеточных схем для уравнения переноса применяется либо метод конечных элементов, либо метод характеристик, либо slice-balance метод. В методе конечных элементов (МКЭ) [2, 5] в ячейках и на гранях искомое решение и правая часть представляются разложением по базисным функциям. При этом коэффициентами разложения служат значения решения и правой части в вершинах ячеек и граней. Сеточные уравнения находятся после подстановки разложений в уравнение (12) и интегрирования получившегося соотношения с различными весовыми функциями. В методе характеристик (МХ) [15, 20] уравнение (12) решается отдельно в каждой ячейке в предположении, что на «входных» гранях значения решения известны. Далее находятся значения построенного решения в вершинах ячейки и «выходных» граней. В slice-balance методе [19] ячейка разбивается на подъячейки плоскостями, параллельными вектору Q.m. Тогда в каждой подъячейке имеется только одна «входная» и одна «выходная» грани. Сеточная схема строится отдельно в каждой подъячейке методом конечных элементов и выполняется интерполяция на всю ячейку решений в подъячейках. Slice-balance метод является гибридом методов конечных элементов и характеристик.
В программном комплексе РАДУГА-Т расчет может быть выполнен на сетке из ячеек-выпуклых многограников, среди которых могут быть тетраэдры и треугольные призмы. В комплексе реализованы схемы метода характеристик для ячеек-тетраэдров, метода конечных элементов для ячеек-тетраэдров и треугольных призм [7, 9]. Для ячеек иных типов в настоящее время используется шаговая схема первого порядка точности [1].
2.2. Итерационные методы решения базовой задачи
Грани ячеек, аппроксимирующие границу dG, будем называть граничными. Остальные грани будем называть внутренними. Для каждого вектора направления скорости Q.m все граничные грани делятся на «входные», где Q.m • n(r*) <0, и «выходные», где Q.m •
п(г*) > 0. Определим следующие вектора.
- Вектор Y содержит значения решения для всех векторов Q.m в ячейках, на внутренних гранях и на «выходных» граничных гранях.
- Вектор содержит значения решения на граничных «входных» гранях.
- Вектор +Ф содержит значения решения на граничных «выходных» гранях.
- Вектор f содержит значения граничного источника на «входных» гранях.
- Вектор F содержит значения правой части в ячейках.
Тогда сеточные уравнения, аппроксимирующие задачу (12), (13), можно записать в форме
LV = SV + F, ~V = f(+V). (14)
Матрицы L и S аппроксимируют операторы L и S в уравнении (12) соответственно.
Для решения системы уравнений (14) можно использовать метод простой итерации
L^+^S^ + F, Wn+1 = f(+Vn). (15)
Здесь n — номер итерации. Это итерации третьего уровня или внутренние итерции.
Одна итерация в методе (15) требует обращения оператора L. Для каждого узла Q.m элементы вектора ^Vn+1 находятся последовательным расчетом во всех пространственных ячейках. Первой рассчитывается прилегающая к границе dG ячейка, для которой все граничные грани являются «входными». Найденные при этом значения решения на «выходных» гранях являются значениями решения на «входных» гранях для соседних ячеек. Подчеркнем, что порядок расчета ячеек зависит от вектора Q.m.
Поскольку простые итерации могут сходиться очень медленно, вместо (15) может использоваться двухшаговый итерационный KP1 метод [3]
L \^п+4/г = § + р, -Wn+1/S =f(+wn), (16)
D L tyn+1/2 = DS tyn+1/2 +DS (Vn - 4n+1/2), ^Vn+1 = Wn+1/2 + P tyn+1/2. (17)
Здесь на втором Р1-шаге (17) к решению, полученному методом простой итерации на первом K-шаге (16), добавляется ускоряющая поправка tyn+1/2. Поправка tyn+1/2 ищется как линейная функция компонент вектора направления скорости Q.m [3]. Уравнение для поправки находится суммированием по индексу m уравнения для ошибки Ф — Фn+1/2 (это суммирование реализует оператор D). В пространство искомого решения Ф поправка фn+1/2 переводится оператором Р.
Также используется DSA (Diffusive Synthetic Acceleration) метод, где поправка tyn+1/2 ищется как функция, не зависящая от вектора Q.m.
Достоинство KP1 и DSA методов заключается в том, что ускоренный процесс (16), (17) сходится к тому же решению, что и неускоренный процесс (15). Более популярен метод грубосеточного ребаланса [2, 4, 15], не обладающий таким качеством. Однако в нем усреднение выполняется еще и по группам пространственных ячеек, что еще больше снижает размерность поправки. В программном комплексе РАДУГА-Т из двухшаговых итерационных методов используется только KP1 метод.
2.3. Распараллеливание вычислений на распределенной памяти
Разобьем пространственную сетку на J непересекающихся подобластей, каждую подобласть присоединим к своему процессору. Метод простой итерации (15) принимает вид
L} Ч^1 = У? + , ~lVjn+1 = K+VP, Г = 1.....J. ?*}), 3=1,J (18)
или
L} Ч^1 = S} Ч*," + F}, "Ч^1 = f( +Ч'"+1, У = 1...../, / Ф ;) j=1,..., J. (19)
Здесь j и _/' — индексы подобластей. Уравнения (18) отвечают BJ (Block-Jacobi) методу, а уравнения (19) — BS (Block-Seidel) методу.
В BJ методе итерации выполняются как по правой части, так и по значениям решения на границе подобласти. В BS методе итерации выполняются только по правой части. Это достигается за счет того, что каждая итерация совершается за несколько проходов, между проходами выполняется обмен данными между соседними подобластями. Фактически в BS методе выполняется сквозной расчет по всей области, как при расчете на одной подобласти. При этом значения плотности потока нейтронов, выходящих из одной подобласти и входящих в другую, пересылаются между соответствующими процессорами.
В BJ методе число итераций растет с ростом числа подобластей. В BS методе число итераций не зависит от числа подобластей, но каждая итерация требует большего времени, чем итерация BJ метода. Задержки возникают вследствие того, что на разных проходах каждый процессор рассчитывает разное число ячеек. Кроме того, число актов обмена очень велико. Таким образом, часть времени процессор проводит в ожидании получения необходимых ему для расчета данных от других процессоров.
После завершения простой итерации (18) или (19) может быть вычислена ускоряющая поправка в соответствии с уравнением (17). Подчеркнем, что уравнение (17) решается во всей области. При применении на первом шаге BJ метода уравнение для поправки меняется при увеличении числа подобластей и поправка ф"+1/2 ускоряет итерации как по правой части, так и по значениям решения на границах подобластей. При использовани на первом шаге BS метода уравнение для поправки не зависит от числа подобластей и поправка ускоряет только итерации по правой части.
Основная трудность при применении BS метода состоит в нахождении такой последовательности расчета ячеек, при котором было бы минимально время ожидания процессором данных от соседних процессоров. В то же время использование ускоряющей поправки позволяет снизить число простых итераций. Поэтому в некоторых комплексах применяется BJ метод [12, 15, 20]. Тем не менее развиваются и методы построения оптимальных последовательностей расчета ячеек для BS метода. Особенно популярен здесь KBA метод, названный по именам его авторов (Koch, Baker, Alcouffe) [13, 16, 17].
В программном комплексе РАДУГА-Т реализованы и BJ метод, и BS метод, причем каждый из них может быть выполнен как при с использованием, так и без использования ускоряющей процедуры (17). Таким образом, возможно применение одного из четырех итерационных процессов: BJ, BJ+P1, BS, BS+P1.
Решение уравнения для поправки в программном комплексе РАДУГА-Т также выполняется параллельно, тогда как в других комплексах ее решение поручается одному процессору [15]. В комплексе РАДУГА-Т для решения системы (17) используется параллельный метод минимальных невязок. Данный метод опирается на построение системы базисных векторов; причем основные совершаемые операции — это вычисление скалярных произведений и умножение матрицы на вектор. В параллельном счете каждый процессор выполняет операции со своими элементами векторов и матриц; далее после межпроцессорных обменов формируются единые для всей области базисные вектора.
2.4. Дополнительное распараллеливание на общей памяти
Поскольку процессоры обычно являются многоядерными, наряду с распараллеливанием на верхнем уровне на распределенной памяти возможно дополнительное распараллеливание на нижнем уровне на общей памяти.
При обращении оператора L, распараллеливание на нижнем уровне возможно и по пространственным ячейкам, и по узлам Пт, и по энергетическим группам. Распараллеливание
по энергетическим группам приводит к увеличению числа итераций в методе Зейделя (8), (9). Оно применяется в [15], где при построении ускоряющей поправки выполняется и усреднение по энергии, и, тем самым, ускорение итераций (8), (9). Распараллеливание по пространственым ячейкам на нижнем уровне привлекается, если используется BS метод сквозного счета [13, 16, 17]. Оптимально распараллеливание по узлам 0.т, так как расчет для всех узлов выполняется независимо. В программном комплексе РАДУГА-Т, также как в [20], распараллеливание на нижнем уровне выполняется только по узлам 0.т.
Распараллеливание нижнего уровня в методе минимальных невязок при решении системы линейных алгебраических уравнений для поправки (17) в программном комплексе РАДУГА-Т выполняется по элементам матриц и векторов. На обоих шагах (К и P1) необходимо выполнять большое число актов обмена. Для выполнения обменов асинхронно выделяется одно процессорное ядро.
3. Программный комплекс РАДУГА-Т
Программный комплекс РАДУГА-Т предназначена для решения краевых задач для уравнения переноса нейтронов (1)—(2) в групповой аппроксимации (3)—(7) сеточным методом на параллельных компьютерах. Перечислим кратко его основные характеристики.
3.1. Основные возможности программного комплекса
В программном комплексе реализована возможность расчета на неструктурированной сетке, состоящей из ячеек-выпуклых многогранников. Единственное требование к сетке — отсутствие «висячих» вершин, то есть вершина каждой грани должна быть вершиной хотя бы одной другой грани. Допустимо присутствие в сетке ячеек разных типов.
В программном комплексе реализованы сеточные схемы метода конечных элементов (для ячеек-тетраэдров и треугольних призм), и сеточные схемы метода характеристик (для ячеек тетраэдров). Если для ячейки какого-либо типа не реализована сеточная схема, то расчет в ней проводится по шаговой схеме первого порядка точности [1].
Для решения системы сеточных уравнений может быть использован один из четырех итерационных методов:
- BJ (Block-Jacobi) — двойные итерации по правой части и граничным значениям решения без ускорения сходимости;
- BJ+P1 — BJ с ускорением сходимости внутренних итераций;
- BS (Block-Seidel) — итерации по правой части без ускорения сходимости;
- BS+P1 — BS с ускорением сходимости внутренних итераций. Внешние итерации выполняются без ускорения.
При вычислениях используется двухуровневый алгоритм распараллеливания, где на верхнем уровне распараллеливание выполняется по пространственным подобластям, на нижнем уровне — по направлениям скорости 0.т.
Разбиение пространственной сетки на подобласти выполняется в самом программном комплексе одним из трех методов: жадным, бисекции или координатным. Время расчета существенно зависит от использованного метода разбиения сетки [6].
Возможен счет с продолжением, счет с ненулевым начальным приближением. Отличие программного комплекса РАДУГА-Т от других комплексов решения уравнения переноса на неструктурированных сетках [2, 5, 11, 12, 15, 20] состоит в том, что программный комплекс РАДУГА-Т включает большой набор алгоритмов и методов: сеточные схемы двух типов, 4 итерационных алгоритма, 3 встроенных алгоритма декомпозиции сетки, возможность расчета на сетках с ячейками разных типов. Сравнение характеристик программного комплекса РАДУГА Т и других комплексов представлено в таблице. Прочерк означает отсутствие данных.
Таблица
Сравнение характеристик программных комплексов для решения уравнения переноса
нейтронов
Программ- Типы Сеточ- Итерационный Средства и переменные
ный ком- ячеек ные метод распараллеливания
плекс схемы
ATTILA [11] тетраэдр МКЭ BJ, BJ+DSA OpenMP (—)
THOR [20] тетраэдр МХ BJ MPI (r, W)
ODETTA [2] тетраэдр МКЭ BJ, BJ+Ребаланс OpenMP (W)
ARES [12] тетраэдр МКЭ BJ —
APOLLO [15] тетраэдр МХ BJ MPI+OpenMP (r, W, E)
PMSNSYS [5] тетраэдр МКЭ BS, BS+Ребаланс MPI (W) + OpenMP (r)
РАДУГА-Т тетраэдр МКЭ BJ, BJ+P1 MPI (r) + OpenMP (W)
треуголь- МХ BS, BS+P1
ная призма
3.2. Особенности программной реализации
Программный комплекс написан на языке FORTRAN, для параллельных вычислений используются команды MPI (верхний уровень) и OpenMP (нижний уровень).
Иные пакеты не используются. Программный комплекс совместим с операционными системами Windows и Unix.
Для экономии оперативной памяти в каждый момент времени расчета в ней находятся только те значения решения, которые относятся к рассчитываемой энергетической группе. После завершения расчета каждой группы происходят обмены с жестким диском.
Все итерационные методы в программном комплексе реализованы в виде единой системы вложенных циклов. Включение пространственных ячеек нового типа требует только подключения модуля обращения оператора L} в ячейках такого типа.
3.3. Подготовка пространственной сетки
Геометрия области расчета в задачах о переносе нейтронов в АЗ и в радиационных защитах реальных ЯЭУ является очень сложной, содержит множество разномасштабных элементов. Эта геометрия может быть задана самыми разными способами:
- с помощью одного из CAD-форматов;
- посредством специального языка в одной из программ моделирования переноса нейтронов методом Монте-Карло (отечественной MCU или зарубежной MCNP);
- с помощью бумажных чертежей.
Для построения сетки используется свободно распространяемый пакет SALOME, в котором геометрия может быть задана в отдельном входном файле на языке Python при помощи библиотеки API. Использование этого подхода позволяет сократить трудозатраты на обработку больших объемов геометрических данных, а также адаптироваться к появлению их новых форматов.
Создана специальная программа CONVERTGEOM, переводящая описание геометрии с языка MCU на язык Python. При задании геометрии в виде чертежей разработан «служебный» формат RAD описания геометрических данных. Файл в формате RAD преобразуется той же программой CONVERTGEOM в файл на языке Python. На основе этого файла с помощью пакета SALOME строится неструктурированная сетка.
Такая сетка не всегда обладает высоким качеством. В частности, сетка может быть излишне сгущена в одной подобласти и слишком разрежена в другой. Для улучшения качества сетки используется пакет ANIMBA [21].
Рис. 1. Схема подготовки пространственной сетки
Следующий этап в подготовке сетки состоит в переводе ее в формат программного комплекса РАДУГА-Т. В этом формате содержится информация о материале каждой ячейки и о соседних с ней ячейках. Такой перевод выполняется Препроцессором программного комплекса РАДУГА-Т. Им же создаются файлы в формате VTK (Visualization ToolKit) для визуализации сетки с помощью пакета PARAVIEW. Схема подготовки пространственной сетки приведена на рис. 1.
РАДУГА-Т
Бинарный файл решения
Файл решения в VTK формате
> 1
PARAVIEW
Постпроцессор
1 f
>
Функционалы решения
Визуализация решения
Рис. 2. Схема обработки решения Постпроцессором
После выполнения расчета уже Постпроцессором может быть выполнен расчет различных функционалов решения, в частности, средних значений решения в заданных телах (параллелепипедах, цилиндрах и др.). Также может быть выполнена интерполяция решения с неструктурированной сетки в точки с заданными координатами. Пакет PARAVIEW также может быть использован для визуализации решения, см. рис. 2.
4. Вычислительные эксперименты
В этом разделе приводятся результаты исследования эффективности алгоритма распараллеливания программного комплекса РАДУГА-Т. Рассматривается только метод BS+P1. Находится коэффициент ке;; в двухгрупповой (Q=2) модели легководного реактора [18], см. рис. 3, с опущенным стержнем. Используется пространственная сетка из 108 тыс. ячеек, 271 тыс. граней. Сетка по направлениям скорости содержит M = 80 узлов.
Используется схема метода характеристик. Полное число неизвестных в уравнении переноса в одной энергетической группе — 74 млн., а в системе уравнений для поправки — 3,7 млн. Для ее решения используется 39 тыс. ортонормированных векторов.
Рис. 3. Геометрия модели легководного реактора
Расчеты выполняются на многопроцессорном компьютере МВС-10П (МСЦ РАН). В расчете каждой подобласти используется 2 ядра. Одно из них выполняет вычисления, другое — межпроцессорные обмены. На одной внешней итерации (8), (9) выполняется 28 внутренних итераций (19).
Каждый расчет состоит из двух шагов: К и Р1. Каждый шаг разделяется на два этапа: вычисления и обмены, всего 4 этапа. Доли временных затрат (%) на каждый из 4-х этапов для различного числа подобластей 3 приведены на рис. 4.
Рис. 4. Распределение времени по этапам расчета
Можно видеть, сколь большую часть времени занимают обмены, причем время обменов увеличивается с ростом числа подобластей 3, особенно на К-шаге. Это связано ростом числа актов обмена, увеличивающегося из-за роста числа проходов, см. рис. 5. На Р1-шаге с ростом 3 число базисных векторов не меняется. Здесь доля времени обменов сначала растет, затем снижается.
Рис. 5. Число проходов и актов обмена на К-шаге в зависимости от числа подобластей
На рис. 6 представлены величины а ускорения вычислений всего расчета и отдельно каждого шага, найденные по формуле
а = Т(2)/ТЦ). (20)
Здесь Г(/) — время расчета на J подобластях. Идеальное ускорение аШеа1 = ]/2 на рис. 6 представлено штрихованной линией. Можно видеть, что эффективное ускорение вычислений на P1-шаге сочетается с малым ускорением на К-шаге. Более того, при увеличении числа подобластей ускорение вычислений на К-шаге прекращается. Это происходит за счет увеличения числа актов обмена данными между подобластями, см. рис. 5. Суммарное ускорение при небольшом числе подобластей определяется ускорением P1-шага. С ростом числа подобластей доля времени выполнения К-шага возрастает, см. рис. 4, и суммарное ускорение вычислений снижается.
7
Рис. 6. Ускорения вычислений в зависимости о числа подобластей
Отметим, что выбранная для тестирования стратегия — использование в каждой подобласти двух ядер, одно из которых выполняет вычисления, а другое межпроцессорные обмены — представляется оптимальной для метода BS. В BS методе на каждом проходе для каждого вектора Пт рассчитывается разное число пространственных ячеек, причем ячейки рассчитываются в строго заданном порядке. Использование для вычисления двух и более ядер в одной подобласти приводит к задержкам на синхронизацию. Для другого итерационного метода оптимальным будет другая стратегия.
Заключение
Представлен программный комплекс РАДУГА-Т, предназначенный для решения краевой задачи для уравнения переноса нейтронов в 3D областях на неструктурированных сетках. Описаны способы задания геометрии области расчета в сеткопостроителе, используемые сеточные схемы и итерационные методы решения систем сеточных уравенений, методы распараллеливания вычислений.
Приведены результаты исследования эффективности распараллеливания одного из четырех итерационных методов на модельной задаче о легководном реакторе. Выбранный метод обладает тем преимуществом, что решаемые системы линейных алгебраических уравнений не зависят от числа пространственных подобластей. Таким образом, от числа подобластей не зависит и число итераций. Такой эффект, однако, достигается за счет роста числа актов межпроцессорных обменов при увеличении числа подобластей. В рас-
сматриваемой задаче время счета с увеличением числа подобластей от 2 до 64 уменьшеа-
ется всего в 14 раз. Для других итерационных методов результаты исследования эффективности распараллеливания будут иными.
Развитие комплекса РАДУГА-Т планируется по следующим направлениям.
- Разработка сеточной схемы для расчета в ячейке-четырехугольной пирамиде; это даст возможность выполнять расчет на комбинированных сетках, содержащих ячейки-тетраэдры и ячейки-призмы.
- Разработка методов ускорения сходимости внешних итераций (8), (9) и (10), (11).
- Повышение эффективности распараллеливания вычислений на К-шаге.
Литература
1. Басс Л.П., Волощенко А.М., Гермогенова Т.А. Методы дискретных ординат в задачах о переносе излучения. Москва: ИПМ, 1986. 231 с.
2. Белоусов В.И., Грушин Н.А., Сычугова Е.П., и др. Некоторые результаты верификации кода ODETTA для неоднородных задач // Вопросы атомной науки и техники. Серия: Физика ядерных реакторов. 2018. № 3. С. 46-53.
3. Коконков Н.И., Николаева О.В. KP1 итерационный метод решения уравнения переноса в трехмерных областях на неструктурированных сетках // Журнал вычислительной математики и математической физики. 2015. Т. 55, № 10. С. 95-108. DOI: 10.7868/S0044466915100154.
4. Николаев А.А., Усенков В.В., Афанасьев П.Б., и др. Современное состояние развития программного обеспечения расчета переноса ионизирующего излучения в активных зонах и радиационной защите реакторных установок со свинцово-висмутовым теплоносителем // Вопросы атомной науки и техники. Серия: Ядерно-реакторные константы. 2017. № 1. С. 129-143.
5. Николаев А. А. Совершенствование геометрических опций SN кода PMSNSYS-II // Вопросы атомной науки и техники. Серия: Физика ядерных реакторов. 2017. № 1. С.143-147.
6. Николаева О.В., Гайфулин С.А., Басс Л.П. О декомпозиции неструктурированной сетки при решении уравнения переноса нейтронов на параллельных компьютерах // Параллельные вычислительные технологии (ПаВТ'2019): Труды международной научной конференции (Калининград, 2-4 апреля 2019 г.). Челябинск: Издательский центр ЮУрГУ, 2019. С. 362-372.
7. Николаева О.В., Гайфулин С.А., Басс Л.П. Детальное моделирование эксперимента IRON 88 на установке ASPIS в (r,z)- и (x,y,z) геометриях // Известия вузов. Ядерная энергетика. 2019. № 3. С. 135-147. DOI: 10.26583/npe.2019.3.12.
8. Николаева О.В., Казанцева А.С. Сравнение свойств сеточных схем для решения уравнения переноса на неструктурированных тетраэдрических сетках // Вопросы атомной науки и техники. Серия: Математическое моделирование физических процессов. 2019. № 1. С. 3-18.
9. Николаева О.В., Казанцева А.С. Точность схем метода конечных элементов для решения уравнения переноса на неструктурированных тетраэдрических и призматических сетках // Вопросы атомной науки и техники. Серия: Математическое моделирование физических процессов. 2020. № 1. С. 3-18.
10. Селезнев Е.Ф., Березнев В.П. Использование диффузионного приближения при расчете реактора с полостями // Известия вузов. Ядерная энергетика. 2018. № 2. С. 6777. DOI: 10.26583/npe.2018.2.07.
11. Vassiliev O.D., Wareing T.A., Davis I.M., et al. Feasibility of a Multigroup Deterministic Solution Method for 3D Radiotherapy Dose Calculations // International Journal of Radiative Oncology, Biology, Physics. 2008. Vol. 72 P. 220-227. DOI: 10.1016/j.ijrobp.2008.04.0572017.
12. Chen Y., Zhang B., Zhang L., et al. ARES: A Parallel Discrete Ordinates Transport Code for Radiation Shielding Applications and Reactor Physics Analysis // Hindawi Science and Technology of Nuclear Installations. 2017. Article ID 2596727. DOI: 10.1155/2017/2596727.
13. Colomer G., Borrell R., Trias F.X., et al. Parallel Algorithms for Sn Transport Sweeps on Unstructured Meshes // Journal of Computational Physics. 2013. Vol. 232. P. 118-135. DOI: 10.1016/j.jcp.2012.07.009.
14. Kim J.W., Lee Y.O. A Deep Penetration Problem Calculation Using AETIUS: An Easy Modeling Discrete Ordinates Transport Code UsIng Unstructured Tetrahedral Mesh, Shared Memory Parallel // EPJ Web of Conferences. 2017. Vol. 153. P. 06025. DOI: 10.1051/epjconf/20171530.
15. Lenain R., Masiello E., Damian F., et al. Domain Decomposition Method for 2D and 3D Transport Calculations Using Hybrid MPI/OPENMP Parallelelizm // Mathematics and Computation (M&C), Supercomputing in Nuclear Applications (SNA) and the Monte Carlo (MC) Method: Proc. of the Joint International Conference (Nashville, Tennessee, April, 19-23, 2015). LaGrange Park, IL, American Nuclear Society, 2015. URL: https://hal-cea.archives-ouvertes.fr/cea-02506817/document (дата обращения: 14.09.2020).
16. Pautz Sh.D. An Algorithm for Parallel Sn Sweeps on Unstructured Meshes // Nuclear Science and Engineering. 2002. Vol. 140, no 2. P. 111-136. DOI: 10.13182/NSE02-1.
17. Plimpton S.J., Hendrickson B., Burns Sh.P., et al. III, Rauchwerger L. Parallel Sn Sweeps on Unstructured Grids: Algorithms for Prioritization, Grid Partitioning, and Cycle Detection // Nuclear Science and Engineering. 2005. Vol. 150. P. 267-283. DOI: 10.13182/NSE150-267.
18. Takeda T., Ikeda H. 3-D Neutron Transport Benchmarks // Journal of Nuclear Science and Technology. 1991. Vol. 28, no 7. P. 656-669. DOI: 10.3327/jnst.28.656.
19. Vega R.M., Adams M.L. Transport Sweeps Using an Improved Slice Balance Approach with LDFE and GPU Acceleration // Mathematics & Computational Methods Applied to Nuclear Science & Engineering: Proc. of the International Conference (Jeju, Korea, April, 16-20, 2017). URL: https://www.kns.org/files/int_paper/paper/MC2017_2017_1/ P056S01-11VegaR.pdf (дата обращения: 14.09.2020).
20. Yessayan R., Azmy Y., Schunert S. Development Of A Parallel Performance Model For The THOR Neutral Particle Transport Code // Mathematics & Computational Methods Applied to Nuclear Science & Engineering: Proc. of the International Conference (Jeju, Korea, April, 16-20, 2017). URL: https://www.osti.gov/servlets/purl/1369430 (дата обращения: 14.09.2020).
21. Adaptive Numerical Instruments 3D. URL: https://sf.net/p/ani3d/ (дата обращения: 01.09.2020).
Николаева Ольга Васильевна, к.ф.-м.н., с.н.с., Институт прикладной математики
им. М.В. Келдыша РАН (Москва, Российская Федерация)
Гайфулин Сергей Андреевич, н.с., Институт прикладной математики им. М.В. Келдыша РАН (Москва, Российская Федерация)
Басс Леонид Петрович, к.ф.-м.н., с.н.с., Институт прикладной математики
им. М.В. Келдыша РАН (Москва, Российская Федерация)
DOI: 10.14529/cmse210106
CODE RADUGA T FOR SIMULATING NEUTRONS FLUXES IN NUCLEAR POWER STATIONS
© 2021 O.V. Nikolaeva, S.A. Gaifulin, L.P. Bass
Keldysh Institute of Applied Mathamatics (Miusskaya Sq. 4, Moscow, 125047 Russia) E-mail: nika@kiam.ru, sagarp@kiam.ru, bass@kiam.ru Received: 14.07.2020
Design and operation of nuclear power stations (NPS) are followed by simulation of neutron propagation in these stations. It is necessary to consider borders of multi-scale structural elements consisting of different materials. It is desirable to use parallel computer because of big size of NPSs. To keep such conditions, algorithms and codes for solving the integro-differential transport equation on unstructured grids are being developed. In the paper such algorithms included into the code RADUGA T are outlined. Grids, grid schemes, iterative methods to solve grid equations are presented. Calculation parallelization methods for hybrid computers are considered, MPI and OpenMP techniques are used. Methods of building, improvement, decomposition and visuzlisation of spatial grids are considered. Software implementation is described. The algorithms of the code RADUGA T are compared with methods of other codes. Computation parallelization efficiency study results are presented. The problem of neutron multiplication factor calculation in a light-water reactor model is solved. The multi-processor computer MVS-10P of the Joint SuperComputernal Center is used. Acceleration of each algorithm being used and summary acceleration are given. Keywords: transport equation, unstructured frids, grid schemes, iteration methods, parallel calculation.
FOR CITATION
Nikolaeva O.V., Gaifulin S.A., Bass L.P. Code RADUGA T for Simulating Neutrons Fluxes in Nuclear Power Stations. Bulletin of the South Ural State University. Series: Computational Mathematics and Software Engineering. 2021. Vol. 10, no. 1. P. 75-89. (in Russian) DOI: 10.14529/cmse210106.
This paper is distributed under the terms of the Creative Commons Attribution-Non Commercial 3.0 License which permits non-commercial use, reproduction and distribution of the work without further permission provided the original work is properly cites.
References
1. Bass L.P., Voloschenko A.M., Germogenova T.A. Discrete Ordinate Method in Radiation Transport Problem. Moscow, KIAM Press, 1986. 231 p. (in Russian)
2. Belousov V.I., Grushin N.A., Sychugova E.P., et al. Some Results of Verification of Code ODETTA for Shielding Calculations. Nuclear Science and Technology: Nuclear Reactor Physics. 2018. No. 3. P. 46-53. (in Russian)
3. Kokonkov N.I., Nikolaeva O.V. An iterative KP1 Method for Solving the Transport Equation in 3D Domains on Unstructured Grids. Computational Mathematics and Mathematical Physics. 2015. Vol. 55, no. 10. P. 1698-1712. DOI: 10.7868/S0044466915100154.
4. Nikolaev A.A., Usenkov V.V., Afanasiev P.B., et al. Actual State of Development of the Software for Radiation Transport Calculatons in Core and Shield of Reactor Facities with Lead-Bismuth Coolant. Nuclear Science and Technology: Nuclear Reactor Physics. 2017. No. 1. P. 143-147. (in Russian)
5. Nikolaev A.A. Improvement of Geometric Options of SN-code PMSNSYS-II. Nuclear Science and Technology: Nuclear Reactor Physics. 2017. No. 1. P. 129-143. (in Russian)
6. Nikolaeva O.V., Gaifulin S.A., Bass L.P. On decomposition of a unstructured grid to solving the neutron transport equation on parallel computers. Parallel Computational Technologies (PCT'2019): Proceedings of the International Scientific Conference (Kaliningrad, Russia, April, 2-4, 2019). Chelyabinsk, Publishing of the South Ural State University, 2019. P. 362-372. (in Russian)
О.В. HHKO^aeBa, C.A. raft^ynHH, ^.n. Bacc
7. Nikolaeva O.V., Gaifulin S.A., Bass L.P. Detailed Simulation of Winfrith IRON 88 Benchmark (ASPIS) in (r, z)- and (x, y, z)-Geometries. Izvestiya vuzov. Yadernaya Energetika. 2019. No. 3. P. 135-147. DOI: 10.26583/npe.2019.3.12. (in Russian)
8. Nikolaeva O.V., Kazantseva A.S. Comparison Between the Properties of Grid Schemes for Solving the Transport Equation on Unstructured Thetrahedral Grids. Nuclear Science and Technology. Series: Mathematical Modeling of Physical Processes. 2019. No. 1. P. 3-18. (in Russian)
9. Nikolaeva O.V., Kazantseva A.S. Accuracy of FEM Schemes for Solving the Transport Equation on Unstructured Thetrahedral and Prismatic Grids. Nuclear Science and Technology. Series: Mathematical Modeling of Physical Processes. 2020. No. 1. P. 3-19. (in Russian)
10. Seleznev E.F., Bereznev V.P. Using the Diffusive Approximation for Reactor with Cavities Calculation. Izvestiya Vuzov. Yadernaya Energetika. 2018. No. 2. P. 66-67. DOI: 10.26583/npe.2018.2.07. (in Russian)
11. Vassiliev O.D., Wareing T.A., Davis I.M., et al. Feasibility of a Multigroup Deterministic Solution Method for 3D Radiotherapy Dose Calculations. International Journal of Radiative Oncology, Biology, Physics. 2008. Vol. 72. P. 220-227. DOI: 10.1016/j.ijrobp.2008.04.0572017.
12. Chen Y., Zhang B., Zhang L., et al. ARES: A Parallel Discrete Ordinates Transport Code for Radiation Shielding Applications and Reactor Physics Analysis. Hindawi Science and Technology of Nuclear Installations. 2017. Article ID 2596727. DOI: 10.1155/2017/2596727.
13. Colomer G., Borrell R., Trias F.X., et al. Parallel Algorithms for Sn Transport Sweeps on Unstructured Meshes. Journal of Computational Physics. 2013. Vol. 232. P. 118-135. DOI: 10.1016/j.jcp.2012.07.009.
14. Kim J.W., Lee Y.O. A Deep Penetration Problem Calculation Using AETIUS: An Easy Modeling Discrete Ordinates Transport Code UsIng Unstructured Tetrahedral Mesh, Shared Memory Parallel. EPJ Web of Conferences. 2017. Vol. 153. P. 06025. DOI: 10.1051/epjconf/20171530.
15. Lenain R., Masiello E., Damian F., Sanchez R. Domain Decomposition Method for 2D and 3D Transport Calculations Using Hybrid MPI/OPENMP Parallelelizm. Mathematics and Computation (M&C), Supercomputing in Nuclear Applications (SNA) and the Monte Carlo (MC) Method: Proceedings of the Joint International Conference (Nashville, Tennessee, April, 19-23, 2015). LaGrange Park, IL, American Nuclear Society, 2015. URL: https://hal-cea.archives-ouvertes.fr/cea-02506817/document (accessed: 14.09.2020).
16. Pautz Sh.D. An Algorithm for Parallel Sn Sweeps on Unstructured Meshes. Nuclear Science and Engineering. 2002. Vol. 140, no. 2. P. 111-136. DOI: 10.13182/NSE02-1.
17. Plimpton S.J., Hendrickson B., Burns Sh.P., et al. Parallel Sn Sweeps on Unstructured Grids: Algorithms for Prioritization, Grid Partitioning, and Cycle Detection. Nuclear Science and Engineering. 2005. Vol. 150. P. 267-283. DOI: 10.13182/NSE150-267.
18. Takeda T., Ikeda H. 3-D Neutron Transport Benchmarks. Journal of Nuclear Science and Technology. 1991. Vol. 28, no. 7. P. 656-669. DOI: 10.3327/jnst.28.656.
19. Vega R.M., Adams M.L. Transport Sweeps Using an Improved Slice Balance Approach with LDFE and GPU Acceleration. Mathematics & Computational Methods Applied to Nuclear Science & Engineering: Procedeengs of the International Conference (Jeju, Korea, April, 16-20, 2017). URL: https://www.kns.org/files/int_paper/paper/MC2017_2017_1/ P056S01-11VegaR.pdf (accessed: 14.09.2020).
20. Yessayan R., Azmy Y., Schunert S. Development Of A Parallel Performance Model For The THOR Neutral Particle Transport Code. Mathematics & Computational Methods Applied to Nuclear Science & Engineering: Procedeengs of the International Conference (Jeju, Korea, April, 16-20, 2017). URL: https://www.osti.gov/servlets/purl/1369430 (accessed: 14.09.2020).
21. Adaptive Numerical Instruments 3D. URL: https://sf.net/p/ani3d/ (accessed: 14.09.2020).
