Pull to refresh
28
0
Anton Vasin @tony-space

Software Engineer

Send message

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

Проиллюстрирую проблему на мысленном эксперименте. Есть три шара A, B и C. Шар A и B до момента разрешения коллизий пересекаются друг с другом. C находится в свободном состоянии. После обнаружения коллизии A и B, и ёё разрешения может оказаться, что теперь A или B сталкивается с C. Хотя изначально этой коллизии не было. И это проблема. Причём результат ещё и зависит от того, в каком порядке мы обходим объекты: A -> B -> C или C -> B -> A.

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

Вообще, есть ну очень хороший и доступный материал от Карнеги-Меллона, рекомендуемый к ознакомлению.
https://www.cs.cmu.edu/~baraff/sigcourse/

То что я описал выше подробно описано в лекциях Differential Equation Basics и Particle  Dynamics.

Там есть и про идеальные связи, которые вы пытаетеесь сделать. Лекция Constrained Dynamics.

Пришла в голову такая мысль.

А что если мы скрестим JPEG с вариационным исчислением?

Шаг 1. JPEG использует дискретное косинусное преобразование (DCT) для вычисления аналитической аппроксимации g(x,y) дискретной функции f(x,y), где значения обоих функций -- это яркость пикселя. Иначе говоря из матрицы пикселей f получаем функцию g из суммы косинусов. Большую часть этих косинусов можно выкинуть за ненадобностью. Хорошая новость в том, что в DCT косинусы не считаются для каждой точки сетки -- только для строк и столбцов (O(2N) вместо O(N^2)).

Шаг 2. Зная аналитически заданную поверхность g(x,y), можно аналитически посчитать кратчайший путь методом вариационного исчисления. Иначе говоря -- найти геодезическую линию на поверхности.

Интересно, что первый шаг -- это просто матричное умножение. Второй шаг чуть сложнее. Для простых кейсов -- это может и пушкой по воробьям, но предположу, что на огромных сетках, где есть плавные перепады высот, это может дать выигрыш. Пусть и ценою некоторой потерей точности. Бонусом также можно накладывать аналитические ограничения.

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

Девиртуализация работает и в контекстах когда компилятор точно не уверен. И механизм прост как дверь: https://youtu.be/w0sz5WbS5AM?t=3088

Я подчёркиваю главную мысль в первом и последнем абзацах своего изначального комментария: автор пытается приложить столько усилий оптимизировать то что тормозить будет в последнюю очередь. Небольшой "потенциальный" оверхед виртуальных функций погоды не сделает, зато код станет портируемым, без UB и понятным другим разработчикам для сопровождения.

Принимая во внимание контекст задачи, где мы пытаемся абстрагироваться от реализации системы рендеринга, то вирутальные функций -- это как раз тот самый инструмент, который нужно здесь использовать. Этот инструмент был придуман ровно для этой задачи и подобным им. И яркий тому пример этого -- контекст в любой версии DirectX
https://docs.microsoft.com/en-us/windows/win32/api/d3d11/nn-d3d11-id3d11devicecontext

Имея абстракный интерфейс на руках:
1. Можно полностью избавиться от define, переложив создание объектов на абстракную фабрику.
2. Подключение реализации фабрики и реализации интерфейса можно разрулить через CMake: если собираем под Windows, подключаем одни CPP файлы, если под Linux, то другие.
3. Нет необходимости городить кастомные аллокаторы / деалокаторы. Если очень уж нужно, то всё это можно переложить на фабрику, и так называемый deleter умных указателей.

Вызов виртуальных функций на современных десктопных и мобильных CPU ничего уже не стоит, потому что branch-prediction, instruction cache секции кода, и оптимизиции компиляторов по девиртуализации.

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

У Visual Studio, к примеру, сборщиком по умолчанию является MSBuild, а у CLion - CMake.

CMake не является системой сборки. Это генератор скриптов для системы сборки. Если вы используете Linux, то CMake по-умолчанию будет генерировать Makefile'ы. Если же вы сидите на Windows и у вас стоит VS, то CMake по-умолчанию генерирует скрипты для MSBuild. На маках так вообще две нативные системы: Makefiles и XCode (не знаю как он под капотом работает, может тоже через Makefiles).

Справедливо и другое: CLion может вызывать MSBuild, а Visual Studio может работать с Ninja.

Короче, гуглите про ключ -G для CMake.

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

void foo (unsigned x=0);
int y = 42;
foo(y); //warning: signed/unsigned mismatch

#include <iostream>

#include <optional>

auto generate(int from, int to)

{

    return [=, cur = from]() mutable -> std::optional<int>

    {

        if(cur == to)

            return std::nullopt;

        return cur++;

    };

}

int main()

{

    auto gen = generate(0, 10);

    while(auto val = gen())

    {

        std::cout << *val << std::endl;

    }

    return 0;

}


Это конечно не труъ-генераторы, но симулировать их поведение в ранних стандартах можно через лямбды (строго говоря, через функторы; лямбы -- сахар над ними)

Автор оригинала забыл опустить наверное самую важную особенность лямбд, введённую с C++14 -- это generic lambda. В сочетании с fold expression (c++ 17) можно теперь не городить лютые шаблоны, а делать всё куда лакончинее и проще:

#include <iostream>

#include <string>

int main()

{

    auto sum = [](auto x, auto y, auto... args)

    {

        return ((x + y) + ... + args);

    };

    //prints: 8

    std::cout << sum(3, 5) << std::endl; 

    //prints: hello world one two 

    std::cout << sum(std::string{"hello "}, "world ", "one ", "two") << std::endl;

    return 0;

}

Фактически, лямбды стали синтаксическим сахаром над шаблонными функторами.

Нашёл вероятную причину. Ваша видеокарта поддерживает Compute Capability 5.0. В настройках проекта стоит поддержка 6.1. Можете попробовать поставить опцию в настройках «compute_50,sm_50»

Такая проблема может быть связана со старыми драйверами. CUDA 10.2 требует, чтобы стояла версия от 441.22.

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

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

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

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

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

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

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


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

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

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

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

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

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

Частицы ведут себя как неупругие шарики, примерно так же, как и в DEM. В частности демпфирование пришлось добавить для улучшения численной стабильности. Можете посмотреть как используется функция SpringDamper: github.com/tony-space/WingSimulator/blob/master/Simulation/libCudaSimulation/CDerivativeSolver.cu
Тут на первой же картинке очевидно, что дело идет о сверхзвуковом обтекании — ясно видна ударная волна. Судя по углу ее конуса (на глаз), там почти 3M. Так что частицам приходится взаимодействовать.

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


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

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


Некоторые модели идеального газа допускают эластичные столкновения между частицами. Главное чтобы не было дальнодействующих сил.
en.wikipedia.org/wiki/Ideal_gas#Microscopic_model
Замечание дельное, исправил. Спасибо.
вы пишите про 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;
}

Корректное моделирование газов не является целью данной симуляции. Целью же является проверка правдоподобности модели подъёмной силы Ньютона. А именно, будут ли наблюдаться какие-либо интересные эффекты, если взять большое количество шариков и профиль крыла.
Корректное моделирование газов не является целью данной симуляции. Она не может быть использована для реальных расчётов. Для этих целей нужно использовать методы из CFD и соответствующие программные пакеты (Solidworks, Autodesk Airflow). Поэтому, модель не пройдёт вышеуказанные тесты.

Целью же является проверка правдоподобности модели подъёмной силы Ньютона. А именно, будут ли наблюдаться какие-либо интересные эффекты, если взять большое количество шариков и профиль крыла.
Именно! Сам впервые столкнулся на практике именно на видеокартах Turing. После того, как отстрелил себе ногу, начал изучать вопрос детально.

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

Два доклада Фёдора Пикуса многое прояснили в этой области. Поэтому я решил ими поделиться.

Information

Rating
Does not participate
Location
Москва, Москва и Московская обл., Россия
Date of birth
Registered
Activity