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

В этом ноутбуке мы продолжим исследовать простой метод для симуляции двумерных несжимаемых жидкостей для визуальных эффектов. Эта работа в основном основана на статье Джоса Стама «Stable Fluids» (SIGGRAPH 1999), а также на туториале Карла Симса.

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

То, что получится в итоге

Адвекция

С математической точки зрения, адвекция описывает общий процесс перемещения чего-либо по скалярному или векторному полю в направлении, заданном другим векторным полем. Зачем нам это? Увидим позже.

Для простоты начнём со скалярных полей.

grid = Table[{0, 0}, {5}, {5}];  grid[[3, 3]] = {0, 1}; u = Table[1., {5}, {5}];  (* наше скалярное поле, которое мы будем перемещать *)

Можно представить u как массу чего-то, что будет переноситься потоком grid.

Наивный подход

Смотрим на границы

Можно интерполировать скорость на границах и затем перемещать массу между ними. Мы можем записать это следующим образом

$\frac{u_{q}^{t} - u_{q}^{t-1}}{\delta t} = - \frac{v_{q} + v_{q+1}}{2} \frac{u_{q} + u_{q+1}}{2} + \frac{v_{q} + v_{q-1}}{2} \frac{u_{q}^{t-1} + u_{q-1}^{t-1}}{2}$

Предупреждение: это неверное выражение

Покажем, почему 😀. Реализация влоб в коде будет выглядеть так

advect[v_, u_, δt_:0.1] := With[{max = Length[v]}, With[{   take = Function[{array, x,y}, If[x >= 1 && x <= max && y >= 1 && y <= max, array[[x,y]], array[[1,1]] 0.]] },   Table[         With[{       v2 =  (*FB[*)((take[v, i, j+1] + take[v, i, j])(*,*)/(*,*)(2.0))(*]FB*).{0,-1},       v4 =  (*FB[*)((take[v, i, j-1] + take[v, i, j])(*,*)/(*,*)(2.0))(*]FB*).{0,1},       org = u[[i,j]]     },        org + (                v4  (*FB[*)((take[u, i, j-1] + org)(*,*)/(*,*)(2.0))(*]FB*) + v2  (*FB[*)((take[u, i, j+1] + org)(*,*)/(*,*)(2.0))(*]FB*)                ) δt     ]        , {i, max}, {j, max}] // Chop  ]]

Теперь, если мы протестируем это на нашем поле u

u = Table[1., {5}, {5}];

arrow = Graphics[Arrow[{{0,0}, {1,0}}], ImagePadding->None, ImageSize->50];  Table[MatrixPlot[Nest[advect[grid, #]&, u, n], ImageSize->200], {n, {0, 1, 100}}]; Riffle[%, arrow] // Row

Видите эти синие области — это отрицательная плотность. Что конечно нефизично. В нашем коде

v4  (*FB[*)((take[u, i, j-1] + org)(*,*)/(*,*)(2.0))(*]FB*) +  v2  (*FB[*)((take[u, i, j+1] + org)(*,*)/(*,*)(2.0))(*]FB*)

Если org (текущая ячейка) равна 0, это не означает, что она не может быть уменьшена ещё больше, если у нас есть вектор скорости к ближайшей ячейке.

Метод донорных ячеек (Donor-Cell advection)

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

$v_{q+1/2} = \frac{v_{q} + v_{q+1}}{2} \\ v_{q-1/2} = \frac{v_{q} + v_{q-1}}{2}$

Используя эти обозначения, можно записать метод донорской ячейки следующим образом

$\frac{u_{q}^{t} - u_{q}^{t-1}}{\delta t} = - v_{q+1/2}\begin{cases} u_{q}^{t-1} ,& v_{q+1/2} > 0 \\ u_{q+1}^{t-1} & \end{cases} + v_{q-1/2}\begin{cases} u_{q-1}^{t-1} ,& v_{q-1/2} > 0 \\ u_{q}^{t-1} & \end{cases}$

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

advect[v_, u_, δt_:0.1] := With[{max = Length[v]}, With[{   take = Function[{array, x,y}, If[x >= 1 && x <= max && y >= 1 && y <= max, array[[x,y]], array[[1,1]] 0.]] },   Table[         With[{       v1 =  (*FB[*)((take[v, i-1, j] + take[v, i, j])(*,*)/(*,*)(2.0))(*]FB*).{1,0},       v2 =  (*FB[*)((take[v, i, j+1] + take[v, i, j])(*,*)/(*,*)(2.0))(*]FB*).{0,-1},       v3 =  (*FB[*)((take[v, i+1, j] + take[v, i, j])(*,*)/(*,*)(2.0))(*]FB*).{-1,0},       v4 =  (*FB[*)((take[v, i, j-1] + take[v, i, j])(*,*)/(*,*)(2.0))(*]FB*).{0,1},       org = u[[i,j]]     },        org + (                v1 (*TB[*)Piecewise[{{(*|*)take[u, i-1, j](*|*),(*|*)v1 > 0(*|*)},{(*|*)org(*|*),(*|*)True(*|*)}}](*|*)(*1:eJxTTMoPSmNkYGAo5gESAZmpyanlmcWpTvkVmUxAAQBzVQdd*)(*]TB*) + v3 (*TB[*)Piecewise[{{(*|*)take[u, i+1, j](*|*),(*|*)v3 > 0(*|*)},{(*|*)org(*|*),(*|*)True(*|*)}}](*|*)(*1:eJxTTMoPSmNkYGAo5gESAZmpyanlmcWpTvkVmUxAAQBzVQdd*)(*]TB*)  +                  v4 (*TB[*)Piecewise[{{(*|*)take[u, i, j-1](*|*),(*|*)v4 > 0(*|*)},{(*|*)org(*|*),(*|*)True(*|*)}}](*|*)(*1:eJxTTMoPSmNkYGAo5gESAZmpyanlmcWpTvkVmUxAAQBzVQdd*)(*]TB*) + v2 (*TB[*)Piecewise[{{(*|*)take[u, i, j+1](*|*),(*|*)v2 > 0(*|*)},{(*|*)org(*|*),(*|*)True(*|*)}}](*|*)(*1:eJxTTMoPSmNkYGAo5gESAZmpyanlmcWpTvkVmUxAAQBzVQdd*)(*]TB*)                ) δt     ]        , {i, max}, {j, max}] // Chop  ]]

В блокноте этот код выглядит чуть более «математичным» за счёт синтаксических сладостей редактора

То же самое, но в редакторе WLJS Notebook

Теперь, если мы запустим тот же тест, результат должен быть корректным

u = Table[1., {5}, {5}]; arrow = Graphics[Arrow[{{0,0}, {1,0}}], ImagePadding->None, ImageSize->50];  Table[MatrixPlot[Nest[advect[grid, #]&, u, n], ImageSize->200], {n, {0, 1, 100}}]; Riffle[%, arrow] // Row

Мы можем убедиться в этом и «численно»

Nest[advect[grid, #]&, u, 100] // Chop // MatrixForm

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

Дифференциальное уравнение переноса

Я хотел бы обратить ваше внимание на дифференциальное уравнение адвекции

$\frac{\partial u}{\partial t} + (\mathbf{v} \cdot \nabla) u = 0$

Если попытаться дискретизировать его на одномерной сетке

$u^{t+1}_{q} - u^{t}_{q} + \frac{\delta t}{\delta q} (f^{t}_{q + 1/2} - f^{t}_{q - 1/2}) = 0$

Где $f_{q\pm\frac{1}{2}}$ — это поток на границе ячейки . Точная форма потока определяется используемым алгоритмом решения этого уравнения, а мы ведь уже знаем одно из них. Тогда для метода донорной ячейки получаем

$f_{q\pm 1/2}^{t} = v_{q\pm 1/2} \begin{cases} u_{q}^{t} ,& v_{q\pm 1/2} > 0 \\ u_{q \pm 1}^{t} & \end{cases}$

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

Импульс жидкости

Второй закон Ньютона в уравнении Навье-Стокса определяется следующим членом

$\frac{\partial \mathbf{v}}{\partial t} + (\mathbf{v} \cdot \nabla) \mathbf{v} + ... = 0$

В трёх точках скрываются явления вязкости, давление и все внешние силы. Пока ограничимся первыми двумя членами уравнения.

Ничего не напоминает? Это фактически адвекция на собственном поле скоростей или self-advection term.

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

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

vgrid = Table[{0.,0.}, {i,15}, {j,15}]; vcolors = Table[1.0, {Length[vgrid]}, {Length[vgrid]}]; vfps = 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}]; (* a very stipid way of drawing lines *)                    ];          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]

подписываемся на начало кадра и пересчитываем все поле

vtime = AbsoluteTime[];  EventHandler["vframe", Function[Null,   vgrid = advect[vgrid, vgrid] // N;    vcolors = Map[Function[val, (*FB[*)((π + 2.0 ToPolarCoordinates[val]// Last)(*,*)/(*,*)(3.0 π))(*]FB*) ], vgrid, {2}];      vfps = (*FB[*)(((vfps + 1 / (AbsoluteTime[] - vtime)))(*,*)/(*,*)(2.0))(*]FB*) // Round;   vtime = AbsoluteTime[];  ]];

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

With[{   win = CurrentWindow[],    currentCell = ResultCell[],   canvas = CreateUUID[] },    EventHandler[win, {"Closed" -> Function[Null,     Delete[currentCell]    ]}];    Graphics[{     Table[With[{i=i, j=j},               Offload[{          Hue[vcolors[[i]][[j]]],         Arrow[{{i,j}, {i,j} +  vgrid[[i]][[j]]}]       }]           ], {i,15}, {j,15}],          EventHandler[Graphics`Canvas[], mouseHandler[vgrid, win, MetaMarker[canvas]]],       AnimationFrameListener[vfps // Offload, "Event"->"vframe"],           MetaMarker[canvas],      Text[vfps // Offload, {0,0}]   },      Controls->False,      ImageSize->500,      PlotRange->{{-0.5,15.5}, {-0.5,15.5}},      ImagePadding->None,      TransitionType->None   ] ]

Как результат вы увидите следующее

Пробные частицы

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

Билинейная интерполяция

В Wolfram Language существуют готовые решения на этот счёт, однако для общности мы сделаем интерполяцию вручную. Это в любом случае будет полезно, так как потом захочется перенести код на GPU.

Возможно, это не самая быстрая версия, но достаточно хороша для учебных целей

bilinearInterpolation[array_, {x0_, y0_}] := Module[   {rows, cols, x = y0, y = x0, x1, x2, y1, y2, fQ11, fQ12, fQ21, fQ22},      {rows, cols} = Take[Dimensions[array], 2];      x = Clip[x, {1, cols}];   y = Clip[y, {1, rows}];      x1 = Floor[x];    x2 = Ceiling[x];   y1 = Floor[y];    y2 = Ceiling[y];      fQ11 = array[[y1, x1]];   fQ12 = array[[y2, x1]];   fQ21 = array[[y1, x2]];   fQ22 = array[[y2, x2]];       If[x2 == x1,     If[y2 == y1,       fQ11,       1/(2 (y2 - y1)) * (         fQ11 (y2 - y) +         fQ21 (y2 - y) +         fQ12 (y - y1) +         fQ22 (y - y1)       )     ],     If[y2 == y1,       1/(2 (x2 - x1)) * (         fQ11 (x2 - x) +         fQ21 (x - x1) +         fQ12 (x2 - x) +         fQ22 (x - x1)       ),       1/((x2 - x1) (y2 - y1)) * (         fQ11 (x2 - x) (y2 - y) +         fQ21 (x - x1) (y2 - y) +         fQ12 (x2 - x) (y - y1) +         fQ22 (x - x1) (y - y1)       )     ]   ] ]

Можно проверить, корректно ли оно работает «на глаз», растянув какие-нибудь мелкие изображения

bmp = Table[{(*SpB[*)Power[Cos[y](*|*),(*|*)2](*]SpB*), (*SpB[*)Power[Sin[x](*|*),(*|*)2](*]SpB*), 0.5} // N, {x, 0, Pi, Pi/5}, {y, 0, Pi, Pi/5}]; ibmp = Table[bilinearInterpolation[bmp, {x,y}],  {x, 1, 6 - 0.1, 0.0165 2}, {y, 1, 6 - 0.1, 0.0165 2}];  {   Image[ImageResize[Image[bmp, "Real32"], Scaled[25], Resampling->"Nearest"], Magnification->2],    Image[ibmp, "Real32", Magnification->2] } // Row

Слева — исходник, справа — интерполяция нашей функцией.

Равномерное распределение

Еще одна мелкая подзадача — как распределить наши пробные частицы равномерно по полю. Однако в стандартной библиотеке Wolfram Language есть готовое решение

RandomPointConfiguration[   HardcorePointProcess[10000, 0.4, 2],   Rectangle[{1, 1}, {5, 5}],   Method -> {"LengthOfRun" -> 10000000} ]["Points"] // Point // Graphics

Двигаем частицы

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

advectParticles[v_, pts_, δt_: 0.02] := Map[Function[p, p + δt*(bilinearInterpolation[v, p])], pts]

Симуляция несжимаемой жидкости в реальном времени

Наконец-то! Теперь мы можем связать решения наших уравнений из первой части и второй воедино

$div \mathbf{v} = 0, \qquad \frac{\partial \mathbf{v}}{\partial t} + (\mathbf{v} \cdot \nabla) \mathbf{v} = f_{ext}$

Откуда f_ext

Здесь, надеюсь, эксперты в гидродинамике могут меня поправить. Под $f_{ext}$ я подразумеваю не только мышь пользователя, но и эффекты давления, которые устраняют дивергенцию методом описанным в первой части. Про сам метод проекции: не совсем понятно, вовлекает ли он в себя элементы временного уравнения Навье-Стокса и если да, то как.

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

Предполагается, что вы уже объявили функцию removeDivergence из первой части, а также объявили advect, bilinearInterpolation и advectParticles а также mouseHandler

fgrid = Table[{0.,0.}, {i,15}, {j,15}]; fcolors = Table[1.0, {Length[fgrid]}, {Length[fgrid]}]; ffps = 0;  particles = RandomPointConfiguration[       HardcorePointProcess[10000       , 0.4, 2],       Rectangle[{1+4,1+4}, {15-4,15-4}], Method -> {"LengthOfRun" -> 10000000}]["Points"];

подписываемся на кадр со следующим алгоритмом

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

ftime = AbsoluteTime[];  fpipeline = Composition[removeDivergence, removeDivergence, advect[#,#, 0.2]&];  EventHandler["fframe", Function[Null,    fgrid = fpipeline[fgrid];      fcolors = Map[Function[val, (*FB[*)((π + 2.0 ToPolarCoordinates[val]// Last)(*,*)/(*,*)(3.0 π))(*]FB*) ], fgrid, {2}];      particles = With[{p = advectParticles[fgrid, particles, 0.3]},     advectParticles[fgrid, p, 0.3]   ];    ffps = (*FB[*)(((ffps + 1 / (AbsoluteTime[] - ftime)))(*,*)/(*,*)(2.0))(*]FB*) // Round;   ftime = AbsoluteTime[];  ]];

а теперь сам виджет

With[{   win = CurrentWindow[],    currentCell = ResultCell[],   canvas = CreateUUID[] },    EventHandler[win, {"Closed" -> Function[Null,     Delete[currentCell]    ]}];    Graphics[{Arrowheads[Medium/2],     Table[With[{i=i, j=j},               Offload[{          Hue[fcolors[[i]][[j]]],         Arrow[{{i,j}, {i,j} +  fgrid[[i]][[j]]}]       }]           ], {i,15}, {j,15}],           EventHandler[Graphics`Canvas[], mouseHandler[fgrid, win, MetaMarker[canvas]]],            AnimationFrameListener[ffps // Offload, "Event"->"fframe"],           MetaMarker[canvas],      PointSize[0.02], Point[particles//Offload],     Text[ffps // Offload, {0,0}]   },      Controls->False,      ImageSize->500,      PlotRange->{{-0.5,15.5}, {-0.5,15.5}},      ImagePadding->None,      TransitionDuration->35    ] ]

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

Скрытый текст

advect = Compile[{{v, _Real, 3}, {u, _Real, 3}, {δt, _Real, 0}} , With[{max = Length[v]}, With[{   },   Table[         With[{       v1 =  (*FB[*)((If[i-1 >= 1, v[[i-1, j]], {0.,0.}] + v[[i, j]])(*,*)/(*,*)(2.0))(*]FB*).{1,0},       v2 =  (*FB[*)((If[j+1 <= max, v[[i, j+1]], {0.,0.}] + v[[i, j]])(*,*)/(*,*)(2.0))(*]FB*).{0,-1},       v3 =  (*FB[*)((If[i+1 <= max, v[[i+1, j]], {0.,0.}] + v[[i, j]])(*,*)/(*,*)(2.0))(*]FB*).{-1,0},       v4 =  (*FB[*)((If[j-1 >= 1, v[[i, j-1]], {0.,0.}] + v[[i, j]])(*,*)/(*,*)(2.0))(*]FB*).{0,1},       org = u[[i,j]]     },        org +  (                v1 If[v1 >0, If[i-1 >= 1, u[[i-1, j]], {0.,0.} ], org]  + v3 If[v3>0,If[i+1 <= max, u[[i+1, j]], {0.,0.} ], org]+                  v4 If[v4 >0, If[j-1 >= 1, u[[i, j-1]], {0.,0.} ], org] + v2 If[v2>0, If[j+1 <= max, u[[i, j+1]], {0.,0.} ], org]                ) δt      ]        , {i, max}, {j, max}] // Chop  ]] ];

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}] ] ];

bilinearInterpolation = Compile[{{array, _Real, 3}, {v, _Real, 1}}, Module[   {rows, cols, x = v[[2]], y = v[[1]], x1, x2, y1, y2, fQ11, fQ12, fQ21, fQ22},      (* Get the dimensions of the array *)   {rows, cols} = {Length[array], Length[array]};      (* Clip points to the boundaries *)   x = Clip[x, {1, cols}];   y = Clip[y, {1, rows}];      (* Find the bounding indices *)   x1 = Floor[x];    x2 = Ceiling[x];   y1 = Floor[y];    y2 = Ceiling[y];      (* Get the values at the four corners *)   fQ11 = array[[y1, x1]];   fQ12 = array[[y2, x1]];   fQ21 = array[[y1, x2]];   fQ22 = array[[y2, x2]];      (* Perform bilinear interpolation *)   If[x2 == x1,     If[y2 == y1,       fQ11,       1/(2 (y2 - y1)) * (         fQ11 (y2 - y) +         fQ21 (y2 - y) +         fQ12 (y - y1) +         fQ22 (y - y1)       )     ],     If[y2 == y1,       1/(2 (x2 - x1)) * (         fQ11 (x2 - x) +         fQ21 (x - x1) +         fQ12 (x2 - x) +         fQ22 (x - x1)       ),       1/((x2 - x1) (y2 - y1)) * (         fQ11 (x2 - x) (y2 - y) +         fQ21 (x - x1) (y2 - y) +         fQ12 (x2 - x) (y - y1) +         fQ22 (x - x1) (y - y1)       )     ]   ] ] ];

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

ссылка на оригинал статьи https://habr.com/ru/articles/856442/