CC BY
98
1. МОДЕЛИРОВАНИЕ НАНОСИСТЕМ И НАНОЭЛЕКТРОНИКА
1.1. ИССЛЕДОВАНИЕ ЭФФЕКТИВНОСТИ ИСПОЛЬЗОВАНИЯ ГРАФИЧЕСКИХ ПРОЦЕССОРОВ ДЛЯ РЕШЕНИЯ ВЫЧИСЛИТЕЛЬНЫХ ЗАДАЧ НАНОТЕХНОЛОГИЙ
Попова Нина Николаевна, доц., доц., канд. физ.-мат. наук, факультет вычислительной математики и кибернетики, Московский Государственный Университет им. М.В. Ломоносова, popova@cs.msu.su
Никишин Николай Глебович, аспирант, факультет вычислительной математики и кибернетики, Московский Государственный Университет им. М.В. Ломоносова, mail@nikishinng.ru
Аннотация: Вычислительные нанотехнологии неотъемлемо связаны с использованием современных высокопроизводительных систем. С быстрым развитием аппаратного и программного обеспечения графических процессоров (GPU) их использование набирает популярность для задач, требовательных к вычислительной мощности. К числу таких задач относятся исследования в области наноматериалов, изучение и создание которых проводятся с использованием масс-спектрометров. Работа посвящена моделированию поведения ионов в ловушках масс-спектрометров на основе ионного циклотронного резонанса и преобразования Фурье. Решение задачи проводится на основе модели частиц в ячейке, которая применяется для прямого моделирования поведения ионов. Вычисления на графическом процессоре проводятся с использованием библиотек для GPU cuFFT и CULA.
В работе показано, что использование GPU-ориентированных библиотек существенно упрощает разработку параллельных алгоритмов для графических процессоров и позволяет достичь хорошую производительность параллельных приложений. Решение задачи выполнено на суперкомпьютере «Ломоносов», установленном в МГУ имени М.В. Ломоносова. В работе показано, что различные стратегии распределения параллельных процессов по узлам многопроцессорной вычислительной системы могут существенно влиять на производительность всей программы из-за одновременного доступа нескольких процессов к графическим ускорителям.
Полученные в работе результаты могут быть полезными для моделирования больших молекулярных структур, решения задач в области вычислительных нанотехнологий на современных высокопроизводительных параллельных вычислительных системах с использованием графических ускорителей
Ключевые слова: высокопроизводительные вычисления, неоднородные вычислительные системы, нанотехнологии, масс-спектроскопия, CUDA, Быстрое преобразование Фурье.
1.1. PERFORMANCE STUDY OF GRAPHIC PROCESSORS USAGE IN COMPUTATIONAL NANOTECHNOLOGY PROBLEMS
Popova Nina N., PhD., associate professor, PhD in Physics and Mathematics, Lomonosov Moscow State University, Department of Computational Mathematics and Cybernetics. popova@cs.msu.su
Nikishin Nikolai G., postgraduate, Computational Mathematics and Cybernetics Department, Lomonosov Moscow State University, mail@nikishinng.ru
Abstract: Computational nanotechnology is strongly connected with modern high-performance computing. With fast evolution of graphical processing units (GPU), general-purpose computing (GPGPU) became a popular choice for computationally demanding tasks. Tasks in nanomaterial researches with mass-spectrometers usage for analysis and material creation are good examples of such kind of tasks. Paper is devoted to ions behavior modeling in traps of mass-spectrometers based on Fourier transform. We use heterogeneous computational systems with GPU inside for calculation. Particle-in-cell model is used for direct modeling of ions behavior. We also use two GPU libraries: CULA and cuFFT.
Paper shows, that some GPU-oriented libraries could significantly ease the development of parallel algorithms for GPU and allow to get good performance of parallel applications. Calculations were performed on several systems including "Lomonosov" supercomputer in Moscow State University. Paper shows, that different strategies of mapping parallel processes to nodes could significantly effects on performance because of parallel access to an every single GPU from multiple processes.
Results obtained in this work could be useful for big molecular structure modeling, for solving computational nanotechnology problems on modern high-performance parallel computational systems with GPU
Index terms: HPC, heterogeneous computing systems, nanotechnology, mass-spectrometry, CUDA, Fast Fourier transform
1. Введение.
Вычислительные нанотехнологии неотъемлемо связаны с использованием современных высокопроизводительных систем. Использование ускорителей, интегрируемых в состав процессоров, является основной тенденцией развития таких систем. С быстрым развитием аппаратного и программного обеспечения графических процессоров (GPU) их использование набирает популярность для задач, требовательных к вычислительной мощности [1]. Системы, основанные на использовании GPU, привлекают особое внимание исследователей [2]. Это ставит вопрос об эффективности применения таких систем для моделирования больших молекулярных структур, в особенности, для разработок в сфере нанотехнологий [3]. Исследования, связанные с изучением и созданием новых наноматериалов, активно проводятся с использованием масс-спектрометров. В работе исследуется эффективность высокопроизводительных вычислений при проведении программного моделирования поведения ионов в масс-спектрометрах, построенных на основе ионного циклотронного резонанса и преобразования Фурье. В основу решения рассматриваемой задачи принята широко распространенная модель частиц в ячейках [4]. Для прямого моделирования изучаемых устройств, а также для исследования физических эффектов, влияющих на точность работы этих устройств, предложены методы и реализован ряд программ [5], [6]. В них расчет парного кулонов-ского взаимодействия частиц осуществляется через вычисление на трехмерной разностной сети коллективного самосогласованного поля на каждом шаге по времени на основе решения задачи для уравнения Пуассона с использованием быстрого преобразования Фурье [5]. Данная работа продолжает исследования в этом направлении. В работе предлагается параллельный алгоритм, основанный на использовании графических ускорителей. Основное внимание в работе уделяется тому, насколько GPU позволяют ускорить вычисления в модели частиц в ячейке. Для решения краевой задачи для уравнения Пуассона предлагается использование библиотека cuFFT для быстрого преобразования Фурье. В работе исследуется эффективность нахождения решения задачи Пуассона для различных параметров и проводится анализ асимптотической сложности решения. Расчет внешних удерживающих полей на разностной сетке от системы электродов производится с помощью параллельных алгоритмов решения пограничных интегральных уравнений [6]. Задача сводится к решению алгебраической системы уравнений с плотной матри-
цей. Решение системы проводится с использованием библиотеки CULA, реализующей основные операции линейной алгебры на основе технологии NVIDIA CUDA. Предложенный параллельный алгоритм решения трехмерной задачи реализуется на основе технологии параллельного программирования MPI. В работе проводится сравнение различных стратегий распараллеливания вычислений с возможностью использования нескольких параллельных потоков на GPU. Практическая реализация задачи выполнена на суперкомпьютере «Ломоносов» НИВЦ МГУ и на ряде вычислительных систем с различными графическими ускорителями.
Полученные в работе результаты и приведенные на их основе рекомендации могут быть полезны при построении численных кодов для решения различных задач, связанных с применением вычислительных нанотехнологий.
2. Задача моделирования эволюции частиц в ячейке.
Рассмотрим модель частиц в ячейке при моделировании поведения многих ионов в ловушках различных установок, в частности, масс-спектрометров [5]. Чтобы определить траектории частиц, необходимо решить N уравнений Ньютона для каждой частицы под действием силы Лоренца.
m dV- = q.E""(r., t) + q [v . dt
NP q,q,(r - r;)
"Z ,,. . ,/) + E""(r,) + F„n(v,) + E""(r)cosY(t),
4nea I r - r I3
где электрическое поле Eопределяется по вычисленным потенциалам: utrap(r,)— потенциалу удерживающего поля ловушки, urf (r,,t) - потенциалу поля возбуждения (заданного по времени) и uc(rht) - куло-новскому потенциалу парного взаимодействия частиц, который необходимо рассчитывать на каждом шаге по времени.
,. wall г
Через u обозначим значение потенциала на стенках электродов. Тогда:
E = -Vu, u = ulap(r) + uIf (r, t) + u c(r, t) + uwan (1)
Наиболее трудозатратными в методе частица-сетка являются процедуры расчета сил кулоновского взаимодействия частиц и вычисления поля ловушки, удерживающего ионы, которое создаётся системой электродов произвольной формы. При этом в задаче моделирования масс-спектрометра приходится моделировать эволюцию в течение нескольких миллионов шагов по времени порядка миллиона ионов. Использование графических ускорителей для выполнения данных процедур представляет основной интерес.
1
3. Решение краевой задачи для уравнения Пуассона с целью определения кулоновского взаимодействия ионов.
Во избежание попарного вычисления взаимодействий между частицами, вводится пространственная трехмерная разностная сетка (хк, у, zm). Потенциал кулоновского поля определяется плотностью заряда р(г,Д), который вычисляется на сетке при помощи интерполяции зарядов частиц qi внутри ячейки, в узлы данной разностной сети:
1 ' p
P(k'l' m) = V0l "^qWk-'m (' у '
(2)
где Vol есть объем ячейки, а WkJ,m - веса интерполяции, величина которых зависит от положения частицы относительно узла (k, I, m). Тогда коллективный
c
потенциал кулоновского поля u находится из решения первой краевой задачи для уравнения Пуассона на пространственной сетке
Au --p(r, t), u 4ouni= u ™\r, t),
(3)
j1=0 j2=0 j3=
f a^-j ] sm(j ] sinj
l N1 ) 1 N2 ) l N3
Выполнив синус-преобразование Фурье левой и правой частей, перейдём к уравнению в пространстве частот:
Ш = К [А ]>
где Хь(к,!,т) - собственные числа разностной задачи Дирихле. Они равны [7]
4
Лк (k, l, m) = ^т-fsi
a=1 ha
kanha
2la
ki = k, k2 = l, k3 = m.
Поделив обе части уравнения на собственные числа, получим образ решения задачи:
й с = FSLPhI
ик = л •
\
Для получения потенциала на сетке остаётся только сделать обратное преобразование.
uh = F-
FP
А
и по найденному потенциалу вычисляется сила, действующая на все ионы внутри ячейки:
F(x,, у,,г,) = Y^i, у,,г,)£grM. (4)
grid
В результате количество действий для учета парного взаимодействия меняется с O(Np• Np) на O(Np- M), где M - число ячеек сети.
Чтобы избежать флуктуаций силы при переходе через границы ячеек, число рассматриваемых частиц должно быть велико по сравнению с размерностью пространственной сетки. Для этого, как правило, берут порядка ста частиц на одну ячейку.
Переформулируем задачу (3) для трехмерной разностной сетки размера N1x N1x N2x N3, равномерно покрывающей прямоугольную область [0, LJ х[0, L2] х[0, L3]. Тогда координаты узлов сетки задаются следующими выражениями: xt = kh{, i = 1,2,3, (5)
где hi = Li / Ni— шаги разностной сетки по каждому из направлений. Постановка задачи при этом примет вид:
Аhuh = ph> uh Lund=0,
где uch и ph - сеточные функции, заданные значениями функций uc и р в узлах сетки, а Дл - разностная аппроксимация трехмерного оператора Лапласа.
Для удовлетворения нулевым граничным условиям в качестве функций базиса разложения следует использовать ортогональную систему произведений синусов, обращающихся в нуль на границе области. Дискретное синус-преобразование Фурье задаётся следующей формулой:
N1-1N2-1N3-1~
Fs[f ](k,l,m)= y y Y
4. Программная реализация решения краевой задачи для уравнения Пуассона с использованием библиотеки cuFFT.
Программная реализация коллективного потенциала кулоновского поля действий с использованием GPU представлена в виде модуля, состоящего из трёх процедур.
Первая процедура включает в себя инициализацию работы с графическим ускорителем, выделение необходимой памяти на GPU, вычисление собственных значений для дальнейшего использования в программе.
Вторая процедура предназначена для решения краевой задачи для уравнения Пуассона в памяти GPU с использованием дискретного преобразования Фурье.
Третья процедура освобождает выделенную память. Чтобы избежать повторных вычислений, инициализация и завершение работы с ускорителем вынесены за пределы временного цикла.
Взаимодействие с основной программой осуществляется через интерфейсный модуль, обеспечивающий возможность вызова процедур, написанных на языке NVIDIA CUDA, из основной программы на языке Fortran. Для гарантии соответствия типов при его описании использовались именованные константы, предоставляемые модулем ISO_C_BINDING, который входит в стандарт языка Fortran 2003. В главной процедуре, реализующей решения уравнения, для выполнения преобразования Фурье используются функции библиотеки cuFFT версии 5.5 из состава программного пакета CUDA Toolkit. Библиотека cuFFT не содержит процедур для синус-преобразований. В ней содержатся только процедуры дискретного преобразования Фурье (ДПФ). Для получения синус-преобразования на интересующей нас области расширим нашу сеточную функцию в два раза по каждому из направлений, продолжив её нечётно относи-
2
0
0
тельно бывших границ. Это преобразование задаётся следующими равенствами:
/(¡1. ¡2. ¡з) = /(¡1. ¡2. ¡з).
/ (Ж. - ¡„ ¡2. ¡3) = / ( 11.2М1 - ¡2. ¡3) = / ( ¡,. ¡2.2^3 - ¡з) = -/ ( ¡,. ¡2. ¡з).
/(Ж. - ¡, .Ж2 - ¡2. ¡з) = /(2#, - ¡,. ¡2.2^3 - ¡з) = /(¡,.Ж2 - ¡2.2^3 - ¡з) = /(¡,. ¡2. ¡3).
/(Ж. - ¡, .Ж2 - ¡2.2^3 - ¡3) = -/( ¡,. ¡2. ¡з). ¡. =1.....N. -1.
/(¡,. ¡2. ¡з) =0. ¡, = 0.N.. 2М,.! = 1.2.3.
В этих предположениях, для любого i = 1,2,3 верно, что
2 N. -1
2п;м/(2 N.)
X f (k, jt ,... jN¡ >= {j = 2 ^ - jt }=
h =Nt
N. N.-1
2n<2 N. - h )m/<2 N.) ^ „ , . -2n¡m/<2N¡.)
= Xf (...,2N¡ - j;,...)e2'
j'¡ =i
' = - X f (•••, h,K)e
j¡=o
Используя эту формулу, покажем, что ДПФ от расширения нашей функции даст нужный результат на интересующей нас части сетки.
2 N. -12^ -12 N. -1Г
2<*1— 12N 2-12N3—1 р
F[f ](k,l,m) = X X X [f (Л'h> Л)'
2j/(2N1)e2n11ll(2N2)e2nj3ml(2N3) ] =
h =0 h =0 h =0
2N -12N -1
= X Xle^2N1)ej2N2)X[f (11,12,1з)(е ^ - - e
j1=0 j2=0 [ 2 N. -12 N--1 Í
j3 =
2j/(2 N1)e2nj1U(2 N2)
j1=0 j2=0 [ j3=o
= X X^e2nj1k'(2N1>e2n1l/<2N2)2¡ X
f (j1, j2, 1з) sin|
= /(jJi.J]sinff^]sin j] = -8«F,[/](*,l,m).
После проведения прямого и обратного преобразований необходимо провести нормировку на общее количество элементов сетки, так как процедуры cuFFT не делают это автоматически.
Решение находится по следующей цепочке действий. Значения сеточной функции ph пересылаются из основной оперативной памяти в память GPU. Затем эти значения переупорядочиваются из массива с индексацией «по столбцам», принятого в языке Fortran, в массив с индексацией «по строкам», естественного для языков семейства C. Сеточная функция ph нечетно продолжается в два раза по каждому из направлений. Далее производится дискретное преобразование Фурье продолженной сеточной функции. Все мнимые части результата (только они не равны нулю) делятся на собственные значения разностной задачи. Затем производится обратное дискретное преобразование Фурье образа решения. После обратного переупорядочивания «по столбцам» массив значений сеточной функции uc пересылается обратно в ОЗУ.
В предположении, что исходная трехмерная сетка содержит Ng узлов по каждому из направлений, количество операций с плавающей точкой для осуществления прямого и обратного преобразования Фурье, а также деления на собственные числа оценивается по следующим формулам
Nop = 2N FFT + (2N g)3, N FFT = 5(2Ng)3 log2 ((2 Ng)3), (6)
где Nop - общее количество операций, а NFFT - количество операций, необходимое для реализации алгоритма быстрого преобразования Фурье над массивом
комплексных чисел. Общая асимптотическая сложность составляет O(Ng3 log(Ng)).
5. Оценки эффективности однопроцессорной программной реализации с использованием различных графических процессоров.
Для оценки технических характеристик графических процессоров рассчитываются следующие величины:
- максимально возможная пропускная способность памяти,
- пиковая производительность на операциях с одинарной точностью,
- пиковая производительность на операциях с двойной точностью.
Данные величины могут быть получены с помощью вызовов CUDA runtime API и зависят только от характеристик графических карт.
_ Memory Clock rate (MHz)-106 • MemoryBusWidth (Bits)/8
2n3m/(2 N3) -2n3m/(2N3) I
BWthe„(GB/S) = -
10'
(7)
Perf^ (GFLOPS) = 2 • Total CUDA Cores • GPU Clock rate (GHz), Perf1dpak (GFLOPS) = 2 • Total CUDA DPU • GPU Clock rate (GHz).
Время решения задачи для уравнения Пуассона замерялось на каждом шаге по времени и затем вычислялось среднее время по 100 шагам. Решение считалось удовлетворительным, если абсолютная ошибка по сравнению с точным аналитическим решением составляла величину ~10-5.
Производительность программ с использованием GPU оценивается с помощью двух метрик [8]: BWeff -достигнутой пропускной способности памяти и Perfeff - достигнутой производительности. Значение BWeff вычисляется по формуле:
BWeff (GB/s) = ^b+Wb, (8)
t •lO9
где RB - количество прочитанных из памяти GPU (в байтах), а WB .количество записанных в память данных. В нашем случае RB = WB = N3, то есть нам требуется чтение сеточной функции ph и запись ответа uh.
N
Perfeff (GFLOPS) = - op
t-10
9 '
(9)
Мы не будем рассматривать сетки с N < 32 с целью достаточной загрузки мультипроцессоров GPU. Время решения задачи в зависимости от числа точек (одинаковое число точек по каждому из направлений) на устройстве GeForce GTX TITAN для одинарной и двойной точности показано в таблице 1.
Таблица 1
Время решения задачи Дирихле на устройстве GeForce GTX TITAN. Усреднение по 100 шагам по времени
Ng Одинарная точность Двойная точность
32 0.001152 0.001475
64 0.004569 0.007384
128 0.026510 0.053721
256 0.259035 0.520362
Аппроксимация в смысле наименьших квадратов кривыми Сх3^2(х), соответствующими асимптотической оценке сложности, представлена на рис. 1. Накладные расходы начинают играть все меньшую
роль с увеличением размера задачи, уступая вычислительной составляющей. Время вычислений для нахождения решения составляет от 50% до 70% от общего времени выполнения процедуры и растет с увеличением размера задачи. Отношение времени, затраченного на вычисление решения, к общему времени выполнения процедуры, включающему пересылку и подготовку данных, для различных размеров задачи показано в таблице 2.
50 100 150 200
Число точек по каждому из направлений
Рисунок 1. Время решения первой краевой задачи для уравнения Пуассона для одинарной (+) и двойной (х) точности. Пунктиром показаны кривые вида Cx3log2(x) наилучшим образом аппроксимирующие числовые данные в смысле наименьших квадратов.
Таблица 2
Доля вычислений от общего времени решения задачи одинарной точности на GeForce GTX TITAN
Ng Решение (с) Общее время (с) Доля решения в % от общего
64 0.002461 0.004569 53.862990
128 0.016318 0.026510 61.554131
256 0.188552 0.259035 72.790163
Сравнение эффективности различных ускорителей было проведено для характерного для задачи числа точек по каждому из направлений Ng = 128. Соответствующие результаты приведены в таблице 3.
Таблица 3
Сравнение времени решения задачи Дирихле, пиковой и достигнутой пропускной способности памяти, пиковой и достигнутой производительности на GeForce GT 640, Tesla C2075 и GeForce GTX TITAN. Усреднение по 100 итерациям для решения задачи с одинарной точностью Ng = 128.
GPU Время, с BWtheor, ГБ/с BWf, ГБ/с Perf sP , Pe peak ТФЛОПС PerfeSP, ТФЛОПС
GeForce GT 640 0.090254 40.08 0.05 0.81 0.0448
Tesla C2075 0.037220 150.34 0.11 1.03 0.1086
GeForce GTX TITAN 0.026510 208.00 0.16 3.54 0.1525
Из таблицы видно сильное превосходство ускорителей GeForce GTX TITAN и Tesla C2075 над GeForce GT 640. Это объясняется большим количеством потоковых мультипроцессоров (14 против 2) и повышенной
пропускной способностью памяти. GeForce GTX TITAN оказывается быстрее Tesla C2075 за счет большего количества ядер на мультипроцессоре.
В исследуемой задаче количество точек сетки выбиралось разумным для серийных расчетов по модели частиц в ячейках и составило порядка 1,7 • 107 точек. Такое количество данных допускается широтой пропускания памяти, но после чтения исходных данных мы каждый раз вынуждены выполнять транспонирование и расширение массива в 8 раз, что вызывает заметные задержки, так как память GPU обладает высокой латентностью. В зависимости от ускорителя достигнутая производительность составляет от 44 до 150 ГФЛОПС, что соответствует 4-10% от пиковой производительности. Описанные выше особенности алгоритма не позволяют в полной мере задействовать все ресурсы GPU, однако решение задачи осуществляется за хорошее время и может использоваться в приложениях.
6. Вычисление потенциалов удерживающего поля от электродов произвольной формы.
Для вычисления потенциала поля от незамкнутых электродов произвольной формы требуется решить уравнение Лапласа
АУ = 0 (10)
с граничными условиями первого рода на электродах, окружающих ловушку:
У iEk = Uk, (11)
где k = 1,2,3,...,K - номер поверхности электрода, и Uk - заданный постоянный потенциал на каждой поверхности.
Мы ищем решение как сумму потенциалов двойного слоя Wk, создаваемого каждой поверхностью:
V (M ) = JWt (M).
k
Неизвестные дипольные моменты vk(P) на каждой поверхности дают нам решение задачи - потенциал V(M). После определения vk(P) мы подставляем дипольные моменты в выражение для Wk:
^^ (M ) = "Я
дп„
(
\
V Rmp J
vk (P)dap.
Здесь Р есть точка на поверхности !к, М - точка наблюдения (точка трехмерной разностной сети, покрывающей объем ловушки),
Кир ^ представляет
собой расстояние между точкой Р и точкой М, пр -нормаль к поверхности в точке Р, а vk(P) является неизвестной поверхностной плотностью дипольного момента, которая различна на каждой поверхности.
Дипольные моменты ук(Р) получаются из граничного условия для потенциала двойного слоя, когда мы помещаем точку М последовательно на все поверхно-(к)
сти, М = Р0 . Так мы получаем систему уравнений для
1
L
k
дипольных моментов vk(P) на каждой поверхности. ^(Р) определяются граничными условиями, а именно значениями ик на каждой поверхности и потенциалами, наведенными другими поверхностями. В каждом уравнении один из интегралов в сумме при i = к и Р0' = Р имеет сингулярность, что приводит к тому, что в численной схеме такие члены должны иметь специальную аппроксимацию.
Внешняя часть поверхности является отталкивающей, а внутренняя - притягивающей.
fe4)- Я (P)d«p = U. + $ Ц
dn
R_.
V (P )d&p.
В левой части мы оставляем интегрирование только по поверхности к (точка Р), при этом точка наблюде-
п (к)
ния потенциала Р0 тоже находится на этой поверхности. В правой части уравнения находятся приложенный потенциал на поверхности к и сумма потенциалов от всех поверхностей которые индуцируют электрическое поле на данной поверхности. Идея итерационной процедуры состоит в фиксации ди-польных моментов для всех поверхностей, кроме к. Правая часть уравнения известна на текущем итерационном шаге и может быть вычислена. По решениям для всех поверхностей на этом шаге можно вычислить новую правую часть. Система уравнений переписывается в форме
( \ 1
р')- Я dt
R
v. (P)dap = Ф.
со следующей правой частью,
( \
1
Ф„
и + $ ff i
i=1, i^k '
R
Vi (Pi )da'.
(12)
(13)
(n)
есть
число сеточных точек по двум направлениям, и^ решение для п-ой поверхности.
Пусть 5 есть номер итерации. Тогда получаем следующую итерационную процедуру:
2ЯМ;
N w2
'О-"1 J V W(,,0,s+1)>/'!0's+1)
-
fkm + £ I I ^kmij
' kmij ij
Первый член правой части представляет собой заданный поверхностный потенциал, а второй найден на предыдущем шаге и представляет собой влияние
других поверхностей на распределение дипольного момента с номером п0.
Каждое уравнение переписывается в форме алгебраической системы
А У = Ф, (14)
где Y - линейный вектор решения, Ф - линейный вектор известной правой части, и А есть матрица, включающая геометрические характеристики электродов. Для решения можно использовать соответствующие стандартные процедуры решения систем линейных алгебраических уравнений с плотной матрицей из библиотеки LAPACK.
Электрический потенциал, индуцированный на каждой поверхности электрода, может быть вычислен параллельно для всех поверхностей с номерами п = 1,2,...,т. На каждой итерации определяются поля Цпс(п5,пс), индуцированные собственной поверхностью процесса на всех других поверхностях с номерами пс. Редукционной операцией суммирования, реализуемой функцией МР1_А11гес1исе, во всех процессах определяется полное поле, индуцированное всеми электродами, и новая правая часть в уравнениях.
На новой итерации используется модифицированное граничное условие для каждой поверхности. В каждом МР1-процессе на графических процессорах решается алгебраическая система уравнений (14) со своей матрицей и правой частью. Схема параллельного алгоритма показана на рис. 2.
7. Решение системы интегральных уравнений на гибридных вычислительных системах c CPU и GPU.
После написания квадратур для интегралов, полная система для набора поверхностей (12) может быть переписана в форме алгебраических уравнений.
Обозначим интегралы для поверхности n после дискретизации как N1 N2
YYF (п) и (n)
¿_t¿_t kmij ij i=0 j =0
где F учитывает квадратурную формулу, N1, N2 -
Рисунок 2. Схема параллельного алгоритма расчета удерживающего поля ловушки. В каждом параллельном процессе решение алгебраической системы уравнений проводится на GPU.
8. Решение алгебраических систем уравнений на GPU для каждой поверхности с помощью библиотеки CULA.
Решение алгебраической системы для каждой из поверхностей электрода выполняется на GPU с использованием библиотеки CULA v.R17. В рассматриваемой задаче параметры выбраны так, что алгебраическая система уравнений (14) имеет матрицу размера N1x N2, N1 = N2 = 35, 45, 55, 65, 75. Расчеты по решению задачи проводились на суперкомпьютере «Ломоносов». Для решения выбирались процессорные узлы, оборудованные двумя 4-ядерными CPU и двумя GPU Tesla X2070, по характеристикам аналогичными Tesla C2070. Исследовались различные стратегии распределения параллельных процессов по
1
е
е
р.- ,р
p: ,р
Е
P ,P
.
N.. N N
О j=0
n=1, i=0 /=0
узлам. Первая стратегия состояла в распределении 20 MPI-процессов по 1 процессу на вычислительное ядро процессора. Доступ к каждому GPU осуществлялся одновременно несколькими процессами. Вторая стратегия подразумевала, что каждый из 20 MPI-процессов запускается на своем вычислительном узле, при этом в вычислениях задействовалось только одно из четырех ядер процессора.
На рис. 3 показаны схемы обеих стратегий.
Таблица 4.
Сравнение среднего времени решения СЛАУ одним процессом при использовании первого и второго вариантов.
N Вариант а) (с) Вариант б) (с) Отношение а) к б)
35 0.3082 0.0491 6.27
45 0.5082 0.0778 6.53
55 0.5855 0.1487 3.94
65 1.2202 0.3128 3.90
75 1.7293 0.5168 3.35
Рисунок 3. Схемы распределения задачи расчета удерживающего поля ловушки по узлам суперкомпьютера «Ломоносов».
Первая стратегия характеризуется тем, что несколько процессов конкурируют за ресурсы GPU. Каждый из ускорителей единовременно может выполнять лишь один CUDA-контекст, соответствующий одному MPI-процессу. Процедура решения СЛАУ в CULA не монолитна и разбивается на некоторое количество подзадач, что позволяет переключать контексты нескольких процессов по мере необходимости. Накладные расходы на такие переключения сокращаются при увеличении размеров обрабатываемых данных.
При использовании второй стратегии каждый процесс получает доступ к GPU сразу же, без ожидания, и выполняет подзадачи одного процесса до получения окончательного ответа.
Средние времена решения СЛАУ при использовании первого и второго вариантов и их отношение показаны в таблице 4. При увеличении числа параллельно обращающихся к одному GPU процессов разброс времени возвращения из процедуры решения также увеличивается. На рис. 4 показано время выполнения процедуры решения СЛАУ в зависимости от условного номера параллельного процесса при характерном числе точек в матрице N = 65. На рисунке также показано время решения при использовании лишь одной видеокарты GeForce GT 740M.
10 12 Номер процесса
Рисунок 4. Время решения линейной системы уравнений при параллельном выполнении вычислений при использовании схем а) (о) и б) (х). Время вычислений показано в зависимости от условного номера параллельного процесса
при числе точек в матрице системы уравнений N = 65. В каждом параллельном процессе на CPU используется GPU
для решения алгебраической системы уравнений. Для сравнения показано также время решения при использовании лишь одной видеокарты GeForce GT 740M (+). 9. Заключение.
В работе представлены результаты реализации модели частиц в ячейке для моделирования работы масс-спектрометров на основе ионного циклотронного резонанса и преобразования Фурье на гибридных вычислительных системах с использованием GPU.
Использование графических ускорителей для выполнения дискретного преобразования Фурье и решения систем линейных алгебраических уравнений позволяет достигать хорошего времени решения возникающих подзадач. К ним относятся определение кулоновского взаимодействия ионов с помощью решения первой краевой задачи для уравнения Пуассона и параллельное вычисление полей от каждой поверхности электрода через решение алгебраических систем. Учёт особенностей архитектуры вычислительной системы при организации параллельных вычислений с использованием GPU позволяет избежать дополнительных накладных расходов на перегрузку контекстов, связанных с переключением потоков.
Представленные в работе результаты могут быть полезны при разработке параллельных программ для моделирования молекулярных структур в сфере разработок в области нанотехнологий в силу аналогич-
ных подходов к использованию исследованных процедур.
Работа выполнена при поддержке грантов РФФИ №1301-12078 офи-м, №14-07-00654 и №14-07-00628.
This work was supported by the Russian Foundation for Basic Research, projects no. 13-01-12078 ofi_m, no. 14-07-00654 and no. 14-07-00628.
Список литературы:
1. Wilt N. The CUDA Handbook: A Comprehensive Guide to GPU Programming. Addison-Wesle, 2013. P. 528
2. Parallel Computing: From Multicores and GPU's to Pet-ascale / Chapman B., Desprez F., Gerhard R., et al. IOS Press, 2010. Vol. 19. P. 739
3. Попов А.М. Вычислительные нанотехнологии: учебное пособие. М.: КНОРУС, 2014. 312 c.
4. Хокни Р., Иствуд Дж. Численное моделирование методом частиц: пер.с англ. М.: Мир, 1987. 640 с.
5. Nikolaev E. N. Heeren R. M. A. Popov A. M. Realistic modeling of ion cloud motion in a Fourier transform ion cyclotron resonance cell by use of a particle-in-cell approach // Rapid Communication in Mass Spectrometry. 2007. Vol. 21(22). P. 3527-3546.
6. Misharin A., Popov A. Parallel numerical code «parTfield» to simulate ion optical elements for any electrode geometry // Proc. 60th ASMS Conference on Mass Spectrometry and Allied Topics, Vancouver, Canada, May 20-24. 2012.
7. Самарский А. А., Николаев Е.С. Методы решения сеточных уравнений. М.: Наука, 1978. 590 с.
8. Malony A.D., et al. Parallel Performance Measurement of Heterogeneous Parallel Systems with GPUs // International Conference on Parallel Processing. ICPP 2011, Taipei, Taiwan, September 13-16. 2011. P. 176-185.
