Очередная реализация fluid sim методом Эйлера, но в блокноте WL. Часть 1

Здесь мы рассмотрим простой метод симуляции несжимаемой жидкости в 2D для визуальных эффектов в интерактивном блокноте (впервые). Основная идея основана на работе Йоса Стама Stable Fluids, представленной на SIGGRAPH 1999, а также на учебном пособии Карла Симса.

Версия на английском и конечный блокнот

Часть 1. Дивергенция и игры с сеткой

Часть 2. Адвекция и первая симуляция

Часть 3. Чернила

Кстати, на Хабрахабре есть отличные статьи на эту тему

• https://habr.com/ru/articles/588372/ (Общий подход)

• https://habr.com/ru/articles/748430/ (Numba + не Эйлер)

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

Игры с сеткой

Создадим для примера такой двумерный массив

grid = Table[Cross[{i,j,0}, {0,0,1}][[;;2]], {i, 5}, {j, 5}]; % // MatrixForm 

Мы можем визуализировать её как векторное поле

grid // ListVectorPlot

При работе в любой численной сеткой для лучшей производительности можно использовать Map или MapIndexed с чистой функцией внутри. Это позволяет компилятору Wolfram Kernel выполнять JIT-компиляцию. К примеру такая операция

Map[Function[vector, ((*GB[*){{0(*|*),(*|*)-1}(*||*),(*||*){1(*|*),(*|*)0}}(*]GB*)) . vector], grid, {2}] ListVectorPlot[%]

Как результат такой операций мы эффективно перпендикуляры к нам векторам

Добавляем сахарку

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

.js core.ArrowGuide = async (args, env) => {   const dir = await interpretate(args[0], env);   const mag = dir.map((el)=>el*el).reduce((c,a) => c+a);    const angle = Math.round(180 * Math.atan2(dir[0], dir[1]) / Math.PI);      console.log(angle);   env.element.style.transform = `rotate(${angle}deg)`;      if (mag < 0.01) { //не показывать, если слишком маленький     env.element.style.opacity = 0.5;     env.element.innerHTML = ".";   } else {     env.element.innerHTML = "&uarr;";   } }

Теперь мы определим стандартную форму отображения для специального символа ArrowGuide[{x,y}], которая будет использовать описанную JavaScript-функцию для отображения стрелки внутри ячейки

ArrowGuide /: MakeBoxes[a_ArrowGuide, StandardForm] :=    ViewBox[a // First, a]

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

• a // First будет напечатано в редакторе (скрыто от глаз)

• a будет использоваться для отображения скрытого от глаз.

К примеру, теперь если написать в ячейке

ArrowGuide[{1,1}]

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

Теперь можно использовать эту обертку для отображения матрицы, применив ее в функции Map к нашему полю векторов, а затем вывести это поле как матрицу

Map[ArrowGuide, grid, {2}] // MatrixForm 

Также определим вспомогательную функцию для удобного отображения матриц

showAsMatrix[l_] := Map[ArrowGuide, l, {2}] // MatrixForm 

Дивергенция

Расходится что-то… Можно думать об этом, как о дифференциальном операторе, и как о характеристике векторного поля в заданной точке. Если вы забыли про эти штуки, закончив второй курс бакалавриата, то в принципе не страшно, так как оно обретет понятный смысл в наших жидкостях.

С сжимаемыми жидкостями работать тяжело, поэтому постулируем, что все несжимаемое

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

Пример

Создадим сетку с нулевыми векторами

grid = Table[{0,0}, {i,5}, {j,5}]; grid // showAsMatrix

Теперь добавим один вектор в центре

grid = ((*GB[*){{(*VB[*)({0, 0})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMYAIAiY0L1A=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0, 0})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMYAIAiY0L1A=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0, 0})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMYAIAiY0L1A=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0, 0})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMYAIAiY0L1A=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0, 0})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMYAIAiY0L1A=="*)(*]VB*)}(*||*),(*||*){(*VB[*)({0, 0})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMYAIAiY0L1A=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0, 0})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMYAIAiY0L1A=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0, 0})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMYAIAiY0L1A=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0, 0})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMYAIAiY0L1A=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0, 0})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMYAIAiY0L1A=="*)(*]VB*)}(*||*),(*||*){(*VB[*)({0, 0})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMYAIAiY0L1A=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0, 0})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMYAIAiY0L1A=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0, 1})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMmSBlAImRC9U="*)(*]VB*)(*|*),(*|*)(*VB[*)({0, 0})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMYAIAiY0L1A=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0, 0})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMYAIAiY0L1A=="*)(*]VB*)}(*||*),(*||*){(*VB[*)({0, 0})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMYAIAiY0L1A=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0, 0})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMYAIAiY0L1A=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0, 0})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMYAIAiY0L1A=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0, 0})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMYAIAiY0L1A=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0, 0})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMYAIAiY0L1A=="*)(*]VB*)}(*||*),(*||*){(*VB[*)({0, 0})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMYAIAiY0L1A=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0, 0})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMYAIAiY0L1A=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0, 0})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMYAIAiY0L1A=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0, 0})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMYAIAiY0L1A=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0, 0})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMYAIAiY0L1A=="*)(*]VB*)}}(*]GB*));

ну или более читаемое для редактора Хабра

grid[[3, 3]] = {0, 1}; grid // showAsMatrix

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

{intX, intY} = {     ListInterpolation[grid[[All,All,1]]],     ListInterpolation[grid[[All,All,2]]] }; (* по умолчанию - бикубическая *)  (* вспомогательная функция, чтобы раскрасить результат *) CoolColor[ z_ ] := RGBColor[z, 0.5, 1 - z];  Div[{intX[x,y], intY[x,y]}, {x,y}];  ContourPlot[%, {x,1,5}, {y,1,5}, ColorFunction->CoolColor]

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

Запрещаем жидкости сжиматься

Этот подход известен как проецирование

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

где — это новое значение вектора скорости в координатах , а — значения векторов из предыдущего состояния в соседних координатах .

Мы как бы добавляем поток от областей высокого давления (конвергентных) к областям низкого давления (дивергентным).

Вот реализация данного алгоритма на WL

removeDivergence[grid_] := With[{   (* вышли за границу сетки - получим нуль *)   take = Function[{array, x,y}, If[x >= 1 && x <= Length[grid] && y >= 1 && y <= Length[grid], array[[x,y]], {0,0}]] },   MapIndexed[Function[{val, i},      val + (1/8.0) (       ((take[grid, i[[1]] - 1, i[[2]] - 1] + take[grid, i[[1]] + 1, i[[2]] + 1]).{1,1}){1,1} +        ((take[grid, i[[1]] - 1, i[[2]] + 1] + take[grid, i[[1]] + 1, i[[2]] - 1]).{1,-1}){1,-1} +        (take[grid, i[[1]]-1, i[[2]]] + take[grid, i[[1]]+1, i[[2]]] - take[grid, i[[1]], i[[2]]-1] - take[grid, i[[1]], i[[2]]+1]){2,-2} + take[grid, i[[1]], i[[2]]] (-4)      )   ], grid, {2}] ]

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

А давайте проверим по шагам 0, 1, 2 итерации

Table[Nest[removeDivergence, grid, i] // showAsMatrix, {i, 0,2}] // Row 

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

Давайте проверим это напрямую вычислив дивергенцию явно

CoolColor[ z_ ] := RGBColor[z, 0.5, 1 - z];  plotDiv[grid_] := Module[{intX, intY, x, y},    {intX, intY} = {       ListInterpolation[grid[[All,All,1]]],       ListInterpolation[grid[[All,All,2]]]   };     With[{div = Div[{intX[x,y], intY[x,y]}, {x,y}]},     ContourPlot[       div, {x,1,5}, {y,1,5},        ColorFunctionScaling -> False,        ColorFunction -> ColorData[{"ThermometerColors", {-1, 1}}]     ]   ] ]

А теперь построим графики

{   {Style["Original", 14], plotDiv[grid]},    {Style["8 iterations", 14], plotDiv[Nest[removeDivergence, grid, 8]]} } // Transpose // Grid 

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

Интерактивный пример

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

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

grid = grid // removeDivergence;  Graphics[{   Table[     Arrow[{       {i,j},        {i,j} + 1.5 grid[[i]][[j]]     }], {i, 5}, {j, 5}] }, PlotRange->{{0,6}, {0,6}}, ImagePadding->None]

Что насчет цвета? Для цвета можно рассчитать Hue на основе полярного угла нашего вектора

Graphics[{   Table[{     Hue[((π + 2 ToPolarCoordinates[grid[[i]][[j]]]// Last)/(3 π))],     Arrow[{       {i,j},        {i,j} + 1.5 grid[[i]][[j]]     }]   }, {i, 5}, {j, 5}] }, PlotRange->{{0,6}, {0,6}}, ImagePadding->None]

Т.е. мы фактически воссоздали функцию VectorPlot с помощью примитивов. Следующим шагом будет реализация динамики и интерактивности.

Реактивность примитивов

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

Нам нужно связать примитив Arrow с конкретным элементом массива grid. Для этого можно использовать обертку Offload

grid = Table[{0.,0.}, {i,5}, {j,5}]; grid[[3,3]] = {0,1.0}; (* зададим источник *)  Graphics[{   Table[With[{i=i, j=j}, (* подставляем не символы а уже числа *)      Arrow[{       {i,j},        {i,j} + 1.5 grid[[i]][[j]]      } // Offload ]      ], {i, 5}, {j, 5}] }, PlotRange->{{0.5,5.5}, {0.5,5.5}}, ImagePadding->None]

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

grid = grid // removeDivergence;

Должно быть как-то так

Собираем все вместе

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

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

(* Все будущие динамические переменные *) bgrid = Table[{0.,0.}, {i,15}, {j,15}]; bcolors = Table[1.0, {Length[bgrid]}, {Length[bgrid]}]; bfps = 0;

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

mouseHandler[bgrid_, win_, canvas_] = Module[{   arrow = False,   dest = {0,0},   start = {0,0} }, {       "mouseup" -> Function[xy,         With[{v = -Normalize[start - xy]},           Do[              With[{p = Clip[#, {1,15}] &/@ Round /@ ((xy - start) a + start)},               bgrid[[p[[1]],p[[2]]]] = {v[[1]], v[[2]]};             ];           , {a, 0, 1,0.1}]; (* очень глупый способ рисования линии *)                    ];          Delete[arrow];          arrow = False;              ],        "mousemove" -> Function[xy,         dest = xy;       ],            "mousedown" -> Function[xy,         start = xy;         dest = xy + {1,1}; (* просто на всякий случай *)                If[arrow =!= False, Delete[arrow]];                  arrow = FrontSubmit[{           AbsoluteThickness[3], Gray,            Arrow[{xy, dest // Offload}]         },            canvas,            "Window"->win,            "Tracking"->True          ];              ]   }];  SetAttributes[mouseHandler, HoldFirst]

Здесь задается набор обработчиков (или списка правил) для событий от мыши. Они модифицируют поле скоростей по символу (а не по значению, из-за атрибута HoldFirst) линейно интерполируя направление скоростей от точки где была нажата клавиша мыши к точке, где была отпущена. Очевидно, что нужна визуализация этого пути стрелкой, а для этого служит функция FrontSubmit, которая фактически «исполняет» примитив Arrow, в окне win, в контексте canvas, ссылка на который передается в аргументе.

Давайте также назначим обработчик на пока неизвестное событие "bframe", который будет обновлять поле, а также считать кадры

btime = AbsoluteTime[];  EventHandler["bframe", Function[Null,   bgrid = removeDivergence[bgrid] // N;   bcolors = Map[Function[val, (*FB[*)((π + 2.0 ToPolarCoordinates[val]// Last)(*,*)/(*,*)(3.0 π))(*]FB*) ], bgrid, {2}];      bfps = (*FB[*)(((bfps + 1 / (AbsoluteTime[] - btime)))(*,*)/(*,*)(2.0))(*]FB*) // Round;   btime = AbsoluteTime[];  ]];

А теперь соберем все в интерактивный виджет

With[{   win = CurrentWindow[],    currentCell = ResultCell[],   canvasId = CreateUUID[] },    EventHandler[win, {"Closed" -> Function[Null,     Delete[currentCell] (* удалим ячейку если блокнот закрыли *)   ]}];    Graphics[{     Table[With[{i=i, j=j},               Offload[{          Hue[bcolors[[i]][[j]]],         Arrow[{{i,j}, {i,j} +  bgrid[[i]][[j]]}]       }]           ], {i,15}, {j,15}],           EventHandler[Graphics`Canvas[], mouseHandler[bgrid, win, MetaMarker[canvasId]]],           AnimationFrameListener[bfps // Offload, "Event"->"bframe"],           MetaMarker[canvasId], (* та самая ссылка внутрь Graphics *)     Text[bfps // Offload, {0,0}]   },      Controls->False,      ImageSize->500,      PlotRange->{{-0.5,15.5}, {-0.5,15.5}},      ImagePadding->None,      TransitionType->None   ] ]

Здесь есть вспомогательные конструкции, которые получают ссылку на текущее окно, а также помечают место внутри Graphics — canvasId и все такое.

Результат может выглядеть вот так

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

Здесь мы лишь удовлетворили свойство несжимаемости жидкости.

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

removeDivergence = Compile[{{grid, _Real, 3}}, With[{   max = grid // Length },   MapIndexed[Function[{val, i},      val + (*FB[*)((1)(*,*)/(*,*)(8.0))(*]FB*) (       (         (           If[i[[1]] - 1 >= 1 && i[[1]] - 1 <= max && i[[2]] - 1 >= 1 && i[[2]] - 1 <= max, grid[[i[[1]] - 1, i[[2]] - 1]], {0.,0.}]                       + If[i[[1]] + 1 >=1 && i[[1]] + 1 <= max && i[[2]] + 1 >= 1 && i[[2]] + 1 <= max, grid[[i[[1]] + 1, i[[2]] + 1]], {0.,0.}]                  ).{1,1}                ){1,1} +        (         (           If[i[[1]] - 1 >= 1 && i[[1]] - 1 <= max && i[[2]] + 1 >= 1 && i[[2]] + 1 <= max, grid[[i[[1]] - 1, i[[2]] + 1]], {0.,0.}]                      + If[i[[1]] + 1 >= 1 && i[[1]] + 1 <= max && i[[2]] - 1 >= 1 && i[[2]] - 1 <= max, grid[[i[[1]] + 1, i[[2]] - 1]], {0.,0.}]                    ).{1,-1}                ){1,-1} +        (         If[i[[1]]-1 >= 1 && i[[1]]-1 <= max, grid[[i[[1]]-1, i[[2]] ]], {0.,0.}]                  + If[i[[1]]+1 >= 1 && i[[1]]+1 <= max, grid[[ i[[1]]+1, i[[2]] ]], {0.,0.}]                  - If[i[[2]]-1 >= 1 && i[[2]]-1 <= max, grid[[ i[[1]], i[[2]]-1 ]], {0.,0.}]                  - If[i[[2]]+1 >= 1 && i[[2]]+1 <= max, grid[[i[[1]], i[[2]]+1]], {0.,0.}]                ){2,-2}                   + grid[[ i[[1]], i[[2]] ]] (-4)      )   ], grid, {2}] ],  "CompilationTarget" -> "C",   "RuntimeOptions" -> "Speed"];

Увидимся в следующей части 🧙