Комментарии 29
Офигенно, спасибо огромное. Те посты, за которые ценится Хабр.
статья классная!
if (blockIdx.x == 0) { md->engElecField = 0; } ... if (!threadIdx.x) // 0th thread { atomicAdd(&md->engElecField, shEngElField); }
Так делать нельзя. Синхронизации между блоками не предусмотрено, и никто не гарантирует, что первый блок выполнится первым. Соответственно обнулиться может уже после того, как какие то другие блоки добавят свое значение. Соответственно сумма скорее всего вычисляется неправильно. Во втором проходе то же самое.
Обнулять нужно в отдельном kernel ну или через cudaMemset() заранее.
В put_periodic():
atomicAdd(&md->specAcBox[type].y, 1);y — это описка? Везде вроде x.
Входной параметр step равен количеству частиц (атомов) обрабатываемых одним потоком.
Очень плохая идея. Если step > 1, получается плохой патерн доступа к памяти в warp'е — потоки обращаются к памяти не последовательно, а с шагом в step, т.е. warp будет читать не один кусок памяти целиком, а несколько маленьких кусочков — по одному на каждый поток. Это может замедлять выполнение в несколько раз.
Ну и это только то, что заметил в первой части :)
Спасибо огромное за комментарий.
согласен, не подумал об этом, поправил в кодах, добавил UPD
спасибо, поправил
Это не было опечаткой, x подразумевался для частиц летящих назад, в отрицательную область по x, а y — для положительного направления. Но поскольку там 3 измерения, пришлось всё равно переделать, выделить разные массивы для разных направлений. Сейчас все немного по-другому
Так делать нельзя… обнулиться может уже после того, как какие то другие блоки добавят свое значение...
согласен, не подумал об этом, поправил в кодах, добавил UPD
Очень плохая идея. Если step > 1, получается плохой патерн доступа к памяти в warp'е — потоки обращаются к памяти не последовательно, а с шагом в step
спасибо, поправил
y — это описка? Везде вроде x.
Это не было опечаткой, x подразумевался для частиц летящих назад, в отрицательную область по x, а y — для положительного направления. Но поскольку там 3 измерения, пришлось всё равно переделать, выделить разные массивы для разных направлений. Сейчас все немного по-другому
Так делать нельзя… обнулиться может уже после того, как какие то другие блоки добавят свое значение…
согласен, не подумал об этом, поправил в кодах, добавил UPD
Все равно остались такие же места в temp_scale() и tstat_nose(), только они в конце, а не в начале:
if (!blockIdx.x)
if (!threadIdx.x)
{
md->engKin = md->teKin;
}
Также никто не гарантирует, что этот код в блоке 0 выполнится после всех остальных блоков, скорее даже наоборот, он выполнится до того, как большинство других блоков начнет исполняться и там вычислится:
k = sqrt(md->teKin / md->engKin);Вообще не стоит ориентироваться на номер блока и выполнять в зависимости от него какой то код — не могу даже придумать случая, где это не было бы ошибкой.
По поводу оптимизаций можно много всего написать. Но мне кажется, здесь имеет смысл оптимизировать только вычисление сил, т.к. оно должно занимать практически все время вычислений. И здесь бы я сделал все таки обработку каждой частицы в отдельном потоке, а ячейки использовал только для поиска соседей. Т.е. в каждом потоке обрабатывается одна частица, но для нее перебираются не все соседи, а только частицы в той ячейке, в которой она находится и в соседних. И избавился бы от atomic'ов — для этого сохранял бы силы только для обрабатываемой частицы — да в этом случае для каждой пары частиц компонент сил будет вычисляться дважды, но это ничто по сравнению с избавлением от atomic'ов. Вычислительные операции на gpu это практически самая дешевая операция — лучше что-то лишний раз посчитать, чем читать из глобальной памяти, а тем более с использованием atomic'ов — одной из наиболее тяжелых операций.
Все равно остались такие же места в temp_scale() и tstat_nose()/поправил, в temp_scale убрал только в начале, в конце — неважно, какой из блоков это сделает.
не все соседи, а только частицы в той ячейке, в которой она находится и в соседних
при использовании алгоритмов cell list так и происходит
да в этом случае для каждой пары частиц компонент сил будет вычисляться дважды, но это ничто по сравнению с избавлением от atomic'ов
думаете, что два раза считать корни, экспоненты и т.д. быстрее, чем один раз, но с атомик?
поправил, в temp_scale убрал только в начале, в конце — неважно, какой из блоков это сделает.
Я бы все таки перенес это присваивание в другое место, например туда, где вычисляется и присваивается md->k. Просто чтобы код был понятнее, но это уже, понятное дело, относится к стилю, а не к корректности кода.
думаете, что два раза считать корни, экспоненты и т.д. быстрее, чем один раз, но с атомик?
По моему опыту — да, гораздо быстрее. Как я уже писал, чистые вычисления на gpu в большинстве случаев не являются узким местом. В большинстве случаев все упирается в доступ к глобальной памяти, ну а глобальные атомики — это еще медленнее.
Как следующий этап оптимизации я хотел предложить не просто обработку каждой частицы в отдельном потоке, но обработку их в порядке их расположения в ячейках. Т.е. частицы, расположенные в одной ячейке должны обрабатываться подряд — в идеале потоки каждого warp'а или даже блока обрабатывают частицы из одной ячейки. Тогда потоки каждого warp'а будут синхронно обращаться к одним и тем же соседям, а не разным и случайным, что должно давать очень существенное ускорение. Именно про это и написано в презентации, ссылку на которую дал pavel_kudinov. И там уже рассматривается как раз способы сортировки частиц, чтобы они шли в порядке принадлежности к ячейкам. И как раз уже время выполнения этой задачи становится основной проблемой.
Да, я посмотрел. Попробую добавить такую сортировку. В CUDA нет готовых функций сортировки? По мануалу не нашёл.
Ну там есть вроде какая то реализация Radix, не знаю в библиотеках она типа npp или просто в примерах, никогда ее не использовал. Да и здесь она будет не очень эффективна. Проще реализовать свою. Тут все достаточно просто.
Как я понимаю информация о частицах в каждой ячейке уже есть в каком то виде. Из нее эта сортировка получается достаточно просто: для начала нужно знать количество частиц в каждой ячейке. Алгоритм такой:
1. Выделяется массив целых чисел размером равным общему количеству частиц. Значение каждого элемента этого массива будет номером частицы.
2. Выделяется массив смещений, являющийся массивом целых чисел размером равным общему количеству ячеек.
3. Массив (2) заполняется следующим образом: нулевой элемент содержит 0, каждый следующий i-й элемент равен сумме i-1'го элемента и количества частиц в i-1'й ячейке. Заполнять проще на обычном cpu, т.к. заполнение должно происходить последовательно, если будет являться узким местом, можно написать специальную редукцию на gpu.
4. Для каждой ячейки записываются в массив (1) последовательно номера принадлежащих ей частиц, начиная со смещения из массива (2) для этой ячейки.
На этом сортировка закончена. После этого при вычислении каждый поток просто должен обрабатывать частицу, номер которой хранится в массиве (1) по индексу, равному глобальному индексу потока.
Как я понимаю информация о частицах в каждой ячейке уже есть в каком то виде. Из нее эта сортировка получается достаточно просто: для начала нужно знать количество частиц в каждой ячейке. Алгоритм такой:
1. Выделяется массив целых чисел размером равным общему количеству частиц. Значение каждого элемента этого массива будет номером частицы.
2. Выделяется массив смещений, являющийся массивом целых чисел размером равным общему количеству ячеек.
3. Массив (2) заполняется следующим образом: нулевой элемент содержит 0, каждый следующий i-й элемент равен сумме i-1'го элемента и количества частиц в i-1'й ячейке. Заполнять проще на обычном cpu, т.к. заполнение должно происходить последовательно, если будет являться узким местом, можно написать специальную редукцию на gpu.
4. Для каждой ячейки записываются в массив (1) последовательно номера принадлежащих ей частиц, начиная со смещения из массива (2) для этой ячейки.
На этом сортировка закончена. После этого при вычислении каждый поток просто должен обрабатывать частицу, номер которой хранится в массиве (1) по индексу, равному глобальному индексу потока.
в CUDA есть thrust в котором и fast_radix_sort и inclusive_scan
возьмите пример Particles из CUDA samples за основу, там сетка сделана на основе fast_radix_sort
суть не в порядке warp'ов и shuffle/отсутствии — эти оптимизации стали не актуальными уже на GTX 400, архитектуру очень давно оптимизировали под обычный код (достаточно померить, чтобы убедиться), потому что с warp/shuffle/coalescing в голове было невозможно писать нормальный понятный код. просто много научных работ по теме сделанных ещё на самых первых древних CUDA, и с них тянется по публикациям этот «синдром поиска глубинного научного смысла»
самое главное в fields_v3 это применение inclusive_scan prefix_sum вместо fast_radix_sort, что и дает ускорение на порядок
сетку можно строить БЕЗ сортировки, на основе алгоритма prefix_sum, и работает это в 10 раз быстрее
самое главное в fields_v3 это применение inclusive_scan prefix_sum вместо fast_radix_sort, что и дает ускорение на порядок
сетку можно строить БЕЗ сортировки, на основе алгоритма prefix_sum, и работает это в 10 раз быстрее
Спасибо огромное за статью, очень круто! Расскажите, что за идея с non-constant force field? Какие проблемы решаете? Это связано с constant pH или polarizable force fields? Есть какие-то практические результаты?
По поводу технической стороны — смотрели ли на OpenMM и AMBER? AMBER разрабатывается при непосредственном участии инженеров Nvidia, по производительности очень быстрый код с около нулевой нагрузкой на CPU. От авторов AMBER есть публикации и по технической части, как там всё устроенно, особенно mixed precision model (SPFP).
По поводу технической стороны — смотрели ли на OpenMM и AMBER? AMBER разрабатывается при непосредственном участии инженеров Nvidia, по производительности очень быстрый код с около нулевой нагрузкой на CPU. От авторов AMBER есть публикации и по технической части, как там всё устроенно, особенно mixed precision model (SPFP).
Благодарю.
Идея простая: в классической молекулярной динамике типы частиц, валентные связи и углы задаются постоянными, на все время моделирования. Классическая МД, например, не умеет рвать связи или образовывать новые. Можно обратится к методу ReaxFF, но он слишком медленный. Ну и принцип там другой. Из результатов: попытка описания механизма электропроводности в ванадатных стеклах.
Идея простая: в классической молекулярной динамике типы частиц, валентные связи и углы задаются постоянными, на все время моделирования. Классическая МД, например, не умеет рвать связи или образовывать новые. Можно обратится к методу ReaxFF, но он слишком медленный. Ну и принцип там другой. Из результатов: попытка описания механизма электропроводности в ванадатных стеклах.
смотрели ли на OpenMM и AMBER?Из исходников я ковырял только DL_POLY и Moldy. Честно говоря, я бы хотел, чтобы распараллеливание мне сделали специалисты в этом деле. И я даже общался с математиками на этот счет и обещал какое-то вознаграждение, но дело как-то не пошло и я решил сделать всё сам. Об общих вычислениях на графических процессорах я уже был наслышан, поэтому сразу стал пробовать CUDA. Если есть ссылки на статьи с исходниками или алгоритмами — буду очень признателен.
механизма электропроводности в ванадатных стеклах.И как, сработало?
Мне кажется, что валентные связи фиксируются в МД не столько потому, что это технически проще, а потому, что это физически легче описать и параметризовать. Не очевидно, как «правильно» с точки зрения хим. физики сделать эти связи динамическими. В данный момент есть несколько подходов для частных случаев: QM/MM, constant pH, в каком-то смысле — polarizable force fields.
Не знаю, на сколько сейчас быстры DL_POLY и Moldy, но из того, что пробовал самые быстрые на GPU — это AMBER (non-free open source) и OpenMM (free & open source). Gromacs медленнее. Статья про модель точности AMBER — вот, но у них и другие (третий слайд) публикации с алгоритмами и тех. решениями есть.
Идея простая: в классической молекулярной динамике типы частиц, валентные связи и углы задаются постоянными, на все время моделирования.
От потенциалов зависит, нет? LAMMPS отлично отрабатывает изменение гибридизации углерода при использовании Tersoff или AREBO. Или LAMMPS это не классическая МД?
LAMMPS это не классическая МД?Скорее «продвинутая» или «расширенная».
Классическая. Потенциал Терсоффа — классический трехчастичный потенциал. Функциональная форма довольно сложная и способна воспроизводить несколько гибридизаций углерода. AREBO, судя по названию, использует формализм порядка связей (bond order). Логическим продолжением этих методов является ReaxFF, где потенциал ещё более усложнен, включает, например, такие вклады как штраф за перекоординацию. Потенциальная энергия всё ещё остается непрерывной функцией координат.
Добрый день! Нет времени сейчас погрузиться в детали, но вы что-то слишком переусложнили с сеткой
наиболее каноничный вариант — это example из стандартного набора CUDA samples из SDK под названием Particles — www.youtube.com/watch?v=RqduA7myZok
наиболее быстрый вариант из всех, что мне удалось найти (подробно исследовал вопрос производительности) — это www.fluids3.com — разница с CUDA sample в применении inclusive_scan prefix_sum вместо fastRadixSort (по производительности в 10 раз!)
PDF по теме — on-demand.gputechconf.com/gtc/2014/presentations/S4117-fast-fixed-radius-nearest-neighbor-gpu.pdf
p.s. исследую тему симуляции частиц на CUDA 8 лет — habr.com/ru/post/458612
наиболее каноничный вариант — это example из стандартного набора CUDA samples из SDK под названием Particles — www.youtube.com/watch?v=RqduA7myZok
наиболее быстрый вариант из всех, что мне удалось найти (подробно исследовал вопрос производительности) — это www.fluids3.com — разница с CUDA sample в применении inclusive_scan prefix_sum вместо fastRadixSort (по производительности в 10 раз!)
PDF по теме — on-demand.gputechconf.com/gtc/2014/presentations/S4117-fast-fixed-radius-nearest-neighbor-gpu.pdf
p.s. исследую тему симуляции частиц на CUDA 8 лет — habr.com/ru/post/458612
я бы еще раз особо обратил внимание читателей на самое главное: эффективность распараллеливания (любого) в классической МД драматически сильно зависит от того, как коротко можно обрезать парный потенциал. Собственно всякие ухищрения и придумывают как бы это сделать с наименьшими потерями в точности расчета — ведь на первый взгляд даже небольшая неточность после 1е9 итераций может сделать ценность вашей симуляции нулевой. Общего решения тащем-та нету, поэтому придумывают специфические трюки для специфических задач. Отсюда еще отдельный вопрос как обосновать то приближение потенциала, которое вы выбрали. А так, конечно, сам по себе метод не хитрый.
Метод, действительно, не сложный, так что и школьник может написать основу МД программы. Хитростей в методе 3:
1. термостаты / баростаты;
2. учет электростатики;
3. ускоряющие техники.
1 и 3 мы немного обсудили здесь, а 2 я как раз ковыряю. Сложность там в том, что электростатический потенциал в периодических граничных условиях так просто не обрезать. Эталонным считается суммирование по Эвальду (Particle Mesh Ewald, PME).
1. термостаты / баростаты;
2. учет электростатики;
3. ускоряющие техники.
1 и 3 мы немного обсудили здесь, а 2 я как раз ковыряю. Сложность там в том, что электростатический потенциал в периодических граничных условиях так просто не обрезать. Эталонным считается суммирование по Эвальду (Particle Mesh Ewald, PME).
Вопрос не по теме. Вам приходилось видеть программы для моделирования радиационной химии? Именно МД? Ищу уже года два.
Забыл, спасибо за статью. Вдохновляет.
Зарегистрируйтесь на Хабре, чтобы оставить комментарий
Перенос молекулярной динамики на CUDA. Часть I: Основы