Филиал ТюмГНГУ в г.Сургуте1, НИУ МАИ филиал «Стрела»2, Сургутский государ- ственный университет3

Abstract

В данной работе рассматривается использование технологий параллельного программирования CUDA и OpenCL в задаче численного моделирования отклонения быстрых заряженных частиц изогнутым кристаллом. Разработан алгоритм расчета непрерывного потенциала и атомного коэффициента диффузии с помощью графического ускорителя. Ключевые слова: моделирование атомного коэффициента диффузии, разложение в тригонометрический ряд Фурье, суперкомпьютер АПК-5 ФГУП «РФЯЦ-ВНИИЭФ», графический ускоритель.

-

Использование параллельных технологий в задаче численного моделирования отклонения быстрых заряженных частиц изогнутым кристаллом* диусы изгиба кристалла в направлении оси ОХ и OY, соответственно, которые изменяются с глубиной z; z    t ; m  m0 ;  1   2 12 − Лоренц-фактор;    c ; с – скорость света;  − скорость заряженной частицы в направлении оси OZ; m0 - масса покоя; U x, y − непрерывный потенциал [ 13 ], который имеет вид  2 g 2  U r  d 3  V (g)  exp    cosg  r  ri , g i  2  (2) где V ( g ) − компонента Фурье изолированного атома; g 2  g x2  g 2y − модуль вектора обратной решетки; g  g x ; g y ; g x  2nx ax ; g y  2ny a y ; nx , n y − целые числа; ax , a y − периодичности расположения атомов кристалла по направлениям осей ОХ, OY, соответственно; 2 − средний квадрат флуктуаций тепловых колебаний атомов кристалла; d − постоянная кристаллической решетки; r  x; y; ri − координаты атомов кристалла.

Многократное рассеяние каналированных частиц описывается с помощью дифференциального уравнения для скорости изменения среднего квадрата флуктуаций поперечной энергии E , построенного в [5] или с помощью системы из десяти дифференциальных уравнений для средних квадратов флуктуаций координат и скоростей, которые были получены в [6].

Важную роль при описании многократного рассеяния играет выбор коэффициента диффузии, который описывается с помощью компонент диффузионной матрицы вида  Dxx x, y Dxy x, y Dx, y    .

 Dyx x, y Dyy x, y Если суммарный коэффициент диффузии вычисляется в приближениях [7,8], то, как было показано в [ 11, 13 ], Dxx x, y  Dyy x, y и Dxy x, y  Dyx x, y  0 . Одна из компонент атомного коэффициента диффузии каналированных частиц в окрестности изолированной атомной цепочки имеет вид [ 9-11 ]

a  V r  V r  Dxy x, y  z    x y

T  V r  x

T V r  y   , (3) T  где r  x 2  y 2 ;  T − усреднение по независимым тепловым колебаниям атомов кристалла; V r  — потенциальная энергия взаимодействия каналированной частицы изолированной атомной цепочки.

Выражение (3) при разложении в двойной тригонометрический ряд Фурье и использовании усреднений [ 13 ] определит одну из компонент атомного коэффициента диффузии каналированных частиц

Dxy  x, y   Dyx  x, y   4 2a

z    nxmyV  g V q   axayd 6 g q i  exp   2 g  q 2 2  exp   2  g 2  q2  2  sin g r  ri   sin q r  ri    , (4) где q  qx , q y ; q 2  q x 2  q y 2 ; qx  2mx ax , q y  2m y a y ; mx, my − целые числа; az − 2 период кристаллической решетки в направлении оси OZ; g  q2   gx  qx 2   g y  qy  . 3. Численный расчет

Программный комплекс TROPICS [4] основан на численном решении системы уравнений (1) с уравнением, полученного в [5] или системой уравнений [6], что эквивалентно решению кинетического уравнения Фоккера−Планка [5,6]. Это позволяет решать задачу отклонения быстрых заряженных частиц изогнутыми кристаллами для одной частицы, поэтому если задать пучек частиц с заданным размером пучка и расходимостью, то в результате задача легко масштабируется для большого числа частиц в пучке. Для решения системы дифференциальных уравнений используется метод Рунге-Кутты 4 порядка, который в программе реализован с использованием векторного представления данных и специального класса valarray языка С++ библиотеки STL. При интегрировании системы дифференциальных уравнений для одной частицы необходимо несколько раз вычислять значения выражений (2), (4) в результате, если использовать центральный процессор время расчета для электрона с энергией 855 МэВ, который отклоняется (111) плоскостным каналом изогнутого кристалла кремния толщиной 30 мкм с шагом 1.8 нм составляет 16 и 28 минут при использовании подходов [5] и [6], соответственно. В связи с этим чтобы минимизировать время интегрирования используется аппроксимация для производных непрерывного потенциала и коэффициента диффузии. В результате анализа скорости аппроксимации выявлено, что большая часть времени затрачивается на расчет значений производных от непрерывного потенциала и коэффициента диффузии. Алгоритм [ 14 ] расчета непрерывного потенциала и атомного коэффициента диффузии с использованием технологий параллельного программирования CUDA и OpenCL был включен в состав программного комплекса TROPICS. Время расчета одного значения непрерывного потенциала и атомного коэффициента диффузии сильно зависят от количества членов разложения ряда Фурье nx , ny , mx , my , а также количества атомов кристалла. Если nx  ny  mx  my  50 , то скорость расчета 100 значений одной из компонент атомного коэффициента диффузии составляет 2062, 180 и 3 секунд на четырех ядерном центральном процессоре Core i5-3450, суперкомпьютере АПК-5 ФГУП «РФЯЦ-ВНИИЭФ» и графическом ускорителе Geforce GTX-660Ti, соответственно. При разработке алгоритма было применено несколько оптимизаций [ 15 ]: одна из которых понижение точности вычислений за счет использования типа данных float; использование тригонометрических функций пониженной точности для sin(x), cos(x), exp(x); а также использование заранее вычисленных контакт и функций. Применение этих оптимизаций позволило улучшить быстродействие, сэкономить память ускорителя, а также получить относительную погрешность на уровне 2.65%. Таким образом с помощью реализации алгоритма вычисления атомного коэффициента диффузии на графическом ускорителе позволило добиться ускорения вычислений в 700 раз по сравнению с вычислениями на центральном процессоре. Стоит отметить что результат сильно зависит атомного форм-фактора изолированного атома [ 16 ]. В таблице 1 представлена зависимость времени расчета одного значения от количества членов разложения на разных вычислительных устройствах. Видно, что технологии CUDA и OpenCL показывают аналогичные результаты, а центральный процессор не позволяет использовать атомный коэффициент диффузии в своих расчетах. Таблица 1. Зависимость времени расчета одного значения атомного коэффициента диффузии от количества членов разложения N.

N 50 100 200 300 tCUDA, сек tOpenCL, сек tCPU, сек 0.04 1 9 41 0.05 1.2 10 43 20 333 5328 Скорость аппроксимации производных непрерывного потенциала в результате использования вычислительных мощностей видеоускорителя занимает до 3 минут и 10 часов для атомного коэффициента диффузии на сетке 256х256. В итоге время решения задачи отклонения одного электрона с энергией 855 МэВ, который отклоняется (111) плоскостным каналом изогнутого кристалла кремния толщиной 30 мкм с шагом 1.8 нм составило меньше секунды.

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

Алгоритм работы программы можно описать следующим образом. Задаются начальные условия: 1) параметры пучка частиц, их количество, масса, энергия и заряд частиц; 2) тип кристалла, его температура, толщина, изгиб и кристаллографическое направление; 3) выбирается потенциал изолированного атома; 3) выбирают подход для моделирования многократного рассеяния [5] или [6]; 4) задают параметры интегрирования и аппроксимации. После определения начальных условий выполняется аппроксимация производных непрерывного потенциала и коэффициента диффузии. Далее выполняется расчет траекторий движения частиц. По завершению расчета можно построить графики траекторий движения частиц, угловые распределения и ионизационные потери энергии.

Выполнено моделирование отклонения быстрых заряженных частиц изогнутым кристаллом с учетом суммарного коэффициента диффузии и атомного коэффициента диффузии с начальными условиями [1]. В компьютерном эксперименте использовался структурный и атомный форм−факторы, а, также, фактор Дебая−Валлера. Средний квадрат флуктуаций тепловых колебаний атомов кристалла рассчитывался по алгоритму [ 17 ] при температуре 294 ᵒК кристалла кремния. Для потенциала изолированного атома использовалось приближение Дойля−Тернера [ 18 ]. В плоскости XOZ лежит как радиус изгиба Rx плоскости (111), так и угол вылета частиц θx, которые измерялись относительно первоначального направления пучка. Угол разориентации между осью <110> и направлением падения частиц в плоскости (111) был взят равным 50 мрад, так как в работе [1] не было указано его значение. Начальные значения точек влета были равномерно распределены в следующих пределах ∆x ∙ ∆y = 0.07 мм ∙ 0.2 мм = 0.014 мм2, а углы влета нормально распределены вокруг среднего значения согласно алгоритму [ 19 ]. Модельное среднеквадратичное отклонение, которое определяет угловую расходимость пучка частиц, было выбрано x=30 и y=70 мкрад. Толщина кристалла 30.5 мкм. При численном решении системы дифференциальных уравнений движения шаг интегрирования равнялся 1.8 нм. Расчеты были произведены для 100000 частиц. Время расчета каждого распределения составило около 40 минут на четырех ядерном процессоре с частотой 3 ГГц. На рис. 1 представлены экспериментальные [1] и расчетные угловые распределения электронов с энергией 855 МэВ, которые отклонялись (111) плоскостным каналом изогнутого кристалла кремния с радиусами изгиба Rx=33.5 мм и Ry=0, соответственно. Рис. 1. Угловое распределение пучка электронов с энергией 855 MeV в (111) плоскостном канале кремния. На рисунке эксперимент [1] обозначен кружками, результаты компьютерного моделирования: сплошная линия с использованием атомного коэффициента диффузии [ 11 ] и крестики – расчет с использованием суммарного коэффициента диффузии в приближении [7,8]. Разрешение детектирующей системы равнялось 30 μrad. Литература

Mazzolari A. et al. Steering of a sub-GeV electron beam through planar channeling enhanced by Rechanneling // Physical Review Letters. 2014. Vol. 112. P. 135503.

Wistisen T.N. et al. A study of Channeling, Volume Reflection and Volume Capture of 3.35 -14.0 GeV Electrons in a bent Silicon Crystal // SLAC-PUB-16435. 2015. URL: http://slac.stanford.edu/pubs/slacpubs/16250/slac-pub-16435.pdf (дата обращения: 16.01.2016). Штанов Ю.Н. Программный комплекс TROPICS. URL: http://wwwinfo.jinr.ru/programs/jinrlib/tropics/index.html (дата обращения: 30.11.2015). 5. Кощеев В.П., Моргун Д.А., Штанов Ю.Н. Моделирование процесса отклонения протонов и π-мезонов кристаллом кремния // Письма в Журнал Технической Физики. 2012. Т.38. Вып.12. C.87-94. 12. Боресков А.В., Харламов А.А. Основы работы с технологией CUDA. М.: ДМК Пресс, 2010.

232 с. 19. Marsaglia G., Bray T.A. A convenient method for generating normal variables // SIAM Review. 1964. Vol. 6. No. 3. P. 260–264. Numerical modeling of deflection of fast charged particles bent crystal with using parallel technologies

Yu.N. Shtanov1, V.P. Koshcheev2, D.A. Morgun3, T.A Panina3, A.S. Fokin3 In this paper the use of technologies of parallel CUDA programming and OpenCL in numerical modeling of deflection fast charged particles bent crystal. An algorithm built for calculating the continuous potential and atomic diffusion coefficient using the graphics accelerator.

Mazzolari A. et al. Steering of a sub-GeV electron beam through planar channeling enhanced by Rechanneling // Physical Review Letters. 2014. Vol. 112. P. 135503.

Wienands U. et al. Observation of deflection of a beam of multi-GeV electrons by a thin crystal // Physical Review Letters. 2015. Vol. 114. P. 074801.

Wistisen T.N. et al. A study of Channeling, Volume Reflection and Volume Capture of 3.35 -14.0 GeV Electrons in a bent Silicon Crystal // SLAC-PUB-16435. 2015. URL: http://slac.stanford.edu/pubs/slacpubs/16250/slac-pub-16435.pdf (accessed: 16.01.2016). Shtanov Yu.N. TROPICS the software package. URL: http://wwwinfo.jinr.ru/programs/jinrlib/tropics/index.html (accessed: 30.11.2015).

Koshcheev V.P., Morgun D.A., Shtanov Yu.N. Modeling the deflection of protons and π–mesons in a bent silicon crystal // Technical Physics Letters. 2012. Vol. 38. No. 6. P. 593–596. Koshcheev V.P., Morgun D.A., Shtanov Yu.N. Modeling the deflection of relativistic particles in axial and planar channels of a silicon crystal // Technical Physics Letters. 2013. Vol. 39. No. 10. P. 924–927.

Kitagawa M., Ohtsuki Y.H. Modified dechanneling theory and diffusion coefficients // Physical Review В. 1973. Vol.8. No.7. P. 3117-3123.

Ol'khovskiy I.I., Sadykov N.M. Klassicheskie kineticheskie uravneniya dlya orienta-tsionnykh effektov s uchetom dvukhchastichnoy korrelyatsionnoy funktsii kristalla [Classic kinetic equations for directional effects, taking into account the two-particle correlation function of the crystal] // Trudy 10 Vsesoyuznogo soveshchaniya po fizike vzaimodeystviya zaryazhennykh chastits s kristallami, ch.1. [Proceedings of 10th All-Union Meeting on the Physics of charged particlescrystals interaction] Moscow: Publishing of Moscow State University. 1980. P. 173-179.

References

Dan. Vid. Selsk. 1965 . 34 . No. 14 . 64 P.

10. Bazylev V.A. , Zhevago N.K. Izluchenie bystrykh chastits v veshchestve i vo vneshnikh polyakh [The emission of fast particles in matter and in external fields] . Moscow: Publishing of Nauka, 1987 . 269 P.

11. Koshcheev V.P. Correlation function and nuclear coefficient of diffusion of channelized particles // Russian Physics Journal . 1997 . Vol. 40 . No. 8. P. 736 - 740 .

12. Boreskov A.V. , Kharlamov A.A. Osnovy raboty s tekhnologiey CUDA [Principles of work with CUDA technology] . Moscow: Publishing of DMK Press, 2010 . 232 P.

13. Koshcheev V.P. , Morgun D.A. , Panina T.A. Stokhasticheskaya dinamika effekta kanalirovaniya v kristallakh i nanotrubkakh [Stochastic dynamics of channeling in crystals and nanotubes] . Khanty-Mansiysk: Publishing of Poligrafist . 2008 . 100 P.

14. Shtanov Yu .N., Morgun D.A. , Fokin A.S. Parallel'noe vychislenie potentsialov dlya osevogo kanalirovaniya na videokartakh [Parallel computing potentials for axial channeling on the cards] // Svidetel'stvo o gosudarstvennoy registratsii programmy dlya EVM № 2014618657 [Certificate of state registration of the computer № 2014618657] . Moscow: Publishing of Rospatent, 2014 .

15. Pryjmak A.V. Optimizatsiya vychisleniy na CUDA pri modelirovanii neustoychivosti lengmyurovskikh voln v plazme [Optimization computation on CUDA in the simulation instability of Langmuir waves in plasma] // Visnyk NTUU «KPI». Informatyka, upravlinnja ta obchysljuval'na tehnika [Bulletin of the National Technical University of Ukraine "Kiev Polytechnic Institute": Informatics, Management and Computer Science]. 2013 . No. 58. P. 125 - 130 .

16. Koshcheev V.P. , Shtanov Yu.N. , Morgun D.A. , Panina T.A. Modelirovanie koeffitsienta diffuzii kanalirovannykh chastits [Modelling of the diffusion coefficient of channeled particles] // Mezhdunarodnyy zhurnal «Vestnik kibernetiki» [ International journal “Proceedings of cybernetics”] . 2015 . № 4 ( 20 ). P. 45 - 51 .

17. Gemmel D.S. Channeling and related effects in the motion of charged particles through crystals // Reviews of Modern Physics. 1974 . Vol. 46 , No. 1. P. 129 - 235 .

18. Doyle P.A , Turner P.S. Relativistic Hartree-Fock X -ray and electron scattering factors // Acta Crystallographica. A. 1968 . Vol. 24 . P. 390 - 397 .

19. Marsaglia G. , Bray T.A. A convenient method for generating normal variables / / SIAM Review. 1964 . Vol. 6 . No. 3. P. 260 - 264 .