ПК обогнал суперкомпьютеры в решении задачи трехчастичного рассеяния

Сложные задачи, для решения которых требуется большая вычислительная мощность, регулярно возникают в разных областях наук. Естественный ответ на это - рост числа и производительности суперкомпьютеров. Однако вычислительное время на них по-прежнему остается дорогим, а конкуренция научных коллективов за него - высокой. Российские физики из НИИЯФ МГУ имени М. В. Ломоносова нашли другой подход к одной из важных ресурсоемких задач квантовой теории рассеяния. Они смогли свести решение уравнений Фаддеева к виду, в котором максимально используются сильные стороны видеокарты обычного персонального компьютера: многопоточное параллельное выполнение одинаковых процедур. Этот результат открывает широкие возможности для решения многих других вычислительных задач.

Суперкомпьютеры и сложные задачи

В настоящее время наблюдается быстрый рост потребностей в высокопроизводительных вычислениях во всех областях естественных наук. Это вызвано как потребностями математического (компьютерного) моделирования сложных многостадийных процессов, характеризующихся наличием большого числа степеней свободы (такие процессы встречаются, например, в гидро- и газодинамике, теории плазмы, астрофизике, ядерной физике, теории поля на решетках и т. д.), так и переходом к точному численному решению многомерных дифференциальных и интегральных уравнений в реалистичной постановке. Как вполне естественный ответ на эти новые запросы сейчас наблюдается стремительный рост числа и установленной производительности суперкомпьютеров по всему миру, причем наиболее быстрый рост пиковой производительности супермашин наблюдается в Китае и США (см. рейтинг TOP500). Например, самый «мощный» суперкомпьютер по состоянию на июль 2016 года - китайский Sunway TaihuLight - имеет производительность порядка 100 петафлопс, что в три раза больше, чем у предыдущего лидера, и в сотни раз больше способностей суперкомпьютеров начала XXI века.

В квантовой теории рассеяния есть тоже класс задач, предъявляющих очень большие требования к вычислительным ресурсам. Речь идет о точном решении (в численном смысле) задач рассеяния с несколькими частицами (так называемые few-body systems) в непрерывном спектре. Напомним читателю, что связанные состояния систем, когда движение частиц ограничено, обладают фиксированными (квантованными) энергиями и тем самым отвечают области дискретного спектра гамильтониана системы. В настоящее время вполне возможно аккуратно рассчитать дискретные спектры систем, состоящих из многих сотен и даже тысяч частиц (таких как электроны и ядра атомов), что, в частности, включает в себя точный расчет структуры атомов и больших молекул (этим занимается квантовая химия). Однако рассеяние частиц соответствует их неограниченному движению, поэтому состояния рассеяния существуют для непрерывного (континуального) набора энергий и таким образом отвечают области непрерывного спектра гамильтониана системы. Также при одной энергии возможны различные варианты (каналы) рассеяния для одной и той же системы. Таким образом, задачи в непрерывном спектре являются качественно другими и значительно более сложными, чем задачи дискретного спектра. Неудивительно поэтому, что решение квантовой задачи рассеяния уже трех частиц вызывает большие трудности, причем как в чисто математическом, так и в вычислительном плане.

Важность точного решения задачи трехчастичного рассеяния (и задач для большего числа частиц) в квантовой механике обусловлена тем, что оно очень чувствительно к внутренней структуре рассматриваемых квантовых систем и позволяет как бы заглянуть внутрь ядер и атомов и найти ключи к ряду важных проблем современной физики (например, лучше понять природу ядерных сил). К тому же многие физические задачи, такие как рассеяние дейтронов на нуклонах и ядрах, рассеяние ряда легких ядер (^6,7Li, ¹³C и т. д.) на других ядрах или рассеяние электронов и позитронов на атомах и молекулах, а также ряд задач в физике твердого тела и т. д., могут быть описаны в рамках квантовой модели трех частиц.

Вспомним также, что математически строгое решение даже ньютоновской задачи трех тел в классической механике, в частности, изучение поведения системы Земля - Луна - Солнце в полном виде не найдено до сих пор, несмотря на фундаментальные исследования многих великих математиков в этой области, включая Лагранжа, Пуанкаре и других.

С математической точки зрения квантовая задача рассеяния трех тел кажется ничуть не проще, а вероятно, даже существенно сложнее классической. Например, до сих пор нет математически строгого решения задачи трех кулоновских частиц (то есть взаимодействующих друг с другом по закону Кулона) при энергиях выше трехчастичного порога и т. д. Тем не менее, огромный прогресс был достигнут в работах знаменитого советского и российского математика Людвига Фаддеева, который построил строгую математическую теорию для квантовой задачи трех тел (Л. Д. Фаддеев, 1966. «Математические вопросы квантовой теории рассеяния для системы трех частиц»), нашел математически правильное уравнение для задачи трехчастичного рассеяния и сформулировал условие разрешимости этих уравнений, которые теперь повсюду в мире известны как уравнения Фаддеева.

Эта во всех смыслах пионерская работа Фаддеева со времени своего появления в начале 60-х годов прошлого века, вызвала огромный резонанс во многих разделах теоретической и ядерной физики и стимулировала большой поток ежегодных конференций во всем мире по так называемым малочастичным системам. (В настоящее время в мире регулярно проводится по крайней мере три серии «малочастичных» конференций: мировая, европейская и Азиатско-Тихоокеанская, каждая с частотой раз в три года. В этом году европейская конференция проходит в середине августа в Дании. Отметим, что буквально вчера на этой конференции было принято решение об учреждении медали имени Людвига Фаддеева, которой будут отмечаться ученые, внесшие выдающийся вклад в квантовую теорию нескольких частиц.) Причем значительная часть исследований, представленных на этих конференциях, особенно в начальные годы после выхода пионерских работ Фаддеева, была посвящена методам решения его уравнений и различным аппроксимациям в задаче рассеяния трех тел. Это было обусловлено тем, что в общем виде уравнения Фаддеева, например, в интегральной формулировке представляют собой большую систему многомерных сингулярных интегральных уравнений, причем интегралы в ядре этих уравнений берутся с переменными пределами, зависящими от других переменных. Это, в частности, приводит к тому, что прямые итерации таких уравнений (когда в правую часть последовательно подставляется решение, найденное на предыдущем шаге итерации) оказываются невозможными, так как неизвестные функции в левой и правой частях уравнения определены для разных аргументов. В итоге при итерациях приходится переинтерполировать текущее решение в каждой точке на огромной многомерной решетке. Это приводит к необходимости многомерных интерполяций текущего решения на каждом шаге итераций, общее число таких переинтерполяций достигает миллионов и даже десятков миллионов. Сюда нужно еще добавить проблемы с регуляризацией сложных сингулярностей в ядре уравнений Фаддеева и очень громоздкие схемы нахождения матриц перевязки большого числа угловых моментов квантовых частиц, что приводит к суммам по промежуточным проекциям угловых моментов частиц высокой кратности (вплоть до 30-кратных сумм) для нахождения каждого матричного элемента матрицы ядра (а таких матричных элементов в реальной задаче может быть много миллионов и даже миллиардов).

В итоге вычислительная сложность решения уравнений Фаддеева для реальных задач в ядерной и атомной физике оказалась такой, что вплоть до появления больших суперкомпьютеров, главным образом в США, в ередине 80-х годов прошлого века, найти точное решение такой задачи за приемлемое время было невозможно. Более того, в конце 80-х годов время решения полной трехчастичной фаддеевской задачи на большом американском суперкомпьютере составляло 2-3 недели (см., например, статью A. Picklesimer, R. A. Rice, and R. Brandenburg, 1991. D degrees of freedom in trinuclei: I. The Hannover one-D model и другие статьи этого коллектива авторов из серии «D degrees of freedom in trinuclei»). Было совершенно очевидно, что никакой обычный компьютер с таким объемом вычислений не справится. Однако очень быстрый прогресс в вычислительной технике, достигнутый в последние десять лет, превзошел самые смелые ожидания и самые дерзкие мечты. И теперь решение полной системы фаддеевских уравнений для трехчастичной системы можно получить на обычной настольном компьютере за 15-30 минут. Статья группы российских ученых с этим результатом опубликована в июльском номере Computer Physics Communications.

Волновые пакеты и дискретизация непрерывного спектра квантовой системы

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

Суть этого подхода в его современной интерпретации состоит в разбиении непрерывного спектра гамильтониана системы на полосы (бины) конечной ширины D. В каждой такой полосе из точных функций непрерывного спектра путем интегрирования по энергии строится волновой пакет (поскольку он не зависит от времени, мы называем его стационарным). В результате, получается дискретный набор пакетированных состояний непрерывного спектра, который в дальнейшем используется как базис для вычислений. Следует отметить, что интегрирования по полосам дискретизации оказывается достаточно, чтобы волновые функции стали нормируемыми подобно функциям дискретного спектра. Отличие таких пакетных состояний от связанных состояний гамильтониана состоит в том, что мы можем строить разбиения непрерывного спектра произвольным образом, при этом энергии пакетных состояний будут задаваться различными сетками дискретизации. В пределе бесконечно малых ширин разбиения, результаты полученные в таких различных базисах стремятся к точному решению для исходного непрерывного спектра

В импульсном пространстве волновой пакет имеет вид одиночного импульса (рис. 3, слева), а пакетированная волновая функция (например, волновая функция дейтрона - связанного состояния нейтрона и протона) имеет вид гистограммы (рис. 2, справа).

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

Например, вместо интегрального уравнения Липпмана-Швингера для t-матрицы, описывающей рассеяние частицы силовым полем V в квантовой механике получим простое матричное уравнение

где все операторы представлены матрицами в пакетном базисе и g_k -элемент матрицы свободной резольвенты, имеющей диагональный вид. Причем, благодаря дополнительному усреднению энергетических зависимостей операторов по полосам дискретизации, которое мы используем, особенности исходного интегрального ядра сглаживаются, и вместо сингулярного интегрального уравнения получаем легко решаемое матричное уравнение.

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

Существенным преимуществом пакетной схемы является «пиксельная форма» матричных уравнений, когда элементы всех матриц вычисляются независимо друг от друга. Это позволяет относительно «несложно» распараллелить очень трудоемкий расчет огромного числа матричных элементов многомерного интегрального ядра, так что время расчета при такой массивно-параллельной реализации можно сильно (то есть в десятки и сотни раз) сократить. Для реализации этого удачным решением является использование графического процессора (GPU), или видеокарты, которой оснащен практически каждый компьютер для работы с компьютерной графикой. Разумеется, для таких массивно-параллельных вычислений подходит далеко не всякая графическая карта, а только специализированная под научные вычисления. Но они теперь уже широко представлены на рынке. Например, компания NVIDIA в последние годы расширила свою основную специализацию с разработки графических карт для игровых приставок на их адаптацию под чисто научные задачи. Вообще, можно сказать, что в связи с широким распространением графических вычислений во всех областях современной науки на наших глазах происходит компьютерная революция в научных вычислениях и широкий переход к массивно-параллельной реализации многих компьютерных программ.

«GPU - компьютерная революция»: что это такое?

Графическая карта - один из основных элементов игровых приставок и большинства компьютеров - позволяет в режиме реального времени генерировать быстро сменяющие друг друга картинки на экране. Для этой цели графическая карта включает программируемые арифметические устройства - пиксельные шейдеры (pixel shaders), которые позволяют вычислить состояние пикселя с координатами (x, y) на экране. В общем случае пиксельный шейдер получает на входе координату и другую дополнительную информацию и должен выдать конечный цвет данного пикселя, включая его освещенность, степень затенения и т. д. При этом, поскольку при работе в масштабе реального времени изображение на мониторе (которое сейчас уже включает миллионы пикселей) должно сменяться много раз за секунду, скорость этих графических устройств должна быть огромной.

В начале 2000-х годов для программистов, использующих 3D-графику, были разработаны специализированные программные средства, OpenGL и DirectX, позволяющие писать инструкции для GPU. Однако почти сразу им пришло в голову, что вместо цвета пикселя шейдер может вычислять любую величину, привязанную к координатам. Если на вход подавать числовые данные, содержащие не цвета, то ничего не мешает написать шейдер, выполняющий с ними произвольные вычисления. Результат можно считать в программу, а GPU безразлично, как именно он интерпретируется. В итоге оказалось, что можно использовать очень высокую скорость арифметических вычислений с помощью графической карты для многих целей, не имеющих ничего общего с графикой на компьютере. Для облегчения программирования даже была создана специальная архитектура CUDA. В результате была создана реальная возможность для сверхбыстрых вычислений с помощью массивно-параллельной реализации.

В настоящее время серийные графические карты, которые свободно продаются на рынке, включают в себя порядка 4 тысяч вычислительных ядер (а некоторые - даже больше!). Сравните это с 4-8 ядрами для обычного процессора. Они позволяют производить параллельные вычисления посредством многих десятков и сотен тысяч потоков - так называемых нитей (threads). И хотя скорость вычисления вдоль каждой такой нити несколько уступает скорости центрального процессора, параллельное использование огромного числа нитей позволяет добиться ускорения на порядки величин.

Однако не все так просто. В силу основных особенностей своей конструкции, GPU-вычисления выполняются в моде SIMD (single instruction on multiple data), что подразумевает, что вычисления вдоль всех нитей должны выполняться по одной и той же инструкции (команде), и, следовательно, тут не допускаются условные операторы и переходы. Есть и другие ограничения.

Поэтому сама возможность сверхбыстрых вычислений на GPU предполагает нетривиальное умение эффективно распараллелить полную задачу с учетом ограничений, налагаемых SIMD-модой. Самый простой и очевидный пример для большого ускорения вычислений посредством GPU - это вычисление элементов матрицы высокой размерности. Обычно для элементов таких матриц имеются аналитические формулы (разной степени сложности), и тогда полная матрица делится на много тысяч или миллионов блоков, причем матричные элементы внутри каждого блока считаются посредством параллельных нитей (по одним и тем же формулам).

Однако часто полная матрица оказывается такой высокой размерности, что она не помещается ни в кэш-память GPU, ни в быструю память (RAM) компьютера, и ее можно записать только на жесткий диск. Но при этом многократное считывание кусков матрицы с диска занимает так много времени, что весь выигрыш от параллельной реализации пропадает. Тем не менее, имеется много способов обойти эти ограничения и добиться высокой степени ускорения GPU-вычислений.

На сегодняшний день ультрабыстрые GPU-вычисления внедрены в большое число научных проектов и приложений, в частности, в квантовой физике многих тел, в квантовой теории поля, химии больших молекул, в геофизике, медицине, финансовой математике и т.д. При этом во многих приложениях использование GPU-вычислений дает совершенно уникальные новые возможности (см., например, книги «Технология CUDA в примерах. Введение в программирование графических процессоров» и «CUDA Fortran для ученых и инженеров»).

Использование GPU-вычислений в квантовой теории многочастичного рассеяния

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

где U - искомая матрица оператора перехода, V - матрица оператора взаимодействия двух частиц (которая имеет простую блочную структуру) и G₁ - матрица резольвенты гамильтониана канала (имеющая диагональный вид). Самым сложным с технической точки зрения элементом этого уравнения является матрица оператора перестановки  , которая имеет полную размерность системы

Чтобы понизить размерность этого уравнения, обычно выполняется разложение всех операторов по сферическим гармоникам, затем получившиеся угловые моменты по каждой (из двух) координат Якоби связываются с переменными спина и изоспина. В итоге для полного момента всей системы J получается огромная система зацепляющихся многомерных интегральных уравнений (в среднем 40-60 зацепляющихся уравнений для каждого значения J). Для решения всего одного из этих уравнений нужно сделать дискретизацию, отвечающую проектированию на базис стационарных волновых пакетов с использованием 50-100 бинов по каждой координате Якоби. Тогда в среднем получается система алгебраических уравнений порядка N = 10⁴-60 = 600 000, с полной размерностью матрицы N-N. Такая матрица состоит из 300-500 миллиардов матричных элементов, так что она не влезает в быстрой памяти компьютера, ее можно разместить только на большом жестком диске. К счастью, основная матрица оператора перестановки представляется в виде ( =OP_0O^T), где O - блочная матрица и ₀ - очень разреженная матрица перекрывания для плосковолновых функций в различных наборах трехчастичных импульсов, которая заполнена едва ли на 1% (или даже меньше). Это позволяет разместить матрицу ₀ в оперативной памяти, однако появляется необходимость предварительной проверки всех ее элементов и отбора среди них отличных от нуля. Такая проверка миллиардов матричных элементов тоже требует много компьютерного времени. И, наконец, следует упомянуть, что вся эта огромная работа должна выполняться для каждого значения полного момента J, которое может принимать полуцелые значения 1/2, 3/2, ..., 17/2!

Оказалось, что распараллеливание на десятки тысяч потоков позволяет выполнить всю эту гигантскую численную работу на обычном персональном компьютере за приемлемое время.

Общая схема решения задачи в пакетном представлении состоит из следующих шагов:
    1) заполнение сеток дискретизации для каждого набора импульсных координат, а также вычисление блочных матриц V, O и диагональной матрицы G₁;
    2) предварительный отбор ненулевых элементов матрицы  ₀;
    3) вычисление ненулевых элементов матрицы  ₀;
    4) решение матричного уравнения методом итераций.

Поскольку самыми трудоемкими в этой схеме являются шаги 2) и 3), мы провели распараллеливание именно для них. При этом скорость вычислений увеличивается на порядки величин и достигает фантастических значений для распараллеленных частей полной программы в случае относительно простого вида двухчастичных взаимодействий. В частности, мы нашли, что полная матрица перекрывания  ₀ содержит примерно (для фиксированного J) 260 миллионов ненулевых элементов, при заполнении каждого из которых требуется численно посчитать двойной интеграл от полиномов Лежандра высокого порядка в комбинации с серией других арифметических операций. Мы с большим удивлением обнаружили, что при массивно-параллельном вычислении этих 260 миллионов сложных интегралов на обычном ПК с графическим процессором GeForce GTX 670 для этого требуется всего 3 секунды, а с учетом их предварительного отбора - всего 8 секунд.

Следует здесь отметить, что процедура предварительного отбора ненулевых элементов матрицы позволяет не только уменьшить количество используемой памяти, но и существенно важна для ускорения процедуры распараллеливания на GPU. Дело в том, что прямое вычисление элементов разреженной матрицы на GPU неэффективно, поскольку параллельные потоки оказываются загружены неравномерно, и большая часть процессоров вынуждена простаивать, пока остальные выполняют вычисления. После же предварительно отбора ненулевых элементов (который нам также удалось распараллелить!), параллельных потоков нужно значительно меньше и они загружены равномерно. Этот пример ясно показывает, каких поразительно высоких скоростей можно достичь при хорошо организованном вычислительном процессе на GPU.

В итоге полное решение огромной системы зацепляющихся интегральных уравнений для всех нужных J (с учетом всех шагов, в том числе и выполняемых последовательно) для полностью реалистических сложных двухчастичных взаимодействий можно выполнить за время порядка 15-30 минут. По контрасту, та же Фаддеевская задача на большом суперкомпьютере Blue Gene может занимать более одного дня и стоит очень дорого.

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

Источник: V. N. Pomerantsev, V. I. Kukulin, O. A. Rubtsova, S. K. Sakhiev. Fast GPU-based calculations in few-body quantum scattering // Computer Physics Communications. 2016. V. 204. . 121-131.

Владимир Кукулин

Источник: elementy.ru

Комментарии: