Высокоскоростные алгоритмы моделирования систем с дальнодействующими силами

Рассмотрим систему N классических частиц с массами m _i и зарядами q _i, динамика которой подчиняется уравнениям движения Ньютона, а взаимодействие зарядов системы с единичным зарядом определяется Кулоновским (дальнодействующим) потенциалом U:

m _i d ² r _i/ dt ² = – q _i∂ U /∂ r _ij

U = ∑ _{j> i} q _j| r _i – r _j|^–1

Простейшим методом ее моделирования является прямой расчет сил суммированием парных взаимодействий всех частиц системы друг с другом, с последующим численным интегрированием уравнений движения. Так как в такой системе каждая частица взаимодействует со всеми остальными, то число операций пропорционально N × (N –1). С помощью 3-го закона Ньютона (действие равно противодействию), можно учесть симметричность вклада каждого парного взаимодействия в суммарные силы, действующие на частицы, и сократить число операций вдвое, однократно рассчитывая этот вклад для каждой пары. Но такая оптимизация все же не устраняет квадратичную зависимость вычислительной нагрузки от количества частиц в системе, т.е. увеличение количества частиц в моделируемой системе в 10 раз все равно приводит к 100-кратному увеличению времени расчетов.

Активные попытки преодоления этого вычислительного барьера привели к появлению большого числа численных методов, в которых с целью экономии числа операций используется замена точного прямого расчета сил различными приближениями, а также всевозможные алгоритмические усовершенствования (см. обзорные работы [17], [18]). Условно их можно подразделить на 3 категории: суммирование Эвальда, сеточные методы и мультипольные или иерархические.

Метод Эвальда [19] позволяет рассчитать взаимодействие частиц некоторой «ячейки» (ограниченной области), с ее отражениями, бесконечно транслированными по одному или нескольким измерениям пространства, разделив медленно сходящийся ряд Кулона такой системы:

U = ∑ _n∑ _j q _j| r _ij+ n L |^–1

на сумму двух быстро сходящихся рядов в реальном и «обратном» пространствах:

U = ∑ _n∑ _j q _j| r _ij+ n L |^–1erfc(α | r _ij+ n L |) +

+ 4π L ^–3∑ _{k≠ 0}∑ _j q _j exp[– i kr _ij – k ²(2α)^–2] – q _i(2α)/π ^1/2

Здесь L – период решетки, α – произвольный, но существенный параметр, отвечающий за сходимость обоих рядов суммы, а последнее слагаемое – «потенциал самодействия», нивелирующий соответствующий вклад от ряда в «обратном» пространстве. Область применение метода строго ограничена периодическими системами, однако он широко используется в исследова­ниях, где важно исключить поверхностные эффекты, неизбежно возни­кающие при моделировании небольших изолированных систем. В работе Перрама и др. [20] показано, что существует оптимальный выбор параметров, при котором вычислительная нагрузка метода Эвальда растет пропор­ционально N^3/2.

Сеточные методы раскладывают заряд частиц по узлам дискретной сетки, затем решают на ней уравнение поля, используя различные быстрые алгоритмы (быстрое преобразование Фурье, метод Гаусса-Зейделя, последо­вательной релаксации), а потом интерполируют силы, рассчитанные в узлах сетки, чтобы получить силы, действующие на отдельные частицы. Основная вычислительная нагрузка связана с решением уравнений поля и для боль­шинства алгоритмов пропорциональна N _Glg N _G, где N _G – количество узлов сетки. Областью эффективного применения сеточных методов так же являются периодические системы, так как для моделирования изолированной системы ее приходится окружать зеркальными копиями, чтобы скомпенсиро­вать вклад транслированных отражений. Подробный анализ многочисленных методов этой категории выходит за рамки данной работы (см., например, две монографии по теме [21] и [22]).

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

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

В 1986 году Барнс и Хат опубликовали алгоритм TreeСode [23], в котором физическое пространство иерархически подразделялось на области для учета относительного положения частиц. Результирующая древесная структура (рис. 7.1), затем использовалась для группировки удаленных частиц в «кластеры», с заменой всех взаимодействий частица–частица одним взаимодействием частица–кластер, существенно уменьшая число операций.

Рис. 7.1. Древесная структура quadtree (двумерный случай).

Основой всех иерархических методов является древовидное представление физического пространства. Корень дерева представляет клетку, которая содержит все частицы системы. Дерево строится рекурсивно – в двухмерном варианте сначала делят корневую клетку на четыре клетки-потомка, а затем дробят полученные клетки, если количество частиц в них больше некоторого фиксированного числа. Полученное таким образом дерево называется деревом квадрантов (quadtree), поскольку у каждой клетки, не являющейся «листовой», по четыре клетки-потомка. Форма дерева соответствует распределению частиц в пространстве, поскольку дерево имеет больше уровней в тех областях, где больше частиц. В трехмерном пространстве родительские клетки де­лятся на восемь потомков, и структура данных называется деревом октантов (octree).

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

r _c = å r _i|q_i| / å |q_i|

M = å q_i

Dx = å q_ir_ix

Qxy = å q_i(3r_ixr_iy − r_i²δ _xy)

Здесь r _c – центр зарядов клетки, r_ix – х-компонента расстояния между i-той частицей и центром зарядов, M – монопольный момент распределения зарядов клетки, Dx и Qxy – соответствующие компоненты дипольного и квадрупольного моментов, δ _xy – символ Кронекера.

Затем проходом «вниз» (от корня к листьям) составляются «списки взаимодействий», когда для каждой частицы данной клетки определяется допустимость расчета ее взаимодействия с частицами «дальней» клетки как с кластером, на основе некоторого «критерия допустимости». Простейшим критерием является представленное в оригинальной работе отношение S/D – длины ребра клетки-кластера S к расстоянию до взаимодействующей частицы D (рис. 7.2).

Рис. 7.2. Отношение S/D для разных уровней дерева.

Если отношение меньше некоторого заданного параметра θ, регулирующего точность метода, то мультипольным разложением этой «дальней» клетки-кластера аппроксимируется действие всех вложенных в нее клеток-потомков на рассматриваемую частицу (R – расстояние от частицы до центра зарядов «дальней» клетки):

U = MR^–1 + R^–3å r_xD_x + R^–5å å r_xr_yQ_xy / 2

иначе по тому же критерию проверяются клетки-потомки. Взаимодействие с частицами клеток данной ветви дерева, которые не удовлетворяют критерию – считается прямым способом. Если задать параметр точности θ равным нулю, все частицы системы будут считаться прямым способом, на практике θ выбирают в диапазоне 0.3 – 1.0, в зависимости от задачи.

На каждом временном шаге дерево перестраивается, поскольку частицы перемещаются, и, следова­тельно, их распределение меняется. Однако, перемещения относительно малы, поэтому дерево можно перестраивать, когда частицы перемещаются за пределы своей клетки. В любом случае центры зарядов и мультипольные разложения всех клеток нужно пересчитывать на каждом временном шаге. Число узлов в дереве октантов равно log₂N^1/3 = lgN/3lg2 ~ lgN, поэтому построение дерева требует порядка NlgN операций. Также было показано [24], что количество взаимодействий выделенной частицы с однородным сферическим облаком пропорционально lgN/θ ², таким образом, вычисление сил в системе из N частиц также требует порядка NlgN операций.

Наличие дополнительных операции в алгоритме TreeСode, отсутствующих в прямом расчете: построение и заполнение дерева и списков взаимодействий частиц, расчет центров зарядов и мультипольных разложений, однозначно говорит о том, что метод будет эффективнее прямого расчета лишь при количествах частиц в системе, превышающих некоторое значение, которое определяется в основном выбором параметра точности θ.

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

В 1987 году, вскоре после появления TreeСode, Грингард и Рохлин опубликовали альтернативный быстрый мультипольный метод – Fast Multipole Method (FMM) [25], который может рассматриваться как элегант­ное усовершенствова­ние алгоритма TreeСode, однако, был разработан не­зависимо. В настоящее время FMM входит в десятку лучших алгоритмов 20-го века и широко используется во многих областях физики, математики и компьютерных наук. Алгоритм FMM основан на том факте, что мультипольное разложение с бесконечным порядком членов содержит полную информацию о распределении зарядов в пространстве. Как и в методе TreeСode взаимодействие между ближайшими соседями вычисляется прямым расчетом, а дальние частицы рассматриваются группами, однако, эти вклады считаются по-другому.

Первое основное отличие FMM – дальнодействие рассматривается как единый общий вклад «дальнего поля», который вычисляется как мультиполь­ное разложение высокого порядка. Метод TreeСode рассматривает только взаимодействия типа частица-частица и частица-клетка, в FMM помимо этого вычисляются взаимодействия типа клетка-клетка, как было показано в оригинальной работе [25], это позволяет полу­чить алгоритм с линейной зависимостью вычислительной нагрузки от числа частиц.

FMM a priori позволяет получить произвольную точность (в пределах машинной), если используется мультипольного разложение достаточного порядка, поэтому существенным является выбор наименее ресурсоемкого численного представления мультипольных разложений, операторов их комбинации и трансляции их центров. Первые варианты FMM были двумерны, ориентированы на однородное распределение зарядов и использовали сложную мат. модель представления полей и потенциалов [25], позднее была сформулирована адаптивная трехмерная реализация, использующая сферические гармоники [26].

Второе основное отличие FMM от TreeСode состоит в способе, по которому управляют точностью аппроксимации сил. Метод TreeСode представляет систему в виде дерева мультипольных разложений в центрах зарядов отдельных клеток, и определяет допустимость аппроксимации отношением расстояния между рассматриваемой частицей и центром заряда удаленной клетки к длине ее ребра. В FMM используется иерархия вложенных сеток, а распределение заряда в клетке представля­ется с помощью мультипольного разложения около геометрического центра клетки, а также используется другой критерий допустимости мультипольного разложения: считается, что две клетки достаточно удалены друг от друга, если расстояние между ними превосходит длину ребра клетки на данном уровне иерархии. Таким образом, аппроксимации в FMM используются намного чаще, чем в методе TreeСode, но это компенсируется более точным описанием распределения зарядов внутри клетки.

Кратко алгоритм FMM устроен следующим образом (см. рис. 7.3):

Рис. 7.3. Упрощенная схема алгоритма FMM.

Построение сеточной иерархии, делением клеток на потомки, начиная с корневой, включающей все частицы системы – проход «вниз» (от корня к листьям).

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

Трансформация мультипольных разложений всех клеток в локальные (разложения в ряд Тейлора около центров клеток), применение локального разложения каждой клетки к клеткам из ее списка взаимодействий. Трансляция (с добавлением) локальных разложений клеток с верхних уровней своим клеткам-потомкам – проход «вниз» (от корня к листьям). Эта дополнительная фаза реализует взаимодействие клетка-клетка, в результате в каждой клетке самого нижнего уровня будет находиться локальное разложение поля всей системы зарядов около геометрического центра данной клетки.

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

На первый взгляд, FMM, с его теоретически предсказываемой линейной зависимостью времени расчетов от числа частиц, превосходит все остальные методы. Однако, этот алгоритм содержит существенную дополнительную вычислительную нагрузку, связанную с расчетом мультипольных и локальных разложений и их трансформацией и трансля­циями, поэтому на практике вопрос ставится иначе: при каком числе частиц стоит использовать FMM вместо TreeCode или прямого расчета? Ответ на него не прост, так как время расчетов зависит не только от числа частиц, но и от требуемой точности вычислений.

В алгоритме TreeСode точность определяется параметром θ, для FMM ту же роль играют число членов L в разложениях и глубина R иерархической структуры, Шмидт и Ли показали [27], что его полиномиальная зависимость от N, L и R такова:

P = aBL² + bNL² + cBL⁴ + dB∙ (g1 + g2(N/B) + g3(N/B)²).

Здесь a, b, c, g1, g2, g3 – машинно-зависимые параметры, а B = 8R количество клеток на нижнем уровне иерархии R. Это выражение показывает, что отношение (N/B) должно быть мало, чтобы квадратичной зависимостью можно было пренебречь в пользу линейной – на практике это означает увеличение глубины сеточной иерархии с увеличением числа частиц.

Отметим также, что одно из предположений оригинального варианта FMM – более-менее однородное распределение зарядов в пространстве; существенная неоднородность потребует одновременно большей глубины иерархии и большего числа членов в разложениях, что приведет к увеличению вычислительной нагрузки. Поэтому FMM в оригинальной форме плохо применим к неоднородным и динамичным системам, где возможны резкие изменения плотности с течением времени. Для решения этого вопроса были разработаны адаптивные методы [26], в большинстве своем являющиеся гибридами TreeCode и FMM.

Приведем график (см. рис. 7.4) сравнения производительности разных алгоритмов: TreeCode, FMM и прямого расчета из работы [17], данные производительности в этой работе были получены на вычислительной системе SunBlade Sparc, на модели с однородным начальным распределением зарядов, там же приводим данные производительности нашей собственной поточно-параллельной версии прямого расчета. Для метода TreeСode, практически независимо от числа частиц в системе, выбор параметра θ равным 0.2 и 0.5 дает порядок ошибки при вычислении сил 0.1% и 1%, соответственно. А точность метода FMM зависит не только от количества членов в разложениях, но и от глубины иерархической структуры, которую необходимо увеличивать с ростом числа частиц в системе, этот эффект более детально исследовал Эзелинк [28]. На графике приводятся данные его работы для разложений 4-го, 7-го и 10-го порядков.

Кривые на рис. 7.4 демонстрирует несколько важных моментов: алгоритм TreeСode превосходит FMM в приложениях, где приемлема низкая точность вычисления сил ~0.1–1%, в которых он обгоняет прямой расчет на системах из 5000 частиц и более, тогда как FMM при той же точности делает это лишь на системах из 10000 частиц; однако если точность играет решающее значение, то FMM подтверждает свое преимущество над остальными методами, правда лишь на системах из 100000 частиц и более. Наша поточно-параллельная версия прямого расчета превосходит по скорости остальные методы вплоть до точки, соответствующей системам размером в 300000 частиц, давая при этом максимальную точность (в пределах машинной погрешности).

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

Рис. 7.4. Сравнение зависимости времени расчетов одного шага ММД от числа частиц в системе для разных методов по данным работы [17] и нашей поточно-параллельной методики.

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

Параллельная версия классического метода Эвальда достаточно легко реализуема: во-первых, непосредственно разделяется расчет рядов в реальном и «обратном» пространстве; во-вторых, вычисление ряда в реальном пространстве может быть распределено между процессорами по аналогии с прямым расчетом; а ряд в «обратном» пространстве может быть преобразован в суммирование по частицам, умноженное на структурный фактор, которое также считается параллельно [29].

Распараллеливание сеточных методов, использующих быстрое преобразование Фурье (FFT) – довольно проблематично, т. к. бинарная рекурсивная природа FFT мало пригодна для реализации на архитектурах с распределенной памятью. Беккером была предложена альтернативная схема для решения уравнений поля на сетке [30], в которой вместо FFT исполь­зуются итеративные методы, пригодные к распределенной параллельной реализации.

Реализация параллельных версий мультипольных алгоритмов требует наибольших дополнительных усилий, но привлекает потенциальной возможностью моделирования систем из 10⁹ частиц на современных архитектурах с терафлоповой производительностью. Параллельная версия неадаптивного, двумерного FMM была предложена самим Грингардом в начале 90-х [31]. Эта схема, построенная на распараллеливании независимых стадий алгоритма, относительно эффективна на параллельных архитектурах с общей памятью, но мало пригодна для распределенных параллельных архитектур.

Из всех рассмотренных алгоритмов TreeCode Барнса-Хата, вероятно наиболее трудоемок в параллельной реализации. На параллельных архитектурах с общей памятью векторизация стадий построения дерева и создания списков взаимодействий позволяет добиться относительно неплохого прироста в скорости. Однако на распределенных архитектурах возникает трудно решаемая проблема – для построения списков взаимодействий необходим доступ ко всей структуре дерева, а поиск удаленных узлов дерева, находящихся в памяти других процессоров с использованием адресации указателями, которая типична для скалярных версий TreeCode, приводит к огромным коммуникационным издержкам. Сэлман и Уоррен [32] для решения этой проблемы изобрели практически новый алгоритм, отказавшись от указателей в пользу хэшированной схемы адресации.

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