Если вам неинтересно долгое скучное введение, то переходите сразу к разделу о методике виртуальных труб. Но меня это немного расстроит.
Я одержим генерацией рельефа, играми на основе сеток, симуляциями и тому подобным. И часто во всём этом присутствует вода, или, по крайней мере, её присутствие кажется естественным.
Допустим, вы генерируете карту для стратегической игры, но не хотите, чтобы границы карты были заполнены непроходимой пустотой (как в олдскульных RTS). Разве не будет здорово, если граница будет заполнена водой, как на этой карте из одного моего заброшенного проекта?
Она создаёт красивую естественную границу, а может, и позволяет добавить связанные с водой механики, например хождение под парусом, рыбалка, торговля и морские бои. Вот 3D-вид похожего острова из того же проекта:
Или, допустим, вы делаете мирную игру-симулятор деревни/городка и хотите, чтобы в городке была река как источник питьевой воды, для рыбалки, транспорта и просто ради красоты. А может, река вам нужна просто как разделитель между областями, например, как в моей первой выпущенной игре:
Надеюсь, я вас убедил, что вода в играх очень полезна.
Но у неё есть проблемы.
Проблемы с водой
В большинстве игр отсутствует модификация рельефа, что вполне логично, ведь не везде она требуется. Но даже если она есть, то разработчики часто реализуют его достаточно просто: например, в игре наподобие Civilization можно вполне насыпать холм (хотя, он может многое поломать с точки зрения геймплея), который никак не будет взаимодействовать с водой.
Однако в моей игре модификация рельефа нужна:
Почему? Хм... хотя бы потому, что так написано в дизайн-документе!
А если серьёзно, то в моём случае на то были две причины:
Некоторые ресурсы в игре в буквальном смысле добываются из почвы, например, земля, песок и глина, и логично, что при их выкапывании слой почвы снимается (в противном случае мы получим бесконечный ресурс, что нас не устраивает)
Аналогично для ресурсов наподобие камней и руд металлов. Я предпочитаю, чтобы при добыче они выкапывались из почвы (а не на земле лежал булыжник «золотой руды», как это бывает во многих играх), и логично, что вырытая почва тоже пропадает
Моя игра не допускает строительство зданий на склонах, поэтому важно иметь возможность выравнивать рельеф перед постройкой, например, как в The Sims
Я хочу дать игроку множество инструментов для творческого самовыражения, и модификация рельефа будет одним из них
Да и будем откровенны, всем нравится модификация рельефа
Ну ладно, а как это всё связано с водой? Допустим, у нас есть озеро или даже лужа и вы хотите прорыть её границу, чтобы открылся проход для выливания воды. Как поведёт себя вода?
Существует множество вариантов решения задачи:
Вода никуда не денется и останется там, где была сгенерирована изначально (разумеется, это скучно и глупо)
Вода не течёт; всё ниже определённой высоты считается водой, как в игре Sapiens
Вода тоже существует ниже определённого уровня, но игрок не может копать настолько глубоко или он может копать только на один уровень вниз, чтобы срывать холмы/горы, как в RimWorld
Течение воды реализовано в виде крайне простой модели, как в Minecraft
Вода течёт на основании хорошей, но всё равно простой модели в стиле клеточных автоматов, как в Dwarf Fortress
Ни одно из этих решений меня не устраивало. Это хорошие альтернативы, если я не найду варианта получше, но в качестве основной модели воды они ощущались слишком скучными. Dwarf Fortress оказалась ближе остальных к тому, что мне нужно, но её модель слишком блочная, к тому же спроектирована для 3D, что мне не требуется (см. следующий раздел).
В течение нескольких лет (да, в буквальном смысле) я пассивно искал модель, отвечающую моим требованиям. Я прочитал кучу литературы о симуляциях жидкостей, в особенности о симуляциях климата и океанов/рек, и это меня очень испугало, потому что большинство моделей было безумно сложным, а их применимость в моей игре казалась неочевидной, потому что научное сообщество обычно не касается дизайна игр. Однажды мне даже удалось вычислить что-то вроде ожидаемого течения потоков вокруг острова, решив уравнение сохранения массы для потока жидкости вида ∇(mass⋅velocity)=0:
Уже кое-что, но мне нужна была динамическая модель, а не вычисление статических потоков. (Но я думаю, она всё равно была бы полезна для любопытной генерации климата на карте.) Кстати, теоретически можно легко превратить динамическую модель в статическую: достаточно добавить уравнения вида dXdt=0 для всех переменных состояния X, приняв, что у нас есть стабильное состояние (то есть идеально стабильное решение, не меняющееся со временем, например медленно текущая река).
Разумеется, можно было просто сдаться. В конце концов, важно помнить, что я не разрабатываю какой-то симулятор реальности, я делаю игру и могу изменять правила, как мне угодно. Однако я не мог избавиться от мысли, что в моём случае может сработать какая-нибудь хитрая комбинация простых формул. И я был прав!
Кстати, если в процессе чтения этого раздела семь раз завопили «Тимбербо-о-орн», то знайте: именно эту модель я и буду описывать.
Если вы знаете другие модели, которые, как вам кажется, подойдут для моего случая, то напишите, мне бы хотелось услышать о них!
Подготовка
Я уже несколько раз повторил «мой случай», но что же я имею в виду? Вот список возможностей, которая должна иметь моя симуляция воды:
Симуляция, вероятно, должна работать на сетке, предпочтительно на той же сетке, которую я использую для рельефа
Средний масштаб симуляции — примерно метр; меня не волнуют крошечные брызги, а километровая симуляция была бы слишком грубой для моего градостроительного симулятора (я с лёгкостью смогу имитировать высокочастотные визуальные детали без необходимости их симуляции)
Меня вполне устроит, если вода будет полем высот над рельефом, то есть она не будет течь вертикально и не будет иметь зазоров в вертикальных сечениях, потому что рельеф моей игры сам по себе представлен в виде поля высот; к тому же это сводит задачу к 2D
Вода должна уметь течь (очевидно)
Вода не должна волшебным образом исчезать из-за ошибок симуляции (такое происходило в некоторых моделях, которые я пробовал раньше)
Симуляция должна быть контролируемо стабильной
Симуляция должна быть достаточно быстрой; в идеале затраты на один такт симуляции должны линейно зависеть от размера симуляции (то есть несколько циклов for для области действия симуляции)
Как можно догадаться, я нашёл такую модель, которой и посвящена эта статья. Но для начала перечислю
Решения, которые мне не подошли
Позвольте описать пару популярных вариантов, не подошедших для моего сценария.
Безумно популярная методика реализации симуляций жидкостей, создающая впечатляющие результаты — это гидродинамика сглаженных частиц (Smoothed Particle Hydrodynamics). Однако эта методика решает совершенно иную задачу! Она создаёт реалистичные высокодетализированные и красивые симуляции, но мне это не требуется. Помните, что выше я говорил о «среднем масштабе»? Частицы воды диаметром 1 метр мне не подойдут, они будут выглядеть, как наполненные водой шары, но снижение размеров частиц слишком сильно влияет на производительность. Мне нужна не высокодетализированная симуляция, а быстрая и убедительная!
Наверно Stable Fluids Джоса Стэма — самая известная работа по симуляции жидкости в сфере компьютерной графики, но она тоже решает другую задачу. Она работает с полным объёмом жидкости, а не со свободной поверхностью, как мне нужно. Это отличие бассейна с водой от воды над рельефом. Кроме того, этот алгоритм далеко не быстрый: некоторые этапы симуляции требуют итеративного решения разреженной линейной системы, что реализуемо, но всё равно затратно.
На самом деле, я когда-то реализовал Stable Fluids, это интересная впечатляющая модель, только не для моей задачи (она решает полные уравнения Навье-Стокса, а не уравнения мелкой воды, которые мне нужны):
Уравнения мелкой воды
Предупреждение: я не физик, поэтому этот раздел может оказаться полной чушью.
Когда мы говорим о нахождении математической модели чего-то, то обычно подразумеваем уравнения, которые нужно решить. В общем случае движение жидкостей описывается уравнениями Навье-Стокса или более простыми уравнениями Эйлера. Однако эти уравнения тоже предназначены не для поверхности свободной воды, а для объёма, полностью заполненного жидкостью.
В этих уравнениях используются такие переменные, как вязкость жидкости, давление, плотность, термодинамическая работа, тензор напряжения и другие аспекты, и чтобы понять, что из этого известно и неизвестно, а что можно вывести из другого, нужно пройти полный курс динамики жидкостей (спустя годы ленивого изучения этой темы я так и не знаю ответа). Можно заметить, что здесь отсутствует понятие наподобие «объёма воды». У нас есть плотность массы, но все уравнения делятся на неё, поэтому не может быть нулевой плотности или границы между жидкостью и воздуха. (Если не ошибаюсь, их с лёгкостью можно реализовать при помощи методик particle-in-cell и marker-and-cell.)
Даже в своём простейшем виде эти уравнения используют давление p, которое известно и меняется со временем, но нет уравнения того, как же конкретно оно изменяется! Например, нет уравнения dpdt=… Вместо этого изменения поля давления от времени косвенно встроены в другие уравнения. Это тесно связано с этапом проецирования в солвере Stable Fluids.
Нам же нужно следующее: взять что-то типа уравнений Навье-Стокса, предположить, что у нас есть слой воды поверх какого-то рельефа (в данном случае обычно называемого дном), и выполнить некое усреднение вдоль вертикали, чтобы мы получили чисто двухмерные уравнения перемещений воды. Именно эту задачу решают уравнения мелкой воды!
Под «мелкой» подразумевается, что типичный вертикальный размер столба воды будет гораздо меньше, чем типичные размеры по горизонтали. Обычно это более-менее истинно, например, в случаях моделирования климата или паводков: высоты воды (допустим, метры или десятки метров в случае реки) гораздо меньше, чем горизонтальные расстояния (километры).
Такие уравнения обычно используются в науке в ситуациях «вода над рельефом», и они решают задачу, которую нам нужно решить. Я не буду приводить здесь сами уравнения, можете найти их в Википедии, потому что пост посвящён не этому. Я опишу саму методику решения и попытаюсь объяснить её принцип.
Неравномерные сетки
Перед тем, как мы приступим к описанию полного решения, нам нужно поговорить о сетке. Обычно при численном решении дифференциальных уравнений мы дискретизируем площадь симуляции в сетку (допустим, в сетку из квадратов) и храним значения в поле (вещи наподобие векторов скорости и давления) в каждой ячейке или вершине сетки.
Например, при решении волнового уравнения (которое лучше работает для звуковых и электромагнитных волн, чем для волн воды) мы храним высоту воды u(i,j) и её производную по времени du(i,j) для каждой ячейки. Код обновления значения выглядит так:
du(i,j) += (u(i+1,j)+u(i,j+1)+u(i-1,j)+u(i,j-1)-4*u(i,j)) * dt/dx/dx;
u(i,j) += du(i,j) * dt;Возможно, всё кажется непонятным, но основной принцип в том, что обе величины хранятся в одной сетке. Это приводит к красивым конечно-разностным уравнениям, которые проще выразить в коде и понять.
Такие сетки называются collocated-сетками, что я всегда читал, как co-located: обе величины находятся в одной сетке. А для динамики жидкостей такие сетки очень неприятны!
Одна из причин заключается в том, что они запутывают конечно-разностные методы. В волновом уравнении есть вторая переменная, которую можно вычислить красивым симметричным образом:
(f(x+1) + f(x-1) - 2*f(x)) / (dx * dx)Но в динамике жидкостей используются первые производные, в вещах наподобие «ускорение жидкости пропорционально разности между объёмами воды слева и справа». Как же вычислить такие конечные разности? (f(x+1) - f(x)) / dx приводит к смещению симуляции вправо, игнорирует некоторые направленные эффекты и в целом приводит к нестабильности. (f(x) - f(x-1)) / dx имеет те же проблемы. (f(x+1) - f(x-1)) / (2*dx) симметрична, но игнорирует значение в текущей ячейке f(x), что тоже приводит к нестабильности. В общем случае, наивная дискретизация уравнений динамики жидкостей обычно крайне нестабильна, и вместо решения уравнения всё просто колеблется дурацким образом.
Ещё один способ понять, что эта сетка неидеальна — рассмотреть следующую ситуацию: предположим, что мы храним в каждой ячейке высоту воды и вектор скорости воды. Теперь рассмотрим ячейку, в которую вода втекает слева и справа, а вытекает в верхнюю и нижнюю ячейки. Вот иллюстрация:
Общая скорость в этой ячейке равна нулю! Но это ерунда, ведь мы чётко видим, что вода активно течёт.
Подобные аспекты заставили исследователей пересмотреть методики на основе сеток и изобрести так называемые неравномерные сетки (staggered grid). Есть различные их вариации, но типичная, используемая для симуляции жидкостей, работает так: мы сохраняем, допустим высоту/плотность/так далее воды в квадратные ячейки, но храним скорость в границах между ячейками. В вертикальных границах (то есть в границах между соседями по горизонтали) хранится скорость по горизонтали и наоборот. Вот изображение:
Синими стрелками обозначены хранящиеся значения, по одному на каждую границу между ячейками. Обратите внимание, что крайние границы тоже присутствуют: они соответствуют пограничным условиям нашей симуляции (мы вернёмся к этому чуть позже).
Итак, если нам нужна в качестве области симуляции квадратная сетка N×N, то нам нужно будет хранить
Массив N×N, допустим, для высоты воды
Массив (N+1)×N для скорости по X
Массив N×(N+1) для скорости по Y
Это может показаться необычным, но чем раньше вы к этому привыкнете, тем более качественные симуляции жидкостей получите.
Методика виртуальных труб
Мы наконец-то добрались до той методики, которую я использую для симуляции воды над рельефом. Эта методика называется «виртуальные трубы» (virtual pipes), потому что она выведена из допущения о том, что ячейки воды соединены друг с другом воображаемыми трубами некоторого радиуса. В качестве справочных материалов для реализации методики я использовал две статьи. Можно заметить, что в первой статье рассматриваются многоуровневые столбы воды и вертикальные соединения, а вторая в основном посвящена гидравлической эрозии; ни то, ни другое я не реализовывал, потому что цель у меня была другая.
Во-первых, ради простоты забудем пока о рельефе, его тривиальным образом можно добавить позже. В неравномерной сетке мы будем, как и в предыдущем разделе, хранить три значения:
Массив N×N
waterдля высоты поверхности воды в этой ячейкеМассив (N+1)×N
flowXдля общего расхода воды между соседними по горизонтали ячейкамиМассив N×(N+1)
flowYдля общего расхода воды между соседними по вертикали ячейками
flowX(i,j) — это горизонтальный расход (также называемый истечением) между ячейками water(i-1,j) и water(i,j), если это не крайняя граница (i == 0 или i == N); аналогично для flowY.
Кстати, ниже я использую запись array(i,j) для обозначения элемента массива i,j. В моём реальном коде на C++ 2D-массивы имеют перегруженный operator() , обеспечивающий доступ к элементу для чтения и записи, так как многомерный operator[] есть только в C++23.
Обратите внимание, что мы храним расход, но не вектор скорости. Чтобы понять, что такое расход, представьте, что мы поместили между двумя соседними ячейками волшебный занавес, способный замерять проходящий через него объём воды. Поток проходящей через этот занавес воды имеет определённую площадь поперечного сечения и движется с некой скоростью. Если перемножить эти два значения, можно получить объём воды за единицу времени, и именно эту величину мы храним.
Расход обычно ведёт себя лучше, чем скорость при отсутствии воды. Для двух соседних пустых ячеек расход, очевидно, будет нулевым, потому что никакая вода не перемещается. А скорость — это нечто наподобие расхода, делённого на площадь поперечного сечения воды, то есть ноль, делённый на ноль, а это всегда вызывает проблемы. Можно сказать, что в этом случае скорость очевидно равна нулю, но что, если скорость и расход очень малы, но ненулевые? Мы попадаем в условия проблем и неоднородностей чисел с плавающей запятой, поэтому нам придётся задавать некие пороговые значения, ниже которых скорость считается нулевой. Всё это очень запутанно, поэтому работа с расходом вместо скоростей изящным образом решает все эти проблемы.
Но хватит болтовни; в методике используются три этапа:
Давайте разберёмся с ними всеми по порядку.
Ускорение расхода
Если у нас есть две соседние ячейки воды и их высоты различаются, то логично будет, что вода будет течь из более высокого в более низкий столб. Именно этот процесс и выполняет ускорение расхода: мы берём все внутренние (то есть не находящиеся на краях массива) границы и ускоряем расход в них на основании разности уровней воды в соответствующих ячейках воды. Мы делаем это для расходов по X:
for (int y = 0; y < N; ++y)
for (int x = 1; x < N; ++x)
flowX(x,y) += (water(x-1,y) - water(x,y)) * g * dt * A / dx;и аналогично по Y:
for (int y = 1; y < N; ++y)
for (int x = 0; x < N; ++x)
flowY(x,y) += (water(x,y-1) - water(x,y)) * g * dt * A / dy;Здесь dx и dy — это горизонтальный и вертикальный размер ячеек (длины соответствующих «труб»), dt — временной такт симуляции (к этому я вернусь ниже), g — гравитация, A — площадь поперечного сечения виртуальной трубы. В статье берётся A=dx*dx, а я взял A=1. В целом, A и g используются только на этом этапе и всегда в виде произведения, поэтому для наших простых целей можно просто игнорировать A и сделать вид, что она слилась с g.
Обычно на этом этапе также добавляется трение. Трение просто снижает расход в каждой итерации, обычно приводя к схождению симуляции к статическому стабильному состоянию. В статье рекомендуется использовать коэффициент pow(friction,dt), чтобы сделать его независимым от dt (что тесно связано со степенным сглаживанием, о котором я рассказываю в другой статье). Таким образом, код симуляции принимает следующий вид:
for (int y = 0; y < N; ++y)
for (int x = 1; x < N; ++x)
flowX(x,y) = flowX(x,y) * pow(friction,dt)
+ (water(x-1,y) - water(x,y)) * g * dt / dx;и аналогично для Y:
for (int y = 1; y < N; ++y)
for (int x = 0; x < N; ++x)
flowY(x,y) = flowY(x,y) * pow(friction,dt)
+ (water(x,y-1) - water(x,y)) * g * dt / dy;friction обычно имеет значение от 0 до 1. friction=0 означает максимальное трение, то есть расход из предыдущего такта симуляции полностью уничтожается. friction=1 означает полное отсутствие трения. Из-за этого я использую формулу pow(1-friction,dt), чтобы сделать значение трения более понятным. Разумеется, мы должны предварительно вычислить коэффициент трения pow(1-friction, dt) в начале такта симуляции, а не вычислять его для каждой ячейки.
С временным тактом dt всё чуть хитрее. Очевидно, что чем он больше, тем быстрее идёт симуляция. Однако большие значения dt также ведут к нестабильности. В общем случае существует известный критерий Куранта — Фридрихса — Леви, общий для всех симуляций жидкостей. Если упростить, он гласит, что коэффициент dx/dt должен быть не меньше максимальной скорости жидкости. На практике это значит, что мы должны уменьшать dt (таким образом снижая значение dx/dt), пока симуляция не станет стабильной. Я использовал значения примерно от 0.001 до 0.01.
То есть этот этап, по сути, ускоряет расход воды, наподобие dvdt=a=Fm из законов Ньютона.
Обновление столбов воды
Да, сначала мы рассмотрим третий этап симуляции, потому что он будет важен для понимания второго.
Наверно, этот этап самый простой. Для каждой ячейки воды мы просто прибавляем или убираем воду, согласно вертикальному и горизонтальному расходу, соседнему с этой ячейкой:
for (int y = 0; y < N; ++y)
for (int x = 0; x < N; ++x)
water(x,y) += (
flowX(x, y) + flowY(x,y )
- flowX(x+1,y) - flowY(x,y+1)
) * dt/dx/dy;flowX(x,y) и flowY(x,y) текут в сторону ячейки water(x,y), поэтому мы прибавляем их. flowX(x+1,y) и flowY(x,y+1) вытекают из ячейки, поэтому мы вычитаем их.
На этом этапе просто происходит перемещение воды между ячейками согласно вычисленному расходу. Наподобие dxdt=v из законов Ньютона.
Масштабирование истечения
Второй этап, наверно, сложнее всего. Дело в том, что на этапе 3 у нас может возникнуть проблема: если исходящий расход достаточно велик, то величина воды в ячейке может стать отрицательной! Это плохо, потому что отрицательной величины воды быть не может (в отличие от, например, электромагнитных волн).
К счастью, существует простое решение под названием «масштабирование истечения» (outflow scaling). Мы смотрим на расходы, соседние для ячейки воды, и учитываем только исходящий расход, то есть только расход, удаляющий, а не прибавляющий воду. Мы получаем общий исходящий расход, суммируя эти расходы, и сравниваем его с истинным объёмом воды в этой ячейке. Если мы обнаружим, что на третьем этапе мы попытаемся удалить больше воды, чем есть в ячейке, то просто уменьшаем объём исходящего расхода, чтобы объём воды оставался положительным (или нулевым). Вот код:
for (int y = 0; y < N; ++y) {
for (int x = 0; x < N; ++x) {
float total_outflow = 0.f;
total_outflow += max(0.f, -flowX(x,y));
total_outflow += max(0.f, -flowY(x,y));
total_outflow += max(0.f, flowX(x+1,y));
total_outflow += max(0.f, flowY(x,y+1));
float max_outflow = water(x, y) * dx*dy/dt;
if (total_outflow > 0.f) {
float scale = min(1.f, max_outflow / total_outflow);
if (flowX(x,y) < 0.f) flowX(x,y) *= scale;
if (flowY(x,y) < 0.f) flowX(x,y) *= scale;
if (flowX(x+1,y) > 0.f) flowX(x+1,y) *= scale;
if (flowX(x,y+1) > 0.f) flowX(x,y+1) *= scale;
}
}
}Немного запутанно из-за того, что нам нужно фильтровать только исходящий расход, зато выполняет нужную нам задачу.
Высота рельефа
Теперь нам нужно, чтобы вода двигалась по рельефу, но где в этих уравнениях рельеф? На самом деле, добавить рельеф очень просто. Представьте две соседних ячейки воды с равными высотами столбов воды: они не пытаются переливаться друг в друга, потому что уровень одинаков. Теперь представьте, что под левой ячейкой есть какая-то ненулевая высота рельефа. Она поднимает высоту воды, и теперь поверхность этой ячейки выше, чем у правой ячейки, поэтому вода начинает двигаться.
То есть при вычислении ускорения расхода нам просто нужно заменить высоту столба воды на высоту поверхности воды, то есть на высоту рельефа плюс высоту столба.
Если terrain(x,y) — это высота рельефа в ячейке, то нам достаточно лишь изменить вычисления ускорения:
for (int y = 0; y < N; ++y)
for (int x = 1; x < N; ++x)
flowX(x,y) = flowX(x,y) * pow(friction,dt)
+ (
water(x-1,y) + terrain(x-1,y)
- water(x,y) - terrain(x,y)
) * g * dt / dx;и аналогично для Y:
for (int y = 1; y < N; ++y)
for (int x = 0; x < N; ++x)
flowY(x,y) = flowY(x,y) * pow(friction,dt)
+ (
water(x,y-1) + terrain(x,y-1)
- water(x,y) - terrain(x,y)
) * g * dt / dy;Пограничные условия
При решении дифференциальных уравнений с частными производными (чем мы втайне здесь и занимаемся!) всегда важно учитывать происходящее на краях симуляции. В случае с методикой виртуальных труб пограничные условия косвенно определяются значениями потока крайних границ, то есть.
Левый край:
flowX(0,y)Правый край:
flowX(N,y)Нижний край:
flowY(x,0)Верхний край:
flowY(N,0)
Мы можем задать их так, как нам удобно:
Если присвоить значение 0, то они будут похожи на стены: волны воды будут просто сталкиваться с ними
Если задать в них входящий расход (положительный для левого/нижнего края, отрицательный для правого/верхнего), то они будут добавлять воду в симуляцию
Если задать в них исходящий расход (отрицательный для левого/нижнего края, положительный для правого/верхнего), то они будут убирать воду из симуляции
Для систем наподобие воды над рельефом кажутся логичными пограничные условия с исходящим расходом (вода на краях карты будет просто исчезать через край). Если некоторые части крайней границы пересекаются с рекой, то, вероятно, они должны создавать входящий расход.
Важно задавать эти значения в начале каждого такта симуляции, потому что масштабирование истечения может их изменять (и тогда крайняя граница с исходящим расходом превратится в границу-стену).
Вязкость
Статья также добавляет в симуляцию вязкость, просто масштабируя расход на некоторый коэффициент, зависящий от высоты воды (то есть отличающийся от трения). Смысл заключается в том, что меньшим слоям воды сложнее двигаться из-за различных внутренних сил, а большие слои воды могут перемещаться более свободно. Это можно просто превратить в умножение на
где H — это уровень воды в ячейке, откуда берётся расход (например, в левой ячейке для положительного расхода по X и в правой ячейке для отрицательного расхода по X), а ν — константа вязкости.
Код выглядит так:
for (int y = 0; y < N; ++y) {
for (int x = 1; x < N; ++x) {
float H = (flowX(x,y) > 0.f) ? water(x-1,y) : water(x,y);
H *= H;
if (H > 0.f)
flowX(x,y) *= H/(H + 3*dt*viscosity);
}
}И аналогично для Y:
for (int y = 1; y < N; ++y) {
for (int x = 0; x < N; ++x) {
float H = (flowY(x,y) > 0.f) ? water(x,y-1) : water(x,y);
H *= H;
if (H > 0.f)
flowY(x,y) *= H/(H + 3*dt*viscosity);
}
}Думаю, это могло бы быть полезно для чего-то наподобие потока магмы, но для воды я эти расчёты не использовал. Влияние вязкости важно в малых масштабах (в статье приводится пример потока крови внутри органов), но в крупномасштабных рельефах вязкость почти никак не влияет.
Полный код симуляции
Итак, вот полный код используемой мной симуляции:
// Инициируем расход на крайних границах
for (int i = 0; i < N; ++i) {
flowX(0,i) = ...; // левая граница
flowX(N,i) = ...; // правая граница
flowY(i,0) = ...; // нижняя граница
flowY(i,N) = ...; // верхняя граница
}
// Предварительно вычисляем коэффициент трения
float frictionFactor = pow(1-friction,dt);
// Ускоряем расход по X
for (int y = 0; y < N; ++y)
for (int x = 1; x < N; ++x)
flowX(x,y) = flowX(x,y) * frictionFactor
+ (
water(x-1,y) + terrain(x-1,y)
- water(x,y) - terrain(x,y)
) * g * dt / dx;
// Ускоряем расход по Y
for (int y = 1; y < N; ++y)
for (int x = 0; x < N; ++x)
flowY(x,y) = flowY(x,y) * frictionFactor
+ (
water(x,y-1) + terrain(x,y-1)
- water(x,y) - terrain(x,y)
) * g * dt / dy;
// Масштабируем истечение, чтобы предотвратить появление отрицательных объёмов воды
for (int y = 0; y < N; ++y) {
for (int x = 0; x < N; ++x) {
float total_outflow = 0.f;
total_outflow += max(0.f, -flowX(x,y));
total_outflow += max(0.f, -flowY(x,y));
total_outflow += max(0.f, flowX(x+1,y));
total_outflow += max(0.f, flowY(x,y+1));
float max_outflow = water(x, y) * dx*dy/dt;
if (total_outflow > 0.f) {
float scale = min(1.f, max_outflow / total_outflow);
if (flowX(x,y) < 0.f) flowX(x,y) *= scale;
if (flowY(x,y) < 0.f) flowX(x,y) *= scale;
if (flowX(x+1,y) > 0.f) flowX(x+1,y) *= scale;
if (flowX(x,y+1) > 0.f) flowX(x,y+1) *= scale;
}
}
}
// Обновляем столбы воды
for (int y = 0; y < N; ++y)
for (int x = 0; x < N; ++x)
water(x,y) += (
flowX(x, y) + flowY(x,y )
- flowX(x+1,y) - flowY(x,y+1)
) * dt/dx/dy;Вот и всё! Поначалу это может показаться пугающим, но после объяснения каждого этапа всё становится довольно логичным и интуитивно понятным. В конечном итоге, основная часть симуляции — это всего лишь четыре цикла for и несколько 2D-массивов с достаточно простыми формулами.
Мы можете посмотреть полный код обновления на C++, хотя там есть много чего другого.
Вот, как работает симуляция:
Видео моего симулятора воды на WebGPU, который я выпустил недавно.
(Частицы в видео используются только для визуализации, в самой симуляции они не участвуют).
После подбора хороших значений параметров (dt и g) симуляция работает чертовски устойчиво, при этом удовлетворяя всем моим требованиям; к тому же она более-менее походит на воду. Ура!
Недостатки модели
Разумеется, эта модель неидеальна. Одна из самых очевидных проблем заключается в отсутствии инерции и диффузии скоростей. Быстрый поток воды, впадающий в озеро, не будет распространяться дальше внутри озера, а потечёт в всех направлениях, игнорируя накопленную инерцию. Возможно существование двух потоков воды, которые движутся в противоположных направлениях и не взаимодействуют друг с другом (если уровни воды одинаковы).
Кроме того, симуляция обычно создаёт такие волны, когда вода впервые попадает в какую-то область. Выглядит немного странно, но вполне приемлемо:
Бонус: сетки шестиугольников/треугольников
В моей игре используется сетка из треугольников. Boris The Brave кратко описал все преимущества таких сеток, поэтому вам стоит прочитать его статью.
Сетки треугольников можно рассматривать как двойники сеток шестиугольников: достаточно соединить центры соседних шестиугольников отрезками, и они образуют сетку из правильных треугольников. Если судить по потрясающей статье Red Blob Games о сетках шестиугольников [перевод на Хабре], я использую осевую систему координат с шестиугольниками, вершины которых направлены вверх.
Здорово, что мы можем хранить такую сетку как обычные 2D-массивы (просто немного смещённые) и обращаться к вершинам такой сетки по координатам X,Y:
Для симуляции воды на такой сетке я буду хранить высоту столбов воды в вершинах сетки, что упростит рендеринг поверхности воды (это просто триангуляция, как и в случае с самим рельефом).
С потоками всё чуть сложнее. Существует расход в направлении X, в направлении Y и... в направлении Z? Давайте назовём это так. Вот иллюстрация:
Простите за мои рисунки, не было времени сделать хорошее векторное изображение
Итак, для каждой вершины у нас есть:
расход X, поступающий слева
расход X поступающий вправо
расход Y, поступающий слева снизу
расход Y, поступающий вправо вверх
расход Z, поступающий справа снизу
расход Z, поступающий влево вверх
Для сетки вершин N×N у нас есть:
Массив (N+1)×N расходов по X
Массив N×(N+1) расходов по Y
Массив (N+1)×(N+1) расходов по Z, в котором левое нижнее и правое верхнее значения не используются
Таким образом,
flowX(x,y)— это расход из (x-1,y) в (x,y)flowY(x,y)— это расход из (x,y-1) в (x,y)flowZ(x,y)— это расход из (x,y-1) в (x-1,y)
Не так много изменений по сравнению со случаем сетки квадратов, нам нужно лишь добавить расход по Z в 1) пограничные условия, 2) ускорение, 3) масштабирование истечения и 4) обновление воды. Самое сложное здесь — не запутаться с индексами.
Код такой симуляции на C++ выложен на bitbucket. Он работает достаточно неплохо и, надеюсь, симуляция чуть более однородна, чем на сетке квадратов.
Я ещё не реализовал эту систему в своей игре, есть множество других важных задач. Но я уже решил, что в моей игре точно будет хотя бы какая-то простейшая симуляция воды, и это меня очень радует. Представьте, что можно будет рыть траншеи, чтобы автоматически поливать большие фермы, или повернуть реку, чтобы она затопила соседнюю деревню. Возможности бесконечны.
