Comments 78
Это всё 2Д? Трехмерное не поддерживается?
Каким алгоритмом перебираются пары частиц? Cell List, Verlet List ? Можно глянуть особенности реализации? Это ваша часть или библиотечная?
Да, это 2д симулятор. Весь код мой, потому что он, в основном - кернелы для CUDA, написанные на Питоне, которые Numb'ой компилятся на лету под GPU - мои вкусы весьма специфичны) Весь остальной код на Питоне - это обертки для вызовов CUDA Kernel'ов и всякая вспомогательная фигня, склеивающая эти вызовы в итоговый продукт.
Я использую что-то вроде hash grid, чтобы обеспечить локальность взаимодействий - в 2д массив записываю индексы частиц в соответствии с их координатами (я это зову полем индексов), а дальше в дело вступает параллельность GPU (ну, точнее, там все расчеты на GPU проводятся, просто особенно важно это именно в вычислении сил, это самая вычислительно-затратная часть, не считая рендера).
Я беру каждую частичку и стартую NxN тредов для нее, где N - это округленный до ближайшего инта максимальный радиус взаимодействия. Его можно поменять, но по дефолту он равен 5, то есть по 25 GPUшных тредов на частицу. Каждый из них смотрит свою ячейку поля индексов, вычисляет силы от частиц из этой точки и плюсует их в переменную, хранящую силы, действующие на частицы.
Посмотреть это можно в репозитории, в начале поста ссылка. Вычисление сил находятся в simba_ps/physics_kernels.py : compute_forces
Каждый из них смотрит свою ячейку поля индексов
свою и соседние?
Нет, только свою. На каждую частицу порождается NxN тредов, смотрящих ее окружение.
тогда получается есть шанс, что не все частицы, которые рядом, провзаимодействуют? если наша частица близко к границе ячейки?
Нет, почему же? Представьте поле из 5х5 клеточек, центрированное на каждой из частиц. В смысле, к каждой частице как бы привяжем такой квадратик из 5х5 клеток, поместив частицу, для которой происходит расчет, в центральную. Круг с радиусом, равным максимальному радиусу, на котором частицы еще взаимодействуют, оказывается вписан в этот квадрат.
Представил. Как перебрать все частицы внутри квадрата? Как получить список таких частиц?
Список и не надо, каждая клетка квадрата - это отдельный поток. Все, что ему нужно сделать - это обратиться к "полю индексов" в соответствии со своей позицией, и посмотреть, есть ли в этом месте частица. Если есть - посчитать, какое от нее будет воздействие на ту, к которой он принадлежит, и добавить результат к переменной-аккумулятору сил. Это же на видюхе все считается, каждая клетка каждого квадратика обсчитывается отдельным тредом.
посчитать, какое от нее будет воздействие на ту, к которой он принадлежит, и добавить результат к переменной-аккумулятору сил. Это же на видюхе все считается, каждая клетка каждого квадратика обсчитывается отдельным тредом.
это всё понятно, вопрос вот где:
обратиться к "полю индексов" в соответствии со своей позицией, и посмотреть, есть ли в этом месте частица
Можно здесь подробнее, как это работает? Что такое "поле индексов"?
есть ли в этом месте частица
одна? может быть в клетке больше одной частицы?
Ну я же выше писал... И в самой статье. В зарубежной литературе это зовут hash grid.
Я использую что-то вроде hash grid, чтобы обеспечить локальность взаимодействий - в 2д массив записываю индексы частиц в соответствии с их координатами (я это зову полем индексов),
Да, там может быть несколько частиц, считайте, что на каждой клеточке стоит колонка, в которой может вообще не быть частиц, а может быть одна или несколько, тред этой клеточки по ним последовательно проходит, просто пока не встретит невалидный индекс.
Так как 
я задаю около 1.0 (в редких случаях, типа симуляции химии с молекулами может быть что-то типа 0.7), система естественным образом стремиться частицы раскидать так, чтоб в клеточке больше одной не было. Они могут туда влететь на короткое время, за счет своей энергии, но ненадолго, и редко когда их там больше штук 3-5.
Да, если по каким-то причинам в симуляции бешеная плотность частиц, то этот "столбик" может заполниться, наверное - для этого у меня там даже есть print, который об этом скажет, но до сих пор такого не было, потому что это все раньше отвалится из-за недостаточного шага интегрирования, либо вовсе до этого не дойдет, если шаг будет достаточным (тогда частицы просто успеют провзаимодействовать раньше и разлететься).
Я почему это все зову по-своему - во-первых, я это дело писал сам, в смысле, чисто решал эту задачу, как мне казалось правильно и эффективно, и только потом узнал, что у этого есть какое-то свое название, а во-вторых - я же это все в машинном обучении использовал, а там все эти индексы пропадают, сетка видит только окружение частиц, как картинку. Для нее это - как будто обсчет небольшого кусочка поля, она выполняет локальный расчет на своем квадратике. Поэтому, по аналогии и "поле индексов" - такое фейковое поле, которое обеспечивает локальность расчетов, когда их производит не сетка, а написанный ручками код.
Жаль, на google Colab запустить не получится.
А ведь такими методами, если всё-таки смоделировать в объемном случае (он нужен в виду природы турбулентности) горение топлив с учётом 20..50 компонентов и верифицировать подобную модель, то с помощью её уже можно будет обучать сети совместимые с RANS моделями турбулентности. Вдруг получим модель на основе сети с точностью и быстродействием превосходящим существующие?
Ну, one step at a time. Как я упоминал в статье, практическая сторона этого вопроса мне тоже очень интересна, поэтому я, собственно, и начал эксперименты - но не стал сразу замахиваться на сверхсложные модели, нужно сначала отработать технологию в 2д) А так - да, очень может быть, это одно из направлений, в котором сетки активно применяют. Я точно знаю, что Дисней разрабатывал технологию, ускоряющую расчет рассеивания света в облаках за счет обученной нейросетки - там как раз такой же принцип - учим на медленной симуляции, потом используем быструю сетку.
Но такое применение довольно очевидно любому, кто хоть чуть-чуть сталкивался с нейросетями. Проблемы начинаются именно в реализации, там хренова тонна вещей, которые можно реализовать по-разному, и, вообще говоря, никто не скажет заранее, какой из вариантов лучше, какой хуже, а какой вообще не заработает. Я вот совсем не сразу подобрал подходящий энкодинг, архитектуру, да и вообще задолбался этот пайплайн строить и отлаживать) Но все-таки сетка обучилась, хотя пока, конечно, медленнее чем мой, довольно оптимизированный, код расчета. Сейчас вот отдохну немного от релиза, и продолжу эксперименты. Мне очень хочется реалтайм дистилляцию получить) Хочу алхимию)
А почему, кстати, не запустить в коллабе? Я этот код удаленно на своей линуксовой машине стартовал вообще без проблем. Только рендер надо отрубить - ну, точнее, не надо врубать, он отдельно от симулятора. И мне кажется должно стартануть.
Вот это да!
Как таковых радиусов у частиц нет, они представлены не шариками в пространстве, а именно точками, у которых есть координаты и скорость.
Я так понимаю "массы" у частиц нет, т.е у летящей частицы скорость просто падает когда она подлетает к покоящейся, с другой стороны оно так и задумано судя по всему.
Кстати вот такая формула для нового значения скорости:
new_v = v + dt * 0.5 * (a + new_a)
а зачем в ней 0.5?
Еще вопрос: dt наверное равен времени кадра в мс?
Нет-нет, масса есть, конечно. Потенциал отвечает за действующие силы, а ускорение из этих сил считается, как обычно, a = F / m.
А 0.5 там потому что так устроено это интегрирование Вереле - по сути этот интегратор говорит "Между дискретными отсчётами будем считать, что всё происходит равноускоренно/равнозамедленно", то есть, под действием постоянной силы. В качестве этого ускорения берется средняя величина между ускорением в прошлом отсчёте и новым, посчитанным из сил на этом шаге. Отсюда и 0.5)
dt - это время между отсчётами, да. Но это не совсем время между кадрами, на каждый кадр идёт N под-шагов. Величина N зависит от энергий в симулируемой системе. Это одна из поразительных вещей, вроде и простая, но не сразу очевидная: виртуальная энергия не "бесплатная", нельзя просто так взять и добавить энергию в столкновение, например - за нее придется платить возросшими затратами реальной энергии, идущей на вычисления, потому что придется делать больше шагов интегрирования. Иначе "бомбанет")
А 0.5 там потому что так устроено это интегрирование Вереле - по сути этот интегратор говорит "Между дискретными отсчётами будем считать, что всё происходит равноускоренно/равнозамедленно",
А для позиции выбран обычный способ:
d_pos = v * dt + (0.5 * dt ** 2) * a
positions[idx, c_idx] = positions[idx, c_idx] + d_pos
а почему бы не использовать метод Верле и там и там?
Нет-нет, масса есть, конечно. Потенциал отвечает за действующие силы, а ускорение из этих сил считается, как обычно, a = F / m.
Ясно, а что происходит с "энергией" которой обладают движущиеся частицы?
она как то учитывается при взаимодействии?
Это и есть "Метод Вереле", все целиком. Тут нет двойного интегрирования, в смысле, это не то что мы сначала применяем этот интегратор к x'' получая x', а потом к x' получая x. Вместо этого позиция считается напрямую, как если бы тело двигалось под действием постоянной силы. Мы сразу получаем точное (ну, в пределах точности метода) приращение координат.
Не очень понял вопрос - в каком смысле, что происходит с энергией? Частицы обладают потенциальной энергией из-за своей позиции в потенциале Леннарда-Джонса, она меняется вместе с изменением их координат. Еще они обладают кинетической энергией, которая половина квадрата скорости помноженная на массу. Она меняется вместе с изменением скорости при движении частицы в направлении минус градиента потенциала - как и положено, они друг в друга перетекают и их сумма для всей системы даже более-менее остается постоянной, за вычетом дрифта из-за неточности интегрирования. Если шагов интегрирования мало, а энергии столкновений большие, то этот дрифт становится больше. Если очень большие, то он начинает расти с каждым отсчетом и симуляция "взрывается".
она как то учитывается при взаимодействии?
учитывается - чтобы приодолеть отталкивание соседней частицы нужно обладать достаточной энергией, приближаясь всё ближе частица замедляется силами отталкивания, при этом потенциальная энергия взаимодействия - растёт. Без термостатов сумма этих энергий - константа. С термостатами будут колебания.
"all-Python" не тоже самое, что и "чистый Python".
А в чем разница? В том, что тут дергается нумба, которая из части этого питона собирает код для GPU?
Ну, я подходил с практической точки зрения - "чистый питон" в том смысле, этот код можно сразу пристыковать к любому питоновскому ML-пайплайну, не городя никаких "переходников", ничего не перекомпиляя и т.п.
Оффтоп: детерминизм симуляции у вас, напомнили мне похожее видео. Идея использования этого для криптографии тоже интересная - по сути это простой алгоритм, но с ооочень длинным ключом (состоянием).
Эх, на такой либе, написать бы продолжение Powder Toy. В свое время не один вечер убил, экспериментируя с разными материалами в попытках добиться чего-то впечатляющего.
Не хватает немного объяснение теории как взаимодействуют частицы в столкновении и как действует гравитация со схемами
Выглядит круто!
Только ссылки бы на рускоязычную Википедию давать (соответствующие статьи в большинстве случаев там есть) - статья на русском всё же.
Мне всегда было интересно, какие знания и как можно дистиллировать в коэффициенты современных сеток. Получится ли обучить ее динамике частиц? А как будет выглядеть выход за установленные пределы по энергии — обычная симуляция просто взрывается, а как себя поведет сеть? А можно ли "на дурачка" промоделить жидкости на уровне частиц, а в сетку дистиллировать знания об их усредненном поведении, чтобы она сама вывела внутри себя аналог уравнений CFD (computational fluid dynamics, уравнений динамики жидкостей), и потом использовать для симуляции уже ее, не симулируя сотню тысяч частиц?
Возможно вы в курсе этой публикации (оригинал), где решается даже боле сложная задача, и обучение производится не на симуляциях, а на эмпирических данных, что логичнее. Задача состоит в оценке размерности сложных динамических систем, которые пока не описываются аналитически, путем обучения на более простых, описание которых уже известно. А эта система сама выводит законы из эмпирических данных. В этой работе исследуются некоторые методические сложности применения обучения нейросетей для вывода физических законов из эмпирическим данным и предлагается вариант их решения.
Что касается обучения на физических симуляциях, то подобное делается давно, например, для задачи трех тел. Эта работа посвящена классификации рассеяния ионов. Использование МО для решения физических задач устойчивый, быстро растущих тренд в физических исследования, а астрофизических носит просто взрывной характер, из-за наблюдательной специфики области.
Возможно также слышали о работах Ванчурина, который вместе с соавторами, вообще вознамерился свести всю физику, включая фундаментальную, и космологию к обучению нейросетей (1, 2, 3, популярно 4, 5, видео). Вот такие дела)
Не исключено будущее именно за такими методами исследования и описания в прикладной и фундаментальной физике, из-за возрастающей сложности объектов исследования, результаты которых невозможно будет свести к традиционным аналитическим моделям. А если и удастся, то для их решение для реальных систем из многих частиц может не хватить никаких вычислительных ресурсов, см., например, лагранжиан СМ. Только приближенные, использующие упрощающие предположения и удовлетворяющие некоторым критериям точности. Но это та область в которой решение задач с помощью обучения нейросетей представляется очень эффективным методом.
Спасибо автору за публикацию и удачи в поисках!
и обучение производится не на симуляциях, а на эмпирических данных, что логичнее
Смотря для какой задачи логичнее) Например, я упоминал в комментариях выше, что подход "обучить на медленной симуляции быструю сеть, чтобы после использовать ее вместо симулятора" тоже очень распространен.
Да, я в курсе, что такие исследования проводятся и это теоретически возможно на самых разных уровнях реализации. Но, как в том фильме - "Винтовок много, но эта - моя") как правило в таких задачах огромное количество "подвижных частей", столько всего можно реализовать по-разному, столько тонкостей, которые влияют на результат и которые в статьях обычно опускаются.
Мне интересно именно посмотреть, что получится на практике конкретно у меня, с моими ресурсами, конкретно на тех задачах, которые интересуют меня.
Про Ванчурина не слышал, спасибо, ознакомлюсь. Да-да, мне тоже кажется, что такие модели и исследования ещё внесут масштабный вклад в наше понимание мира, рад, что есть люди, которые этим занимаются на таком уровне, целыми НИИ.
Когда я работал над обучением сетки динамике частиц, меня прямо поразило то, как некоторые известные из физики вещи как будто бы сами естественным образом вылезают из этой задачи - ну, например, касаемо лагранжиана.
Я подбирал лосс-функцию исходя чисто из желания "чтобы симуляция не бомбила и выглядела более-мне правдоподобно", и пришел к тому, что надо оптимизировать дельту энергий состояний, помноженную на дельту времени - в смысле, что на практике было бы идеально совершать максимально длинные шаги, при этом вызывающие минимальное отклонение в энергии. Ну, или если шаг фиксированный, то просто минимизировать эту величину. Величина получается в джоуль-секундах. "Где-то я это видел", подумал я) И вспомнил, что это, вообще говоря, action, который в The Principle Of Stationary Action, фундаментальнейшем принципе физики)
Добавим сюда же следующее наблюдение: в идеале надо бы шаг интегрирования сделать не общим для всех частиц, а локальным - ну то есть, вот в этой симуляции столкновения, допустим, приходится делать достаточно много под-шагов, чтобы не "бомбануло", из-за высокой энергии столкновения. Но ведь вдали от точки столкновения в течение долгого времени частицы весьма "холодные", они себя нормально ведут даже при очень маленьком количестве шагов интегрирования! Если бы удалось реализовать увеличение количества под-шагов (уменьшения шага интегрирования) в зависимости от плотности энергии... Но опять я это где-то видел - это же time dilation!)
Вообще, если уж так фантазировать, то это выглядит так, будто вселенная старается оптимизировать вычислительные ресурсы при помощи принципа стационарного действия, замедления времени и второго закона термодинамики (стараясь максимально рассеять энергию в пространстве, повышая энтропию и снижая плотность энергии, а значит, и количество необходимых вычислительных мощностей) - все ресурсы идут на вычисление всяких крутых высокоэнергетичных штук типа нейтронных звёзд и черных дыр, а остальное достается менее вычислительно-затратная областям)
Но, конечно, это запросто может быть "поиск глубинного смысла" там, где его нет. В любом случае, у меня пока даже не получилось сделать шаг локальным - симуляция взрывается при взаимодействии двух областей с разными шагом интегрирования) Возможно, с этим можно что-то сделать, но пока я этот вопрос отложил. В любом случае, это все очень интересно)
Смотря для какой задачи логичнее
Логичнее для тех задач, для которых имеется необходимый объем эмпирических данных, нет мат. моделей, или расчет по ним весьма затратен. Например, мы хотим предсказать динамику пламени, которую зафиксировали на видео. Для этого подаем на вход сети кадр и пытаемся предсказать следующий. После обучения она должна с некоторой точность предсказывать эту динамику для других условий. Но есть задачи для которых такой набор эмпирических данных получить затруднительно. Для той же задачи 3-х тел. Где в природе найти такую систему и долго наблюдать и фиксировать ее поведение, которое затем использовать для обучения? Желательно для разных случаев с разными параметрами? В этом случае можно использовать только вычисленные траектории. Однако в этом случае нет гарантий, что поведение реальных систем будут именно таким, т.к. обучение происходило по симуляциям, а они подвержены не только накоплению вычислительных ошибок, но ограничений исходных формальных моделей. Например, модель ньютоновская, а в реальности, в случае сильных потенциалов, тела могут испытывать взаимодействия описываемые ОТО. Или вообще неизвестных эффектов. Конечно этому подвержено и обучение на конечном наборе эмпирических данных, но в меньшей степени. Оно как раз имеет смысл для учета необычного поведения, не поддающегося формальному описанию. Вероятно, в каких-то случаях, полезно и комбинированное обучение, как по эмпирическим данным, так и симуляциям.
Но, конечно, это запросто может быть "поиск глубинного смысла" там, где его нет.
Да, конечно. В реальности никакого пламени и его динамики нет, да и частиц с примитивным законом взаимодействия также. Это приближенные модели, не важно, на уровне восприятия, или абстрактного описания. Мне, как физику по образованию занятому в психофизиологических исследованиях, трудно ассоциировать эти модели с реальностью. Скорее это отражает устройство мозга, его нейросетей, которые экономят энергию на функционирование, не только по тому что это дорогой для организма ресурс, но и опасный из-за чрезмерного выделении тепла при метаболизме. Разработчикам микропроцессор этот эффект хорошо известен, и является одной из основных проблем в их разработке. Мозг эволюционно настроен на экономию описания на всех уровнях начиная от восприятия и завершая высшими формами мышления у человека. Что там в реальности, что там во Вселенной достоверно не известно. Физики для описания поведения объектов используют разнообразные вариационные принципы, которые не могут обосновать, и они носят статус эвристических, методологических. А в нейрофизиологии управление поведением, восприятие, мышление, и др. ментальные феномены регулируются принципом свободной энергии, носящего фундаментальный характер, и могущего внести ясность в источник происхождения физических вариационных принципов. Возможно эта та причина благодаря которой подход предлагаемый Ванчуриным имеет под собой некоторые основания — не сама Вселенная является обучаемой нейросетью, а ее ментальный образ в обученных нейросетях мозга.
Это понятно, я скорее о том, что есть чисто практически задачи - ускорить расчет рассеивания света в облаках, например - как это сделали исследователи из RnD отдела Диснея. У них была модель, которая давала подходящие результаты, но вычисление было медленным. Они использовали ее для обучения нейросети, при помощи которой получили реалтайм расчет - он, может, и расходится с реальностью на каком-то знаке после запятой, но так как это не для фундаментальных исследований, а для реалтайм генерации реалистичных изображений, эти расхождения не так важны.
По поводу второго тезиса - ну, "вселенная как нейросеть" это, конечно, может быть чересчур, но сам подход "мы ищем модель, которая будет согласоваться с наблюдениями и при этом минимизировать вычислительные затраты", как мне кажется, может оказаться полезным.
Я имею в виду, что минимизация экшена (того, что джоуль-секудны) и замедление времени в присутствии высокой плотности энергии, вместе с возрастающей энтропией (то есть, с все возрастающим рассеиванием этой энергии всем возможным количеством способов), как будто бы одновременно являются механизмами, которые позволяют сократить вычислительные затраты.
Разумеется, я отдаю себе отчёт в том, что настоящие физические модели, которые описывают квантовые поля - это совсем не то же самое, что точечные частички в потенциале) Но принципы, которые я упомянул, мне кажется, в той или иной мере могут быть экстраполированы.
Когда делаете шаг локальным, он кратен? Т.е отличается у соседних ячеек в целое число раз?
Нет, не думаю. Я по-разному пробовал, но если я правильно помню, последний проверенный вариант выглядел примерно так: есть dt между кадрами, делю его, скажем, на 20 под-шагов для всех частиц. Но силу, действующую на частицу я пересчитываю не на каждом шаге, а в зависимости от плотности энергии.
То есть, если частица в "холодном" окружении, сила будет рассчитана, грубо говоря, пару раз за эти 20 шагов, остальное время она будет двигаться под действием постоянной силы. Если она в высокоэнергетическом "горячем" окружении, то пересчет будет выполнен чаще, вплоть до расчёта на каждом шаге.
А энергию частицы передают только через столкновения, тепловое излучение не учитывается?
Ну, конечно не учитывается - это же нужно, получается, параллельно моделировать электромагнитные эффекты - совсем другая задача. Максимум, что я в этом плане делал - "подкладывал" под симулируемые частицы вторую симуляцию, распространения волн в плоскости - это не особо "физично" было, в том смысле, что я не стремился этим промоделировать какие-то реальные эффекты, скорее было просто интересно, что из этого получится (я так симулировал звук от своего парового двигателя на этих частичках) - ну вот можно в это поле энергию сливать) Будет poor man's "излучение")
Вот так это выглядит. Тут это "поле" специально только в нижней части, но вообще я его подкладывал и прямо под частицы. В моем коде частицы вызывали изменения амплитуды в этом поле, а само оно на них никак не влияло, но думаю, можно при желании закодить и так, чтобы энергия перетекала в него и обратно.
При большом количестве частиц они ведут себя как жидкость или газ. Не пробовали решать уравнения Навье-Стокса? Интересно было бы посмотреть на поведение разных математических моделей в одинаковых условиях.
Да, ведут, вон я их кипячу и дистиллирую)
Модели, которые описывают усреднённую динамику частиц я пока не использовал - собственно, это и есть одна из причин, по которым я начал этот проект - мне интересно, можно ли получить нейросеть, которая будет повторять динамику этих моделей (типа того же Навье-Стокса), но при этом выведет ее сама из данных, полученных при симуляции жидкостей на уровне частиц.
Очень интересно. По идее что-то получиться должно. Мне кажется главный вопрос жизни, вселенной и всего такого - сможет ли нейросеть стать более эффективной с точки зрения точности и вычислительной эффективности по сравнению с классической системой уравнений Навье-Стокса. На стороне классики - отсутствие лишних вычислений связанных с ненужными коэффициентами нейросети. На стороне нейросети - способность обобщать, соответственно, не нужно использовать сплошную мелку сетку.
Оказывается в последнее время исследования на тему моделирования жидкости и газа при помощи нейросетей ведутся довольно активно.
Вот короткое видео, в котором показано что-то вроде MVP. Есть ссылки на статьи, проект на Гитхаб
https://www.youtube.com/watch?v=ISp-hq6AH3Q
Отличная презентация того как устроено моделирование жидкости при помощи Навье-Стокса и нейросетей.
https://www.youtube.com/watch?v=IXMSOSEj14Q
На самом деле, информации очень много, хотелось бы привести больше ссылок, но лучше просто погуглить "neural network navier stokes".
Надеюсь, будет полезно.
Del — Ошибся веткой.
Есть такая игра
https://noitagame.com/
Автор очень хорошо реализовал физику частиц в 2D. Вам будет интересно взглянуть и может быть почитать devlog.
Я знаю и люблю эту игру. Но физика частиц там очень специфическая, в силу ограничений, накладываемых самой задачей - им нужна очень быстрая реалтайм модель, в которой не так важна физическая достоверность и разнообразие взаимодействий)
Поэтому тамошние частички очень упрощённо взаимодействуют друг с другом и не образуют макро-структур. Вы не увидите в Нойте результатов столкновений типа того, что в заголовке поста, не будет слипающихся друг с другом "осколков", не будет настоящего кипения и т.п. Это я не к тому, что "а-ха-ха, смотрите, а у меня лучше", а к тому, что это разные модели для разных задач.
Если у меня получится ускорить вычисления при помощи нейросети, попробую собрать на этом какой-нибудь интерактив)
Кстати, то количество частиц, что в примере impact, в целом, более-менее реалтаймово тянется. ФПС на 20 вместо 50, но всё-таки. Но все красивые штуки, в основном, требуют большую плотность энергии, а значит - увеличение количества под-шагов интегрирования. Дистилляция, например, тоже даёт примерно 30 ФПС, но она при этом в 10 раз медленнее идёт, чем в гифке и на видео - в игре это не особо прикольно было бы, ждать 10 минут, пока у тебя зелье сварится.
Хэй! оч круто!
Делал такое 12 лет назад на с+куда+опенгл на дряхлой 8400GS... оч нравилось видеть результат симуляции.
Также делал на дипломе молекулярную динамику разделительной центрифуги для двухкомпонентной смеси и даже наблюдался коэффициент разделения, несмотря на скромное количество частиц.
Вот даже историческое демо осталось: https://youtu.be/QEkC1k0gdnI
Красивый проект.
А вы пробовали задавать расстояние, в котором частицы взаимодействуют, в зависимости от силы взаимодействия. То есть отсекать все, сила взаимодействия с которыми меньше epsilon.
А зачем? Чтобы узнать силу взаимодействия и увидеть, что она меньше, чем эпсилон, надо сначала ее посчитать, то есть посчитать взаимодействие со всеми частицами из окружения, что сводит на нет смысл такой "оптимизации".
Другое дело, что я сходу не считаю взаимодействия на расстоянии больше некоторого предельного радиуса - это получается само собой за счёт того, что все взаимодействия тут локальные, частица взаимодействует только со своим непосредственным окружением как раз в этом предельном радиусе. Это автоматически ведёт к тому, что она не взаимодействует с частицами на расстояниях, на которых сила пренебрежимо мала.
Очень красивые симуляции. За САЕ и МВК отдельный респект: очень люблю вторую часть Гоблинов, и эта любовь даже передалась моим детям.
По сути: не смотрели в сторону квантовомеханических симуляций в потенциале ЛД? Это сейчас хот-топик в экспериментальной атомной физике и ультрахолодной химии. Там минимум этого потенциала формирует резонанс Фешбаха, которым управляют с помощью небольшого магнитного поля. См, например, здесь
https://nplus1.ru/material/2023/02/08/supercool-four
и далее по ссылкам
Ну и в целом в квантовой физике дефицит наглядных представлений. Оно и понятно: волновая функция - это функция трех координат, поэтому даже типичные атомные орбитали выглядят на картинках учебников условно и не совсем так, как их стоило бы воспринимать на самом деле.
Я обычно своим студентам привожу в пример игру Quantum Moves, в которую можно поиграть в браузере. Эта игра предлагает решить проблему перемещения атома оптическим пинцетом по нужной траектории, "не расплескав" его волновую функцию. Там, конечно, тоже двумерия, но во всяком случае оно дает студентам интуитивное представление о квантовом поведении объектов с помощью численного решения динамического уравнения Шрёдингера.
Если я правильно помню, то волновая функция еще и растет размерностью вместе с количеством частиц, то есть если одна частица может двигаться в трех измерениях, ее волновая функция будет трехмерная, а если таких частиц две - то уже в шести)
В смысле, пространство-то и остается трехмерным, но именно изменения в конфигурации волновой функции путешествуют уже в их совместных шести.
Нет, в сторону квантовых симуляций я пока не смотрел. Это тоже, конечно, интересно) Но я все-таки хотел бы сначала попробовать ускорить расчеты в 2д при помощи нейросетки. Не знаю, получится или нет, я пока только в самом начале исследований.
Все верно) это мы ещё спин не учитываем)
Очень интересно, что у вас выйдет! Если нейросети смогут предсказывать исходы физических процессов быстрее, чем прямой расчёт, это может быть революцией. Я думаю о том, как это сможет ускорить работу физических движков, чтобы запускать реалистичную графику на слабых компьютерах
Я бы не назвал это революцией, потому что это уже делают. Но это точно повысит привлекательность этого кода до коммерческой)
Если нейросети смогут предсказывать исходы физических процессов быстрее, чем прямой расчёт, это может быть революцией.
Есть уже довольно старая работа по предсказанию ВФ молекулярной электронной структуры материалов, и соответственно их свойств, в квантовой химии. Ценна также списком ссылок на нее.
Визуализация для уроков физики хороша. Желаю вам успешного движения глубже в квантовый мир!
Вы баффнули Хабр тортом на несколько месяцев.
Хочу сразу разъяснить читателям. В реальной жизни межатомные взаимодействия описываются именно такими кривыми, как приведенная в статье Леннарда-Джонса или потенциал Морзе (электростатические взаимодействия для атомов в составе одной молекулы).
Атомы можно ОЧЕНЬ ГРУБО представить в виде отталкивающихся друг от друга магнитов, помещенных в липкие слаймы. Магниты - это компактные ядра (которые отталкиваются друг от друга в реальности не из-за магнитных взаимодействий, это ради наглядности!), а слаймы вокруг - электронные облака (которые слипаются - аналог спаривания электронов с образованием ковалентной химической связи).
Так вот, если вы возьмете две такие штуки, и будете вжимать друг в друга, они будут яростно сопротивляться - чем ближе, тем сильнее отталкиваются магниты, поэтому потенциал слева на графике идет в бесконечность. Потом, если отпустить их, система придет в равновесие: отталкивание магнитных сердцевин уравновесится тягучестью слаймовых оболочек, будет два слипшихся шарика, потенциальная энергия этой системы придет к минимуму (та самая потенциальная яма). А если растащите еще дальше, то слаймы будут сопротивляться, пока липучести хватает, затем они оторвутся - взаимодействие прекратится вовсе. Поэтому справа график приближается, а затем становится равен нулю в зависимости от расстояния между объектами.
Еще раз: у реальных атомов между ядрами будут не магнитные, а кулоновские и квантовые силы отталкивания, а на больших расстояниях их вместе удерживают ковалентные, электростатические ионные, Ван-дер-Ваальсовы силы взаимодействия. Но надеюсь, аналогия "отталкивающе-липких шариков" получилась наглядной.
Это насчет межАТОМНЫХ взаимодействий. С межМолекулярными чуть по-другому: если атомы в молекуле стабильны, существует только отталкивание. К нестабильным (готовым образовать химическую связь) атомам могут присоединяться другие, начиная с определенных энергий столкновения. Это и продемонстрировано автором в его экзотермической реакции.
Осталось ответить на главный вопрос.
Зачем все это?
Будучи студентом химико-технологического вуза в конце нулевых, я бы поклонился вам в ноги за такие визуализации. И вообще, для популяризации и наглядности я бы очень советовал вам (ну или кто еще этот код подхватит) продолжать.
Куда можно продолжать?
физико-химические процессы: симуляция идеального газа, термодинамических изопроцессов, диффузии, адсорбции;
химические процессы: кинетика реакций 0, 1, 2 порядка, радикальные цепные процессы;
технологические процессы: реология твердых частиц, процессы в псевдоожиженном слое, осмос, ректификация;
Даже в упрощенном виде, для демонстрации, а не реальных расчетов, такие симуляции ОЧЕНЬ полезны.
В радикальных цепных процессах было бы интересно посмотреть на формированиие ударной волны
Спасибо на добром слове) У меня на этот год намечены несколько проектов, если никакой очередной пиздец не наступит, попробую их довести до значимых майлстоунов. Один проект прямо сильно на базе его кода, другой - скажем так, на опыте, который я получил из реализации симулятора, с некоторым заимствованием кода из него, и третий - это вот исследования, связанные с обучением нейросети этой динамике.
Вот бы автор статьи, Onigiri и Техношаман52 объединили усилия и создали супер клеточный автомат с эволюцией, деревьями, светом, ресурсами и агентами.
С какой скоростью идут расчеты и при каком количестве частиц? Есть планы добавить третье измерение? Очевидно, что вычисления идут в кернелах, но все же питон не дает оверхеда?
Оверхед от питона очень незначительный, большую часть времени идут расчеты сил. На втором месте рендер, но он дергается намного реже.
Примерно так это выглядит)
Ну вот пример из заголовка, с достаточно высокоэнергетическим столкновением при 10 под-шагах выдает у меня около 50 FPS, реалтайм. При ~30к частиц. У моей подруги на 2060 дает 60 FPS на тех же настройках. Все зависит от количества под-шагов, а оно зависит от предельной энергии, которую может выдержать симуляция, не развалившись) В этом случае 10 хватает, чтобы не бомбануло.
В более затратных, типа дистилляции, конечно, либо сам процесс будет идти медленно - у меня тоже примерно 40 FPS выдает, но сама дистилляция длится 10 минут (в видео ускорено в 10 раз) - либо надо повышать энергию частиц вместе с количеством под-шагов.
65к частиц примера vortex, где они с большими энергиями крутятся, у меня идет на 15 фпс примерно, на ноуте.
Нет, 3д я пока не планирую, там нужно много чего менять будет прямо в самом сердце алгоритма, в расчете сил. Сейчас он базируется на 2д hash-grid (или, как я это зову, "поле индексов"), и этот код я достаточно долго оптимизировал - собственно, эти графики - это как раз и есть процесс оптимизации, я пробовал разные штуки и замерял скорость. В 3д нужно уже что-то другое придумывать и снова оптимизировать.
Я, скорее, вместо 3д буду пытаться 2д ускорить, очень уж хочу дистилляцию быструю и в реалтайме.
Очень крутая статья, спасибо! Если в школе по физики была бы подобная визуализация, я бы реже прогуливал уроки. :)
Творческих успехов, буду ждать более подробнее про то, как реализованы хотя бы несколько примеров!
Я правильно понимаю, что вы пишете аналог симулятора LAMMPS, только на Python? Чем не устраивает этот симулятор на C++? Он вроде под открытой лицензией и часто используется в научных исследованиях для симуляции частиц и химических реакций.
Ну, я на это уже отвечал. "Таких винтовок много, но эта - моя")
Хотя на самом деле - не то чтобы много. Даже если опустить тот момент, что мне захотелось разработать эту штуку самому, она еще и на других технологиях. Этот симулятор полностью на питоне, что позволяет его очень легко стыковать с пайплайном машинного обучения. Я его сделал для не для исследований химических реакций, а для исследований в области машинного обучения, прежде всего.
И мне копаться в своем питоновском коде (зная каждую его строку, каждый компромисс, каждую оптимизацию) куда проще, чем прикручивать к питоновскому пайплайну с пайторчем эту махину на C++ и перекомпилять ее каждый раз, когда мне потребуется там что-то отрезать или прикрутить)
Случайно промахнулся и отклонил чей-то коммент, не успел запомнить автора. Автор, прошу прощения, напиши еще раз) Там было про кулоновский потенциал и про то, с ним не получится считать воздействия только ближайшего окружения. Почему, кстати? Он же тоже падает до нуля с расстоянием.
И, кстати, почему некоторые комменты не сразу постятся, а ожидают от меня подтверждения? Это при низкой карме?
Почему, кстати? Он же тоже падает до нуля с расстоянием.
Падает не так быстро, как хотелось бы. Поэтому накапливаются слишком большие ошибки, если поступать с ним как с короткодействующими потенциалами. Обычно электростатику считают методом Эвальда (PME, Particle mesh Ewald), у меня есть небольшая статья про это. Но в одной из работ показано, что можно вместо Эвальда использовать хитрый короткодействующий потенциал. Если надо, найду эту работу и формулу.
ПС, начальный коммент не мой, не знаю, что там было
Лет 10 назад (ВолГУ, 2013) писал дипломную по ускорению алгоритма Барнс-Хата на видеокартах NVIDIA. Изначально был старый код на Фортране и входные данные для него. На курсовых 2010-2012 переписывал одномерную и двумерную версию с использованием PGI CUDA Fortran Compiler, а в дипломной попытался в 3D, но в одиночку не осилил, и нашёл в итоге treecode/bonsai на C++/CUDA, который смог скомпилить под винду и адаптировать входные параметры в нужный формат. Далее всё это запускалось на рабочих станциях с Tesla C2050 под виндой. И визуализатор там тоже был, выглядит как-то так: https://www.youtube.com/watch?v=aPgzo9Mvk6o
Немного не по теме, т.к. BH в основном используется для гравитационных симуляций, но наверное тоже можно использовать для обычных симуляций с частицами?
Да, но тут ситуация немного другая - за пределами нескольких 
потенциал почти нулевой, в основном тут важно именно действие соседних частиц. BH, если я правильно понимаю, скорее для случая, когда есть силы, тянущиеся на огромные расстояния - слипающиеся под гравитацией частички для других, удаленных от них, частиц, выглядят как одна большая гравитирующая масса, поэтому ее нельзя просто так отбросить.
Да, наверное оптимизировать расчёт разбивкой на динамическую сетку здесь не получится.
Хотя, если у нас в симуляции 10,000 частиц, то сложность алгоритма будет O(N*K), где K — максимальное количество соседних частиц? По-идее сейчас количество CUDA-ядер и shared memory сильно больше, чем в 2013, и если попробовать разбить симуляцию на квадраты/кубы (2D/3D) одинакового размера, то каждое ядро сможет параллельно обсчитывать десятки тысяч взаимодействий на цикл. Только вопрос будет в граничных условиях, можно ли обойтись без копирования частиц из соседних блоков памяти.
Поправьте, если не так.
Так он и сейчас оптимизирован, потому что обсчитывается только ближайшее окружение. Что-то типа BH вводить, я так понимаю, имеет смысл, если что-то за пределами этого окружения вносит большой вклад, который нельзя отбросить - не для скорости, а для точности. Типа того, что выше сказали про кулоновское притяжение/отталкивание.
В том коде, который сейчас, как мне кажется, надо смотреть в сторону какой-то хитрой магии с количеством шагов интегрирования. Не дает мне покоя тот факт, что при том же столкновении большинство частиц не требует такого большого количества под-шагов (маленького шага интегрирования), потому что до них далеко не сразу докатывается волна, да и после того, как она схлынет, они тоже вполне могут быть просчитаны с очень небольшим количеством под-шагов. Много шагов интегрирования нужно только там, где высокая плотность энергии, а она совсем не одинаковая в пространстве, да еще и постоянно падает (солдат спит - энтропия растет!). Но пока у меня не получилось интегрировать с разным шагом, симуляция взрывается) Это как раз то, на что идет большая часть вычислительных ресурсов, и частичкам, которые движутся совместно (даже с большой общей скоростью) не особо нужны эти промежуточные точки, но так как шаг общий для всех, на это тратится куча ресурсов.
попробовать разбить симуляцию на квадраты/кубы (2D/3D) одинакового размера
так и делают, это называется cell list
Только вопрос будет в граничных условиях, можно ли обойтись без копирования частиц из соседних блоков памяти.
Все приходится держать в global memory
Напоминает чем-то клеточные автоматы.
Интересная работа. Сам много лет, как хобби делаю подобные симуляторы время от времени. Но в этом году делаю уже ЭМ взаимодействие заряженных частиц с учетом ТО методом учета потенциалов Лиенара-Вихерта и как следствие в 3d.
Проблеме интегрирования (и взрыва) посвятил много вечеров, в конечном итоге пришел к формуле, которая конечно усложняет вычисление одной пары, но значительно уменьшает ошибку.
Суть ее в том, что силу между точками интегрирования, нужно проинтегрировать и получить среднее.
Скажем для кулоновской силы F = k / (R[t] * R[t+1]), действующая сила будет зависеть от будущего положения :) но у этого парадокса таки есть решение, которым поделюсь если найду силы и вдохновение опубликоваться на хабре.
В вашем случае потенциалов "6-12" даже проще, вам нужно учесть, что суммы энергий не изменились до и после просчета для каждой пары: Vi+Ui +Vj+Uj = const . Кстати такая поправка уменьшает проблему порядка расчета и проблему float32 ошибки. Кстати эти размышления приводят к необходимости введения квантования энергии))
Для контроля ошибки можно применять такие показатели:
% отклонения энергии всей системы
средне квадратичное отклонение между расчетами с разными шагами интегрирования (например, 1 секунда за 10 шагов, против 1 секунда за 1000 шагов) Предлагаю создать чат группу Telegram - любителей симуляции частиц))
нужно учесть, что суммы энергий не изменились до и после просчета для каждой пары: Vi+Ui +Vj+Uj = const
Но это не похоже на правду. Или я не так понял это предложение. Если во взаимодействии участвуют несколько частиц, то не изменилась сумма энергий их всех, а не каждой пары в отдельности. А в этом случае не очень понимаю, как мне поможет то, что я знаю локальную сумму энергий, распределиться по частцам-то она могла как угодно. Можно попробовать собирать эту локальную плотность энергии и при ее помощи попытаться локально же масштабировать скорости, но это надо учитывать все частицы, которые были во взаимодействии в t0 и t1 (включая те, что влетели в радиус взаимодействия и вылетели оттуда).
Суть ее в том, что силу между точками интегрирования, нужно проинтегрировать и получить среднее.
Ну так это вычислительно примерно то же самое, что и интегратор более высокого порядка получается же? Хотелось бы не просто сложность вычислений повысить ради точности, а увеличить скорость расчетов.
Вот про метрики для ошибки - совершенно согласен, сам к таким же пришел. Более того, если предположить, что делаем динамический шаг, и нам нужно, чтобы система совершила шаг, максимальный по времени и минимальный по отклонению, получается, что нужно оптимизировать dE*dt - что внезапно совпадает с экшеном (тем, что джоуль-секудны) из принципа стационарного действия!
На мой взгляд, больше всего лишних вычислительных ресурсов тратится на под-шаги интегрирования там, где этого не требуется - я выше про это писал, в примере со столкновением, например, большая часть частиц могла бы быть просчитана с 1-2 шагами на кадр, только в точке удара (и на гребне распространяющейся волны) требуется больше. Поэтому нужно бы что-то вроде time dilation, замедления времени в присутствии большой плотности энергии, сделать количество под-шагов локальным и зависящим от локальной плотности энергии, чем больше - тем больше шагов, тем сильнее "замедляется время". Но на практике у меня пока не получилось это реализовать, взрывается)
Ну так схема Верле как раз берет среднюю силу между текущим и следующим шагом.
Спасибо за проделанную работу и интересную публикацию.
Я сам писал 2D симулятор бильярдных шариков, но когда включал для них гравитацию, у меня, видимо из-за ошибок интегрирования быстро росла энергия всей системы и симуляция разваливалась. Можете подробнее рассказать, натыкались ли на эту проблему, как решали? Можете дать совет, как решать проблему роста энергии?
Да, разваливается из-за ошибок интегрирования. Чем выше энергия столкновений, тем меньше нужен шаг интегрирования, это реальная плата за виртуальную энергию)
Можно добавить больше под-шагов на кадр (уменьшив тем самым шаг интегрирования), можно использовать интегратор более высокого порядка - но это все решения одного плана, за счет повышения вычислительной сложности задачи.
Можно добавить термостат или вовсе жестко ограничивать максимальную силу, действующую на частицу или максимальную энергию частицы - это сделает симуляцию чуть менее реалистичной с точки зрения физики, но в ряде случаев это приемлемо, если правильно подобрать параметры, то визуально будет по-прежнему смотреться нормально.
Можно подобрать параметры симуляции, постаравшись снизить энергию столкновений - например, частички, которые падают "в вакууме" под действием гравитации на зафиксированные препятствия, вероятнее "взорвут" симуляцию, чем если они падают в разреженную массу таких же частичек, потому что во втором случае они успеют часть энергии передать окружению. Сюда же - понизить гравитацию или подобрать высоту, с которой они падают так, чтобы они не успели приобрести высокую кинетическую энергию.
Можно реализовать динамический шаг интегрирования, чтобы при превышении частицей какого-то порога энергии, добавлялись под-шаги, разбивая интервал интегрирования, скажем, на два.
Не описан существующий софт/игры симуляторов и почему вам необходимо было создать свой.
Потрясающе! Пропустил летом вашу статью, сейчас случайно наткнулся, искал формулы Леннарда-Джонса
делаю с 2012 года похожий проект с GPU 2D particels симуляцией только на JavaScript, компилирую самописный DSL в CUDA и OpenCL
Только у меня цель - виртуальная эволюция существ, нейросети тоже есть но ориентированы на симуляцию поведения. о возможности предсказания свёрточной нейросетью симуляций и построения на этой основе масштабируемой симуляции я тоже думал, даже делал зарисовки, очень круто что вам удалось это воплотить!
на Хабре писал об опыте создания: https://habr.com/ru/articles/458612/
примеры физики:
https://www.youtube.com/watch?v=0EVQaWqUQOg
https://www.youtube.com/watch?v=UYw2kWtZsYQ
p.s. ещё знаю alien-project.org, тоже 2D CUDA open source симуляция с прицелом на эволюцию, пример его физики: https://www.youtube.com/watch?v=iUA-w5GF0Ck
Simba: Симуляция десятков тысяч частиц в потенциале Леннарда-Джонса на чистом Python с GPU-ускорением