Симуляция подъёмной силы Ньютона методом частиц на CUDA


https://www.youtube.com/playlist?list=PLwr8DnSlIMg0KABru36pg4CvbfkhBofAi


Как-то на Хабре мне попалась довольно любопытная статья “Научно-технические мифы, часть 1. Почему летают самолёты?”. Статья довольно подробно описывает, какие проблемы возникают при попытке объяснить подъёмную силу крыльев через закон Бернулли или модель подъёмной силы Ньютона (Newtonian lift). И хотя статья предлагает другие объяснения, мне бы всё же хотелось остановиться на модели Ньютона подробнее. Да, модель Ньютона не полна и имеет допущения, но она даёт более точное и интуитивное описание явлений, чем закон Бернулли.


Основной недостаток этой модели — это отсутствие взаимодействия частиц газа друг с другом. Из-за этого при нормальных условиях она даёт некорректные результаты, хотя всё ещё может применяться для экстремальных условий, где взаимодействием можно пренебречь.


Я же решил проверить, что же произойдёт в модели Ньютона если её улучшить. Что если добавить в неё недостающий элемент межатомного взаимодействия? Исходный код и бинарники получившегося симулятора доступны на GitHub.


Computational Fluid Dynamics


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


Если вас пугает внешний вид этих уравнений, но вы хотели бы разобраться, то для вас есть очень хорошее объяснение в 7-м томе лекций Ричарда Фейнмана. Загляните в главы 38 (“Упругость”), 40 (“Течение ‘сухой’ воды”) и 41 (“Течение ‘мокрой’ воды”). Если совсем кратко и на пальцах, то система уравнений Навье-Стокса — это векторное уравнение второго закона Ньютона. Эта система определяет равнодействующую всех сил (давления, вязкости и гравитации) для всех направлений. Дополнительно может быть задано второе векторное уравнение, обеспечивающее условие несжимаемости, если нужно описать жидкости.


Принципиально систему уравнений можно решать двумя подходами: методами Эйлера или Лагранжа. Эйлеров подход рассматривает среду как поля векторных и скалярных величин. Подход Лагранжа рассматривает отдельно каждую частицу среды.


Один из многих способов численно решить систему уравнений — это применить метод конечных элементов в связке с адаптивным сеточным методом в случае подхода Эйлера. Зачем нужна адаптивность сетки и как её можно реализовать, подробно и доступно рассказали ребята из SpaceX в своём докладе. Эйлеров подход как правило, но не обязательно, применяется для моделирования замкнутых объемов неразрывных сред, т. е. тех, в которых отсутствуют пустые места и включения других сред. Для иных сред чаще всего реализуют подход Лагранжа, через метод сглаженных частиц (Smoothed-particle hydrodynamics, SPH). Метод активно применяется для моделирования воды с полным набором явлений: брызги, капли, лужи, смачивание поверхностей и т. д. Можно даже сымитировать пену или пузыри, если включить в модель частицы воздуха. Реконструкцию поверхности, а точнее изоповерхности, можно произвести любым интересующим вас способом (screen-space meshes, dual contouring, marching tetrahedra, metaballs). Если вы знаете другие интересные подходы, добро пожаловать в комментарии.


Discrete Element Method


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


Метод дискретного элемента (Discrete Element Method, DEM) показался мне очень привлекательным, потому что он решает практически схожую задачу. Однако он в первую очередь предназначен для сыпучих и гранулированных сред, где у каждой частицы в трёхмерном пространстве есть позиция, ориентация, произвольная форма и масса.



Чтобы упростить вычисления, я решил оставить у частиц только две степени свободы (координаты X и Y), одинаковую массу и радиус. При построении своей модели в угоду производительности я хотел отбросить параметры и факторы, которые порождают эффекты второго порядка. Однако при моделировании сложных систем они могут быть очень существенны. Один из показательных примеров — это использование NASA модели идеального газа вместо реального при проектировании космических челноков. В результате, во время миссии STS-1 проявились различные аномалии при входе в атмосферу. Подробнее в разделе Mission Anomalies.


Тем не менее у DEM есть одна важная особенность — это обнаружение столкновений постфактум (Discrete Collision Detection). Разрешение столкновений происходит простым силовым воздействием по закону Гука.


В противоположность этому подходу существует априорный метод Continuous Collision Detection (CCD), который рассчитывает, когда столкновение произойдёт в будущем. Зная точное время контакта, можно скорректировать временной шаг, и избежать неприятных физических артефактов. Метод активно применяется в современных играх. Для игр CCD очень важен чтобы объекты не туннелировали друг через друга, не проваливались друг в друга и не «застревали в текстурах». Метод поддерживается современными движками, в Unity и в Unreal точно.



Подробный доклад о методе Continuous Collision Detection


Но у CCD есть две большие проблемы. Первая — это при всей своей точности он непропорционально ресурсоёмок. А для большого количества частиц он сделает расчёты совершенно непрактичными. Вторая — при увеличении плотности частиц, время свободного пробега каждой уменьшается, вынуждая симулятор делать всё меньшие и меньшие шаги, чтобы не пропустить столкновения. В одном из моих проектов очень хорошо видно, как временной шаг вырождается практически в ноль по мере увеличения плотности “планетарных сгустков”



Поэтому, из практических соображений я решил воспользоваться апостериорной методикой обнаружения столкновений из DEM.


Искомый эффект пушечного ядра


Разработав первую версию симуляции, я обнаружил, что частицы, сталкивающиеся с фронтальной частью крыла, ведут себя очень знакомо. Передняя кромка крыла за счёт своей каплевидной формы при столкновении с набегающим потоком воздуха ударяется об него, расталкивая молекулы. Отраженные частицы создают ударную волну, которая отражает следующие. А те в свою очередь воздействуют на следующий слой и так далее.


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

Ударная волна создаёт щит (изображение National Geographic)


В фильме этот эффект сравнили с пушечным ядром. Летать далеко — не задача ядер. Их задача — ударять. Им не нужна аэродинамическая форма. Именно поэтому тупой нос челнока так хорошо защищает весь аппарат.


Подробнее по ссылке с таймкодом https://youtu.be/cx8XbaQNnxw?t=2206


Удивительно, но именно этот эффект отлично наблюдается в симуляции. Вероятно, именно это взаимодействие слоёв воздуха в сочетании с углом атаки даёт ту самую подъёмную силу. На тепловой карте (heatmap) давления ниже отчётливо видны области высокого (красный и оранжевый цвета) и низкого (зелёный и синий) давления. Наличие взаимодействия между частицами сделало модель Ньютона интереснее.


Картинки ниже кликабельны и доступны в высоком разрешении.



Тепловая карта давления (скалярная сумма модулей сил)



Та же карта, только при большем масштабе



Тепловая карта ускорений частиц


Архитектура симулятора


Давайте теперь рассмотрим, как устроен симулятор. Приложение состоит из двух независимых модулей:


  1. Рендерер текущего состояния.
  2. Модуль симуляции на CPU или CUDA.

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


Фазовое пространство


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


Такое представление очень удобно, потому что практически любую симуляцию можно выразить в таком виде. Движение точки состояния как правило выражается через обыкновенное дифференциальное уравнение первого порядка (Ordinary Differential Equation, ODE). Это уравнение имеет вид $dx/dt = f(x, t)$, где $x$ — это позиция точки в фазовом пространстве, а $f$ — чёрный ящик, способный определить скорость изменения состояния. Зная $x_0$ и $dx/dt$, можно посчитать следующее значение $x_1 = x_0 + \frac{dx}{dt}dt = x_0 + dx$.


Вне зависимости от того, как работает f, симуляция — это всего лишь процесс численного интегрирования.


Подробнее по теме фазового пространства можно ознакомиться в разделах 'Differential Equation Basics' and 'Particle Dynamics' курса https://www.cs.cmu.edu/~baraff/sigcourse/


На канале 3Blue1Brown также доступны отличные материалы:
https://www.youtube.com/playlist?list=PLZHQObOWTQDNPOjrT6KVlfJuKtYTftqH6


Интегратор


После различных экспериментов я решил остановиться на самом грубом, но в то же время самом простом — методе Эйлера (Forward Euler). Я пробовал использовать метод Рунге-Кутты 4-го порядка (RK4), в том числе и с адаптивным шагом, но для конкретно этого сценария больше подошёл метод Эйлера. Преимущество RK4 в том, что он позволяет делать огромные временные шаги ценой четырёхкратного увеличения вычислений, что в некоторых сценариях оправданно. В моём же случае оказалось, что я привязан к малым временным шагам, из-за необходимости избегать туннелирования частиц друг через друга. Кстати, как работают интеграторы с адаптивным временным шагом опираясь на ошибку, можно почитать в 'Differential Equation Basics' lecture notes, section 3, 'Adaptive Stepsizes'.


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


CPU-версия основной функции симулятора. GPU-версия имеет незначительные отличия.
float CSimulationCpu::ComputeMinDeltaTime(float requestedDt) const
{
    auto rad = m_state.particleRad;
    auto velBegin = m_curOdeState.cbegin() + m_state.particles;
    auto velEnd = m_curOdeState.cend();

    return std::transform_reduce(std::execution::par_unseq, velBegin, velEnd, requestedDt, [](const auto t1, const auto t2)
    {
        return std::min(t1, t2);
    }, [&](const auto& v)
    {
        auto vel = glm::length(v);
        auto radDt = rad / vel;
        return radDt;
    });
}

float CSimulationCpu::Update(float dt)
{
    dt = ComputeMinDeltaTime(dt);
    m_odeSolver->NextState(m_curOdeState, dt, m_nextOdeState);
    ColorParticles(dt);
    m_nextOdeState.swap(m_curOdeState);
    return dt;
}

Вычисление производной состояния


Теперь перейдём к сердцу симулятора — определению той самой функции f, упомянутой в параграфе “Фазовое пространство”. Ниже приведён высокоуровневый код солверов производной для CPU и CUDA версий. Стоит отметить, что CPU версия исторически появилась раньше, так как на ней было проще отладить математику. В CUDA версии появились некоторые улучшения и оптимизации, но суть осталась та же. Отличие состоит в переупорядочивании частиц. Подробнее в разделе “Реордеринг частиц”.


Высокоуровневый алгоритм расчёта производной состояния
//CPU-версия 
void CDerivativeSolver::Derive(const OdeState_t& curState, OdeState_t& outDerivative) 
{
    ResetForces();
    BuildParticlesTree(curState);
    ResolveParticleParticleCollisions(curState);
    ResolveParticleWingCollisions(curState);
    ParticleToWall(curState); 
    ApplyGravity();
    BuildDerivative(curState, outDerivative);
} 

//CUDA-версия 
void CDerivativeSolver::Derive(const OdeState_t& curState, OdeState_t& outDerivative) 
{ 
    BuildParticlesTree(curState);
    ReorderParticles(curState);
    ResetParticlesState();
    ResolveParticleParticleCollisions();
    ResolveParticleWingCollisions();
    ParticleToWall();
    ApplyGravity();
    BuildDerivative(curState, outDerivative);
}

Поиск столкновений между частицами


Как вы уже, наверное, могли понять, поиск столкновений — это краеугольный камень всех расчётов. В видео в самом начале статьи участвует 2’097’152 частиц. Среди всего этого количества нужно каким-то быстрым образом найти все сталкивающиеся пары. Интересно, но у этой проблемы нет однозначно правильного решения. Любой способ — это набор компромиссов и допущений.


Один из возможных вариантов — это использование Uniform Grid, то есть однородной сетки из ячеек на подобии шахматной доски. Одна из реализаций для GPU описана в статье “Chapter 32. Broad-Phase Collision Detection with CUDA”.



Каждая ячейка пространства содержит в себе список объектов (изображение Tero Karras, NVIDIA Corporation)


В этом случае, поиск столкновений в среднем будет занимать порядка $O(1)$. Каждой частице нужно обойти списки в 9 (3x3) или 27 (3x3x3) ячейках для 2D или 3D случая соответственно. Ещё один приятный плюс структуры — это относительная простота распараллеливания её построения. Память под списки можно выделить либо заранее в виде массива, и вычислять выходной индекс через атомарный инкремент, либо строить классический RCU lock-free односвязный список. Nvidia в своих видеокартах уже давно добавила поддержку кучи, поэтому можно вызвать malloc()/free() прямо в device коде, выделяя и освобождая элементы списков.



CppCon 2017: Fedor Pikus “Read, Copy, Update, then what? RCU for non-kernel programmers”


Однако, у этой структуры есть следующий ряд фундаментальных ограничений:


  1. Множество значений координат ограничено размером самой сетки.
  2. Близкие ячейки в евклидовом пространстве как правило расположены далеко в адресном пространстве RAM/VRAM, не разделяя единую кэш-линию, что создаёт дополнительную нагрузку на шину памяти.
  3. При низкой плотности объектов или малом их количестве структура данных начинает потреблять больше памяти, чем сами данные.
  4. Возможно появление чрезмерно длинных списков при большой плотности объектов.
  5. В связи с аппаратными особенностями планирования потоков на GPU, некоторые lock-free структуры не способны работать корректно (https://youtu.be/86seb-iZCnI?t=2311, ссылка с таймкодом).

Другой вариант, который я решил использовать в этой симуляции — это BVH-дерево на основе Z-кривой. Я наткнулся на эту довольно любопытную структуру данных, когда искал альтернативы однородной сетке.


Первая важная особенность этой структуры данных — в её основе лежит фрактальная Z-кривая, она же Кривая Мортона.



Фрактальная Z-кривая (изображение Wikipedia)



Принцип вычисления индекса на кривой — чередование битов координат
(изображение Wikipedia)


Задача этой кривой, как и любой другой space-filling curve, состоит в том, чтобы упаковать пространства высших размерностей в одномерное пространство. Если присвоить каждому объекту в 2D/3D пространстве индекс на любой такой кривой, а затем отсортировать все объекты по этому индексу, то мы увидим, что объекты, расположенные близко в геометрическом пространстве, как правило будут лежать близко и в одномерном пространстве. Это свойство позволяет существенно снизить нагрузку на шину памяти. Кстати, если вам нужно обрабатывать изображения, выполняя различные свёрточные операции и применяя фильтры, возможно, вам стоит хранить пиксели в виде одной из такой кривых, а не в виде матрицы.


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


Детально алгоритм описан в статье “Maximizing Parallelism in the Construction of BVHs, Octrees, and k-d Trees”.
Краткое изложение:


  1. Обход дерева: https://developer.nvidia.com/blog/thinking-parallel-part-ii-tree-traversal-gpu/
  2. Построение дерева: https://developer.nvidia.com/blog/thinking-parallel-part-iii-tree-construction-gpu/

В сухом остатке, алгоритм следующий. Для каждого из N объектов запускается отдельный поток, который вычисляет bounding box и, на основе его центра, вычисляет код Мортона (индекс на Z-кривой). После этого этапа боксы сортируются в порядке возрастания кода.

Формирование кодов Мортона (изображение Tero Karras, NVIDIA Corporation)


Затем инициализируются узлы дерева. В частности, ещё до построения самого дерева, известно, что для N листьев будет существовать N-1 промежуточных узлов. Соответственно, необходимые аллокации и первичные инициализации осуществляются на этом шаге. Далее наступает самый хитроумный этап. Алгоритм ищет различия в кодах, двигаясь от старших бит к младшим. Наличие разницы сигнализирует о том, что нужно сформировать промежуточный узел. Ниже на рисунках представлены однопоточная версия алгоритма и его распараллеленная версия.



Последовательный алгоритм построения префиксного дерева (изображение Tero Karras, NVIDIA Corporation)



Параллельный алгоритм построения префиксного дерева (изображение Tero Karras, NVIDIA Corporation)


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

BVH-структура (изображение Tero Karras, NVIDIA Corporation)


N потоков стартуют с листьев и, поднимаясь к корню, обновляют боксы промежуточных узлов. Так как не определено, какой из детей придёт к родителю первым, то в промежуточных узлах хранится специальный флаг, изначально установленный в ноль. Оба ребёнка с помощью атомарной функции atomicExch() устанавливают флаг в 1. Функция возвращает старое значение, которое было до модификации. Если ребёнку функция вернула 0, то значит он первый. Это также означает, что текущему потоку нельзя модифицировать бокс родителя, потому что бокс его сиблинга может быть ещё не готов. На этом этапе поток завершает своё исполнение. Если же ребёнку функция вернула 1, то можно смело модифицировать родительский бокс, объединяя боксы обоих сиблингов, и снова повторить процесс.


После этого этапа дерево готово к осуществлению запросов.


Реакция на столкновения


В симуляции существует два типа столкновений “частица-частица” и “частица-сегмент профиля”.


Реакция “частица-частица” использует факт того, все объекты уже сохранены в дереве, поэтому существует частная процедура “рефлексивного” обхода, когда листья ищут столкновения друг с другом. Эта оптимизация была предложена Tero Karras. Особенность процедуры в том, что она распознаёт столкновения A-B и B-A как одно и то же столкновение, поэтому оно детектируется только один раз. Для этого при построении дерева вводится дополнительная информация. В промежуточных узлах хранится индекс самого правого листа (rightmost leaf), до которого можно добраться. Например, на рисунке выше rightmost(N2) = 4, а rightmost(N3) = 8. Когда поток, связанный с листом, скажем, O6, будет опускаться от корня, он обратится к промежуточному узлу N2. Благодаря переменной rightmost, он увидит, что лист O6 недостижим из поддерева N2. В этом случае поток O6 должен проигнорировать всё поддерево N2. Однако, потоки, связанные с листьями из поддерева N2, будут проверять поддерево N3. В конечном итоге, если столкновение с O6 и существует, то об этом сообщит только один поток, и он будет из поддерева N2.


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


template<typename TDeviceCollisionResponseSolver, size_t kTreeStackSize>
void CMortonTree::TraverseReflexive(const TDeviceCollisionResponseSolver& solver);

Для случая “частица-сегмент профиля” используется универсальная версия:


template<typename TDeviceCollisionResponseSolver, size_t kTreeStackSize>
void CMortonTree::Traverse(const thrust::device_vector<SBoundingBox>& objects, const TDeviceCollisionResponseSolver& solver);

Здесь TDeviceCollisionResponseSolver — это объект, который должен реализовать следующий интерфейс:


struct Solver
{
    struct SDeviceSideSolver
    {
        ... 
        __device__ SDeviceSideSolver(...);
        __device__ void OnPreTraversal(TIndex curId);
        __device__ void OnCollisionDetected(TIndex leafId);
        __device__ void OnPostTraversal();
    };
    Solver(...);
    __device__ SDeviceSideSolver Create();
}; 

Для каждого тестируемого на столкновение объекта, или листа в случае рефлексивного подхода, создаётся отдельный поток. Каждый поток создаёт свой солвер через фабричную функцию Create(). Далее вызывается метод OnPreTraversal, куда передаётся индекс тестируемого объекта. Если бокс текущего тестируемого объекта перекрыл бокс какого-то листа, вызывается функция OnCollisionDetected с индексом листа. Эта функция отвечает за расчёт физики. После обхода дерева вызывается OnPostTraversal.


Такой формат разрешения коллизий появился неслучайно. С самого начала я реализовал его по-другому. Я разделил обход дерева и вычисление физики на две различные стадии, как это сделал Tero Karras. Однако я столкнулся с проблемой построения списков найденных столкновений. Я попробовал сохранять информацию о коллизиях в виде матрицы NxO, где N — количество тестируемых объектов, O — максимальный размер списка. Но я отказался от этой идеи, потому что при определенных сценариях быстро заканчивалось место в списках. А это в свою очередь создавало различные физические артефакты. К тому же я обратил внимание, что профилировщик сигнализировал о неэффективной работе с памятью (coalesced memory access). Поэтому я решил попробовать подход без списков, который был описан выше. К моему удивлению, способ оказался немного быстрее и без артефактов.


Код солвера “частица-частица”
struct SParticleParticleCollisionSolver
{
    struct SDeviceSideSolver
    {
        CDerivativeSolver::SIntermediateSimState& simState;
        TIndex curIdx;
        float2 pos1;
        float2 vel1;
        float2 totalForce;
        float totalPressure;

        __device__ SDeviceSideSolver(CDerivativeSolver::SIntermediateSimState& state) : simState(state)
        {

        }

        __device__ void OnPreTraversal(TIndex curLeafIdx)
        {
            curIdx = curLeafIdx;
            pos1 = simState.pos[curLeafIdx];
            vel1 = simState.vel[curLeafIdx];
            totalForce = make_float2(0.0f);
            totalPressure = 0.0f;
        }

        __device__ void OnCollisionDetected(TIndex anotherLeafIdx)
        {
            const auto pos2 = simState.pos[anotherLeafIdx];
            const auto deltaPos = pos2 - pos1;
            const auto distanceSq = dot(deltaPos, deltaPos);
            if (distanceSq > simState.diameterSq || distanceSq < 1e-8f)
                return;

            const auto vel2 = simState.vel[anotherLeafIdx];

            auto dist = sqrtf(distanceSq);
            auto dir = deltaPos / dist;
            auto springLen = simState.diameter - dist;

            auto force = SpringDamper(dir, vel1, vel2, springLen);
            auto pressure = length(force);
            totalForce += force;
            totalPressure += pressure;

            atomicAdd(&simState.force[anotherLeafIdx].x, -force.x);
            atomicAdd(&simState.force[anotherLeafIdx].y, -force.y);
            atomicAdd(&simState.pressure[anotherLeafIdx], pressure);
        }

        __device__ void OnPostTraversal()
        {
            atomicAdd(&simState.force[curIdx].x, totalForce.x);
            atomicAdd(&simState.force[curIdx].y, totalForce.y);
            atomicAdd(&simState.pressure[curIdx], totalPressure);
        }
    };

    CDerivativeSolver::SIntermediateSimState simState;
    SParticleParticleCollisionSolver(const CDerivativeSolver::SIntermediateSimState& state) : simState(state)
    {
    }
    __device__ SDeviceSideSolver Create()
    {
        return SDeviceSideSolver(simState);
    }
};

void CDerivativeSolver::ResolveParticleParticleCollisions()
{
    m_particlesTree.TraverseReflexive<SParticleParticleCollisionSolver, 24>(SParticleParticleCollisionSolver(m_particles.GetSimState()));
    CudaCheckError();
}

Во время отладки я обратил внимание, что при высокой плотности частиц, функция OnCollistionDetected как правило вызывается для одних и тех же аргументов среди потоков одного варпа. Типовой сценарий был следующий: если в какой-то области пространства есть частицы A, B, C и D, которые в указанном порядке расположены на Z кривой, то приблизительно происходило вот что:


lock-step # Thread #1 Thread #2 Thread #3
1 OnCollisionDetected
A <-> C
OnCollisionDetected
B <-> C
OnCollisionDetected
C <-> D
2 OnCollisionDetected
A <-> D
OnCollisionDetected
B <-> D
INACTIVE
3 OnPostTraversal(1) OnPostTraversal(2) OnPostTraversal(3)

Как видно из таблицы, на шаге 1 и 2 потоки #1 и #2 выполняли атомарные обращения atomicAdd с одними и тем же частицам C и D в процессе работы функции OnCollistionDetected. Это создаёт дополнительную нагрузку на atomic транзакции.


Начиная с архитектуры Volta, Nvidia добавила в чипы поддержку новых warp-vote инструкций. С помощью инструкции match_any поток может опросить весь warp, получив битовую маску потоков, у которых значение запрашиваемой переменной имеет такое же значение.

Результат работы match_any и match_all для двух кооперативных групп


Коммуникация внутри варпа тоже стала удобнее, потому что появились новые warp shuffle функции с поддержкой маски потоков.

Warp-wide редукция с помощью старых функций без маски


Благодаря этим функциям, потоки перед обращением в глобальную память могут сгруппироваться по признаку общего выходного адреса. Далее эта группа должна выполнить суммирование на уровне регистров SM и уже после этого только один поток обращается в глобальную память. К сожалению, на моём домашнем Pascal (1080 Ti) таких инструкций нет, поэтому я решил попробовать их проэмулировать. Увы, никакого прироста, как и замедления это не дало. Профилировка показала, что хоть нагрузка на atomic транзакции и упала в несколько раз, существенно возросла нагрузка на Arithmetic Workload и увеличилось количество регистров на поток. Заняться разработкой на чипах с Volta или Turing пока не представилось возможным. Хотя, мне всё же удалось протестировать симуляцию на RTX 2060 и найти баг, связанной с atomic операцией. Об этом в разделе “Барьер памяти”.


Другие оптимизации и дополнения


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


SoA


Structure of Arrays — одна из техник, которая позволяет ускорить доступ к памяти в определённых ситуациях.



При работе с деревом на любом из этапов как правило не требуется полный набор атрибутов. А значит вместо того, чтобы хранить каждый узел в виде структуры, всё дерево хранится в виде SoA:


typedef uint32_t TIndex; 

struct STreeNodeSoA 
{
    const TIndex leafs;

    int* __restrict__ atomicVisits; 
    TIndex* __restrict__ parents; 
    TIndex* __restrict__ lefts; 
    TIndex* __restrict__ rights; 
    TIndex* __restrict__ rightmosts; 
    SBoundingBox* __restrict__ boxes; 
    const uint32_t* __restrict__ sortedMortonCodes; 
};

То же самое касается и внутреннего состояния солвера производной состояния:


struct SIntermediateSimState 
{ 
    const TIndex particles; 
    const float particleRad; 
    const float diameter; 
    const float diameterSq; 

    float2* __restrict__ pos; 
    float2* __restrict__ vel; 
    float2* __restrict__ force; 
    float* __restrict__ pressure; 
}; 

В то же время, массив bounding box’ов нет смысла хранить в SoA стиле, потому что во всех сценариях необходимы обе точки. Поэтому боксы хранятся в виде Array of Structures (AoS):


struct SBoundingBox 
{ 
    float2 min; 
    float2 max; 
}; 

struct SBoundingBoxesAoS 
{ 
    const size_t  count; 
    SBoundingBox* __restrict__ boxes; 
}; 

Реордеринг частиц


Так как текущая реализация не строит списки столкновений, а разрешает коллизии прямо на месте, то возникает следующая проблема. После присвоения кодов Мортона центрам боксов, выполняется сортировка самих боксов. Однако остальные параметры частиц остаются неотсортированными. Если в процессе обхода дерева продолжить обращаться к данным в исходном порядке, то мы получаем uncoalesced memory access.


Такой паттерн доступа очень медленно работает на GPU. Для восстановления coalesced memory access, позиции и скорости частиц тоже упорядочиваются вдоль кривой. А после выполнения всех расчётов, силы и давления как выходные величины возвращаются к исходному порядку. Идея не нова и была позаимствована из уже упомянутого доклада SpaceX: https://youtu.be/vYA0f6R5KAI?t=1939 (ссылка с таймкодом).



Восстановление объединённого доступа к памяти (изображение SpaceX)


Такая оптимизация даёт 50% прироста производительности: с 8 FPS до 12 FPS для двух миллионов частиц.


UPD: во избежание недопониманий уточню, что между кадрами происходит 16 шагов симуляции, то есть после этой оптимизации пиковая производительность достигла 192 шагов в секунду. Зачем нужно было так делать, я развёрнуто ответил в комментарии.


Стек в Shared Memory


Оригинальная статья приводит пример реализации, где стек для обхода дерева реализуется в виде локального массива в скоупе функции. Проблема этого подхода в том, что задействуется локальная память потока — область в глобальной памяти. А значит SM начинает выполнять долгие транзакции на чтение и запись, которые ко всему прочему могут оказаться ещё uncoalesced. Суть данной оптимизации в том, чтобы использовать сверхбыструю Shared Memory на кристалле самого Streaming Multiprocessor’а.


Оригинальный код
__device__ void traverseIterative( CollisionList& list,
                                   BVH& bvh, 
                                   AABB& queryAABB, 
                                   int queryObjectIdx)
{
    // Allocate traversal stack from thread-local memory,
    // and push NULL to indicate that there are no postponed nodes.
    NodePtr stack[64]; //<---------------------------- Проблемное место
    NodePtr* stackPtr = stack;
    *stackPtr++ = NULL; // push

    // Traverse nodes starting from the root.
    NodePtr node = bvh.getRoot();
    do
    {
        // Check each child node for overlap.
        NodePtr childL = bvh.getLeftChild(node);
        NodePtr childR = bvh.getRightChild(node);
        bool overlapL = ( checkOverlap(queryAABB, 
                                       bvh.getAABB(childL)) );
        bool overlapR = ( checkOverlap(queryAABB, 
                                       bvh.getAABB(childR)) );

        // Query overlaps a leaf node => report collision.
        if (overlapL && bvh.isLeaf(childL))
            list.add(queryObjectIdx, bvh.getObjectIdx(childL));

        if (overlapR && bvh.isLeaf(childR))
            list.add(queryObjectIdx, bvh.getObjectIdx(childR));

        // Query overlaps an internal node => traverse.
        bool traverseL = (overlapL && !bvh.isLeaf(childL));
        bool traverseR = (overlapR && !bvh.isLeaf(childR));

        if (!traverseL && !traverseR)
            node = *--stackPtr; // pop
        else
        {
            node = (traverseL) ? childL : childR;
            if (traverseL && traverseR)
                *stackPtr++ = childR; // push
        }
    }
    while (node != NULL);
}

Стек в Shared Memory
template<typename TDeviceCollisionResponseSolver, size_t kTreeStackSize, size_t kWarpSize = 32>
__global__ void TraverseMortonTreeKernel(const CMortonTree::STreeNodeSoA treeInfo, const SBoundingBoxesAoS externalObjects, TDeviceCollisionResponseSolver solver)
{
    const auto threadId = blockIdx.x * blockDim.x + threadIdx.x;
    if (threadId >= externalObjects.count)
        return;

    const auto objectBox = externalObjects.boxes[threadId];
    const auto internalCount = treeInfo.leafs - 1;

    __shared__ TIndex stack[kTreeStackSize][kWarpSize]; //Тот самый стек

    unsigned top = 0;
    stack[top][threadIdx.x] = 0;

    auto deviceSideSolver = solver.Create();
    deviceSideSolver.OnPreTraversal(threadId);

    while (top < kTreeStackSize) //top == -1 also covered
    {
        auto cur = stack[top--][threadIdx.x];

        if (!treeInfo.boxes[cur].Overlaps(objectBox))
            continue;

        if (cur < internalCount)
        {
            stack[++top][threadIdx.x] = treeInfo.lefts[cur];
            if (top + 1 < kTreeStackSize)
                stack[++top][threadIdx.x] = treeInfo.rights[cur];
            else
                printf("stack size exceeded\n");

            continue;
        }

        deviceSideSolver.OnCollisionDetected(cur - internalCount);
    }

    deviceSideSolver.OnPostTraversal();
}

Использование Shared Memory позволяет достичь прироста на 43%: с 14 FPS до 20 FPS. Подробнее о доступных типах памяти можно почитать в официальной документации:


https://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#device-memory-accesses


Барьер памяти


Когда я разрабатывал симулятор, я пользовался исключительно одной видеокартой — 1080 Ti поколения Pascal. Пользоваться другими в целях разработки, к сожалению, возможности не было. Но у меня была возможность просить трёх моих знакомых запустить приложение на их игровых ноутбуках с новой на тот момент 20-й серией чипов. Все три ноутбука выдавали изображения с вот такими артефактами.



Артефакт на 20-й RTX серии. Позиции и размер артефактов каждый шаг менялись.


Сначала я думал, что это проблема визуализации, но нигде не мог найти ошибку. Проверка кода самой симуляции тоже не дала результатов. Осознание пришло через полгода после просмотра этого доклада:



Доклад об атомиках и барьерах памяти.


Половина доклада посвящена идее барьеров памяти и почему они важны при работе с atomic-операциями и lock-free структурами. Дело в том, что процессоры имеют тенденцию переупорядочивать выполнение инструкций (Out-of-order execution), но при этом отслеживая и сохраняя зависимости между ними, чтобы гарантировать корректность. В случае с lock-free структурами для процессоров зависимость не очевидна. Поэтому, нужны барьеры памяти, которые подсказывают процессору, что инструкции не могут быть переупорядочены через барьер. Каждая платформа реализует барьеры по-своему, но, к счастью, разработчики стандарта C++ смогли построить наиболее общую модель. Подробное описание каждой семантики барьеров доступно в документации по std::memory_order.


__device__ void CMortonTree::STreeNodeSoA::BottomToTopInitialization(size_t leafId)
{
    auto cur = leafs - 1 + leafId;
    auto curBox = boxes[cur];

    while (cur != 0)
    {
        auto parent = parents[cur];
        //1. Опасная atomic операция
        auto visited = atomicExch(&atomicVisits[parent], 1);
        if (!visited)
            return;

        TIndex siblingIndex;
        SBoundingBox siblingBox;

        TIndex rightmostIndex;
        TIndex rightmostChild;

        auto leftParentChild = lefts[parent];
        if (leftParentChild == cur)
        {
            auto rightParentChild = rights[parent];
            siblingIndex = rightParentChild;
            rightmostIndex = rightParentChild;
        }
        else
        {
            siblingIndex = leftParentChild;
            rightmostIndex = cur;
        }

        siblingBox = boxes[siblingIndex];
        rightmostChild = rightmosts[rightmostIndex];

        SBoundingBox parentBox = SBoundingBox::ExpandBox(curBox, siblingBox);
        boxes[parent] = parentBox;
        rightmosts[parent] = rightmostChild;

        cur = parent;
        curBox = parentBox;

        //2. Спасительный барьер памяти. 
        //Следующая итерация гарантированно увидит результаты всех записей 
        __threadfence();
    }
}

Моя ошибка была в том, что я не использовал никаких барьеров памяти в коде, который строит BVH дерево, но при этом активно использует атомарный флаг. Интересно, что оригинальная статья также не использует никаких барьеров. Скорее всего, помимо новой SIMT модели (разделы Volta SIMT Model и Starvation-Free Algorithms) Nvidia добавила в новые архитектуры начиная с Volta более агрессивную реализацию уже упомянутой Out-of-order execution. Как следствие, операции, которые должны были выполняться до atomicExch(), т.е. ещё на предыдущей итерации цикла, на Turing исполняются уже после. В результате такого агрессивного реордеринга инструкций, второй ребёнок, приходя к родителю думает, что его сиблинг уже вычислил и сохранил бокс в общую память, а на самом деле нет. В результате получается data race.


thrust::device_vector


Я слишком поздно заметил, что thurst::device_vector работает в синхронном режиме. Этот контейнер в своём конструкторе и в методе resize() выполняет полную синхронизацию с GPU через cudaDeviceSynchronize(). Видимо, разработчики руководствовались следующими рассуждениями. Раз вектор на видеокарте, то и конструкторы элементов нужно тоже вызывать на видеокарте. Но так как конструкторы могут кидать исключения, нужно дождаться их исполнения, чтобы словить все исключения. Единственный доступный способ для них — полная синхронизация. Ещё одна из обнаруженных проблем — редукция (reduce, sum, min, max) тоже синхронная, так как возвращает результат на хост. Библиотека Cuda UnBound (CUB) в этом плане куда продуманнее. Кстати, недавно она вошла в состав thrust как бэкенд, хотя раньше её нужно было скачивать отдельно.


Результаты профилировки


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



Картинка кликабельна, можно посмотреть в высоком разрешении


Заключение


Когда я брался за этот проект, я всего лишь хотел использовать GPU как мини-суперкомпьютер, чтобы проверить жизнеспособность модели Ньютона. В итоге задача оказалась куда интереснее и плодотворнее, чем ожидалось. Симуляция показала “эффект пушечного ядра”, а сама работа над проектом вылилась в исследование и долгие часы работы в APOD режиме.


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


Если вы хотели бы начать изучать CUDA, но не знаете, с чего начать, на Youtube есть отличный курс от Udacity “Intro to Parallel Programming”. Рекомендую к ознакомлению.
https://www.youtube.com/playlist?list=PLAwxTw4SYaPnFKojVQrmyOGFCqHTxfdv2


Напоследок, ещё несколько видео симуляций:



CPU-версия, 8 потоков, 131'072 частиц



CUDA-версия, 4.2М частиц, 30 минут симуляции

AdBlock has stolen the banner, but banners are not teeth — they will be back

More
Ads

Comments 46

    +3

    С барьерами памяти есть одна небольшая подстава: на каждой архитектуре они работают по-разному. Так, для x86 справедливо следующее:


    Loads and stores are not reordered with locked instructions

    То есть любая атомарная операция на x86, кроме load/store, которые представляют собой обычный MOV — это неявный барьер. Отсюда получается, что в многопоточном коде под x86 барьеры практически никогда не используются, а компилятор C++ генерирует одинаковый код вне зависимости от используемой модели. Бывают, конечно, редкие исключения: запись в память, минуя кэш, или же синхронизация чисто на load/store без блокировки шины. Но это всё экзотика.


    К чему я клоню: под x86 можно годами писать параллельный код, вообще не зная о существовании барьеров. А вот с GPU, увы, уже так не получится.

      +1
      Именно! Сам впервые столкнулся на практике именно на видеокартах Turing. После того, как отстрелил себе ногу, начал изучать вопрос детально.

      Выяснил, что на ARM'ах, с которыми приходится работать уже на основной работе, тоже есть опасность подобных ситуаций. К счастью, стандарт C++ автоматически защищает несведущих людей, используя по-умолчанию медленный, но безопасный std::memory_order_seq_cst.

      Два доклада Фёдора Пикуса многое прояснили в этой области. Поэтому я решил ими поделиться.
        +1
        К счастью, стандарт C++ автоматически защищает несведущих людей, используя по-умолчанию медленный, но безопасный std::memory_order_seq_cst.

        Он защищает только при использовании типов std::atomic<>. Берём любую самодеятельность: да хоть atomicExch или же interlocked API над интовыми указателями — и всё ломается.

      +3
      Простой вопрос — а вы вообще проверили, что ваш код решает задачу газовой динамики корректно? Скорости звука там правильные, зависимости от направления нет, вот-это-все? Прогоните хотя бы через набор тестов Athena — www.astro.princeton.edu/~jstone/Athena/tests — будет гораздо яснее, на каком свете вы находитесь.
        0
        Корректное моделирование газов не является целью данной симуляции. Она не может быть использована для реальных расчётов. Для этих целей нужно использовать методы из CFD и соответствующие программные пакеты (Solidworks, Autodesk Airflow). Поэтому, модель не пройдёт вышеуказанные тесты.

        Целью же является проверка правдоподобности модели подъёмной силы Ньютона. А именно, будут ли наблюдаться какие-либо интересные эффекты, если взять большое количество шариков и профиль крыла.
          +2
          youtu.be/l-oRMQwQjFw
          на основе шариков и плавают, примерно на тех же принципах как Ваше крыло летит )
          Красота картинки важна.
            0
            Стоило бы добавить на симуляцию вектор подъемной силы и попробовать изменять угол атаки для большего эффекта. Неламинарный поток вы же как то симулировали.
              0
              О какой правдоподобности можно говорить, если вы моделируете ненастоящую физику? При чем тут тег «физика» тогда вообще, если вы про производительность и графику?
              Способность нарисовать картинку, похожую на реальность, она похвальна и осмысленна — в тех же компьютерных играх, к примеру. Но физика, вроде бы, все же претендует на описание настоящей реальности…
            +2
            Глаза болят от такого неграмотного описания вычислительной гидрогазодинамики.
            Температурная карта давления
            Это финиш… Она температурная или давления? Какая такая карта, у вас просто картина распределения в пространстве. Нормально это называть можно так: «Контуры давления», «Распределение поля давления», «Распеределение давления», «Поле давления». Легенды где?
            Один из многих способов численно решить систему уравнений — это применить Метод Конечных Элементов в связке с Адаптивным Сеточным Методом в случае подхода Эйлера.
            Не вводите читателей в заблуждение, основным, самым используемым способо решения является метод конечных объемов. SpaceX конкретно не раскрыли что они использовали, но судя по всему это был метод конечеых объемов или конечно разностный метод.
            Эта система определяет равнодействующую всех сил (давления, вязкости и гравитации) для всех направлений. Дополнительно может быть задано второе векторное уравнение, обеспечивающее условие несжимаемости, если нужно описать жидкости

            Только одно уравнение из уравненИЙ Навье-Стокса это уравнение баланса импульса. «Дополнительно», дополнительно блин, уравнение неразрывности с каких пор у нас стало дополнительным?! При чем здесь описание несжимаемости?! Как вы замыкаете эти два уравнения? Где уравнение состояния жидкости?
            Для иных сред, чаще всего реализуют подход Лагранжа, через метод сглаженных частиц (Smoothed-particle hydrodynamics, SPH).
            Да, чаще всего для рисования красивых картинок, не более того, у него огромные проблемы с постановкой корректных ГУ и с сохранением массы. Все эти качественные красивые эффекты (брызги пена и тд) туфта, когда они появляются не четко тогда, там и таким образом как в реальном процессе.
            Для приемлимых результатов с Лагранжевым подходом используют метод FPM — метод конечных точек.

            И самое интересное, почему в явно сверхзвуковом режиме течения вокруг крылового профиля у вас есть возмущения выше по потоку в расчете? Какие скорости течения вы использовали, как замыкали уравнения навье стокса? Какой жидкости соответствует ваша модель?
            Пока на физику это не похоже, на красивые мультики да.
            Самый главный вопрос, где подъемная сила?!

              +2
              Корректное моделирование газов не является целью данной симуляции. Целью же является проверка правдоподобности модели подъёмной силы Ньютона. А именно, будут ли наблюдаться какие-либо интересные эффекты, если взять большое количество шариков и профиль крыла.
                +2
                Так где же правдоподобность подъемной силы? Где расчет подъемной силы? Где сравнение для доказательства правдоподобности? Где вообще подъемная сила? В заголовке?
                Несоответствие заголовка и содержания. Я зашел почитать про расчет подъемной силы!
                Назовите статью тогда «Рисование красивых картинок навеянных физикой методом частиц на CUDA»
                «Интересные эффекты», визуальные эффекты?
                  0

                  А какие вы параметры брали? Соотношение скорости звука и скорости крыла и размера крыла точно реалистичные? Потому что сверх звуковой полёт это уже другое.

                  +7
                  «Температурная карта» — это вроде устойчивое название такой палитры… ОНа может быть чего угодно — «температурная карта глубин» — тоже не имеет никакого отношения к температуре — просто чем ближе, тем как бы теплее цвета. В данном случае чем больше давление — тем «теплее» цвет… В чем проблема?
                    +9

                    Проблема в том, что устойчивое название такой палитры — "тепловая карта", а не "температурная карта".

                      +1
                      Замечание дельное, исправил. Спасибо.
                        –2
                        Устойчивое название такой палитры «Rainbow», «радуга», в любом модуле ANSYS, COMSOL, SolidWorks, Numeca, LS-Dyna список можно продолжать. Палитры где цвет «теплее» совершенно другие. Палитра вообще не должна определять название таких картинок, хоть градиентом от черного к белому заливай.
                        Проблема только в том, что за общепринятый стандарт начали выдавать какой-то бред, вылезший из каких-то людей видевших то ли изображения с тепловизора то ли расчеты теплообмена и решивших, что это всё про теплоту, даже не температуру.
                          +2
                          Вам всё правильно объяснили, тепловая карта или heatmap — устойчивое выражение, в статье оно употребляется корректно.
                          А про «радугу» впервые слышу, в ваших же перечисленных пакетах, также используется понятие heatmap: www.ansys.com/blog/heat-maps-determine-build-orientations-additive-manufacturing
                            +1
                            Человек выше написал
                            «Температурная карта» — это вроде устойчивое название такой палитры…
                            «Палитры», понимаете палитры, палитра называется «радуга».
                            Высосаный из пальца источник, где heat map применяется для описания отклика от двух параметрво угла наклона. В практике ANSYS никогда не используется heat map для описания пространственных распределений параметров, нате вам ссылку тоже: www.ansys.com/products/fluids/ansys-cfx просто посмотрите как называется распределение давления на картинке.
                            Я против размывания вполне конкретных устоявшихся терминов обеспечивающих максимальную корректность. Я вот работаю в сфере гидрогазодинамики и термодинамики, и если мои коллеги по цеху начнут писать «heat flux heat map», то я очень расстроюсь, каша в словах — каша в головах. Автор, очевидно, не специалист в CFD, и не знаком с физикой, а вы его защищаете, зачем?
                              0
                              Потому что «тепловая карта» — плохой перевод термина «heatmap», но понятно, что имел ввиду автор. Остальные ваши комментарии по делу, но вот с heatmap вы придираетесь.
                                0
                                Ну может быть, накипело. У меня ощущение, что программеры лезут со своей терминологией и поверхностным пониманием во все дыры уже, я уже вакансии физика и инженера не могу на linkedin найти, потому что программеры уже стали и физиками и инженерами и архитекторами и черте ещё чем.
                    +4
                    Однако, у этой структуры есть следующий ряд фундаментальных ограничений:

                    1. Множество значений координат ограничено размером самой сетки.


                    нет, не ограничено, для бесконечного пространства координата приводится к ячейке как остаток от деления на ширину/высоту сетки, играя роль хеша, а при обработке столкновений дополнительно проверяется реальное расстояние между парой частиц-кандидатов

                    насчет остальных аргументов — большинство также вызывают скепсис:

                    интересно сравнить конкретные реализации по производительности

                    вы пишите про 8-12 FPS на 1080 Ti на 2млн частиц, у меня получатся считать на ~50-60 FPS это количество на таком железе, можете попробовать запустить демо-версию, там есть практически 1:1 такая же простая симуляция-коллайдер:

                    store.steampowered.com/app/1196080/Space_Simulation_Toolkit

                    работает на основе inclusive_scan, без fast_radix_sort, с обычной сеткой, это самый быстрый способ из всех, что я исследовал:

                    on-demand.gputechconf.com/gtc/2014/presentations/S4117-fast-fixed-radius-nearest-neighbor-gpu.pdf
                      +1
                      да, насчет пересортировки ради coalescing — при хорошей оптимизации не окупится:

                      применение вообще любой сортировки на каждом кадре симуляции делает построение сетки самой дорогой операцией (примерно x10 по отношению ко всему остальному коду симуляции)

                      именно поэтому fluids_v3 алгоритм такой крутой, сетку можно строить вообще без единой сортировки, за счет prefix_sum (inclusive_scan)
                        +2
                        www.youtube.com/watch?v=ISluqewh9jQ

                        вот 400K частиц на 780 Ti, даже без VBO_INTEROP оптимизации рендера
                          +3
                          вы пишите про 8-12 FPS на 1080 Ti на 2млн частиц, у меня получатся считать на ~50-60 FPS это количество на таком железе


                          Тут довольно тонкий момент. Процитирую один фрагмент:
                          Модули слабо связаны друг с другом и напрямую не зависят. Существует возможность проводить симуляцию в оффлайн режиме и рендерить результаты позднее. В CUDA-версии демонстрации рендеринг, что на видео в начале статьи, происходит каждые 16 шагов.


                          То есть эти 12 FPS нужно расценивать как 192 шага в секунду. Причина по которой я рисую не каждый шаг, а через 16 шагов в том, что передача данных между контекстами CUDA и OpenGL занимает существенное время. Я тем самым уплотняю вычисления и минимизирую задержки на передачу данных через PCI-E шину.

                          Мне известно о наличии OpenGL Interoperability в CUDA API, но конкретно в этом проекте я хотел получить возможность рендерить на другом видеоадаптере. Поэтому модуль, занимающийся расчётами, хоть и использует видеоадаптер для вычислений, никак не использует его для отрисовки.

                          Если вдруг интересно, предлагаю ознакомиться с вот этим файлом:
                          github.com/tony-space/WingSimulator/blob/master/Simulation/libCudaSimulation/CSimulationCuda.cu

                          А именно, вот этот фрагмент:
                          float CSimulationCuda::Update(float dt)
                          {
                          	for (int i = 0; i < 16; ++i)
                          	{
                          		m_odeSolver->NextState(ComputeMinDeltaTime(dt), m_curOdeState, m_nextOdeState);
                          		m_nextOdeState.swap(m_curOdeState);
                          	}
                          	
                          	ColorParticles2();
                          	return dt;
                          }
                          

                          +2
                          Основной недостаток этой модели — это отсутствие взаимодействия частиц газа друг с другом

                          Это потому что межатомное взаимодействие проявляется на расстояниях до 6-8 ангстрем, а при комнатной температуре и атмосферном давлении среднее расстояние между частицами газа порядка 30 ангстрем.

                          У меня некоторый диссонанс остался, вначале вы пишите о симуляции из отдельных атомов, а потом говорите о конечных элементах и прочих вещах, где единицей является какая-то более обобщенная структура. Но за статью спасибо, узнал кое-что интересное, например о match_any, может быть удастся применить с пользой.
                            +1
                            Тут на первой же картинке очевидно, что дело идет о сверхзвуковом обтекании — ясно видна ударная волна. Судя по углу ее конуса (на глаз), там почти 3M. Так что частицам приходится взаимодействовать.

                            А как выглядит этим методом дозвуковой режим?

                            Фраза
                            Основной недостаток этой модели — это отсутствие взаимодействия частиц газа друг с другом. Из-за этого при нормальных условиях она даёт некорректные результаты, хотя всё ещё может применяться для экстремальных условий, где взаимодействием можно пренебречь.

                            прямо пугает. На мой взгляд, модель идеального газа — это и есть нормальные условия, ничего экстремального.
                              0
                              Тут на первой же картинке очевидно, что дело идет о сверхзвуковом обтекании — ясно видна ударная волна. Судя по углу ее конуса (на глаз), там почти 3M. Так что частицам приходится взаимодействовать.

                              А как выглядит этим методом дозвуковой режим?


                              Не стоит воспринимать модель подъёмной силы Ньютона (Newtonian Lift) как адекватную модель для точных газодинамических расчётов. На сайте NASA, ссылка в аннотации перед катом, подробно описано почему она некорректна. Это всего лишь грубая модель взаимодействия газа с профилем крыла.

                              На мой взгляд, модель идеального газа — это и есть нормальные условия, ничего экстремального.


                              Некоторые модели идеального газа допускают эластичные столкновения между частицами. Главное чтобы не было дальнодействующих сил.
                              en.wikipedia.org/wiki/Ideal_gas#Microscopic_model
                            0
                            Спасибо за статью!
                            Сам думаю добраться до лагранжевых методов, но пока всё руки не доходят.
                            Позволю пару придирок.
                            «Метод Конечных Элементов» — вряд ли заслуживает больших букв — и даже Википедия на это намекает. Большие буквы в научно-популярной статье сразу пугают.
                            По поводу тепловых/температурных карт — что-то ни разу не слышал такого жаргонного употребления (а на конференциях по гидродинамике я регулярно бываю). Все говорят — «поле скорости» или «распределение давления». Возможно это некий внутренний жаргон? Но может я просто из другой песочницы — век живи — век учись.
                            По поводу уравнений Навье-Стокса. Эта система определяет равнодействующую всех сил (давления, вязкости и гравитации) для всех направлений.
                            Вообще-то это уравнение сохранения импульса для жидкого объема. Гравитация там слегка лишняя — в том смысле, что там сумма сил куда может входить всё что угодно, например электромагнитные силы.
                            По поводу МКЭ. Еще в 2007 году, Брайен Сполдинг призывал всех отказаться от МКЭ в пользу МКО. Что в общем-то и происходит — пакеты которые популярны используют МКО. Научная среда — другая история, там у каждого свой код.
                            По поводу подъемной силы:
                            «На тепловой карте (heatmap) давления ниже отчётливо видны области высокого и низкого давления.» Ну как бы не очень. Видны цветные области — что они значат — «догадайся мол сама». Но не все читатели регулярно видят такие картинки — желательно пояснить.
                            Как у вас взаимодействуют частицы? Что происходит при столкновении? Твердые сферы? Возможно я проглядел этот момент.
                            За урок по Куде — отдельное спасибо.
                              +1
                              Всё у вас в порядке с научным языком и мнение адекватное состоянию реального CFD, не забивайте голову.
                              Из Лагранжевых методов советую finite point method.
                                0

                                Кстати, а какая сейчас ситуация с «гибридными» методами вроде метода крупных частиц Белоцерковского и Давыдова? Используют сейчас или не прижились?

                                  +1
                                  Это что-то очень древнее. Не, ничего такого не используют, может кто-то где-то в какой-то лабе и то вряд ли. Не в индустрии точно. Вы мне скажите, какая задача, я скажу, как обычно или редко считают.
                                  Развитие методов с Лагранжевым подходом идет в задачах со сложной реологией жидкостей (и не совсем жидкостей), где вязкости в одной области могут отличаться на многие порядки, стекло, пластики, паста.
                                  Гибрид Лагранжева и Эйлерова подхода применяется для многофазных течений, для гомогенных сред нет.
                                    0

                                    Ну, меня не столько метод крупных частиц интересовал, сколько сам гибридный подход. Просто я со времён университета этим не занимаюсь, и он первым на ум пришёл. Да и даже тогда мы его только в учебных целях программировали, а для реальных расчётов, я помню, брали метод Годунова и его модификации вроде PPM. Но то были астрофизические расчёты, там не было каких-то сложных граничных условий (зато было много других проблем).

                                0
                                «Метод Конечных Элементов» — вряд ли заслуживает больших букв — и даже Википедия на это намекает. Большие буквы в научно-популярной статье сразу пугают.

                                Замечание дельное, исправил. Спасибо.

                                По поводу тепловых/температурных карт — что-то ни разу не слышал такого жаргонного употребления (а на конференциях по гидродинамике я регулярно бываю). Все говорят — «поле скорости» или «распределение давления». Возможно это некий внутренний жаргон? Но может я просто из другой песочницы — век живи — век учись.

                                Я знаком с термином прежде всего из англоязычных источников (Heat map). Не исключаю, что в узком профессиональном сообществе приняты другие названия.

                                По поводу уравнений Навье-Стокса. Эта система определяет равнодействующую всех сил (давления, вязкости и гравитации) для всех направлений.
                                Вообще-то это уравнение сохранения импульса для жидкого объема. Гравитация там слегка лишняя — в том смысле, что там сумма сил куда может входить всё что угодно, например электромагнитные силы.


                                Опять же, знаком с уравнениями из англоязычных источников, где они записывается в виде:
                                image. Допускаю, что в русскоязычной литературе принята другая форма. Из того, что я в своё время понял, это просто сумма сил, действующих на элементарный объём dV массой dm. Собственно поэтому у нас и фигурирует плотность как отношение этих двух бесконечно малых величин. Если обе части мы домножим на dV, мы получим силу в Ньютонах, правда бесконечно малую. Поправьте, если я что-то не правильно понял.

                                По поводу МКЭ. Еще в 2007 году, Брайен Сполдинг призывал всех отказаться от МКЭ в пользу МКО. Что в общем-то и происходит — пакеты которые популярны используют МКО. Научная среда — другая история, там у каждого свой код.

                                Замечание разумное. В сеточные методы не углублялся, так как больше интересуюсь бессеточными, поэтому вполне могу обладать устаревшей информацией. Просто интересно, вырождается ли в двумерном случае метод конечных объёмов в МКЭ?

                                «На тепловой карте (heatmap) давления ниже отчётливо видны области высокого и низкого давления.» Ну как бы не очень. Видны цветные области — что они значат — «догадайся мол сама». Но не все читатели регулярно видят такие картинки — желательно пояснить.

                                Добавил текстовую расшифровку в указанное предложение. Цветовую легенду всё-же не стал прикладывать, потому что в достаточно грубой модели Newtonian Lift значения давления и ускорений являются крайне не точными. Труъ-физикам они покажутся чересчур неадекватными.

                                Как у вас взаимодействуют частицы? Что происходит при столкновении? Твердые сферы? Возможно я проглядел этот момент.

                                Частицы ведут себя как неупругие шарики, примерно так же, как и в DEM. В частности демпфирование пришлось добавить для улучшения численной стабильности. Можете посмотреть как используется функция SpringDamper: github.com/tony-space/WingSimulator/blob/master/Simulation/libCudaSimulation/CDerivativeSolver.cu
                                  0
                                  Ага.
                                  Кажется я понял в чем дело.
                                  «Тепловая карта» — относительно новый термин.(С) Википедия.
                                  В исследованиях, гидродинамика и теплообмен часто идут в одном флаконе. И на соседних картинках может быть и распределение температуры и скажем распределение давления. И если называть это одним словом то… председатель секции (или рецензент) будет недоволен.
                                  Насчет Навье-Стокса — суть-то не в силах, а в скорости которую можно извлечь из ускорения. В каком нибудь Навье-Стоксе редуцированном до уравнения Эйлера вообще с силами грустно. Ну я думаю, вы всё правильно понимаете — другой вопрос, что написать популярную статью подобрав точные фразы для всего вовсе не просто.
                                  Так, а что такое "(скалярная сумма модулей сил)"?. Давление это сила делённая на площадь. Или сила — градиент давления. У вас должна быть средняя сила за какой-то период осреднения отнесенная к элементарной площадке — что касается распределения давления по крыловому профилю.
                                  Если про диссипацию и потенциал взаимодействия напишете подробнее — будет интересно. В той программе я вряд ли разберусь. Собственно, вся картина и определяется же потенциалом взаимодействия и количеством частиц.
                                    +1
                                    Как ни странно, суть и не в том, что вы говорите. Суть в сохранении импульса, собственно из названия ясно, а выражено это в том, что работа сил вязкости полностью компенсируется работой давления (статического давления). По простому с оговорками высокое давление жидкость толкает в сторону низкого, а вязкость этому сопротивляется. Вроде достаточно простые фразы для популярной статьи, не так ли? Написать хорошую, простую к пониманию, популярную статью может только человек который очень хорошо разбирается в сложном.
                                    +1
                                    Опять же, знаком с уравнениями из англоязычных источников, где они записывается в виде:
                                    image. Допускаю, что в русскоязычной литературе принята другая форма.

                                    Это шутка такая что-ли? В любой, абсолютно любой (русско-, англо-, немецко-, китайско- итд язычной) литературе минимальный набор уравнений Навье-Стокса представлен двумя дифференциальными уравнениями сохранения импульса и сохранения массы или неразрывности (без разложения по направлениям)! Вы тут пишите одно. У меня вопрос, если вы не способны понять английский текст даже в той же википедии, может вам по русски всё же читать с начала?
                                  0
                                  Вылетает ваша программка, на line 179, CSimulationCuda.cu, вызове m_nextOdeState.resize(particles * 2);
                                  Ошибка «Run-Time Check Failure #3 — The variable 'result' is being used without being initialized.»
                                  CUDA версии 11.0, VS 2019. Сэмплы CUDA собираются и работают «искаропки».
                                  Stack trace:
                                  > libCudaSimulation.dll!thrust::detail::allocator_traits_detail::system<thrust::device_allocator>(thrust::device_allocator & __formal) Line 251 C++
                                  libCudaSimulation.dll!thrust::detail::allocator_system<thrust::device_allocator>::get(thrust::device_allocator & a) Line 360 C++
                                  libCudaSimulation.dll!thrust::detail::allocator_traits_detail::copy_construct_range<thrust::cuda_cub::tag,thrust::device_allocator,thrust::detail::normal_iterator<thrust::device_ptr>,thrust::device_ptr>(thrust::execution_policy<thrust::cuda_cub::tag> & from_system, thrust::device_allocator & a, thrust::detail::normal_iterator<thrust::device_ptr> first, thrust::detail::normal_iterator<thrust::device_ptr> last, thrust::device_ptr result) Line 218 C++
                                  libCudaSimulation.dll!thrust::detail::copy_construct_range<thrust::cuda_cub::tag,thrust::device_allocator,thrust::detail::normal_iterator<thrust::device_ptr>,thrust::device_ptr>(thrust::execution_policy<thrust::cuda_cub::tag> & from_system, thrust::device_allocator & a, thrust::detail::normal_iterator<thrust::device_ptr> first, thrust::detail::normal_iterator<thrust::device_ptr> last, thrust::device_ptr result) Line 291 C++
                                  libCudaSimulation.dll!thrust::detail::contiguous_storage<float2,thrust::device_allocator>::uninitialized_copy<thrust::detail::normal_iterator<thrust::device_ptr>>(thrust::detail::normal_iterator<thrust::device_ptr> first, thrust::detail::normal_iterator<thrust::device_ptr> last, thrust::detail::normal_iterator<thrust::device_ptr> result) Line 284 C++
                                  libCudaSimulation.dll!thrust::detail::vector_base<float2,thrust::device_allocator>::append(unsigned __int64 n) Line 877 C++
                                  libCudaSimulation.dll!thrust::detail::vector_base<float2,thrust::device_allocator>::resize(unsigned __int64 new_size) Line 327 C++
                                  libCudaSimulation.dll!wing2d::simulation::cuda::CSimulationCuda::CopyToGPU() Line 179 C++
                                  libCudaSimulation.dll!wing2d::simulation::cuda::CSimulationCuda::ResetState(const wing2d::simulation::SimulationState & state) Line 113 C++
                                  WingSimulator.exe!SetupState(wing2d::simulation::ISimulation * simulation, std::vector<std::tuple<float,float>,std::allocator<std::tuple<float,float>>> && airfoil) Line 117 C++
                                  WingSimulator.exe!main(int argc, char * * argv) Line 132 C++

                                    0
                                    Недавно сам пытался поставить новую версию CUDA 11.0. Столкнулся со схожей проблемой, и с проблемой профилировки на Pascal. Пришлось делать даунгрейд. Если есть желание, можете попробовать собрать старым SDK: developer.nvidia.com/cuda-10.2-download-archive

                                    Библиотека thrust в очередной раз расстраивает.

                                    P.S. можете зайти в раздел релизов на гитхабе и попробовать запустить бинарник через заранее заготовленные BAT файлы.
                                      0
                                      Ага, спасибо, скачал. Через CPU работает, через CUDA — нет, падает с ошибкой «parallel_for failed: cudaErrorNoKernelImageForDevice: no kernel image is available for execution on the device». Запускал на Thinkpad P71, с двумя графическими адаптерами, Intel HD Graphics 530 и NVIDIA Quadro M3000M.
                                        0
                                        Такая проблема может быть связана со старыми драйверами. CUDA 10.2 требует, чтобы стояла версия от 441.22.

                                        Проверить текущую версию можно в «Панель управления NVIDIA» -> «Справка» -> «Информация о системе»:
                                          0
                                          Да фиг с ним, видео я глянул — красиво! Непонятно, все-таки, почему Ваш код такой «привередливый»; другая красивая open source демка с частицами, базирующаяся на CUDA github.com/ramakarl/fluids3, работает без проблем.

                                          P.S. А драйвер у меня референсный, версии 27.21.14.5182 от 7/19/2020.
                                            0
                                            Нашёл вероятную причину. Ваша видеокарта поддерживает Compute Capability 5.0. В настройках проекта стоит поддержка 6.1. Можете попробовать поставить опцию в настройках «compute_50,sm_50»

                                              0
                                              Попробовал, не помогло (сделал rebuild solution, для верности). Впрочем, у меня падает на ошибке из моего первого сообщения, так что это не решение однозначно. Можно, конечно, попробовать переинсталлировать CUDA. Ну, или, пересоберите у себя с этой опцией и выложите отдельный релиз, может, заработает (обещаю проверить).

                                              P.S. Если хотите, можем перейти на гитхаб для обсуждения, тут не очень подходящее место.
                                    0
                                    интересно написано, увлекло. Подписался даже
                                      0
                                      Такой бы харей, да в powder поиграть
                                        0
                                        В статье про то как летают самолёты, есть видео потоков воздуха вокруг крыла и там поток по верхней кромке идёт вниз и он ускоряется помимо всего прочего, в вашей же симуляции частицы, в конечном итоге, идут вверх.
                                        Насколько я понял статью про подъёмную силу крыла, то там присутствуют такие эффекты:
                                        1. Отрицательного давления не бывает, поэтому верхняя кромка крыла не может тянуть вверх, а только не препятствовать этому, соотвественно, подъёмной силой обладает только нижняя плоскость и как следствие — угол атаки, это решающий фактор.
                                        2. На нижнюю кромку давят частицы не один раз, как в Ньютоновской модели, а переотражаются многократно между крылом и набегающим потоком и фактически, почти все потоки, кроме кромочного проходят мимо, сталкиваясь с подушкой, а не крылом, как было бы в Ньютоновской модели. Но куммулятивный эффект равнозначен движению в более плотной среде, это очень сильно увеличивает эффект подъёмной силы.
                                        3. Если происходит срыв потока, то через переднюю кромку выдавливается против движения подушка в верхнюю часть и переотражения нарушаются, что приводит к резкому падению подъёмной силы. Особой формой выгнутости крыла, можно заполнить область пониженного давления и отсрочить тем самым предавливание потока, при незначительном падении подъёмной силы (если форма верхней кромки будет совпадать по форме с областью низкого давления, то это давление не особо возрастёт), поэтому добиться больших углов атаки без данного эффекта. Поэтому же самолёты легко летают вверх ногами, просто больше рискуют свалиться в штопор.

                                        Only users with full accounts can post comments. Log in, please.