В этой статье мы применим некоторые оптимизации к коду, увеличим разрешение и перейдем к software render нашей жидкости.

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

Измерение и оптимизация ⏱️

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

Версии из предыдущей статьи

ClearAll[advect]; ClearAll[removeDivergence]; ClearAll[bilinearInterpolation];

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

 removeDivergence[grid_] := With[{
  (*BB[*)(*safety checks, which enforce closed boundaries*)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*)
  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 + (*FB[*)((1)(*,*)/(*,*)(8.0))(*]FB*) (
      ((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}]
]

bilinearInterpolation[array_, {x0_, y0_}] := Module[
  {rows, cols, x = y0, y = x0, x1, x2, y1, y2, fQ11, fQ12, fQ21, fQ22},
  
  (* Get the dimensions of the array *)
  {rows, cols} = Take[Dimensions[array], 2];
  
  (* 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)
      )
    ]
  ]
]

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

Теперь запустим бенчмарк

runTest[title_] := With[{},
  testGrid = Table[{0.,0.}, {i,15}, {j,15}];
  testParticles = Table[RandomReal[{1,15}, 2], {i,1000}];

  timing = {0., 0., 0.};

  timing[[1]] = -AbsoluteTime[];
  testGrid = advect[testGrid, testGrid, 0.2];
  timing[[1]] += AbsoluteTime[];

  timing[[2]] = -AbsoluteTime[];
  testGrid = removeDivergence[testGrid];
  timing[[2]] += AbsoluteTime[];

  timing[[3]] = -AbsoluteTime[];
  testParticles = advectParticles[testGrid, testParticles, 0.2];
  timing[[3]] += AbsoluteTime[];

  {
    Style[title, Italic, 12],

    Flatten /@ {
      {Style["time, s", Italic], Round[timing, 0.001]},
      {Style["Max FPS", Italic], Round[1 / (timing // Total), 1]}
    } // TableForm 
  } // Column
];

runTest["Uncompiled"]

Компилируем в байт-код WL

Нужно всего-лишь упростить конструкции, удалив объявление функции внутри Module и заменив Piecewise на обычное условие If. Это делает наш код менее читаемым и убирает всю магию WL, но результат того стоит. Больший пинок производительности идет также из-за строгой типизации переменных в Compile

Код

advect = Compile[{{v, _Real, 3}, {u, _Real, 3}, {δt, _Real, 0}} , With[{max = Length[v]}, With[{
 
},
  Table[ 
  
    With[{
      (*BB[*)(* here we add a lot of manual check for boundary condition *)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*)
      v1 = (*BB[*)( (*FB[*)((If[i-1 >= 1, v[[i-1, j]], {0.,0.}] + v[[i, j]])(*,*)/(*,*)(2.0))(*]FB*))(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeB5AILqnMSXXKr0hjgskHleakFnMBGU6JydnpRfmleSlpzDDlQe5Ozvk5+UVFDGDwwR6dwcAAAAHdFiw="*)(*]BB*).{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 +  (
        (*BB[*)(* here we add a lot of manual check for boundary condition *)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*)
      
        v1 (*BB[*)(If[v1 >0, If[i-1 >= 1, u[[i-1, j]], {0.,0.} ], org])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeB5AILqnMSXXKr0hjgskHleakFnMBGU6JydnpRfmleSlpzDDlQe5Ozvk5+UVFDGDwwR6dwcAAAAHdFiw="*)(*]BB*)  
        
        + v3 (*BB[*)(If[v3>0,If[i+1 <= max, u[[i+1, j]], {0.,0.} ], org])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeB5AILqnMSXXKr0hjgskHleakFnMBGU6JydnpRfmleSlpzDDlQe5Ozvk5+UVFDGDwwR6dwcAAAAHdFiw="*)(*]BB*)+
        
        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*) (
      (
        (
          (*BB[*)(* here we add a lot of manual check for boundary condition *)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*)
          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} +

      (
        (
          (*BB[*)(* here we add a lot of manual check for boundary condition *)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*)
          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} +

      (
        (*BB[*)(If[i[[1]]-1 >= 1 && i[[1]]-1 <= max, grid[[i[[1]]-1, i[[2]] ]], {0.,0.}])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeB5AILqnMSXXKr0hjgskHleakFnMBGU6JydnpRfmleSlpzDDlQe5Ozvk5+UVFDGDwwR6dwcAAAAHdFiw="*)(*]BB*)
        
        + 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}}, 
(*BB[*)(* no big changes here *)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*)
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)
      )
    ]
  ]
]];

Тестируем

Получили ускорение почти на порядок - в 5 раз! Можно - лучше.

Компилируем в псевдо Си

Есть вариант, когда байт-код WL может быть преобразован в Си, затем скомпилирован доступным в системе gcc/clang. Для этого нужно лишь указать одну опцию

"CompilationTarget" -> "C"

Код

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
 ]] , CompilationOptions -> {"InlineExternalDefinitions" -> True}, 
    "CompilationTarget" -> "C", "RuntimeOptions" -> "Speed"];

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}]
], CompilationOptions -> {"InlineExternalDefinitions" -> True}, 
    "CompilationTarget" -> "C",   "RuntimeOptions" -> "Speed"];


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)
      )
    ]
  ]
] , CompilationOptions -> {"InlineExternalDefinitions" -> True}, 
    "CompilationTarget" -> "C","RuntimeOptions" -> "Speed"];

Тестируем

Ура, теперь честные 6 раз!

Основы программного рендеринга

Различают два вида взаимодействия с графическим API: retained mode и immediate mode. Все время до этого мы пользовались первым, когда наша программа пересылает вызовы отрисовки примитивов на плечи браузера и SVG. Если мы увеличим сетку жидкости в 10 или 100 раз, то очевидно, что гораздо выгоднее становится отображать поле скоростей или что-то связанное с ним напрямую как растровое изображение. Это означает, что мы будем передавать готовый кадр браузеру, нарисованный средствами нашей программы.

Начнем с простого примера

Table[Clip[x^2 + y^2 - 10^2, {0, 1}], {x,-10,10, 0.5}, {y,-10,10, 0.5}];
Image[%, Magnification -> 5, Antialiasing->False]

Здесь мы использовали так называемый метод SDF для рисования графики. Мы проходим по всем позициям пикселей в прямоугольнике, вычисляем значения яркости и затем отправляем их в Image, который копирует их в текстуру на GPU.

Мы можем использовать все 4 цветовых канала (RGBA) для рисования нашей картинки

Table[{
  Clip[(x^2+y^2) - 10^2, {0, 1}],
  Clip[(x^2+y^2) - 5^2, {0, 1}],
  Clip[(x^2+y^2) - 3^2, {0, 1}],
  Clip[(x^2+y^2) - 1^2, {0, 1}]
}, {x,-10,10, 0.5}, {y,-10,10, 0.5}];
Image[%, Magnification -> 5, Antialiasing->False]

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

Table[{
  Clip[(x^2+y^2) - 10^2, {0, 1}],
  Clip[(x^2+y^2) - 5^2, {0, 1}],
  Clip[(x^2+y^2) - 3^2, {0, 1}],
  Clip[(x^2+y^2) - 1^2, {0, 1}]
}, {x,-10,10, 0.5}, {y,-10,10, 0.5}];
NumericArray[%, "Real32"];
Image[%, "Real32", Magnification -> 5, Antialiasing->False]

Динамическая картинка

Есть стандартный способ обновления чего угодно в WLJS Notebook - используя технику Offload, как мы делали раньше. И Image его поддерживает

Совет
Для лучшей производительности лучше паковать все данные в UnsignedInteger8. Это даст серьезный выигрыш на передаче байтов и копировании в текстуру GPU.

serialize[arr_List] := NumericArray[arr, "UnsignedInteger8", "ClipAndRound"];

buffer = Table[255.0 {
  Clip[(x^2 + y^2) - 10^2, {0, 1}],
  Clip[(x^2 + y^2) - 5^2, {0, 1}],
  Clip[(x^2 + y^2) - 3^2, {0, 1}],
  Clip[(x^2 + y^2) - 1^2, {0, 1}]
}, {x, -10, 10, 0.5}, {y, -10, 10, 0.5}] // serialize;

Image[buffer // Offload, "Byte", Magnification -> 5, Antialiasing -> False]

Теперь мы можем легко обновлять наше изображение, записывая новые значения в buffer

t := AbsoluteTime[];

task = SetInterval[
  buffer = Table[255.0 {
    Clip[((x + Cos[t])^2 + (y - Sin[t])^2) - 3^2, {0, 1}],
    Clip[((x - Cos[t])^2 + (y - Sin[t])^2) - 5^2, {0, 1}],
    Clip[((x + Cos[t])^2 + (y + Sin[t])^2) - 10^2, {0, 1}],
    Clip[((x - Cos[t])^2 + (y + Sin[t])^2) - 1^2, {0, 1}]
  }, {x, -10, 10, 0.5}, {y, -10, 10, 0.5}] // serialize;
, 100];

SetTimeout[TaskRemove[task], 10000];

Результат верхней ячейки будет меняться со временем как-то так

Если вы дошли до этого места, поздравляю ⭐️ Теперь вы узнали, как писать программный рендер в блокноте WLJS Notebook!

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

Виртуальные чернила ✍️

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

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

И мы уже знаем, как эффективно его решать! Так как функция advect может принимать вектор в качестве поля, которое переносится, мы используем не один, а два вида чернил.

ink = Table[{0.,0.}, {i,50}, {j,50}];

nink = NumericArray[Map[255.0 {#[[1]], 0., #[[2]]} &, ink, {2}], "UnsignedInteger8", "ClipAndRound"];

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

vgrid = Table[{0.,0.}, {i,50}, {j,50}];

ink = 0. ink;
nink = NumericArray[Map[255.0 {#[[1]], 0., #[[2]]} &, ink, {2}], "UnsignedInteger8", "ClipAndRound"];

dest = {0,0};
cink = {1.0,0.2};
vfps = 0;

vtime = AbsoluteTime[];

EventHandler["vframe", Function[Null,
  
  vgrid = advect[vgrid,vgrid, 1.0];
  vgrid = removeDivergence[vgrid];
  vgrid = removeDivergence[vgrid];

  ink = With[{a = advect[vgrid, ink, 1.0]}, advect[vgrid, a, 1.0]];
  nink = NumericArray[Map[255.0 {1.0 - #[[2]], 1.0- #[[1]], 1.0 - #[[1]], 1.0} &, ink, {2}], "UnsignedInteger8", "ClipAndRound"];
  
  vfps = (*FB[*)(((vfps + 1 / (AbsoluteTime[] - vtime)))(*,*)/(*,*)(2.0))(*]FB*) // Round;
  vtime = AbsoluteTime[]; 
]];


With[{
  win = CurrentWindow[], 
  currentCell = ResultCell[]
},

  EventHandler[win, {"Closed" -> Function[Null,
    Delete[currentCell] 

  ]}];

  Graphics[{


    EventHandler[Graphics`Canvas[], {
      "click" -> Function[Null,
        cink = cink // Reverse;
      ],

      "mousemove" -> Function[pos, With[{
          xy = {50.0 - pos[[2]], pos[[1]]}
        }, 
          With[{p = Round[xy]},
            If[p[[1]] <= 50-1 && p[[1]] >=2 && p[[2]] <=50-1 && p[[2]] >=2,

                vgrid[[p[[1]],p[[2]]]] = Normalize[(xy - dest)];

                ink[[p[[1]],p[[2]]]] = cink;
                ink[[p[[1]]+1,p[[2]]]] = cink;
                ink[[p[[1]]-1,p[[2]]]] = cink;
                ink[[p[[1]],p[[2]]+1]] = cink;
                ink[[p[[1]],p[[2]]+1]] = cink;
            ];

          ];
        
          dest = xy;
        ] ]
    }], 

    AnimationFrameListener[vfps // Offload, "Event"->"vframe"], 

    Inset[
    
      Image[nink // Offload, "Byte", Magnification->10]
    
    , {0,0}, {0,0}, {500,500}
    ],
    
    Text[vfps // Offload, {0,0}]

    
  }, 
    Controls->False, 
    ImageSize->500, 
    PlotRange->{{0,50}, {0,50}}, 
    ImagePadding->None
  ]
]

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

Итак, результат

Следующим шагом логично будет перенести расчёты на GPU, используя вероятнее всего OpenCL, так как он пользуется большей поддержкой на совершенно разном железе (буквально Metal API в Apple Silicon).

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

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