Содержание к диссертации
Введение
2. Математическая модель 23
2.1. Моделирование методами молекулярной динамики 23
2.2. Расчет макроскопических свойств 25
2.3. Ансамбли 26
2.4. Начальные условия 27
2.5. Граничные условия 27
2.6. Численная схема 29
2.7. Модель неполярных молекул 31
2.8. Модель полярных молекул на примере молекул воды 33
2.9. Модель металлических включений на примере платины 35
2.10. Модель вода-платина 36
2.11. Метод моделирования гетерогенной кавитации 37
2.12. Термостатирование 38
2.13. Обезразмеривание 40
2.14. Общий алгоритм моделирования методом молекулярной динамики 41
3. Методы ускорения расчётов 43
3.1. Ускорение при помощи архитектурных решений 46
3.2. Иерархическая структура данных 51
3.3. Использование иерархической структуры данных для расчета ближнего взаимодействия 63
3.4. Быстрый метод мультиполей 69
3.5. Быстрый метод мультиполей для гибридных архитектур 79
3.6. Одноуровневая структура данных для ближнего взаимодействия 86
3.7. Краткое описание комплекса программ 104
4. Тестовые и практические расчеты 106
4.1. Динамика системы пар-жидкость 106
4.2. Уравнение состояния и вычисление давления 107
4.3. Функция радиального распределения для молекул воды 110
4.4. Выбор радиуса обрезки потенциала Леннард—Джонса 113
4.5. Исследование растекания капли воды по поверхности платины 116
4.6. Исследование масштабируемости в задаче гетерогенной кавитации 118
4.7. Многокомпонентная нуклеация вблизи подложки 120
Заключение 125
Список условных сокращений 127
Литература 128
- Начальные условия
- Метод моделирования гетерогенной кавитации
- Использование иерархической структуры данных для расчета ближнего взаимодействия
- Исследование растекания капли воды по поверхности платины
Начальные условия
Методы молекулярной динамики (МД) могут применяться для исследования динамических свойств газов, жидкостей, твердых тел и взаимодействий между ними. Кроме методов МД, для моделирования на молекулярном уровне применяются метод твердых сфер [62], метод Монте-Карло [7; 20], методы основанные на квантовой теории (MM/QM) [17]. Выбор метода моделирования основан на решаемой задаче.
Первые попытки моделирования методами МД были осуществлены с появлением первых компьютеров в конце 50-х — начале 60-х годов XX века [19; 53; 98]. Одна из первых таких работ под авторством B.J. Alder и Т.Е. Wainwright вышла в 1957 году [19]. Целью этой работы было исследование фазовой диаграммы системы молекул, представленных в виде твердых сфер, и, в частности, области твердого тела и жидкости. В работе моделировалась динамика 32 и 108 молекул, взаимодействующих как бильярдные шары. В 1960 г. вышла работа J.B. Gibson и др. [53], ставшая первым примером моделирования с непрерывным потенциалом. В работе 1964 г. A. Rahman [98] исследовал свойства жидкого аргона, используя потенциал Леннард—Джонса.
В методе МД для описания движения атомов применяется классическая механика и силы межатомного взаимодействия представляются в форме потенциальных сил. Параметры потенциала могут быть получены из теоретических предпосылок, полуэмпирических или из физических экспериментов. Потенциал описывает два типа молекулярных свойств: 1) связанное взаимодействие (силы, связывающие атомы внутри одной молекул), которое характеризует растягивание межпарных связей, изгиб валентных углов, вращение двугранных углов; 2) несвязанное взаимодействие, которое характеризует дисперсию, электроста тическое взаимодействие и т.д. Также потенциал может включать воздействие внешних сил, например, внешних стенок. Именно расчет несвязанного взаимодействия является самым ресурсоемким, расчет связанного взаимодействия имеет линейную вычислительную сложность. Поэтому в данной работе рассматривается только несвязанное взаимодействие, как наиболее вычислительно сложное, и, соответственно, используются модели жестких молекул. В случае моделирования жестких молекул основным является парное взаимодействие [24]. Потенциалы, применяемые в методе МД для моделирования несвязанного взаимодействия, могут описывать близкодействующие взаимодействия (на каждый атом существенно влияние лишь его ближайшего окружения), к таким относится, например, потенциал Леннард—Джонса, и дальнодействую-щие, например, потенциал Кулона, описывающий электростатическое взаимодействие.
За последние полстолетия использования методов МД с развитием вычислительных технологий и разработкой новых алгоритмов количество частиц, которые возможно промоделировать выросло на несколько порядков: с десятков атомов до триллионов. Так, например, в 1964 в работе A. Rahman [98] численно исследовалась система из 864 атомов аргона взаимодействующих при помощи потенциала Леннард—Джонса с использованием периодических граничных условий. В 2008 году Т.С. Germann и К. Kadau [52] провели демонстрационное моделирование триллиона леннард-джонсовских частиц (с радиусом обрезки 2.5 т) на суперкомпьютере BlueGene/L. Однако расчеты, подобные последнему, носят чисто демонстрационный характер — использование и обработка такого большого числа частиц при проведении широкомасштабных исследований пока не представляются возможными.
Рассмотрим подробнее способы ускорения МД расчетов: алгоритмические и с использованием высокопроизводительных вычислительных технологий.
Существуют разные методы алгоритмического ускорения расчета ближнего взаимодействия. Все они основаны на быстроубывающей природе описывающего его потенциала и, следовательно, малом радиусе взаимодействия. В 1967 г. L. Verlet [122] предложил метод («список Верле»), суть которого заключается в следующем: каждые п шагов для каждого атома АІ создается список всех атомов, находящихся внутри сферы с центром в атоме АІ И радиусом много большим радиуса обрезки потенциала взаимодействия; последующие п—\ шагов при расчете взаимодействия на атом АІ учитываются лишь атомы, находящиеся в этом списке, с проверкой по радиусу обрезки. Радиус сферы и количество шагов выбираются так, чтобы атомы не успели перелететь из сферы, ограниченной радиусом обрезки, в область вне большей сферы и наоборот за это количество шагов. Данный метод позволяет ускорить расчет ближнего взаимодействия в десятки раз, однако, сложность создания такого списка 0(7V2), где N — число атомов, также этот список требует большого объема памяти для его хранения. В 1987 г. М.Р. Allen и D.J. Tildesley [20] предложили метод списка по ячейкам, смысл которого заключается в разделении всей моделируемой области на непересекающиеся ячейки, распределении частиц по этим ячейкам согласно их координатам и использовании связанных списков (список частиц, находящихся в боксе) для определения взаимодействия частиц друг с другом.
Алгоритмическое ускорение расчета дальнего взаимодействия является более сложной проблемой. При использовании алгоритма прямого суммирования (расчет всех парных взаимодействий напрямую), вычислительная сложность алгоритма оценивается как 0(N2), где N — число частиц. Одним из самых первых подходов, примененных для уменьшения его вычислительной сложности, было использование обрезки потенциала на некотором радиусе подобно ближнему взаимодействию. Однако ввиду медленно убывающей природы потенциалов дальнего взаимодействия такого рода метод является очень неточным. Другим методом является прямое суммирование с использованием оценочных точек [117]: в моделируемой области распределяются оценочные точки (например, по прямоугольной сетке); влияние каждого атома в системе суммируется в каждой оценочной точке (потенциальная энергия рассчитывается в дискретных координатах сетки); вместо расчета взаимодействия каждого атома с каждым, рассчитывается влияние оценочных точек на атомы. Вычислительная сложность такого метода 0(NM): где N — число атомов, М — число оценочных точек, однако, этот метод непрактичен, так как с ростом N значительно растет и число оценочных точек.
Метод моделирования гетерогенной кавитации
Ввиду быстро убывающей природы потенциала Леннард—Джонса, при расчетах он обрезается по радиусу rcutoff- Потенциал Кулона в свою очередь убывает медленнее, и обрезание кулоновского взаимодействия может сильно сказаться на результатах моделирования. Поэтому необходим расчет попарного электростатического взаимодействия всех заряженных частицы друг с другом.
Отметим, что предлагаемые в работе подходы и методы ускорения расчетов могут быть применимы для моделирования других типов молекул, в которых взаимодействие моделируется потенциалом Кулона и потенциалом ближнего взаимодействия, например, потенциалом Леннард—Джонса.
Для моделирования твердых включений (дисперсные частицы, подложка) была выбрана модель платины. Атомы платины располагаются согласно кристаллической решетке fee (111). На рис. 2.5 (а) показана схема размещения атомов согласно данной кристаллической решетке: атомы располагаются с учетом периодичности, причем вектор нормали к наклонным (диагональным) плоскостям равен п= (1,1,1). Параметр кристаллической решетки равен 3.92 А. В качестве потенциала взаимодействия между атомами платины используется потенциал Леннард—Джонса (2.8). Взаимодействие между атомами платины и атомами неполярных молекул также описывается при помощи потенциала Леннард—Джонса с параметрами сг, є, вычисленными по правилу Лоренца— Бертло (2.9).
Потенциал взаимодействия воды и платины должен удовлетворять как минимум двум условиям. Во-первых, корректно описывать абсорбционные свойства атомов платины. Во-вторых, расстояние между атомами платины и атомами кислорода должно быть меньше, чем между атомами платины и атомами водорода. Существует несколько потенциалов, которые удовлетворяют вышеописанным требованиям. В работе используется модель, разработанная S.-B. Zhu и M.R. Philpott в 1994 [127]. Потенциал взаимодействия состоит из потенциала отражения (UH2o-cond) и короткодействующих анизотропного (Uan) и изотропного (Uisr) потенциалов:
Процесс кавитации (образование пузырька в объеме жидкости) возникает в ряде областей гидродинамики и имеет важное значение в различных технологических процессах. Кавитация происходит при пониженном давлении в системе, находящейся в метастабильном или нестабильном состоянии.
В работе [81] предлагается метод моделирования гомогенной кавитации путем приложения дополнительного поля репульсивных внешних сил в центр области для образования зародыша с контролем количества молекул в области. Однако такой метод не очень удобен в применении. Другой методикой моделирования процесса кавитации является моделирование растянутой жидкости. Исследованию гомогенной кавитации в растянутой леннард-джонсовской жидкости посвящена работа [76]. В работе [75] исследуется кавитация на гидрофобных частицах в растянутой жидкости (воде), находящейся в свободном объеме. Однако кавитация в метастабильном состоянии носит случайный характер и время, необходимое для прослеживания зарождения в системе пузырька, может значительно превышать возможные времена моделирования методами МД. Поэтому для моделирования гомогенной кавитации методами МД рядом исследователей [1; 74] применяется техника постепенного понижения в системе давления за счет растяжения области моделирования и приведения системы к спинодали. Так в работе [97] исследуется кавитация в бинарной леннард-джонсовской смеси при помощи этого метода с переводом системы из стабильного в метастабильное состояние. В большинстве случаев кавитация происходит не в чистой жидкости, а в присутствии инородных включений (гетерогенная кавитация). Предлагается следующий метод моделирования гетерогенной кавитации в системах неполярных молекул, основанный на подходе к моделированию гомогенной кавитации. Моделируемая область равномерно заполняется молекулами жидкости. В центр области вместо молекул жидкости помещается твердое включение, вырезанное из кристаллической решетки соответствующего вещества. Для моделирования эффекта кавитации после достижения системой равновесия давление в системе равномерно понижается за счет постепенного (раз в несколько десятков тысяч шагов по времени) растяжения области (масштабирования ребер моделируемой ячейки). Таким образом система из стабильного состояния приводится в метастабильное ближе к спинодали. При растяжении положения атомов флюида смещаются соответственно, атомы твердого включения не сдвигаются. Важным является моделирование смачиваемости включения, так как от него зависит место образования пузырька: на включении или в объеме. Смачиваемость включения и подложки моделируется путем изменения леннард-джонсовского параметра энергии взаимодействия атомов флюида с атомами включения ssi. Параметр ег = єаі/єц — параметр лиофильности (ец — параметр потенциала Леннард—Джонса взаимодействия атомов флюида друг с другом). С ростом єг лиофильность включения/подложки увеличивается.
Рассматривается канонический NVT—ансамбль, где температура выступает в качестве контролирующего параметра. Полная энергия при моделировании канонического ансамбля не является постоянной, постоянной во времени остается температура. Для поддержания постоянной температуры, а следовательно кинетической энергии, используются термостаты, которые корректируют скорость атомов, чтобы обеспечить постоянство усредненной температуры всей системы. Для простоты приведем термостаты для случая одноатомных молекул. В случае рассмотрения жестких многоатомных молекул термостаты вводятся аналогичным образом на основе суммарной кинетической энергии поступательного и вращательного движений [102].
Термостат Берендсена. В основе термостата Берендсена [30] лежит масштабирование скоростей частиц. Кинетическая энергия системы, состоящей из одноатомных частиц, принудительно задается равной 3NkbTo/2, где къ — кон станта Больцмана. Корректируется отклонение температуры Т системы от заданной температуры То при помощи масштабирования скоростей частиц на некий коэффициент A: v = AVJ , чтобы скорость изменения температуры была пропорциональная разнице температур: где тт — параметр времени связи системы и термостатирующей ванны, определяющий характерное время, за которое достигается необходимая температура. Это достигается при где At — шаг численного интегрирования. При значениях тт близких к значению At термостат Берендсена обращается в простое масштабирование скоростей на каждом шаге по времени. Как правило, при At = 10 15 с параметр тт выбирается порядка
Термостат Нозе—Гувера. Данный термостат был предложен Нозе [87-89] и Гувером [69] и относится к методам с обратной связью. В случае термостата Нозе—Гувера имеется возможность добавления степени свободы , с помощью которой имитируется взаимодействие частиц с термостатом. Совместные уравнения, связывающие движение частиц системы и влияние термостата, имеют следующий вид: где r — подбираемый параметр, описывающий прочность связи системы с термостатом (имеет размерность времени, поэтому применим термин «время релаксации»). Чем меньше значение т, тем сильнее связь с термостатом. Значение г подбирается в зависимости от моделируемого процесса. При задании начальных условий таких, что скорости молекул соответствуют требуемой температуре, в начале полагается равной нулю.
Использование иерархической структуры данных для расчета ближнего взаимодействия
Выбор глубины восьмеричного дерева 1тах основан на балансе между временем вычисления потенциалов и временем генерации ИСД с учетом гетерогенности алгоритма, о которой будет говориться в следующем разделе. Также структура данных может быть эффективно применена для уменьшения времени расчета ближнего взаимодействия. Это применение подобно расчету разреженного матрично-векторного произведения для потенциала дальнего взаимодействия — рассчитывается взаимодействие только с частицами, находящимися в Е ііті I max) (n — мортоновский индекс бокса, в котором находится рассматриваемая частица), и попадающими в сферу с радиусом обрезки для ближнего взаимодействия (центр сферы в рассматриваемой частице). Вычислительная сложность уменьшается за счет того, что отпадает необходимость в проверке взаимодействия всех частиц в области моделирования. В этом случае выбор Imax зависит и от величины радиуса обрезки rcutoff используемого потенциала ближнего взаимодействия: размер бокса на нижнем уровне (1тах) выбирается большим по сравнению с rCMt0//, для того чтобы сфера взаимодействия радиуса rcutoff для рассматриваемой частицы вмещалась в область, ограниченную самим боксом, где частица находится, и ближайшими соседями этого бокса (область І2(п,Ітах), где п — мортоновский индекс бокса, в котором находится рассматриваемая частица). Увеличение /тоаж, ввиду увеличения общего числа боксов, ведет к повышению затрат на генерацию ИСД, затрат на трансляции в методе FMM и на копирование информации между GPU и CPU о структуре данных и о коэффициентах мультипольного разложения для метода FMM. С другой стороны это уменьшает время, необходимое для расчета разреженного матрично-векторного произведения, в том числе ближнего взаимодействия. Поэтому выбор Ітах является нетривиальной задачей и также зависит от параметров моделирования (например, плотности распределения частиц).
В тестах на производительность генерации ИСД используется равномерное распределение частиц по области, рассматривается жидкий аргон плотностью р = 1000 кг/м3. Радиус обрезки ближнего взаимодействия rcut0ff = 5 7. Результаты приведены для вычислений с использованием чисел с плавающей точкой одинарной точности. Использовалась NVIDIA CUDA версии 3.5.
Рисунок 3.7 иллюстрирует зависимость времени генерации ИСД от числа частиц для нескольких глубин (Ітах) восьмеричного дерева. Результаты представлены для чисел с плавающей точкой одинарной точности, для чисел двойной точности характер результатов (вид графиков) остается таким же. Сплошными линиями показано время генерации ИСД на CPU, пунктирными — на GPU. Видно, что для каждого 1тах линии имеют одинаковые наклоны. Плавный рост, который виден в конце линий и превышение времени генерации ИСД на GPU для 1тах = 4 над временем для 1тах = 5 и на 1тах = 5 над 1тах = 6, связаны -о-в
Время генерации ИСД в зависимости от числа частиц для нескольких выбранных максимальных глубин 1тах восьмеричного дерева на CPU (сплошные линии) и на GPU (пунктирные линии) с тем, что с ростом общего числа частиц растет и среднее количество частиц на бокс (т.к. 1тах не меняется), и при достижении на уровне 1тах примерно 256 частиц на бокс число потоков на блок, равное 64, перестает быть оптимальным. Увеличение числа потоков на блок не приводит к существенному улучшению ситуации ввиду специфики архитектуры GPU (например, распределение регистровой и других видов памяти между потоками на одном потоковом мультипроцессоре). В данном случае, для ускорения генерации ИСД, необходимо увеличивать максимальную глубину 1тах восьмеричного дерева.
Из табл. 3.2 видно, что ускорение построения ИСД на GPU по сравнению с CPU растет с увеличением числа частиц. Это говорит о том, что загрузка (утилизация) GPU увеличивается с ростом числа частиц. Использование GPU позволяет ускорить генерацию ИСД до 40 раз для чисел с плавающей точкой одинарной точности. Для чисел с плавающей точкой двойной точности ускорение достигает порядка 30-32 раз. Ускорение генерации ИСД ограничено наличием сложных, плохо объединяемых обращений к памяти. Благодаря использованию предложенного алгоритма, время генерации ИСД для 50 миллионов частиц с использованием чисел одинарной точности на GPU возможно снизить до 1 секунды. Полученные результаты показывают, что, благодаря вышеописанному Таблица 3.2. Время генерации ИСД в зависимости от числа частиц N и максимальной глубины 1тах восьмеричного дерева для чисел с плавающей точкой одинарной точности
Использование иерархической структуры данных для расчета ближнего взаимодействия В работе предлагается использование ИСД не только при реализации быстрого метода мультиполей, но и для ускорения вычисления ближнего взаимодействия по аналогии с расчетом разреженной части в FMM. В случае рассмотрения задачи, в которой участвует как дальнее, так и ближнее взаимодействия, использование одной структуры данных для обоих взаимодействий позволяет существенно сократить объем вычислительных операций.
Ближнее взаимодействие характеризуется тем, что ввиду его затухающей структуры учет его влияния может быть отсечен на каком-то расстоянии от объекта, на который это влияние рассчитывается. Радиус на котором происходит такое отсечение называется радиусом обрезки (rcutoff)- Таким образом, для учета влияния на частицу А необходимо рассмотреть влияние всех частиц, находящихся внутри некоторой сферы с центром в Л и радиусом rcut0ff- Применительно к ИСД это означает, что необходимо рассмотреть взаимодействие на уровне 1тах со всеми частицами, находящимися в боксах, являющихся соседними с боксом, в котором располагается частица А: с учетом проверки попадания в сферу с радиусом rcutoff- От разбиения пространства на боксы (глубины 1тах) зависит сколько соседних боксов попадает в эту сферу. Для рассматриваемой задачи моделирования методами МД наиболее подходящим оказывается разбиение пространства таким образом, что линейный размер бокса (длина ребра) на уровне 1тах будет больше радиуса обрезки rcutoff- Исходя из этого требования, выбирается 1тах (см. раздел 3.2.5). Это необходимо для того, чтобы сфера взаимодействия для рассматриваемой частицы А: находящейся в боксе п, содержалась в области, ограниченной {п тах)- Такой выбор является оптимальным благодаря плотности рассматриваемых систем и архитектуре GPU: в среднем на бокс при разбиении исходя из rcut0ff приходится порядка 20-200 частиц в зависимости от плотности системы.
Благодаря предложенному подходу, расчет ближнего взаимодействия становится подобным расчету разреженного матрично-векторного произведения в FMM (см. раздел 3.4), причем с использованием одной и той же структуры данных. Это позволяет объединить их расчеты при необходимости в одну процедуру и, тем самым, уменьшить число проходов по боксам и обращений к памяти.
В случае исследования систем, в которых рассматривается только ближнее взаимодействие, генерация всей многоуровневой ИСД является излишней ввиду необходимости проведения вычислений по всей иерархической структуре и хранения информации о ней. Для этого случая алгоритм генерации ИСД может быть модифицирован в алгоритм построения одноуровневой ИСД (структуры данных только на уровне 1тах) следующим образом: производится расчет разбиения только на уровне 1тах без расчета всей иерархической структуры и, соответственно, хранится информация только по уровню 1тах- Это позволяет уменьшить размеры используемой для хранения памяти и ускорить алгоритм.
Исследование растекания капли воды по поверхности платины
Расчеты проводились с использованием чисел с плавающей точкой двойной точности. Использовалась NVIDIA CUDA версии 5.0. Все результаты усреднены по нескольким тысячам (2-10 тысяч) расчетных шагов. В табл. 3.9 представлена зависимость времени проведения расчетов от разбиения области на боксы и, следовательно, от среднего числа атомов в боксе. Результаты представлены для числа частиц N = 106. Моделирование проводилось на одном GPU. Размер блока на GPU выставлен в коде равным 64 потокам, исходя из оптимальности. От разбиения области зависят только операции генерации ОСД и расчета ближнего взаимодействия, поэтому представлены времена только по этим операциям и общее время расчета одного временного шага. При малом числе боксов по отношению к числу атомов увеличиваются затраты на генерацию ОСД, что связано со спецификой операций построения гистограммы, например, с необходимостью последовательного обращения (атомарная операция сложения) в участки памяти, соответствующие боксам. Максимальному Таблица 3.10. Время построения ОСД, вычисления ближнего взаимодействия (ЛД) и общее время расчета одного временного шага разработанного на одном GPU алгоритма для числа частиц N = 106 в зависимости от плотности моделируемой системы при одинаковом разбиении на боксы разбиению на боксы при выбранном rcut0ff = 3.5а соответствует разбиение на 34x34x34 бокса. Однако время расчета ближнего взаимодействия для данного числа боксов не является минимальным, так как усредненное число атомов на бокс примерно в два раза меньше выставленного размера блока и даже меньше размера пучка потоков, и GPU задействуется не полностью (плохая заполняе-мость). При уменьшении числа боксов и, следовательно, увеличении среднего числа атомов на бокс увеличиваются затраты на проверку взаимодействия атомов по радиусу обрезки, что приводит к увеличению времени расчета ближнего взаимодействия. Из результатов видно, что наиболее оптимальными по времени расчета ближнего взаимодействия и общему времени расчета является такое разбиение области, при котором усредненное число атомом на бокс не превышает выставленного размера блока, но и не сильно меньше размера блока. При таком разбиении GPU загружается более полно. Поэтому наиболее оптимальным является разбиение области не согласно плотности и радиусу обрезки, а согласно числу частиц на бокс (с учетом того, что размер бокса должен быть не меньше радиуса обрезки). Дальнейшие результаты приводятся для разбиение области согласно оптимальному времени расчета.
В табл. 3.10 показана зависимость времени расчета от плотности моделируемой системы. Результаты представлены для числа частиц N = 106. Моделирование проводилось на одном GPU. Размер блока на GPU выставлен равным 64 потокам, исходя из оптимальности. Как показано выше, выбор разбиения области согласно числу атомов на бокс является более оптимальным. Поэтому время генерации ОСД не зависит от плотности системы в случае, когда общее
Для данного числа атомов разбиение осуществляется на 28x28x28 боксов. Видно, что с увеличением плотности увеличивается и время расчета ближнего взаимодействия, так как все больше атомов оказывается в зоне взаимодействия, ограниченной радиусом обрезки. На такую же величину увеличивается и полное время расчета одного временного шага.
В табл. 3.11 представлена зависимость времени расчета на одном GPU от числа атомов в системе. Время расчета ближнего взаимодействия составляет в среднем примерно 88-90% от общего времени расчета одного временного шага, а время генерации ОСД составляет в среднем примерно 5% в зависимости от разбиения области, т.е. построение ОСД является очень быстрой частью
На рис. 3.22 представлены времена расчетов одного шага по времени в зависимости от числа атомов в системе с использованием разных методов. Из графиков видно, что использование ОСД позволяет уменьшить вычислительную сложность с квадратичной до линейной. Также видна хорошая масштабируемость алгоритма на несколько GPU.
Проведено сравнение производительности с пакетом LAMMPS [77]. Данный симулятор выбран, так как в нём реализован метод, позволяющий снизить вычислительную сложность всего алгоритма (список соседей по ячейкам), и реализована возможность использования GPU. Сравнение было проведено при одинаковых физических параметрах и равномерном распределении молекул по области. Сравнение проводилось на идентичных вычислительных станциях, описанных в начале главы. Получено, что производительность разработанного алгоритма до 3 раз выше, чем производительность LAMMPS на одном GPU (см. рис. 3.22).
В табл. 3.12 показано время расчета ближнего взаимодействия, время коммуникаций между GPU и между CPU и GPU, время расчета одного времен 100 ного шага на нескольких GPU на одном вычислительном узле в зависимости от размера системы. Видно, что расчет сил взаимодействия распараллеливается линейно между GPU. Остальные вычисления производятся только на одном GPU, без распараллеливания между GPU. Время коммуникаций для двух GPU увеличивается почти линейно. Для четырех GPU время коммуникаций всего в 1.5 раз больше, чем для двух GPU (начиная с систем, состоящих из нескольких миллионов частиц), что достигается за счет разработанного алгоритма копирования данных между GPU с использованием асинхронных функций копирования. Также за счет их использования вычисления и коммуникации частично осуществляются параллельно друг другу.
В табл. 3.12 также показано ускорение расчета за счет использования нескольких GPU на одном вычислительном узле по сравнению с использованием одного GPU. Из таблицы видно, что использование нескольких GPU позволяет добиться достаточно хорошего ускорения с учетом того, что примерно 10-12% расчетов (в зависимости от разбиения области на боксы) нераспарал-лелены на несколько GPU и выполняются на мастер-GPU. По закону Амдала, который показывает ограничение роста производительности вычислительной системы с увеличением количества вычислителей, для 10% последовательного кода максимальное распараллеливание, которого можно достичь на 2-х вычислителях, составляет 1.8 раза, на 4-х — 3.08; для 12% максимальное распараллеливание, которого можно достичь на 2-х вычислителях, составляет 1.79 раза, на 4-х — 2.94. Таким образом, ускорение на 2-х GPU достигает почти максимально возможного значения за счет того, что большую часть коммуникаций удается скрыть при помощи использования асинхронных функций копирования. Однако на 4-х GPU не удается скрыть все коммуникации за счет асинхронности, что приводит к недостаточному ускорению.