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

Исходный код к этой статье выложен на GitHub.

Огонь как явление очень интересен с точки зрения компьютерной графики. Раньше огонь было принято имитировать. Например, в фильме «Властелин колец» Питера Джексона для изображения огня использовались спрайты с огромным количеством дыма (на тот момент симуляция флюидов обходилась слишком дорого, даже при бюджете блокбастера). Когда требовалось моделировать огонь в режиме реального времени, например в видеоиграх, использовались почти исключительно нефизические подходы.

Но в последние 10 лет, благодаря развитию GPU, быстрая симуляция флюидов значительно упростилась. Базовые алгоритмы, описывающие динамику флюидов, не составляет труда реализовать на GPU. В 2009 году компания Industrial Light & Magic воспользовалась этими методами для моделирования и рендеринга огня в фильме о Гарри Поттере. В 2014 году компания NVIDIA выпустила FlameWorks — целую систему для генерации эффектов «огня» и «дыма» в компьютерных играх.

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

Симуляция флюидов

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

1.1 Базовая динамика флюидов

Допустим, ? — это область в пространстве, заполненная флюидом. В любой точке x∈D и в момент времени t флюид имеет скорость u(x,t). С вычислительной точки зрения можно представить двумерное поле скорости u с сеткой N×N, где равномерно распределённые по сетке точки позволяют узнать, каково значение поля скорости в конкретной точке пространства.

Пример: Сетка 16×16 представляет u=(y,−x)

Что произойдёт, если добавить во флюид каплю краски?

Определим скалярное поле ψ(x,t) как плотность красителя в любой точке пространства и в любой момент времени. Перенос величин, например, ψ, помноженный на скорость флюида, характеризуется величиной, которая называется адвекция. Имея некоторое поле скорости флюида и исходное поле плотности для вашего красителя, то было бы интересно рассмотреть, как плотность красителя развивается во всём интересующем нас объёме. Для этого можно было бы построить симуляцию адвекции красителя (его распространения) во флюиде.

Упрощённое представление адвекции

Вычислить адвекцию можно было бы, например, так: взять каждую из точек в сетке и продвинуть её в таком направлении и на такое расстояние, которое преодолела бы частица со скоростью, действующей в данной точке сетки и при такте симуляции Δt, после чего обновить тот узел сетки, который расположен ближе всего к точке «приземления» частицы:

ψ(x+u(x,t),t+Δt)=ψ(x,t)

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

Дифференциальные уравнения в частных производных, описывающие адвекцию

До сих пор в рамках всей этой статьи мы решали дифференциальное уравнение в частных производных! Если нам потребуется вывести стабильный метод для описания адвекции, сначала нам нужно получить явное выражение для данного ДУЧП. Начнём с основ.

Рассмотрим фиксированную область пространства W (т.е., W со временем не меняется). Общая масса красителя в пределах W равна  ∫W​ψdV. С течением времени масса будет меняться так:

Теперь обозначим через S поверхность W, а n будет вектор внешней нормали, определённый вдоль поверхности. В таком случае мы сможем проверить, каков будет массовый расход через поверхность W. В частности, наблюдается, что объёмный расход – то есть объём жидкости, протекающий в секунду через S на единичную площадь равен u⋅n, а массовый расход на единичную площадь равен ψu⋅n.

Так мы получаем закон сохранения массы в интегральной форме:

Можем ли мы избавиться от интегралов и сформулировать что-либо подобное для точек? По теореме о дивергенции, вышеприведённое выражение эквивалентно следующему:  

Далее, для единичной подобласти W=dV, справедливо, что

В результате получаем явное уравнение в частных производных, которое мы должны решить для ψ!

И... на этом можно было бы остановиться, но, может быть, нам удастся ещё сильнее упростить это выражение. В частности, представляется, что можно было бы вынести член ∇⋅u, который обнуляется из-за несжимаемости.

Применив наше ограничение несжимаемости ∇⋅u=0 в конце, получаем скалярное уравнение в частных производных, а далее — и первое из наших уравнений для адвекции несжимаемого потока:

Уравнение 2 для адвекции векторного поля v через поле скорости можно вывести примерно так же, как и скалярное уравнение для адвекции.

Стабильный метод адвекции

Давайте подробнее рассмотрим уравнение (1):

Обратите внимание: в правой части уравнения расположена производная по направлению −u. В результате получаем чудесный новый метод, позволяющий узнать адвекцию ψ для несжимаемой жидкости: начиная с точки x в сетке, проследим скорость флюида вспять, заменив значение из нашей исходной точки на значение той точки, в которой мы «сядем» (если мы «сядем» между точками сетки, то можем применить интерполяцию):

На псевдокоде для GPU:

global Vec2Field u;
          global FloatField density;
          global float dt;

          // В нашей скалярной сетке выполняем для всех наших точек, которые мы хотим обновить
          float advectPoint(vec2 x) {
            vec2 coord = x - dt * getVec2At(u, x);
            return getFloatAt(density, coord);
          }

Такой метод называется полулагранжевой моделью адвекции, в 1999 году его изобрёл Джос Стэм. Эта модель подобна эйлеровской, которая является её аналогом первого порядка точности, но обладает именно теми дополнительными свойствами, которые нам требуются:

1.      Её исключительно легко распараллелить, так как любая точка в сетке обновляется всего один раз за итерацию.

2.      Она безусловно устойчива. Почему? Обратите внимание: в любой точке сетки максимальное значение, до которого может дойти обновление, одновременно является максимальным для всех точек сетки.

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

1.2 Уравнения Навье-Стокса

Итак, нам удалось найти модель, описывающую, как с течением времени развиваются скалярные свойства флюида, при условии, что поток фиксирован. Что же можно сказать о потоке флюида как такового — как поле скорости u влияет само на себя с течением времени?

Уравнения Навье-Стокса для несжимаемого потока жидкости определяют, как с течением времени развивается скорость в любой точке движущегося флюида:

Здесь константа μ означает вязкость жидкости, а F — внешние силы. Но ранее мы также предположили, что наш флюид является невязким, так что μ=0, и пока мы можем просто не учитывать внешние силы. Итак, у нас остаётся всего два члена: самоадвекция и давление.

Если на каждом такте симуляции вычислять эти члены и складывать их, так мы и сможем смоделировать наш флюид! Вот что получится в псевдокоде:

          
let u = createVectorGrid();
          let density = createScalarGrid();
          let p = createScalarGrid();

          while (true) {
            // Решить для следующего поля скорости
            u = advect(u, u);

            p = computePressure(...);
            u = u - gradient(p);

            // Адвекция красителя для нового поля скорости
            density = advect(u, density);
          }          

Давайте подробнее рассмотрим каждый из этих компонентов.

Самоадвекция

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

Откуда берутся другие члены этого уравнения? В результате умножения адвекции u на себя получаем новое поле скорости u′ (вычислимое по вышеприведённому полулагранжевому алгоритму поиска с возвратом):

Давление

Неизвестно, удовлетворяет ли это новое поле условию несжимаемости (например, характеризуется ли нулевым расхождением). Этот момент нужно как-то скорректировать при помощи члена p, означающего давление:

Переформулируем и получим:

Уравнения такого типа называются уравнениями Пуассона, в левой части которого находится оператор Лапласа (лапласиан) неизвестного скалярного поля, а в правой части находится известный скаляр. Решение этого уравнения Пуассона — это, на самом деле, наиболее медленный шаг в процессе вычислений, связанных с динамикой флюидов. Почему — обсудим ниже.

Решаем задачи с давлением

Как же нам решить именно это уравнение в частных производных для p? Что ж, мы знаем, значение интересующего нас поля скорости u′ во всех точках нашей сетки, поэтому можем приблизительно вычислить правую часть уравнения Пуассона, применив повсюду дискретную версию функции расхождения:

где u′=(u,v) в 2 измерениях.

Затем можно воспользоваться дискретной версией лапласиана:

чтобы преобразовать всё это уравнение в линейное с пятью неизвестными.

Но на самом деле мы решаем вышеприведённое уравнение Пуассона для всего пространства сразу. Поэтому, работая с сеткой N×N, получаем систему из N² линейных уравнений, в которых ровно N² неизвестных! Таким образом, задача сводится к старому доброму уравнению

Ax=b

где A — это матрица, применяющая оператор Лапласа ко всей сетке, а b — это вектор, в котором содержится дивергенция поля скорости во всех точках сетки.

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

Эффективное решение задачи с давлением

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

В самом деле, невозможно добиться эффективного и одновременно точного решения системы линейных уравнений. Но при этом нужно отметить, что система линейных уравнений и так является лишь приближением уравнения Пуассона. Причём можно добиться произвольно точного приближения, если подходить к нему итеративными методами. Поэтому достаточно выбрать итеративный алгоритм. Итеративный алгоритм начинается с оценки решения и постепенного его уточнения на каждой итерации. Таким образом, можно выбрать такой алгоритм и прогонять его до тех пор, пока у нас не получится «достаточно хороший» результат. К примеру, есть очень простой и при этом лёгкий в реализации итеративный алгоритм для решения линейных уравнений — он называется метод Якоби.

Начнём с самого первого уравнения в системе:

​На k-й итерации, имея определённую догадку x_k о том, каково может быть решение x, у нас накопится определённая ошибка. Этой ошибкой можно воспользоваться, чтобы уточнить нашу догадку для x₁​ следующим образом:

При работе по методу Якоби имеем «догадки» по поводу всех элементов x, причём вычисляем их параллельно. Такая задача отлично подходит для реализации на GPU. Вот как эта реализация выглядит на псевдокоде:

global FloatField divergence;
          global FloatField pressure;
          global float texelSize;

          // Выполнить для каждой точки в нашей сетке значений давления, которые мы хотим обновить 
          float iterateJacobi(vec2 x) {
            float div = getFloatAt(divergence, x);
            float L = getFloatAt(pressure, x + vec2(-texelSize, 0.));
            float R = getFloatAt(pressure, x + vec2(texelSize, 0.));
            float T = getFloatAt(pressure, x + vec2(0., texelSize));
            float B = getFloatAt(pressure, x + vec2(0., -texelSize));
            return (div - L - R - T - B) / -4.;
          }

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

Резюме: симуляция Навье-Стокса

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

let u = createVectorGrid();
        let density = createScalarGrid();
        let div = createScalarGrid();
        let p = createScalarGrid();

        while (true) {
          // Решить для следующего поля скорости
          u = advect(u, u);

          // Навязать несжимаемость, проецируя давление.
          div = divergence(u);
          for (let i = 0; i < JACOBI_ITERATIONS; i++) {
            p = updatePressure(p, div);
          }
          u = u - gradient(p);

          // Показать адвекцию красителя в контексте нового поля скорости.
          density = advect(u, density);
        }

1.3 Ограничение завихрённости

Нам будет очень удобно хранить в сетке поле скорости, но в результате всякий раз, когда приходится интерполировать значения между точками сетки, возникает нежелательное численное сглаживание. Вкупе с относительно грубым приближением, характерным для полулагранжевой схемы адвекции первого порядка, это приводит к рассеиванию турбулентных вихрей в нашем потоке. Физика явления такова: поле скорости теряет энергию, в результате получается крайне сглаженный, «невыразительный» поток флюида.

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

Мельчайшие турбулентные элементы, различимые в нашей модели — это вихри, центром каждого из которых является одна из точек на пересечениях линий сетки. Можно измерить интенсивность этих вихрей (их завихрённость), взяв в каждой точке ротор от u и амплифицировав их. Математически завихрённость определяется так:

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

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

Здесь постоянная удержания ϵ>0 — это параметр, от которого зависит, в каком количестве мелкие детали возвращаются обратно в поток. Даже при сравнительно низком уровне удержания (около 0-15) разница бывает огромной, особенно в тех симуляциях, где используются полулагранжевые схемы адвекции. При сравнительно высоких уровнях удержания могут получаться сильно стилизованные потоки, накатывающиеся волнами.

 

Что будет, если добавить каплю краски в вихревые потоки

На GPU можно рассчитать ротор и удержание следующим образом:

  global Vec2Field u;
          global float texelSize;

          // Выполнить, чтобы получить ротор для каждой точки в сетке
          float computeCurl(vec2 x) {
            float L = getVec2At(u, x + vec2(-texelSize, 0.)).y;
            float R = getVec2At(u, x + vec2(texelSize, 0.)).y;
            float T = getVec2At(u, x + vec2(0., texelSize)).x;
            float B = getVec2At(u, x + vec2(0., -texelSize)).x;
            return (R - L) - (T - B);
          }

          global Vec2Field curl;
          global float confinement;

          // Выполнить, чтобы определить силу удержания в каждой точке сетки
          vec2 confinementForce(vec2 x) {
            float L = getFloatAt(curl, x + vec2(-texelSize, 0.));
            float R = getFloatAt(curl, x + vec2(texelSize, 0.));
            float T = getFloatAt(curl, x + vec2(0., texelSize));
            float B = getFloatAt(curl, x + vec2(0., -texelSize));
            float C = getFloatAt(curl, x);

            vec2 N = vec2(abs(T) - abs(B), abs(R) - abs(L));
            N = N / length(N);
            return confinement * C;
          }

Полная симуляция с учётом турбулентности:

          let u = createVectorGrid();
          let density = createScalarGrid();
          let div = createScalarGrid();
          let p = createScalarGrid();
          let curl = createVectorGrid();

          while (true) {
            // Решить для следующего поля скорости
            u = advect(u, u);

            // Использовать ограничение завихрённости для амплификации турбулентности в поле скорости.
            curl = computeCurl(u);
            u = u + confinementForce(curl, CONFINEMENT);

            // Навязать несжимаемость с проекцией давления.
            div = divergence(u);
            for (let i = 0; i < JACOBI_ITERATIONS; i++) {
              p = updatePressure(p, div);
            }
            u = u - gradient(p);

            // Впрыскивание красителя в новом поле скорости.
            density = advect(u, density);
          }

Вихревой шум и турбулентность

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

А затем вычислить конечное поле завихрённости по формуле

Быстродвижущиеся сильно турбулентные флюиды, в частности, дым и огонь, особенно выигрывают от ограничения завихрённости в сочетании с вихревым шумом, а на практике поле вихревого шума ϕ одновременно развивается с течением времени, а также испытывает адвекцию со стороны потока флюида.

Симуляция огня

Если вы дочитали до сих пор — можете себя похвалить! При помощи методов, описанных в предыдущем разделе, можно эффективно и точно симулировать жидкости с переменными физическими параметрами (масло, воду, мёд), исходя при этом из того, что жидкость занимает фиксированное пространство в ёмкости. Если вас интересуют другие (смежные) области, например, хочется узнать, как взаимодействуют флюиды, занимающие разные участки сетки, то можно сделать такую симуляцию сетки, в которой будут учитываться различные граничные условия.

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

2.1 Базовая модель сгорания

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

Сформулируем точнее: определим скалярное поле ?, где 0≤ρ≤1 — это плотность топлива, а другое скалярное поле T>0 представляет температуру флюида, распределённого по всей системе. На каждом такте симуляции из-за сгорания топлива повышается температура в системе, а в самом топливе горение продолжается лишь при условии, что поддерживается температура, необходимая для горения этого топлива:

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

Кроме того, разогретые молекулы сбрасывают избыточную температуру, излучая свет, по закону Стефана-Больцмана этот процесс описывается уравнением пятой степени:

Здесь параметр σ_cool​ характеризует скорость остывания. Чтобы симуляция получилась физически корректной, мы установим в качестве этого параметра постоянную Стефана-Больцмана, но в случае с графической симуляцией целесообразнее, чтобы художник мог сам контролировать скорость остывания.

Завершая модель сгорания, отметим, что наше топливо горит постоянно (это представимо, так как частицы ионизированного газа отдают тепловую энергию и переходят в более низкое энергетическое состояние), поэтому на каждом такте симуляции горючее расходуется с некоторой заданной скоростью сгорания γ_fuel​:

2.2 Термическая плавучесть

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

Здесь β — это заданная положительная постоянная плавучести, а j — направленный вверх единичный вектор.

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

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

Код симуляции выстраивается на основе базовых процедур, описывающих флюид:

let u = createVectorGrid();
          let density = createScalarGrid();
          let div = createScalarGrid();
          let p = createScalarGrid();
          let curl = createVectorGrid();
          let fuel = createScalarGrid();
          let temp = createScalarGrid();

          while (true) {
            // Решить для следующего поля скорости
            u = advect(u, u);

            // Этап сгорания
            temp = combust(temp, fuel);

            // Использовать ограничение завихрённости для амплификации турбулентности в поле скорости.
            curl = computeCurl(u);
            u = u + confineVorticity(curl, CONFINEMENT);

            // Добавить термическую плавучесть
            u = u + buoyancy(temp);

            // Навязать несжимаемость с проекцией давления.
            div = divergence(u);
            for (let i = 0; i < JACOBI_ITERATIONS; i++) {
              p = updatePressure(p, div);
            }
            u = u - gradient(p);

            // Адвекция плотности, топливо и температура в контексте нового поля скорости.
            density = advect(u, density);
            temp = advect(u, temp);
            fuel = advect(u, fuel);
          }          

2.3 Отображение огня

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

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

Закон Планка описывает спектральную плотность света, излучаемого абсолютно чёрным телом при заданной температуре T:

где

​а h, c и k — это, соответственно, постоянная Планка, скорость света и постоянная Больцмана.

Реализовав рендеринг излучения абсолютно чёрного тела при помощи фрагментных шейдеров, мы полностью смоделировали огонь!

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