автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.18, диссертация на тему:Методы ускорения расчетов математических моделей молекулярной динамики на гибридных вычислительных системах

На правах рукописи

МАРЬИН Дмитрий Фагимович

МЕТОДЫ УСКОРЕНИЯ РАСЧЕТОВ МАТЕМАТИЧЕСКИХ МОДЕЛЕЙ МОЛЕКУЛЯРНОЙ ДИНАМИКИ НА ГИБРИДНЫХ ВЫЧИСЛИТЕЛЬНЫХ СИСТЕМАХ

Специальность:

05.13.18 — Математическое моделирование, численные методы и комплексы программ

АВТОРЕФЕРАТ

диссертации на соискание ученой степени кандидата физико-математических наук

Уфа — 2015

- 2ГГН 2015

005561948

Работа выполнена в Центре «Микро- и наномасштабная динамика дисперсных систем» ФГБОУ ВПО «Башкирский государственный университет» и в Федеральном государственном бюджетном учреждении науки Институт механики им P.P. Мавлютова Уфимского научного центра Российской академии наук.

Научный руководитель:

Научный консультант:

Официальные оппоненты:

доктор физико-математических наук Гумеров Наиль Асгатович кандидат физико-математических наук, доцент Михайленко Константин Иванович

доктор физико-математических наук, профессор Карпенко Анатолий Павлович ФГБОУ ВПО «Московский государственный технический университет им. Н.Э. Баумана» заведующий кафедрой систем автоматизированного проектирования

кандидат физико-математических наук, доцент Жданов Эдуард Рифович ФГБОУ ВПО «Башкирский государственный педагогический университет им. М. Акмуллы» декан физико-математического факультета

Федеральное государственное бюджетное образовательное учреждение высшего образования «Московский государственный университет имени М.В. Ломоносова», г. Москва

Защита диссертации состоится «16» сентября 2015 г. в 1400 на заседании диссертационного совета Д-212.288.06 на базе ФГБОУ ВПО «Уфимский государственный авиационный технический университет» по адресу: 450000, г. Уфа, ул. К. Маркса, 12.

С диссертацией можно ознакомиться в библиотеке ФГБОУ ВПО «Уфимский государственный авиационный технический университет» и на сайте www.ugatu.su.

Автореферат разослан «•/-/» ссЯ^и^ТЬ- 2015 года.

Ведущая организация:

Ученый секретарь диссертационного совета, д.ф.-м.н., профессор

Г. Т. Булгакова

Общая характеристика работы

Актуальность работы. Необходимость моделирования динамики N тел возникает во многих областях, например, при расчете молекулярно-динамического, гравитационного взаимодействий, и в ряде методов вычислительной гидродинамики, например, в методах граничных элементов и частиц-в-ячейках. В данной работе рассматривается моделирование методами молекулярной динамики (МД). Вычислительные эксперименты с использованием методов МД позволяют описывать и измерять мельчайшие детали процессов, протекающих в наномасштабах. Однако исследование методами МД реальных физических процессов является достаточно ресурсоемкой задачей, так как при моделировании большого числа частиц и учете многопараметрично-сти и стохастичности серьезно возрастают требования к производительности используемого программного кода и вычислительной системы в целом. Так, вычислительная сложность метода прямого расчета парного взаимодействия N частиц пропорциональна числу всех парных взаимодействий 0(N2). Все это накладывает серьезные ограничения на размеры моделируемых систем. Поэтому важной и актуальной является задача ускорения МД расчетов. Оно может быть достигнуто за счет применения двух способов. Первый способ заключается в разработке и применении методов и алгоритмов, позволяющих снизить вычислительную сложность до 0(N) или до O(NlogN). Второй способ — в повышении производительности вычислений за счет использования высокопроизводительных вычислительных систем. В настоящее время наиболее эффективными и доступными для задач динамики N тел являются гибридные вычислительные системы (ГВС), состоящие из многоядерных центральных и графических процессоров (CPU и GPU, соответственно). Однако эффективное использование ГВС требует разработки новых методов и алгоритмов или значительной модификации существующих.

Степень научной проработанности темы. Проблемам ускорения расчета ближнего взаимодействия в задачах МД посвящены работы авторов Allen, Stone, Verlet и др. Эти работы в значительной мере способствовали эффективному моделированию методами МД, однако, не все предложенные ими методы были оптимизированы под ГВС и они не учитывали проблему оптимизации учета периодических граничных условий. Вопросам ускорения расчета дальнего взаимодействия посвящены работы авторов Barnes, Essmann и др. Однако предложенные ими методы не позволяют уменьшить вычислительную сложность расчетов до линейной с контролируемой точно-

стью, в отличие от быстрого метода мультиполей (FMM), разработанного Greengard и Rokhlin. Проблемам ускорения расчетов при помощи FMM на ГВС посвящены труды авторов Barba, Gumerov, Yokota и др. Однако в них не рассматриваются вопросы применения FMM к задачам МД.

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

Для достижения цели были поставлены следующие основные задачи:

1) Разработка методов и алгоритмов ускорения моделирования взаимодействия частиц на гибридных вычислительных системах.

2) Реализация комплекса программ для моделирования методами МД с использованием высокоэффективных методов и гибридных вычислительных систем. Оценка эффективности предлагаемых методов и алгоритмов. Валидация комплекса программ на тестовых задачах.

3) Разработка математического метода моделирования процесса гетерогенной кавитации в неполярной жидкости на инородном включении.

4) Моделирование ряда прикладных задач: исследование влияния величины радиуса обрезки потенциала ближнего взаимодействия на рассчитываемые макропараметры; растекание капли воды по подложке; гетерогенная кавитация на инородном включении в неполярной жидкости; образование пузырька вблизи подложки в бинарной одноатомной смеси.

Основные положения, выносимые на защиту:

1) Метод ускорения моделирования неполярных систем на основе предложенной одноуровневой структуры данных. Оригинальные алгоритмы применения FMM, построения иерархической и одноуровневой структур данных, обладающие линейной вычислительной сложностью и ориентированные на моделирование на ГВС методами МД жестких молекул, взаимодействующих согласно потенциалам Леннард-Джонса и Кулона.

2) Программный комплекс для моделирования методами МД, созданный на основе предложенных методов и ориентированный на ГВС.

3) Математический метод моделирования процесса гетерогенной кавитации в системе неполярных молекул на инородном включении с различной лиофильностью.

4) Анализ влияния радиуса обрезки потенциала Леннард-Джонса на рассчитываемые макропараметры, влияния размеров системы в нанометро-вом диапазоне на величину давления на разрыв в процессе гетерогенной кавитации на инородном включении и на величину контактного угла в задаче растекания капли воды по подложке.

Научная новизна работы заключается в следующем:

1) Разработан уникальный метод ускорения моделирования взаимодействия с усеченным потенциалом на основе предложенной одноуровневой структуры данных, который позволяет избавиться от проверки периодических граничных условий и имеет линейную вычислительную сложность. Предложены новые вычислительные алгоритмы построения иерархической структуры данных и применения FMM на ГВС, направленные на моделирование методам МД полярных систем, взаимодействующих согласно потенциалам Леннард-Джонса и Кулона.

2) Создал оригинальный модульный программный комплекс для моделирования методами МД на ГВС, основанный на использовании FMM, иерархической и одноуровневой структур данных.

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

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

Теоретическая и практическая значимость. Разработанные методы и гетерогенные алгоритмы имеют существенное значение для методов МД, астрофизики, физики плазмы, вычислительной гидродинамики и других областей. Благодаря гибкой модульной структуре, реализованный в работе программный комплекс применим для решения задач из тех областей, где возможно использование FMM, иерархической и одноуровневой структур данных. Разработанные алгоритмы могут быть использованы не только в применении к GPU, но и к другим многоядерным системам.

Достоверность результатов обусловлена применением математически обоснованных методов, сравнением с известными результатами теоретических, экспериментальных и численных исследований.

Апробация работы. Основные результаты, представленные в работе, докладывались на: международных конференциях «Параллельные вычислительные технологии (ПаВТ)» (Новосибирск, 2012 и Ростов-на-Дону, 2014); V Российской конференции с международным участием «Многофазные системы: теория и приложения» (Уфа, 2012); международной школе-конференции «Динамика дисперсных систем: экспериментальные и численные исследования на нано-, микро-, мезо- и макромасштабах» (Уфа, 2014); международной конференции «Science of the Future» (Санкт-Петербург, 2014); международных школах-конференциях «Фундаментальная математика и ее приложения в естествознании» (Уфа, 2011 и 2012); XIV и XV Всероссийских конференциях-школах молодых исследователей «Современные проблемы математического моделирования» (Абрау-Дюрсо, 2011 и 2013); VI Всероссийской конференции «Актуальные проблемы прикладной математики и механики» (Абрау-Дюрсо, 2012); международных конгрессах «ASME International Mechanical Engineering Congress & Exposition» (США, Хьюстон, 2012; США, Сан-Диего, 2013; Канада, Монреаль, 2014); международной конференции «International Conference on Numerical Methods in Multiphase Flows» (США, Стейт Колледж, 2012); семинарах в БашГУ, Имех УНЦ РАН, Мэрилендском университете, УГАТУ и МГУ.

Исследования, результаты которых представлены в диссертации, проводились при поддержке грантов Министерства образования и науки Российской Федерации (11.G34.31.0040), РФФИ (12-01-31083-мол_а) и У.М.Н.И.К.

Автор благодарит своего научного руководителя H.A. Гумерова, научного консультанта К.И. Михайленко, коллектив Центра в лице И.Ш. Ахатова, В.Л. Малышева, Е.Ф. Моисеевой, а также С.Ф. Урманчеева за помощь при подготовке диссертации.

Публикации. Основные результаты по теме диссертации опубликованы в 19 научных работах, из них 6 в рецензируемых журналах, рекомендованных ВАК. Получено 2 свидетельства о государственной регистрации программ для ЭВМ.

Объем и структура работы. Диссертация состоит из введения, четырех глав, приложения, заключения и списка литературы. Полный объем диссертации изложен на 146 страницах и содержит 42 рисунка и 18 таблиц.

Список литературы содержит 127 наименований.

4

Краткое содержание работы

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

В первой главе представлен обзор существующих подходов к ускорению моделирования задач динамики N тел, в частности, задач МД: методы алгоритмического ускорения, направленные на уменьшение вычислительной сложности расчетов ближнего и дальнего взаимодействий; современные подходы к аппаратному ускорению, в особенности, с использованием ГВС.

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

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

П = г = МУ, (1)

М М ГП{

где N — число тел; г*, V* — положение и скорость центра масс г-го тела; тщ — его масса; Г; — сила, действующая на г-ое тело; Щ — суммарный потенциал взаимодействия г-го тела со всеми остальными, вычисляемый по формуле = и(Гц), где г,у = |г,—г_,|. Функция С/(гу) определяет физические

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

(РЧг 1

¿¡5-= 2^" ' (2)

^ = I-1 (М< - о,, х (Ь,,)) , («,, 0)г = 20^,

где q¡ — кватернион г-го тела; <3; — его матричное представление; ъг{ — четырехкомпонентный вектор, введенный для удобства обозначений; —

вектор угловой скорости г го тела: I тензор инерции; М* момент сил.

5

Таким образом, система из N тел в ¿-мерном случае описывается при помощи системы из 2йИ обыкновенных дифференциальных уравнений 1-го порядка (1). При учете вращательного движения число уравнений увеличивается вдвое. Уравнения (1)-(2) решаются численно согласно выбранной численной схеме (в работе используется метод скоростей Верле). Шаг моделирования составляет порядка 1 фс (1 • Ю-15 с).

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

Рассматривается канонический Л^Т-ансамбль (фиксированы число молекул N, объем V и температура Т). Начальные положения молекул и их скорости выбираются в соответствии с геометрией области, плотностью и температурой системы. Температура поддерживается постоянной при помощи термостатов Нозе-Гувера или Берендсена. В качестве граничных условий (ГУ) рассматриваются отражение и периодические ГУ, с помощью последних возможно описать квази-бесконечное пространство.

Потенциалы взаимодействия можно разделить на близкодействующие (на атом существенно влияние лишь его ближайшего окружения) и дальнодействующие. Для моделирования парного взаимодействия неполярных молекул широко применяется близкодействующий потенциал Леннард— Джонса. Для пары атомов в точках с координатами г, и г;-, он равен: = 4е ((с/гу)12 — (ег/гу)6^, где гу = |г; —е—-глубина потенциальной ямы; а — расстояние, на котором энергия взаимодействия становится равной нулю. Параметры е и ег являются характеристиками молекул. При моделировании многокомпонентных систем параметры е и а взаимодействия атомов разных типов рассчитываются согласно правилу Лоренца-Бертло. Для ускорения расчетов взаимодействие согласно потенциалу Леннард-Джонса обрывается на некотором расстоянии гс, называемым радиусом обрезки, и смещается на величину {/и(гс), чтобы избежать нарушения непрерывности:

тт г \ \иьМ])-иь]{гс), Гу^Гс,

ии^лтц) = <

[О, Гу>Гс.

Для моделирования полярных веществ широко используется комбинация потенциалов Леннард-Джонса и Кулона (дальнодействующий потенциал). Потенциал Кулона взаимодействия пары атомов с зарядами Ц), и щ и координатами г; и Г] описывается формулой: £7с(гу) = где гу = |г,—г7|.

Взаимодействие полярных молекул рассматривается на примере молекул воды модели TIP4P, которая имеет плоскую конфигурацию и содержит 4 жестко связанных частицы. Взаимодействие между двумя молекулами воды представляет собой комбинацию парных взаимодействий частиц этих молекул: три заряженные частицы взаимодействуют согласно потенциалу Кулона, дополнительная частица — согласно потенциалу Леннард-Джонса.

Для моделирования твердых включений и подложки используется модель платины, атомы которой располагаются согласно кристаллической решетке fee (111). Взаимодействие атомов платины друг с другом и с неполярными молекулами осуществляется согласно потенциалу Леннард-Джонса. Для моделирования взаимодействия воды и платины используется потенциал, предложенный Zhu S.-B. и Philpott M.R..

Процесс кавитации (образование пузырька в объеме жидкости) возникает в ряде областей гидродинамики и имеет важное значение в различных технологических процессах. Кавитация происходит при пониженном давлении в системе, находящейся в метастабильном или нестабильном состоянии. Одним из основных методов моделирования процесса кавитации является моделирование растянутой жидкости (метастабильное состояние). Исследованию гомогенной кавитации в растянутой леннард-джонсовской жидкости посвящены работы Нормана Г.Э., Куксина А.Ю. и др. Существенный интерес представляет изучение процесса кавитации не в чистой жидкости, а в присутствии инородных включений (гетерогенная кавитация). Кавитация на гидрофобных частицах в растянутой жидкости, находящейся в свободном объеме, исследовалась Koishi Т. и др. Однако кавитация в метастабильном состоянии носит случайный характер, и время, необходимое для прослеживания зарождения пузырька, может значительно превышать возможные времена моделирования методами МД. Поэтому для моделирования гомогенной кавитации методами МД рядом исследователей (Байдаков В.Г., Проценко С.П., Kinjo Т., Matsumoto М.) применяется техника постепенного понижения в системе давления за счет растяжения области моделирования. В работе предложен метод моделирования гетерогенной кавитации в неполярной жидкости на инородном включении с различной лиофильностью, который основан на последнем подходе. Область заполняется молекулами моделируемой жидкости. В центр вместо молекул жидкости помещается твердое включение, вырезанное из кристаллической решетки соответствующего вещества. После достижения системой равновесия давление медленно понижается за счет постепенного растяжения области, при этом положения атомов жидкости смещаются

7

соответственно, атомы твердого включения не сдвигаются. Таким образом, система переводится из стабильного состояния в метастабильное ближе к спи-нодали, что позволяет инициировать процесс кавитации.

Смачиваемость включения и подложки моделируется путем изменения леннард-джонсовского параметра энергии взаимодействия атомов жидкости с атомами включения/подложки esi = ет£ц согласно параметру лиофильно-сти ег (с ростом £г лиофильность увеличивается), где г и — параметр потенциала Леннард-Джонса взаимодействия жидкость-жидкость.

В третьей главе предлагаются методы численного моделирования, позволяющие ускорить расчеты задач N тел и, в частности, МД. Приведено описание предлагаемых структур данных и FMM, а также разработанных алгоритмов.

Вычисления производились с использованием гибридной рабочей станции, оборудованной двумя 6-ядерными CPU Intel Xeon 5660 2.8 ГГц, 32 ГБ RAM и четырьмя GPU NVIDIA Tesla С2075 (6 ГБ RAM и 448 ядер каждая). Для программирования на GPU использовалась технология CUDA.

В работе рассмотрена иерархическая структура данных (ИСД) и предложен метод ее построения, оптимизированный под GPU и обладающий линейной вычислительной сложностью. В ИСД область разбивается с использованием 2d-flepeBa (восьмеричного дерева) иерархически до уровня lmax (глубина дерева). В работе рассмотрен вопрос выбора глубины восьмеричного дерева. Алгоритм генерации ИСД на GPU основан на использовании мортоновской кривой заполняющей пространство для быстрого прохода по 2й-дереву, гистограммы размещения частиц по боксам, бинарной сортировке и параллельном сканировании. Взаимодействие частиц рассчитывается с использованием вспомогательного массива «маркеров», который указывает, где в памяти хранится первая частица, находящаяся в боксе. Разработанный алгоритм позволил достичь ускорения построения ИСД на GPU по сравнению с CPU до 40 раз для чисел с плавающей точкой одинарной точности и до 32 раз для чисел двойной точности. Показана хорошая масштабйруемость.

Предложен метод использования ИСД для ускорения расчета усеченного взаимодействия на ГВС на примере расчета потенциала Леннард-Джонса. На рис. 1 показано время вычисления ближнего взаимодействия на основе потенциала Леннард- Джонса для чисел с плавающей точкой одинарной точности. В ИСД выбирался оптимальный lmax свой для CPU и для GPU. Моделировалась система жидкого аргона плотностью р =

= 1000 кг/м3 при гс = 5а. Получено, что применение GPU позволяет

8

ускорить вычисление методом прямого суммирования (проверка взаимодействия каждой частицы с каждой) до 520 раз при использовании чисел с плавающей точкой одинарной точности и до 280 раз при использовании двойной точности по сравнению с однопоточным кодом на CPU. Асимптотика (рис. 1) под- 1Г

X t-bN2.-" > Ч-9 t-cH

10 10° number of atoms

тверждает, что использование ИСД уменьшает вычислительную сложность алгоритма расчета усеченного взаимодействия с квадратичной до линейной, что является очень важным результатом

при моделировании больших рис Время вычисления ближнего взаимодействия систем. Так, например, для на одном временном шаге методом прямого сумми-125000 атомов время расчета роваш1Я (_.) „ с использованием ИСД (HDS) ( ) в одного временного шага на зависимости от числа частиц CPU без распараллеливания и

без использования ИСД составляет порядка 600 с. Использование ИСД на CPU позволяет уменьшить это время до 0.25 с. Из рис. 1 видна хорошая масштабируемость предложенного метода. При использовании одинакового 1тах для GPU ц для CPU достигается ускорение порядка тысячи раз. При выборе оптимального 1тах для каждой из архитектур ускорение ограничено ускорением ИСД.

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

Предложены метод использования FMM+ИСД для моделирования систем полярных молекул взаимодействующих согласно потенциалам Кулона и Леннард-Джонса и алгоритм применения FMM+ИСД на ГВС (рис. 2). Разработанный алгоритм имеет ряд

GPU

Генерация ИСД к

| Решатель ОДУ

Разложение помупь-типспьному базису

Вычисление

усеченного взаимодействия и ближнего поля дали 'С го взаимодействия

Вычисление локального разложения. Суммирование цикл повремени I

Г

CPU

Пфаллельно _jraGPUnCPU

Трансляционная структура данных

Трансляции коэффициентов разложения из мультмпольного базиса в локальный

Параллельно на GPU и CPU

Рис. 2. Блок-схема гетерогенного алгоритма FMM 9

преимуществ: CPU не простаивает пока на GPU производятся вычисления; CPU и GPU загружены работой, которую они могут выполнять наиболее эффективно; их загрузка сбалансирована; обмен данными между ними минимизирован; копирование данных осуществляется асинхронно.

На рис. 3 показано время вычисления одного временного шага для системы молекул воды для чисел с плавающей точкой двойной точности. Видно что время расчета методом прямого суммирования растет пропорционально квадрату числа частиц на CPU и на GPU. Применение FMM и ИСД позволяют получить линейную масштабируемость алгоритма в зависимости от числа частиц как на CPU, так и на GPU. Оценка эффективности моделирования систем полярных молекул на одном GPU показала, что разработанный программный комплекс в 3-4 раза быстрее, чем симулятор LAMMPS с использованием метода Р3М (Particle-Particle Particle-Mesh) для ускорения расчетов. Одним из существенных недостатков LAMMPS является ограничение на число моделируемых атомов из-за размеров памяти GPU. Разработанный метод позволяет моделировать системы, содержащие в 10 раз больше атомов: например, на GPU NVIDIA Tesla С2050 с 3 ГБ RAM вплоть до 12 • 106 атомов.

Предложены одноуровневая структура данных (ОСД) и метод ее применения для уменьшения вычислительной сложности расчета усеченного (ближнего) взаимодействия. Разработан алгоритм ее построения на GPU, обладающий линейной вычислительной Рис. 3. Время вычисления одного временного шага в сложностью. Область раз- зависимости от числа атомов на CPU (распараллеле-бивается на боксы, размер н0 при помощи ОрепМР на 24 потока) и на GPU с наименьшей стороны которых превосходит гс. Информация по частицам располагается в памяти согласно сквозному проходу по боксам. Это позволяет сгруппировать обращения к памяти при расчете взаимодействий и в других процедурах. Для построения ОСД используется псевдосортировка атомов, основанная на блочной сортировке, как и в случае построения ИСД. Взаимодействие атомов рассчитывается с использованием массива «маркеров» (аналогично с ИСД). За счет удобного прохода

10

использованием метода прямого суммирования (BF) и с использованием FMM+HQII,(HDS)

по боксам ОСД позволяет использовать при вычислении взаимодействия вместо глобальных координат частиц их локальные координаты (относительно центра бокса, в котором они располагаются). Такой подход позволяет избавиться от проверки периодичности при расчете взаимодействия между частицами по кратчайшему радиусу и, следовательно, избежать лишних ветвлений алгоритма. В работе исследовано влияние разбиения области на боксы на производительность. Для ускорения моделирования систем с усеченным взаимодействием разработан алгоритм, позволяющий проводить расчеты на нескольких GPU на каждом вычислительном узле (multi-GPU) при помощи технологии ОрепМР и с использованием нескольких вычислительных узлов, каждая из которых может содержать несколько GPU, при помощи технологии MPI. Из рис. 4 видно, что метод с использованием ОСД позволяет уменьшить вычислительную сложность расчета с квадратичной до линейной, применение GPU позволяет значительно ускорить расчеты. Также видна хорошая линейная масштабируемость алгоритма как для одного, так и для нескольких GPU внутри одного узла. Получено ускорение на одном узле на двух GPU в 1.8 раза, на четырех — в 2.8 раза. Оценка эффективности показала, что на 1 GPU разработанный программный комплекс до 3 раз быстрее, чем симулятор LAMMPS с использованием списка соседей по ячейкам для ускорения расчетов.

Разработанный программный комплекс включает следующие модули: МД процедур (генерация начальных данных, расчет взаимодействия с использованием для ускорения нижеперечисленных модулей, численная ioJ ю" схема, термостатирование,

ю 10:

и 10 И

¡10' Е

" 10 10 10

t-btf* 2 А V t-aN! * J' у * A BF, CPU V ВГ, GPU -----LAÏMPS —Q-SDS, 1 GPU -в— SOS, 2 GPU -1-SDS, 4 GPU

У Г • t-CN * У V У ? -¿¿Ф^ У v.'' " V

10 10° number of atoms

10'

шага с использованием ОСД (SDS), пакета LAMMPS на 1 GPU, метода прямого суммирования (BF) на CPU и на GPU

граничные условия, расчет Рис. 4. Время р^чета одного временного макропараметров и др.); процедур построения ИСД и ОСД; процедур FMM; функций-оболочек для функций библиотеки CUDA для удобного доступа из других модулей; пост-обработки данных численных экспериментов; программы по моделированию методами МД и обеспечения коммуникаций между CPU, GPU и вычислительными узлами с использовав

нием функций вышеперечисленных модулей. Все модули, кроме последнего,

11

реализованы в виде отдельных библиотек и могут быть применены для решения других задач. Модули написаны с использованием языков программирования C/C++, Fortran, Matlab, технологий CUDA, ОрепМР и MPI, и реализованы на CPU и на GPU.

В четвертой главе приведены результаты тестирования разработанного программного комплекса: показано удовлетворительное согласование результатов моделирования тестовых задач с теоретическими и экспериментальными результатами других исследователей. Также в этой главе приведены результаты моделирования ряда практических задач МД.

Исследовано влияния радиуса обрезки гс потенциала Леннард-Джонса на точность моделирования для задачи определения толщины межфазного слоя на границе пузырька в жидком аргоне. Разработан метод численного определения эффективной толщины межфазного слоя для пузырька. С увеличением гс число рассчитываемых взаимодействий растет кубически, поэтому его выбор сильно влияет на скорость моделирования. Получено, что малый радиус обрезки величиной 2.5а + 4а (0.85 + 1.4 им) вносит существенную ошибку в оценку ряда макропараметров (вплоть до 30-40% для случая 0.85 нм). Например, малые значения гс приводят к недооценке плотности жидкой фазы до 10%, переоценке плотности пара в пузырьке и эффективной толщины межфазного слоя, и неправильной оценке его местоположения. Выбор гс, начиная примерно с 5а 1.7 нм), позволяет получить корректные оценки макропараметров.

30,-■-т-

Проведено исследование влия- ш25 | , j £ mSsLuL?iôn°°5:

ния размера системы, как ограни- его т

чивающего параметра математиче- и15

га 10

ской модели, на результат определения контактного угла между каплей воды и гидрофильной подложкой. Моделировалась система, содержащая 5832 молекулы воды и 5184 молекул платины. Получено удовле-

о 5-

300 350 400 450 500 Т, К

Рис. 5. Зависимость контактного угла на границе вода — платина от температуры

творительное согласование (рис. 5) с работой Shi В., Sinha S., Dhir V.K., в которой проводилось подобное численное моделирование, но для в разы меньшей системы и другой модели (1690 молекул воды модели SPC/E). Установлено, что с ростом размеров капли в нанометровом диапазоне величина контактного угла не претерпевает существенных изменений.

Таблица 1. Зависимость давления на разрыв от размеров системы: Ь длина стороны кубической области; N — число атомов; р — усредненное давление; в — стандартное отклонение

Размер системы Ь =14.6 нм (Я = 64016) Ь =21.9 нм (Я = 216016) Ь =29.2 нм ^ = 512016) Ь =36.5 нм № = 1000016)

р ± в, атм -351±11.2 -355±11.1 -355±7.9 -355±4.5

Проведено моделирование гетерогенной кавитации в системе с твердым включением с использованием предложенного в работе метода. Исследовано влияние размера моделируемой системы, как ограничивающего параметра математической модели, на величину давления, при котором происходит нуклеация. Рассматривалась кубическая область, заполненная жидким аргоном (р = 1350 кг/м3, Т = 85 К), в центр которой помещено твердое лиофоб-ное шарообразное включение радиусом 5 А, вырезанное из кристаллической решетки платины. В табл. 1 показана зависимость усредненного по десяти численным расчетам давления на разрыв р (максимальное по модулю отрицательное давление в системе, при котором происходит разрыв и образуется зародыш пузырька) от размеров области. Получено, что с изменением размеров области в нанометровом масштабе среднее давление на разрыв меняется слабо. Таким образом, для исследования давления на разрыв возможно использовать предложенный подход и рассматривать системы, состоящие не из сотен тысяч атомов, а меньшие, что позволяет уменьшить ресурсопотребление.

Проведено моделирование задачи образования пузырька в газожидкостной смеси (неон+аргон) вблизи подложки с плотностью смеси р = 1300 кг/м3 и температурой среды Т = 85 К в зависимости от различных модельных параметров: размеров области, концентрации растворенного газа, лиофильности подложки £г и радиуса обрезки гс потенциала Леннард-Джонса. Установлено, что с ростом концентрации растворенного газа зарождение пузырька происходит раньше, а с увеличением лиофильности подложки — позже (табл. 2). Обнаружено, что с увеличением гс уменьшается время до образования пузырька на подложке и порог по параметру лиофильности ег, при котором

13

Таблица 2. Время, когда на поверхности образуется пузырек, в зависимости от ет для случаев с различными гс (гс = 2 нм, Ь = 11 нм; гс = 3 нм, Ь = 12 нм)_

£т гс = 2 нм гс = 3 нм

0.350 0.47+0.53 не 0.05+0.20 не

0.408 0.9+1.0 не 0.45+0.55 не

0.490 10.3+10.5 не 1.7+1.8 не

пузырек образуется именно на подложке, а не в объеме жидкости, сдвигается вверх.

В приложении приведены алгоритмы генерации ИСД и ОСД. В заключении сформулированы основные результаты и выводы:

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

2) Разработан модульный программный комплекс для моделирования методами МД жестких молекул, взаимодействующих согласно потенциалам Леннард-Джонса и Кулона, на гибридных вычислительных системах. Показана линейная масштабируемость разработанных алгоритмов до десятков миллионов частиц. Показано, что разработанный программный комплекс в 3-4 раза быстрее по сравнению с симулятором ЬАММРБ.

3) Предложен метод моделирования процесса гетерогенной кавитации в неполярной жидкости на инородном включении с различной лиофиль-ностью.

4) Проведены комплексные исследования на основе математического моделирования ряда практических задач с использованием разработанного программного комплекса.

• Исследование влияния радиуса обрезки в математической модели ближнего взаимодействия показало, что он оказывает существенное влияние на исследуемые макропараметры. Получение корректных результатов возможно лишь при условии использования радиуса обрезки больше 1.7 нм. Установлено, что увеличение радиуса обрезки ближнего взаимодействия приводит к уменьшению времени до образования пузырька на подложке.

• Показано, что в задаче растекания капли воды по подложке величины контактных углов не зависят от размеров системы в нанометро-вом диапазоне.

• Выявлено, что при исследовании гетерогенной кавитации в системе жидкого аргона с инородным включением изменение размеров си-

стемы в диапазоне десятков нанометров оказывает несущественное влияние (менее 1%) на величину давления на разрыв.

Список основных публикаций по теме диссертации Публикации в рецензируемых журналах, рекомендованных ВАК

1) Марьин Д.Ф., Малышев В.Л., Моисеева Е.Ф., Гумеров Н.А., Ахатов И.Ш., Михайленко К.И. Ускорение молекулярно-динамических расчетов с помощью быстрого метода мультиполей и графических процессоров // Вычислительные методы и программирование. — 2013. — Т. 14. — С. 483-495.

2) Малышев В.Л., Марьин Д.Ф., Моисеева Е.Ф., Гумеров Н.А., Ахатов И.Ш. Ускорение молекулярно-динамического моделирования неполярных молекул при помощи GPU // Вестник Нижегородского университета им. Н.И. Лобачевского. - 2014. - Вып. 3. - С. 126-133.

3) Моисеева Е.Ф., Марьин Д.Ф., Малышев В.Л., Гумеров Н.А., Ахатов И.Ш. Исследование контактного угла и объема поверхностного нанопу-зырька методами молекулярной динамики // Математическое моделирование. - 2015. - Т. 27, № 4. - С. 115-126.

4) Малышев В.Л., Марьин Д.Ф., Моисеева Е.Ф., Гумеров Н.А., Ахатов И.Ш. Исследование прочности жидкости на разрыв методами молекулярной динамики // Теплофизика высоких температур. — 2015. — Т. 53, № 3. - С. 423-429.

5) Moiseeva E.F., Mikhaylenko СЛ., Malyshev V.L., Maryin D.F., Gumerov N.A. FMM/GPU accelerated molecular dynamics simulation of phase transitions in water-nitrogen-metal systems // Proc. of the «ASME 2012 International Mechanical Engineering Congress & Exposition», November 9-15, 2012. - Houston, USA. - 10 p. - Paper No. IMECE2012-86246.

6) Moiseeva E.F., Malyshev V.L., Marin D.F., Gumerov N.A., Akhatov I.Sh. Molecular dynamics simulations of nanobubbles formation near the substrate in a liquid with dissolved gas // Proceedings of the «ASME 2014 International Mechanical Engineering Congress & Exposition», November 14 20, 2014. Montreal, Canada. 8 p. Paper No. IMECE2014-37050.

Свидетельства о регистрации программы для ЭВМ

7) Марьин Д.Ф., Малышев В.Л., Моисеева Е.Ф., Михайленко К.И., Гумеров Н.А. MDS-W высокопроизводительная библиотека для молекулярно-динамического моделирования воды / 2013 / Свидетельство о государственной регистрации программы для ЭВМ. Per. №2013612088

15

8) Марьин Д.Ф., Малышев В.Л., Моисеева Е.Ф., Гумеров H.A. MDS-А — молекулярно-динамическое моделирование неполярных одноатомных молекул / 2014 / Свидетельство о государственной регистрации программы для ЭВМ. Per. №2014611173

Основные публикации в других изданиях

9) Марьин Д.Ф. Использование графических процессоров при решении задач молекулярной динамики // Труды Института механики УНЦ РАН. Вып. 8 / Под ред. С.Ф. Урманчеева — Уфа: Нефтегазовое дело, 2011 — С. 182 188.

10) Марьин Д.Ф. Ускорение расчета процессов молекулярной динамики при помощи GPU // Параллельные вычислительные технологии (ПаВТ'2012): Труды международной научной конференции [Эл. ресурс] — Челябинск: Издательский центр ЮУрГУ, 2012. - С. 606-611.

11) Марьин Д.Ф. Ускорение молекулярно-динамического моделирования многофазных систем при помощи GPU // Труды Института механики УНЦ РАН. Вып. 9, Ч. 2 / Под ред. С.Ф. Урманчеева - Уфа: Нефтегазовое дело, 2012. - С. 76-79.

12) Малышев В.Л., Марьин Д.Ф., Моисеева Е.Ф. Новая структура данных для расчета ближнего взаимодействия в методах молекулярной динамики // Сб. трудов XV Всероссийской конференции-школы молодых исследователей. — Ростов-на-Дону: Изд-во ЮФУ, 2013. — С. 155-159.

13) Малышев В.Л., Марьин Д.Ф., Моисеева Е.Ф., Гумеров H.A., Ахатов И.Ш. Определение поверхностного натяжения методами молекулярной динамики для одноатомных веществ // Труды Института механики УНЦ РАН. Вып. 10 / Под ред. С.Ф. Урманчеева — Уфа: Нефтегазовое дело, 2013 - 5 С.

14) Малышев В.Л., Марьин Д.Ф., Моисеева Е.Ф., Гумеров H.A., Ахатов И.Ш. Ускорение молекулярно-динамического моделирования неполярных молекул при помощи GPU // Параллельные вычислительные технологии (ПаВТ'2014): Труды международной научной конференции [Эл. ресурс] — Челябинск: Издательский центр ЮУрГУ, 2014. — С. 140-149.

Диссертант

Д. Ф. Марьин

МАРЬИН Дмитрий Фагимович

МЕТОДЫ УСКОРЕНИЯ РАСЧЕТОВ МАТЕМАТИЧЕСКИХ МОДЕЛЕЙ МОЛЕКУЛЯРНОЙ ДИНАМИКИ НА ГИБРИДНЫХ ВЫЧИСЛИТЕЛЬНЫХ СИСТЕМАХ

Специальность: 05.13.18 - Математическое моделирование, численные методы и комплексы программ

АВТОРЕФЕРАТ диссертации на соискание ученой степени кандидата физико-математических наук

Лицензия на издательскую деятельность ЛР № 021319 от 05.01.99 г.

Подписано в печать 21.05.2015 г. Формат 60x84/16. Усл. печ. л. 1,04. Уч.-изд. л. 1,08. Тираж 100 экз. Заказ 206.

Редакционно-издательский центр Башкирского государственного университета 450074, РБ, г. Уфа, ул. Заки Валиди, 32.

Отпечатано на множительном участке Башкирского государственного университета 450074, РБ, г. Уфа, ул. Заки Валиди, 32.