CC BY
48
УДК 519.624.3
Параллельный алгоритм и MPI реализация численного исследования фазовых переходов на основе 3D модели
термического пика
И. В. Амирханов, Е. В. Земляная, Н. Р. Саркар, И. С. Сархадов, З. К. Тухлиев, З. А. Шарипов
Лаборатория информационных технологий Объединённый институт ядерных исследований ул. Жолио-Кюри, д. 6, Дубна, Московская область, Россия, 141980
Мы представляем вычислительную схему и параллельную компьютерную реализацию для численного исследования эволюции температурных полей и фазовых переходов в материалах под действием облучения тяжёлыми ионами высоких энергий. Используется модифицированная модель термического пика, которая описывается системой двух связанных уравнений теплопроводности для температуры электронного газа и температуры ионной кристаллической решётки облучаемого материала. Численное решение этой системы осуществляется на основе условно-устойчивой явно-неявной конечно-разностной схемы в цилиндрической системе координат с использованием разложения функции источника по сферическим гармоникам для учёта нарушения аксиальной симметрии в моделируемой системе (3D). Моделирование динамики фазовых переходов реализовано на основе энтальпийного подхода. Представлена математическая постановка задачи, описана вычислительная схема, приведены особенности параллельной компьютерной реализации на базе технологии MPI (Message Passing Interface). Представлены результаты методических расчётов, проведённых на многопроцессорном кластере К100 (ИПМ РАН, Москва) с различным числом узлов конечно-разностной сетки и с разным числом параллельных процессоров, демонстрирующие эффективность разработанной параллельной CH—b/MPI-программы.
Ключевые слова: моделирование, численные методы, фазовый переход, модель термического пика, параллельный алгоритм.
1. Мотивация и постановка задачи
Радиационная физика представляет собой одно из быстро развивающихся направлений применения математического моделирования. Активные исследования в этой области показали, что радиационно-стимулированные процессы приводят к качественному изменению свойств облучаемых материалов. В основе этих изменений лежат различные физические процессы (тепловые, диффузионные и пр.). Проведение натурных экспериментов в этих областях сопряжено с большими трудностями, поскольку исследуемые процессы происходят за очень малые времена t ~ 10-7 — 10-15с. Поэтому особенно важную роль приобретает развитие методов компьютерного моделирования, в том числе разработка параллельных алгоритмов для повышения эффективности численного исследования. В данной работе представлена параллельная реализация алгоритма для численного решения системы уравнений, описывающей модифицированную модель термического пика (МТП) для исследования тепловых процессов в облучаемых тяжелыми ионами материалах в аксиально несимметричном (3D) случае с учётом фазовых переходов.
МТП описывается системой двух уравнений теплопроводности для электронного газа и кристаллической решётки. Эта система в цилиндрической системе координат имеет вид [1]:
Статья поступила в редакцию 22 сентября 2013 г. Работа выполнена при финансовой поддержке Минобрнауки России в рамках государственного контракта №07.524.12.4019 от 17.05.2012.
Амирханов И. В. и др. Параллельный алгоритм и MPI реализация . . . 207
- д(Те)(Те -Тг) + Ае(г,p,z, t), (1)
™)^ = (1 £ ж) + Тг£ ^) +1 §)) +
+ 9(Тг)(Те -Т.) + А.(г,р,г, V). (2)
Здесь ось г направлена вдоль пучка и перпендикулярна облучаемой поверхности; Те(г, р, г, V), Т.(г, р, г, V) — температурные поля электронного газа и решётки облучаемого образца; Се, Сг и Ае, А. — соответственно, удельные теплоёмкости и коэффициенты теплопроводности электронного газа и решётки; д — коэффициент электрон-фононного взаимодействия электронной подсистемы с решёткой. Ае(г, р, г, V) и А.(г, р, г, V) — объёмные плотности энергии, вносимые ионом в электронную и решёточную подсистемы, имеющие вид:
Ае,г (г ,P,Z, ^ = /1(г) ¡2(р) /з(^5,1пв1,рЬ(г).
Нормированные на единицу функции /¿, (г = 1,2, 3) выбраны следующим образом:
fi(r) = -1 exp(-r/го), ¡2(р) = 2" ( 1 + cosmp + c2m sinmp)
h (t) 1
V2"s t
/ M \
1 + (cim cos mp + C2m sin mp) V m=1 )
2 \ ^max
ex^- ^ 2 Г ) И / ( S'mel(z) + SPh(z))dz = Ео.
Здесь го, V о, ^ — физические параметры [2], Е0 — энергия падающего иона, функции 5;пв1(г) и 5Рь(2) вычисляются программой 8ШМ-2012 (http://www.srim.org) и определяют соответственно энергетические потери иона на возбуждение электронной и фононной подсистем в зависимости от глубины мишени. Предполагая, что параметры Се, С., Ае, А.,д постоянные, т.е. не зависят от Те., решение системы (1)—(2) будем искать в виде:
м
Те,г(г, р, г, V) = То,е,г(г, г, V) + ^ Тт,е(г, г, 1)(с\т совтр + с2т ъ'ттр).
Ш=1
Тогда для функции Тт,е,г(г, г, г), (т = 0,1, 2,..., М) получаем систему уравнений:
г. 9Тт,е . fldi дТте\ m2
Се^7~ = \е Г—---^- +
д д
е д е д д 2 д 2
- д(Тт,е - Тт,г) + Ат,е(г, Z, t), (3)
£, дТт^ /1 д [дТтЛ m2 Тт,,г + д Тт,,г \ + д д д - 2 д 2
+ д(Тт,е - Тт,i ) + Am,i (г ,z, t). (4)
Система (3)—(4) решается независимо для каждого значения т с соответствующими начальными и граничными условиями (см. [2]).
2. Параллельная реализация и обсуждение результатов
Моделирование динамики фазовых переходов типа «плавление - затвердевание» осуществляется на основе задачи Стефана в рамках энтальпийного подхода [3]. Для численного решения полученной системы уравнений использована условно устойчивая явно-неявная конечно-разностная схема, подробно описанная в [4]. Соответствующие разностные уравнения являются неявными только по переменной г и остаются явными по переменной z. Поэтому численное решение сводится к применению L одномерных прогонок относительно г (где L — число узлов дискретной сетки по z). Выбор данной схемы обусловлен её экономичностью и простотой компьютерной, в том числе параллельной реализации.
Компьютерная программа, реализующая указанную схему для модели МТП, написана на языке СН—Н с использованием технологии MPI для организации параллельных вычислений. Алгоритм распределения вычислений по MPI-процессам следующий.
Пусть Np — количество MPI-процессов, одновременно участвующих в решении задачи. Как и выше, M — число гармоник, L — число узлов равномерной сетки по переменной z.
При малом числе параллельных процессов Np ^ M распараллеливание осуществляется только по индексу т: значения m = 0,..., M — 1 распределяются между процессами. Каждый 1-й процесс вычисляет для назначенных ему значений m полные матрицы Те и Ti температур электронной и ионной подсистем в узлах дискретной сетки по z,r. В этом случае отсутствует обмен данными между процессами во время счёта. Взаимодействие процессов происходит лишь на стадии сборки и суммирования по m, в рамках процедуры сохранения результатов. Оптимальное распределение нагрузки на процессы обеспечивается, если M кратно Np.
Если Np > M, MPI-процессы распределяются сначала между значениями m = 0,..., M — 1 так, что каждому значению m назначается число процессов N(mK Далее процессам, предоставленным каждому m, назначаются интервалы
Di = (IL/N(m) + 1, (l + 1)L/N(m)) с номерами 1(1 = 0,1, 2,..., N(m) — 1). Каждый процесс Р ведёт расчёты в своих интервалах I для значений jm-m(P) < j < jmax(P), j = 0,... L — 1. На каждом временном слое для корректного продолжения счёта необходимо обновление граничных значений фрагментов матриц Te i, для чего организован обмен соответствующими данными между соседними процессами: каждый Р-й процесс передаёт и получает обновлённые значения матриц Те,г, соответствующие граничным значениям интервалов Di, от соседних процессоров Р — 1 и Р +1. Сборка и суммирование по m осуществляются в конце счёта и на некоторых, заранее заданных, промежуточных временных слоях.
В табл. 1 представлены результаты тестовых расчётов с разными значениями L и Np, проведённые на кластере K100 (ИПМ РАН, Москва) и демонстрирующие уменьшение времени счёта с ростом числа процессов. Здесь M = 5, а число узлов дискретной сетки по переменной г составляет N = 500. В табл. 2 представлены аналогичные результаты для разных N и Np при L = 5000.
Таблица 1
Время работы (в минутах) MPI-программы, реализующей МТП, при M = 5, N = 500, т = 0.01 х 10-13c в зависимости от числа MPI-процессов Np и числа
узлов L по z.
L Np = 1 Np = 5 Np = 10 Np = 20 Np = 30 Np = 40
10000 72.6 19.2 10.0 7.2 5.6 5.1
15000 109.3 26.1 16.1 9.3 8.3 7.2
20000 161.3 38.0 21.3 13.0 11.0 9.5
Амирханов И. В. и др. Параллельный алгоритм и MPI реализация .
209
Таблица 2
Время работы (в минутах) MPI-программы, реализующей МТП, при М = 5, N = 5000, т = 0.01 х 10-13е в зависимости от числа MPI-процессов Мр и числа
узлов Ь по г.
N Np = 1 Np = 5 Np = 10 Np = 20
500 37.22 9.47 4.2 3.32
1000 82.31 15.71 11.31 8.37
2000 162.23 34.03 22.57 14.26
4000 247.12 96.12 74.28 44.58
В целом расчёты подтверждают эффективность разработанной С++/МР1-программы. Видно, что наибольшая эффективность с точки зрения соотношения «число процессов - время счёта» при указанных параметрах достигается при Мр ~ 20. Далее ускорение вычислений замедляется из-за возрастающего объёма пересылаемых данных при сборке конечных результатов.
Отметим, что тестирование без учёта затрат на конечную сборку и сохранение результатов показывает, что затраты на обмен данными между соседними МР1-процессами на каждом временном слое незначительны по сравнением с затратами на собственно вычисления и практически не влияют на динамику ускорения счёта.
Отметим в заключение, что в [5] были проведены, согласно представленной здесь схеме, оценки размеров областей на поверхности мишени, где под действием ионного облучения температура мишени превышает температуру плавления материала мишени. Наличие таких областей принято интерпретировать в рамках МТП как образование треков. Результаты наших расчётов согласуются с известными экспериментальным оценками, причём согласие улучшается при учёте фазовых переходов. Это свидетельствует об адекватности используемой модели и подтверждает корректность и эффективность вычислительной схемы и компьютерной С++/МР1-реализации.
Литература
1. Каганов М. И., Лифшиц И. М, Танатаров Л. В. Релаксация между электронами и решёткой // ЖЭТФ. — 1956. — № 2(8). — С. 232-237. [Kaganov M.I., Lifschitz I.M., Tanatarov L.V. Relaxation Between Electrons and the Lattice // JETF. — 1956. — No 2(8). — Pp. 232-237. ]
2. Распыление твердых тел под действием тяжелых ионов и температурные эффекты в электронной и решеточной подсистемах / И. В. Амирханов, А. Ю. Ди-дык, И. В. Пузынин и др. // Физика элементарных частиц и атомного ядра. — 2006. — Т. 37, № 6. — С. 1592-1644. [Sputtering of Solids by Heavy Ions and Temperature Effects in the Electronic and Lattice Subsystems / I. V. Amikhanov, A. Yu. Didyk, I.V. Puzynin et al. // PEPAN. — 2006. — Vol.37, No6. — Pp. 15921644.]
3. Численное моделирование динамики температурных полей на плоских мишенях при нестационарном интенсивном лазерном воздействии. / М. П. Галанин, И. С. Ерхов, Е. Ю. Локтионов и др. // Препринт ИПМ им. М. В. Келдыша РАН. — 2008. — № 61. [Numerical Modeling of Temperature Fields on a Flat Target at Unsteady Intense Laser Pulses / M.P. Galanin, I. S. Erhov, E. Yu. Loktionov et al. // Preprint IPM im. M.V. Keldisha RAN. — 2008. — No 61. ]
4. MPI реализация алгоритмов для 2D и 3D моделирования фазовых переходов в материалах, облучаемых тяжелыми ионами, в рамках модели термического пика / И. В. Амирханов, Е. В. Земляная, Н. Р. Саркар и др. // Вестник РУДН. Серия «Математика. Информатика. Физика». — 2013. — № 4. — С. 80-94. [MPI
Implementation to the 2D and 3D Simulation of Phase Transitions in Materials Irradiated by Heavy Ion Beams within the Thermal Spike Model / I. V. Amikhanov, E. V. Zemlyanaya, N. R. Sarker et al. // Bulletin of Peoples' Friendship University of Russia. Series "Mathematics. Information Sciences. Physics" . — 2013. — No 4. — Pp. 80-94. ]
5. Microstructural Modifications Induced by Swift Ions in the NiTi Intermetallic Compound / A. Barbu, A. Dunlop, A. Hardouin et al. // Nucl. Instrum. Meth. — 1998. — No 145.
UDC 519.624.3
Parallel Algorithm and MPI Implementation of Numerical Study of Phase Transition in the 3D Thermal Spike Model
I.V. Amirkhanov, E. V. Zemlyanaya, N.R. Sarker, I. S. Sarkhadov, Z.K. Tukhliev, Z. A. Sharipov
Laboratory of Information Technologies Joint Institute for Nuclear Research 6, Joliot-Curie str., Dubna, Moscow region, Russia, 141980
We present an algorithm and parallel computer code for numerical investigation of the thermal processes and phase transitions in materials irradiated by the high energy heavy ion beams. We employ the modified thermal spike model based on the coupled heat conductivity equations for the electron gas and the ion lattice subsystems in the target sample. This system of equations is numerically solved in the cylindrical coordinate system in axially non-symmetric (3D) case. We utilize an expansion of the source function in spherical harmonics, a finite difference approximation and semi-explicit numerical scheme. The dynamics of phase transitions is implemented on the basis of the enthalpy approach. The mathematical formulation of the problem is given; a numerical scheme is described; a parallel algorithm is presented on the basis of the MPI technique (Message Passing Interface). The test calculations on the K100 multi-processor cluster (KIAM RAS, Moscow) with various dimension of the finite-difference mesh and with different number of parallel processors demonstrate efficiency of the C++/MPI code.
Key words and phrases: modeling, numerical methods, phase transitions, thermal spike model, parallel algorithm.
