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

Код проекта можно найти на моём Github.


Ресурсы, посвящённые симуляции жидкостей, могут быть очень пугающими. Помню, что когда я впервые прочитал статью о ней и увидел уравнения Навье-Стокса, то был в ужасе. Со временем я понял, что сама тема совершенно не сложна. На самом деле, если бы вам дали задачу написать собственный симулятор жидкостей, то вы, вероятно, написали бы что-то подобное на основании интуитивного понимания работы жидкостей. Если, посмотрев на уравнения Навье-Стокса, вы подумали «ага, понятно», то вы, возможно, быстрее бы реализовали симуляцию жидкости, прочитав работы, перечисленные в конце этой статьи. Я же попытаюсь не использовать сжатое выражение и объяснять всё как можно медленнее.

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

Свойства жидкостей

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

Вы можете заметить три следующих аспекта:

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


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

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

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

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

Кроме того, я буду моделировать всё в 2D, потому что так проще для объяснения. Поняв базовую идею, будет несложно перенести её в 3D.

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


Рисунок 1 — Задаём водоём как последовательность полей в ячейке

На основании наших трёх наблюдений мы найдём величину интересующего нас поля для каждой ячейки по следующим трём факторам:

Количества чего-то в этой ячейке = что диффундировало из соседних ячеек + что было принесено откуда-то ещё из-за движения воды + что было добавлено в ячейку пользователем в этом кадре

Иными словами:

FluidStateForFieldX = Diffusion_of_X + Advection_of_X + X_Added_by_the_user

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

Наивная реализация, диффузия

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


Рисунок 2 — Сетка A; поле помеченное крестом — текущая обрабатываемая ячейка

Пусть дана сетка (A). Количество краски в этом кадре d при количестве краски в предыдущем кадре d0 вычисляется как:

d_X= d0_X + diffusionFactor * deltaTime (d0_01 + d0_02+ d0_03 + d0_04 -4*d0_X)

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

Эта система достаточно проста для реализации на GPU или CPU.

Наивная реализация, адвекция

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

Field_amount_for_target[cellPosition + current_cell_velocity * timestep] = field_amount_for_current_cell

Рисунок 3 — Адвекция

Реализация, пользовательский ввод

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

Amount_of_field_X += user_input_for_this_cell

Что наивного в этой реализации?

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

Изъян первый — сжимаемая жидкость?

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

Рассмотрим следующий сценарий. В показанной ниже сетке все соседние ячейки имеют векторы скорости, направленные в одну ячейку. Если все эти молекулы перенесутся в следующем кадре в эту ячейку, то как они там поместятся? Они сожмутся?


Рисунок 4 — Сжимающиеся молекулы воды

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


Рисунок 5 — Вещество, постоянно переносимое из точки

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

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


Рисунок 6 — Жидкость движется туда, куда может, потому что она несжимаема

Изъян второй: неудобство для реализации на GPU

Давайте снова взглянем на вычисление адвекции. Каждая ячейка смотрит на вектор скорости, который имеет в данный момент, а затем распространяет величину какого-то из полей на основании этой скорости в то место, которого достигнет вектор скорости после шага времени. Иными словами, каждая ячейка выполняет запись не в собственную память, а в некое неизвестное место. Такой тип паттерна доступа называется scatter operation.

Если задуматься о многопоточности (на CPU или GPU), то у нас возникнет проблема. Если мы привяжем каждый поток к ячейке, то потоки будут выполнять запись в участки памяти, в которые потенциально должны выполнять запись и другие потоки. Это создаёт условие гонки.

В идеале нам нужна схема, при которой поток привязывается к ячейке и каждый поток выполняет запись в память ячейки, которая ему назначена, а считывание выполняет только из соседних ячеек. Это будет примером gather operation. GPU отлично справляются с gather operations.

Изъян третий: нестабильность

Этот изъян сложно объяснить и полное его рассмотрение находится вне рамок нашей статьи.

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

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

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


Рисунок 7 — Слишком большая дельта времени приводит к взрыву значений. Каждый шаг имеет локальную погрешность, которая со временем складывается в очень большие значения

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

Проецирование

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

Для начала давайте выработаем интуитивное понимание того, как это можно сделать. Вещество не может сжиматься, не может иметь областей с повышенным и пониженным давлением. Когда мы вталкиваем воду в область так, что это повышает мгновенное давление (см. пример на Рисунке 4), вытолкнутые частицы пытаются переместиться из областей повышенного давления в области с пониженным давлением так, чтобы избежать сжимания. Частицы продолжат двигаться из областей с высоким в области с низким давлением, пока общая плотность частиц не станет одинаковой. Если мы сможем вычислить вектор скорости, вызванной этой разницей давлений (вызванной движением водоёма) и прибавить его к имеющемуся вектору скорости, то мы получим поле векторов скорости, не содержащее областей с увеличением давления. То есть мы получим поле векторов скорости, в котором не происходит сжатия. Можно посмотреть на это и иначе: мы компенсируем области, в которые нельзя втолкнуть воду и из которых нельзя её вытолкнуть, изменяя вектор скорости поля.

velocity_field + velocity_due_to_pressure_difference = new_velocity_with_no_pressure_difference

velocity_field = new_velocity_with_no_pressure_difference - velocity_due_to_pressure_difference

Ниже (см. Рисунок 11) мы рассмотрим связь между разностью давления и вызванным ею вектором скорости:

Difference_In_Pressure = -velocity_due_to_pressure_difference

Учитывая это выражение, окончательная формулировка будет такой:

velocity_field = new_velocity_with_no_pressure_difference + Difference_In_Pressure

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

velocity_field = Difference_In_Pressure + divergence_free_velocity
divergence_free_velocity = velocity_field -Difference_In_Pressure

Но как нам вычислить накопившееся в ячейке давление? Рассмотрим следующую схему.


Рисунок 8 — Различные сценарии и то, как они приводят к давлению или его отсутствию

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

Для этого вычисления нам нужно рассмотреть соседние ячейки во всех измерениях, то есть не только по горизонтали, но и по вертикали. На Рисунке 8 посмотрите на схему в правом нижнем углу. Хотя горизонтально давление нарастает (оба вектора направлены в A), вертикально давление спадает (векторы направлены от A). То есть в сумме здесь не происходит нарастание давления.

Если мы вычтем компонент X векторов скорости ячеек справа и слева от интересующей нас ячейки, то получим число, сообщающее нам, движутся ли эти два вектора скорости в одном направлении по оси X, то есть вызывают ли соседние ячейки в этом измерении нарастание давления. Можно сделать то же самое для верхней и нижней ячеек (но с компонентом вектора скорости Y) и сложив два скалярных значения со знаком (вклад в нарастание давления вдоль горизонтальной и вертикальной осей), мы получим одно скалярное число, представляющее количество воды, сходящееся к центру ячейки, или расходящееся от него. Это количество и является дивергенцией поля векторов скоростей текущей ячейки.


Рисунок 9 — Графическая схема и объяснение оператора дивергенции

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

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

Обозначим давление как . Как и вектор скорости, также является полем (каждая ячейка имеет собственное значение для давления), но в отличие от вектора скорости, давление — это только скаляр. На данный момент мы не знаем значение для каждой ячейки, однако давайте посмотрим, сможем ли мы создать уравнение, где значение давления для каждой ячейки будет неизвестным параметром. Посмотрите на Рисунок 9 и задайтесь вопросом: какой должна быть связь давления в центральной ячейке и её соседях, чтобы поле не имело дивергенции?


Рисунок 10 — Поле давления вокруг ячейки 11

Давление каждого поля влияет на 4 окружающих соседа, получая при этом четверть вклада от каждого из окружающих соседей. Поэтому чтобы в поле [11] возникло равновесие, давление в поле [11] должно быть равно сумме всех соседних давлений с соотношением 1 к 4. Можно посмотреть на это и иначе: давление в центральном поле должно быть равно среднему окружающих его давлений.

1/4 * (p[10] + p[01] + p[12] + p[21]) = p[11]
p[10] + p[01] + p[12] + p[21] = 4 * p[11]

Что происходит, когда эти два параметра неравны? Давайте вычислим разность между левой и правой частями:

(p[10] + p[01] + p[12] + p[21])- 4 * p[11]

В случае, когда давление в центре выше, чем вокруг, вещество будет выталкиваться из ячейки [11] и создастся отрицательное значение. Если окружающее давление выше, то вещество будет вталкиваться в [11] и значение будет положительным. Вы уже могли заметить, что это и есть дивергенция, и это значение, которое вычисляет наш оператор дивергенции, поэтому можно записать:

(p[10] + p[01] + p[12] + p[21])- 4 * p[11] = divergence(velocity_field)

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

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

difference_in_p => pressure_difference_induced_velocity

Чтобы преобразовать эту корреляцию в какое-то уравнение, вероятно, потребуется измерение массы или плотности (вспомним f=ma) или какая-то константа.

Как измерить изменение давления поля вдоль основной оси? Это достаточно простая операция, для ячейки [11] (Рисунок 10) мы вычитаем значение ячейки справа, ячейки слева, и результат становится компонентом X разности давлений, а затем можно то же самое сделать в вертикальном направлении, и это значение становится компонентом Y. Посмотрите для примера на Рисунок 11.

Приведённое ниже вычисление тоже является оператором, который можно записать как функцию. Чтобы объяснить функцию, я на секунду сделаю вид, что знаю значения поля давлений. Этот оператор называется оператором градиента, так как мы вычисляем градиент давления вдоль каждой из осей. Обратите внимание, что мы берём скалярное поле и получаем векторное поле. Возможно, вы уже заметили кое-что, связанное с направлением результирующего вектора и вектора скорости, который должна вызывать разность давлений. Если слева присутствует высокое давление, а справа — низкое, то оператор градиента создаёт вектор, направленный справа налево (как показано на Рисунке 11). Однако, как сказано ранее, частицы будут течь из областей высокого давления к низкому давлению, но никак не наоборот, то есть окончательная формулировка выглядит так:

difference_in_p = -pressure_difference_induced_velocity

Рисунок 11 — Оператор градиента, вычисляющий градиент/разность давления в каждом измерении

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

velocity_field = Difference_In_Pressure + divergence_free_velocity

Теперь мы можем записать его так:

velocity_field = Gradient(Pressure) + divergence_free_velocity

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

Divergence(velocity_field) = Divergence(Gradient(Pressure)) + Divergence(divergence_free_velocity)

Дивергенция (divergence_free_velocity) по определению равна нулю, потому что она обнуляется, и теперь у нас осталось одно уравнение и одна неизвестная! Итак, у нас осталось следующее:

Divergence(Gradient(Pressure)) = Divergence(velocity_field)

(p[10] + p[01] + p[12] + p[21])- 4 * p[11] = Divergence(velocity_field)

И это именно то же, что и выведенная ранее формулировка, если выполнить вычисления. Опять-таки, здесь сделано много допущений, например, разделение Divergence(A+B) на Divergence(A) + Divergence(B). Любопытный факт: дивергенция градиента — это тоже оператор, имеющий собственное название. Он называется оператором Лапласа и он является второй произвольной поля. Мы встретим его в других областях, например, в обработке геометрии, в контурных эффектах, компьютерном зрении, выявлении признаков и т.д.

Итак, теперь у нас есть уравнение для вычисления, где все значения (давления) неизвестны. Давайте выясним, как это сделать!

(p[10] + p[01] + p[12] + p[21])- 4 * p[11] = Divergence(velocity_field)

Солверы

На первый взгляд кажется, что уравнение дивергенции давления невозможно решить, так как у нас есть 5 неизвестных и одно уравнение. На самом деле, у насесть уравнение для каждой ячейки в нашем поле. То есть если поле имеет размер N на N, то у нас есть N*N неизвестных и такое же количество уравнений. По сути, это система линейных уравнений, которую вы можете помнить из старшей школы!

Напомним, что система уравнений может выглядеть так:

Уравнение 1) y = x + 1
Уравнение 2) y = 5x — 2

Каждое из этих уравнений имеет собственный график и геометрически точка или точки (если они есть) пересечения этих графиков являются решением уравнения. Другими словами, парами значений x и y, при которых оба уравнения верны.

Из старшей школы вы можете помнить, что один из способов решения этой простой системы является подстановка: подставляем уравнение 1 во второе, вычисляем один параметр, а затем подставляем значение параметра в любое из уравнений, чтобы получить второй параметр. Примерно так:

Уравнение 2) y = (x + 1) = 5x — 2 <= так как y также является x +1
=> x — 5x = -2 -1 = -3
=> x = -3/-4 = 3/4

Подставляем в Уравнение 1) y = 3/4 + 1 = 7/4

X = 3/4 = 0,75, y = 7/4 = 1,75

Однако наша система совершенно непохожа на это уравнение. При сетке 128 на 128 у нас получится 16 тысяч неизвестных и уравнений. Нам совсем не хочется программировать что-то подобное на CPU или GPU. Кроме того, в нашей системе есть нечто уникальное: значение давления для каждой ячейки зависит только от 4 соседей и никак не связано с остальными 16 тысячами неизвестных. Так что если мы запишем одно из этих 16 тысяч уравнений, то увидим, что 99% неизвестных в каждом из них равны нулю. Поэтому мы и не хотим записывать их как матрицу и обращать её, так как большинство элементов будет равно нулю. Есть более быстрый способ!

Прежде чем я объясню этот быстрый способ, давайте совершим с нашей простой системой фокус. Преобразуем два уравнения так, чтобы в левой части у одного был X, а у другого Y.

Уравнение 1) y = x + 1

Уравнение 2) x = y/5 + 2/5

Теперь возьмём любое начальное значение для x и y (например, ноль) и вставим его в правую часть. Получим:

Уравнение 1) y = 0 + 1 = 1

Уравнение 2) x = 0/5 + 2/5 = 2/5 = 0,4

Теперь у нас есть новые значения X и Y, 0,4 and 1. Давайте повторим этот шаг ещё раз, но вместо нулей подставим новое значение.

Уравнение 1) y = 0.4 + 1 = 1,4

Уравнение 2) x = 1/5 + 2/5 = 3/5 = 0,6

Мы можем повторять эти шаги множество раз, на каждом этапе беря предыдущее значение X и Y и снова вставляя его в уравнение. Давайте посмотрим, что произойдёт!

Шаг 0) Input (X, y) = (0, 0), output = (0.4, 1)

Шаг 1) In(0.4, 1), out(0.6, 1.4)

Шаг 2) In(0.6, 1.4), out((1.4 + 2)/5, 0.6 + 1 ) = (0.68, 1.6)

Шаг 3) In(0.68, 1.6), out((1.6 + 2)/5, 0.68 + 1 ) = (0.72, 1.68)

Шаг 4) In(0.72, 1.68), out((1.68+ 2)/5, 0.72 + 1 ) = (0.736, 1.72)

Магия! Даже через четыре шага мы уже почти получили решение. Можно повторить это ещё много раз и рано или поздно приблизиться к точному решению. Этот способ называется итеративным солвером, точнее, методом (солвером) Якоби. Удобно то, что это решение также применимо и к нашей системе, и его очень легко реализовать на GPU. Нам достаточно просто сохранять значение давления в каждой ячейке (как текстуру или вычислительный буфер) на каждом этапе нашей итерации и подставлять его на следующем этапе.

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

Однако остаётся один вопрос: почему это работает? Наиболее точным способом объяснения солверов является использование их матричной формы (диагональной матрицы Якоби для метода Якоби, нижний треугольник для метода Гаусса-Зейделя, специальным образом созданный треугольник для метода релаксации и мультисетка для многосеточного метода), при помощи которой мы также можем запросто объяснить стабильность, точность и т.п. Но я бы хотел избежать этого, чтобы не отпугнуть читателей, желающих обрести базовое понимание того, что они кодят.

Для своего объяснения я бы хотел использовать метод Гаусса-Зейделя вместо солвера Якоби, который мы используем в своей реализации, потому что он достигает результата вдвое быстрее. Разница между двумя этими методами проста, в методе Гаусса-Зейделя мы мгновенно используем значения, вычисленные для наших неизвестных в вычислениях следующих неизвестных, не ожидая целый этап. То есть в приведённом выше примере вместо подстановки 0 и 0 вместо X и Y мы подставим 0 вместо X, вычислим, что Y равен 1, и вместо подстановки 0 вместо Y в уравнение 2 мы сразу же используем значение 1, вычисленное на основе X, что даст нам 0,6. Затем мы переходим на этап 2 и подставляем 0,6 вместо X, чтобы найти новый Y (1,6), и так далее. Заметьте, что мы, по сути, доходим до того же решения, но за половину шагов. Фантастика, правда? Но существует проблема с GPU и многопоточностью, так как мы решаем все ячейки одновременно и не можем использовать эту технику без чего-то вроде красно-чёрного паттерна доступа, который уменьшая вдвое количество этапов, удваивает количество вызовов отрисовки. Поэтому выигрыш минимален.

Вернёмся к нашей простой системе. Ниже показан график двух уравнений. Помните, что решение находится в точке их пересечения. Мы уже знаем эти значения из аналитического решения, это X= 0,75 и Y=1,75

Уравнение 1) y = x + 1

Уравнение 2) y = 5x — 2

Что же на самом деле выполняет наш солвер? В качестве начальных значений давайте возьмём 3 вместо x (ранее мы брали вместо начального значения ноль, но выбрать можно любое). На шаге 1 мы находим значение Y, которое имеет наш график при X = 3. Затем мы смотрим на графике 2, какое значение X соответствует только что найденному Y. Повторяя этот процесс, мы можем заметить, что выполняем проецирование графиков друг на друга, в ту и другую сторону, и поскольку эти графики сходятся к одной точке, мы продолжаем приближаться к точке пересечения как мяч, катящийся вниз по склону.


Рисунок 12 — Графическое представление солверов.

Тогда на N-ном шаге мы получим нечто подобное. Красивая игра в пинг-понг всегда будет приближаться к истинному решению, но никогда не достигнет его и не перескочит его.


Рисунок 13 — Солверы на N-ном шаге, всегда приближаются к ответу, но никогда его не достигают

Разумеется, это работает не для всех уравнений. Два графика могут не сойтись или иметь несколько решений, и т.д. Но в нашем случае это работает.

Чтобы использовать солвер, мы перепишем уравнение таким образом, чтобы поле посередине находилось слева. Тогда для любой ячейки с координатами i и j мы получим:

Solve(Velocity_Divergence)
{
p_new[i, j] =((p[i+1,j] + p[i,j+1] + p[i-1,j] + p[i,j-1]) - Velocity_Divergence[i, j])/4
}

После запуска этой программы мы можем сделать следующее:

CopyOldToNew()
{
p[i, j] = p_new[i, j]
}

И повторять вычисления заново с новыми значениями p. Мы повторим это примерно 30 шагов, пока не получим решение для текущего кадра. Для первого шага мы берём p равным нулю.

То есть в общем наше проецирование будет выглядеть так

Velocity_Div_field = Divergence(velocity_field); 
for(int i = 0; i < 30; i++) 
{
   Solve(Velocity_Div_field);
   CopyOldToNew();
}
deltaP = Gradient(p);
new_divergent_free_velocity = velocity_field - deltaP;

Разумеется, это псевдокод, не стоит писать конвейер вычислительного шейдера таким образом. Вместо этого стоит создать буферы для всего этого, задать источники/получатели и вызывать ядра.

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

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

Реализация диффузии

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

Напомню, что при любом d поля, где d0 было значением поля для каждой ячейки в предыдущем кадре, явная формулировка выглядела так:

d_X= d0_X + diffusionFactor * deltaTime *(d0_01 + d0_02+ d0_03 + d0_04 -4*d0_X)

В случае с диффузией легко увидеть источник нестабильности. При большой дельте времени или показателях диффузии если количество плотности, поступающее от соседних ячеек, меньше того, что покидает центральную ячейку, и diffusionFactor*deltaTime больше 1, тогда значение плотности для нового кадра станет отрицательным, что не имеет смысла.

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

d0_X = d_X - diffusionFactor * deltaTime * (d_01 + d_02+ d_03 + d_04 -4*d_X)

Обратите внимание, что мы выполняем диффузию не d0, а d — значения следующего кадра. Чтобы понять это уравнение интуитивно, нужно разобраться, что же делает диффузия. Диффузия «смотрит», как количество вещества в центральной ячейке соотносится со средним в окружающих её ячейках. Если в центре вещества больше, то она забирает его часть, а если в центре вещества меньше среднего окружающих ячеек, то она добавляет часть. Это приводит к постепенному уравновешиванию значений и достижению однородного состояния. Так как мы предполагаем, что с течением времени всё уравновешивает друг друга, если посмотреть на процесс в обратную сторону (двигаться назад во времени), должно быть истинно противоположное. Происходит антидиффузия, при которой области с плотностью выше средней притягивают больше частиц. Противоположностью функции диффузии является вторая часть приведённого выше уравнения.

Но почему оно более стабильно или является лучшей формулировкой? Если переформулировать это уравнение (что нам всё равно предстоит сделать) так, чтобы новое значение d для центральной ячейки (X) единственное находилось в левой части, то мы получим следующее:

d_X = (d0_X + diffusionFactor * deltaTime * (d_01 + d_02+ d_03 + d_04)) / (1 + 4 * diffusionFactor * deltaTime)

Обратите внимание, что в отличие от явной формулировки, здесь есть только сложение. Поэтому нет никакой вероятности того, что новая плотность будет отрицательной, вне зависимости от дельты времени или diffusionFactor. Следовательно, d_x всегда будет оставаться положительным. Здесь также важно ещё раз повторить, что стабильность и точность — не одно и то же. При большой дельте времени или коэффициенте диффузии мы не получим значимых результатов. Но по крайней мере, они не будут взрываться!

Возможно, вы уже заметили, что вышеприведённая формулировка очень похожа на ту, которую мы получили для проецирования. Вместо d_X мы вычисляли давление в центральной ячейке. Значения с d_01 по d_04 были значениями давления в соседних ячейках, а d0_x было divergence(velocity_field).

Как и в случае с проецированием, на первый взгляд кажется, что это уравнение решить невозможно (одно уравнение и пять неизвестных, d_x и с d_01 по d_04), но это ещё и система линейных уравнений. По одному уравнению для каждой ячейки, и всего N*N. Хорошо то, что мы можем решить её уже имеющимся солвером.

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


Изменив _rDiagonal и _centerFactor, и превратив _b_buffer в divergnece(velocity) в одном и в d_x0 в другом, я могу использовать одинаковый солвер для обоих уравнений

Реализация адвекции

Далее нам нужно исправить нашу адвекцию, чтобы в ней использовалась gather operation, а не scatter operation.

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

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

field_amount_new_frame[cellPosition] += field_amount_old_frame[cellPosion-current_cell_velocity * timestep ]

Рисунок 14 — Обратите внимание, как это соотносится с рисунком для scatter operation. Величины всегда прибавляются к ячейке, которая вычисляется, и не переносятся в какое-то другое место.

Последние примечания

И на этом базовая симуляция жидкости завершена. Имея эту информацию, код (хорошо задокументированный) и документацию на Github, вы сможете получить симуляцию, которую действительно понимаете.

Для простой симуляции код будет выглядеть так (я много раз выполняю проецирование, потому что так выглядит лучше):


Готовый код для простой 2D-симуляции

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

Границы и произвольные границы

Заданные нами операторы (дивергенция, градиент и т.д.) будут иметь проблемы на краях границ, потому что соседние пиксели находятся вне сетки симуляции. Чтобы решить эту проблему, мы выполним другую программу для четырёх углов сетки.


Рисунок 15- GPU Gems, глава 38, внутренние и граничные ячейки

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


Обработка границ в четырёх углах.

С произвольными границами работать немного сложнее.


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


Рисунок 16 — Схема границ с произвольными нормалями из An Improved Study of Real-Time Fluid Simulation on GPU, Enhua Wu et al.

Я запекаю эту карту при помощи вычислительных шейдеров и массива, содержащего все возможные пермутации показанного выше и их нормалей. Индексами этого массива является своего рода битовая маска (не буду вдаваться в подробности, потому что для этого потребуется отдельная статья). Битовая маска позволяет очень быстро обновлять буфер нормалей при изменении динамических границ без ветвления в шейдере. Вот код для создания массива поиска нормалей на CPU и его использования для генерации буфера нормалей границ в вычислительном шейдере, а также код, выполняемый по разу для каждой ячейки, обрабатывающий различные поля в симуляции жидкости, находящиеся в пределах или на краях произвольных границ.





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

Фальшивая 3D-симуляция

В моей демо-сцене с персидским садом видно, что симуляция воды выглядит трёхмерной. Однако это не 3D-симуляция, внутри она всё равно работает как 2D. Я использую буфер давлений (с тональной коррекцией для подходящего диапазона) для определения того, куда нужно выталкивать вершины поверхности воды — вверх или вниз. С физической точки зрения это имеет смысл, так как если бы это была 3D-симуляция, часть воды двигалась бы вверх при высокой разности давлений, и в меньшей степени в стороны. Однако это всё-таки хак, хотя и красиво выглядящий, и гораздо менее затратный, чем выполнение настоящей 3D-симуляции воды.

Освещение

Я реализовал для всех объектов в демо-сцене costume lighting. Это вариация стандартного Cook-Torrance BRDF, о которой подробнее можно прочитать здесь. Сама вода состоит из следующих частей:

1. Зеркального освещения от направленного источника света
2. Преломления (для которого я просто беру текстуру под водой и отображаю её с искажениями/смещениями на основании движения поверхности воды)
3. Отражения

Например, здесь вы видите зеркальную часть освещения от направленного источника света:


Также я добавил немного рассеянного освещения на случай, когда вода будет грязной. Однако этот эффект незавершён, потому что для него требуется SSS и объёмные тени/столбы света. Проект на будущее!

Отражение

Отражение в демо состоит из двух частей. Первая — это планарное отражение:


Однако одного планарного отражения недостаточно для случаев с предельным смещением. Как показано на схеме ниже, оно работает хорошо в случае 1, но для пограничных случаев мне нужно скомбинировать его с отражением в экранном пространстве, поскольку решение имеет их оба.




Рисунок 18 — Одного планарного отражения недостаточно для предельных смещений

Каустика

Мой подход к каустике достаточно стандартен. Существует два типа каустики: 1) подводная, 2) отражённая от воды.

Оба работают с проецированием. Я проецирую текстуру давлений из позиции источника света на землю и беру первую производную (фильтр Собеля) буфера давлений с тональной коррекцией в качестве источника моей каустики. С точки зрения физики это противоположно тому, где должна появляться каустика, но поскольку это выглядело лучше, я выбрал этот вариант. Чтобы улучшить внешний вид преломления, зависящего от длины волны, каустика, применённая к поверхности дна, имеет другие смещения считывания RGB-каналов буфера давлений, что образует плавную границу каустики с небольшим радужным эффектом.


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

Анимация рыб

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

Микродетали

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

Дальнейшая разработка?

Если вы прочитали статью и вам она понравилась, то что ещё можно почитать об этой теме? Хороши представленные ниже ресурсы, кроме того, у Стэма (автора динамики жидкости для игр) есть книга, в которой он раскрывает тему подробнее. Также вы сами можете добавить многое к моей реализации. Самый большой выигрыш, мне кажется, можно получить от перехода к многосеточному солверу. Кроме того, можно оптимизировать систему, использовав вместо структурированных буферов текстуры (чтобы использовать текстурные кэши), или применив для солвера групповую общую память. К тому же мой движок жидкости написан для объяснения системы. Иногда я использую текстуры большего необходимого размера или буферы со слишком большой точностью. Некоторые буферы дублируются в памяти ради простоты использования, да и весь конвейер в целом можно оптимизировать. Ещё одной областью, требующей дополнительных исследований и настроек является выбор более подходящего алгоритма тональной коррекции для преобразования буфера давлений в карту смещений, а также комбинирование традиционных «фальшивых» техник течения воды с симуляцией для добавления микродеталей.

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










Красивые снимки демо-сцены