Новые горизонты суперкомпьютерного моделирования

Разработан параллельный программный код для решения трехмерного нестационарного уравнения Шредингера. Программный код базируется на методе конечных разностей и использует явную численную схему. Простота используемой численной схемы обеспечивает эффективное распараллеливание и высокую производительность программного кода при работе на графических вычислителях. Например, расчет 106 шагов по времени на сетке 1000·1000·1000 (109 точек) занимает всего 16 часов на 16 вычислителях Tesla M2090 суперкомпьютера Ломоносов. Сравнение с другими программами для решения подобных задач, показывает, что производительность разработанного программного кода в 3 раза превосходит существующие аналоги при решении задач одинаковой сложности и эквивалентной стоимости вычислительных ресурсов. (I. K. Gainullin, M. A. Sonkin. Computer physics communications, 188 (2015) 68-75)

1 Зачем это нужно

В данной работе представлен и апробирован параллельный программный код для решения трехмерного (одноэлектронного) нестационарного уравнения Шредингера. Численные схемы и подходы к решению нестационарного уравнения Шредингера достаточно хорошо проработаны. Но, прямое решение трехмерной задачи весьма трудоемко с точки зрения численной сложности (для трехмерного случая). Поэтому для проведения серийных вычислений используют альтернативные (приближенные) способы, например разложение по базисным функциям. В тоже время за последние годы инфраструктура компьютерных вычислений сделала существенный шаг вперед. Технологии суперкомпьютеров позволяют выполнять расчеты в параллельном режиме, что делает трудоемкие задачи более масштабируемыми. Одной из перспективных технологий является использование графических вычислителей. Это требует специальной модификации программного кода, но увеличивает производительность вычислений в 30-100 раз в расчете на 1 процессор. Подобный прирост производительность позволяет моделировать не только «одиночные» трехмерные задачи, но и производить серии расчетов для обобщения и анализа данных.

Нестационарное уравнение Шредингера является основой для решения задачи зарядового (электронного) обмена между атомными частицами и наносистемами. Зарядовый обмен играет важную роль в анализе структуры поверхности твердых тел и методах ее диагностики. Ведь именно анализ заряженных, а не нейтральных частиц лежит в основе большинства методов диагностики поверхности с помощью ионных пучков. На сегодняшний день существуют адекватные методики расчета интегральных параметров зарядового обмена для случая массивных образцов. В тоже время, прямое «первопринципное» моделирование зарядового обмена атомной частицы с атомами поверхности долгое время оставалось нерешенной задачей силу ее численной сложности. Аналитическими методами задача не решается, поэтому необходимо использовать компьютерное моделирование.

2 Реализация параллельных вычислений

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

Для реализации параллельных расчетов на графических вычислителях использовалась технология MPI (Message Parsing Interface). Реализация программного кода в технологии MPI гораздо более трудоемко по сравнению с технологией OpenMP, но позволяет достичь линейного роста производительности при увеличении количества вычислителей.

Трехмерная расчетная сетка была разбита по оси X на срезы с перекрытием (halo) равным двум плоскостям по оси X. На каждом шаге по времени каждый вычислитель проводил расчеты в своей области расчетной сетки, а затем производилась синхронизация данных (передача значений волновой функции в перекрывающихся областях).

Использование графических вычислителей (Graphical Processing Units — GPU) увеличивает выигрыш от распараллеливания программы, но требует существенных дополнительных модификаций программного кода для достижения наилучшего результата. Для расчета непосредственно на графических вычислителях использовался язык программирования Cuda + C. Для оптимизации производительности был применен весь спектр рекомендованных приемов, включая использование общей памяти (shared memory), выравнивание расчетной сетки по размеру расчетного блока для одного мультипроцессора, одновременная асинхронная передача граничных данных и расчет для внутренней области.

3 Производительность расчетов

Производительность программы (ГФлопс) измерялась как количество полезных операций (т. е. операций необходимых для реализации численной схемы, но не вычисления служебных переменных) в единицу времени.

Производительность программы составила 60 ГФлопс на один вычислитель Tesla M 2090 (стоимость порядка 2500 USD) и 110 ГФлопс на один вычислитель Tesla K20 (стоимость порядка 3500 USD). Данные показатели производительности не следует сравнивать с цифрами пиковой производительности, приведенными в маркетинговых материалах компании NVidia (2-4 ТФлопс), т. к. последние не учитывают затрат на передачу данных и вычисление служебных переменных. Также программа демонстрирует линейную масштабируемость по количеству используемых вычислителей, т. е. наблюдается линейный рост производительности при увеличении количества вычислителей. Для примера, расчет 1 миллиона шагов по времени на сетке из 1 миллиарда точек занимает около 16 часов при использовании 16 вычислителей Tesla M2090.

Для сравнения, производительность аналогичного программного кода на «традиционном» процессоре (CPU) Intel Xeon E5-2670 (стоимость порядка 1500 USD) составляет около 10 ГФлопс. Таким образом, технология расчетов на графических вычислителях дает экономию почти на порядок.

Для сравнения с другими аналогичными программами (описанными в научной литературе) была предложена и обоснована величина «нормированное время расчета» (Normalized Calculation Time — NCT), которая учитывает численную сложность задачи и стоимость вычислительных ресурсов. Для примера, значение NCT = 1 сек. означает, что задачу один шаг по времени для системы из одного миллиарда точек программа просчитывает за 1 секунду на вычислительных ресурсах стоимостью 1000 USD. Сопоставление с опубликованными материалами, показало, что разработанный программный код (GPU TDSE Solver) по производительности расчетов в 3-6 раз превосходит имеющиеся аналоги.

Версия программного кода	Используемый компьютер; количество и стоимость задействованных процессоров	Нормированное время расчета, сек.
TDSE GPU Solver	Lomonosov MSU, 16 Tesla M2090 * 1600 USD	1.51
TDSE GPU Solver	GSRV MSU, 1 Tesla K20m * 2700 USD	1.47
TDSE Solver (версия для CPU)	GSRV MSU, 1 Xeon E5-2670 8C * 1500 USD	9.60
Y.-M. Lee, J.-S. Wu, T.-F. Jiang, Y.-S. Chen, Phys. Rev. A 77 (2008) 013414	Taiwan Center for HPC, 128 Itanium-2, 1C, 1.5 GHz * 900 USD	5.47
Y.-M. Lee, J.-S. Wu, T.-F. Jiang, Y.-S. Chen, Phys. Rev. A 77 (2008) 013414	Taiwan Center for HPC, 32 Itanium-2, 1C, 1.5 GHz * 900 USD	9.26
S. X. Hu, L. A. Collins, B. I. Schneider, Phys. Rev. A 80 (2009) 023426	Coyote supercomputer 480 Opteron 2.6 GHz, Infiniband * 1000 USD	5.23
B. I. Schneider, L. A. Collins, S. X. Hu, Phys. Rev. E 73 (2006) 036708; 1D MPI decomposition	LANL QSC, 128 EV68 CB, 1.25GHz * 650 USD	7.09
B. I. Schneider, L. A. Collins, S. X. Hu, Phys. Rev. E 73 (2006) 036708; 2D MPI decomposition	LANL Flash Supercomputer 256 AMD, 2.0 GHz * 250 USD	4.50

4 Физические результаты

В качестве иллюстрации применения разработанного программного кода приведем задачу туннелирования электрона с отрицательного иона водорода в тонкую металлическую островковую пленку (на диэлектрической подложке). Геометрия системы представлена на рисунке снизу. Для простоты рассматриваем цилиндрически-симметричную модельную задачу в приближении закрепленной частицы (ион находится на оси симметрии пленки на заданном расстоянии от ее поверхности).

Описанная модельная задача иллюстрирует электронный обмен в наносистемах. В отличие от макроскопических систем, в наносистемах электронная плотность носит дискретный характер. Также в наносистемах наблюдаются квантово-размерные эффекты, заключающиеся в немонотонной зависимости ключевых параметров от размера системы.

Следующий рисунок иллюстрирует эволюцию электронной плотности для вышеописанной системы. На рисунке показаны изоповерхности электронной плотности в последовательные моменты времени (100, 500 и 2250 атомных единиц времени). Сферическая составляющая соответствует локализации электрона на отрицательном ионе, цилиндрическая локализации внутри тонкой островковой пленки.

Мы видим, что на начальном этапе (0-200 ат. ед.) электрон туннелирует в пленку перпендикулярно поверхности. При этом образуются максимумы электронной плотности (дискретность распределения электрона) только по нормальной координате (т. е. электрон чувствует конечный размер островковой пленки только вдоль нормали к поверхности). На втором этапе (200-2000 ат. ед.) электрон внутри пленки распространяется параллельно поверхности и формируется дискретная структура (кольца), содержащая максимумы как по нормальной, так и по радиальной координате. Наконец, на третьем этапе (свыше 2000 ат. ед.) проявляются гармоники и по азимутальной координате. Данный результат трехмерных расчетов представляется весьма интересным, т. к. система, изначально находившаяся в цилиндрически-симметричном состоянии, эволюционирует в полностью трехмерное состояние.

Доцент кафедры физической электроники Гайнуллин И. К.