Real-time Fluid Simulation Part 3

· 5 min read

Maintainer

Using Wolfram Language and WLJS

In this notebook we will apply some optimizations to the code, expand the resolution and switch to immediate mode of graphics rendering.

Download original notebook

Measure and optimize ⏱️

Let us try to estimate the time we need to do our normal calculations on divergence, advection and bilinear interpolation

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]

For a single run we have

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]}
    } // TableView 
  } // Column
];

runTest["Uncompiled"]

(*GB[*){{(*BB[*)("Uncompiled")(*,*)(*"1:eJxTTMoPSmNiYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCRYg4ZGfkwJRxgkkgkuKMvPSnfIritmAPM+SxJzM5EweIBOiBKQhqDQnNZgNrhYsFlJUmgoAykcZ4w=="*)(*]BB*)}(*||*),(*||*){(*VB[*)(FrontEndRef["bc61c228-82bc-43b5-a78b-b2829082984b"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKJyWbGSYbGVnoWhglJeuaGCeZ6iaaWyTpJhlZGFkaALGFSRIAgx0VRQ=="*)(*]VB*)}}(*]GB*)

Compile

Most probably for pure function JIT compiler should kick in, however, not all expressions are compilable in our case. We can make them by removing function declaration within the Module and replacing Piecewise with just a normal If statement.

It makes our code looks less readable ~~and removes all magic of WL~~, but the result worth it

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

Let us test it again!

runTest["Compiled"]

(*GB[*){{(*BB[*)("Compiled")(*,*)(*"1:eJxTTMoPSmNiYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCRYg4ZGfkwJRxgkkgkuKMvPSnfIritmAPM+SxJzM5EweIBOiBKQhqDQnNZgNrhYsFlJUmgoAykcZ4w=="*)(*]BB*)}(*||*),(*||*){(*VB[*)(FrontEndRef["695430be-0583-497f-8116-2c88bdb6b06d"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKm1mamhgbJKXqGphaGOuaWJqn6VoYGprpGiVbWCSlJJklGZilAAByCRUe"*)(*]VB*)}}(*]GB*)

x5 time faster

However, one should note, that this is still not a machine code, but rather byte-code of Wolfram Kernel internal command representation or something similar to that.

One can go futher and force WL to compile it to C and then link automatically. This will probably require to have gcc or clang installed on your system

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

Now we have new results

runTest["Compiled + C"]

(*GB[*){{(*BB[*)("Compiled + C")(*,*)(*"1:eJxTTMoPSmNiYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCRYg4ZGfkwJRxgkkgkuKMvPSnfIritmAPM+SxJzM5EweIBOiBKQhqDQnNZgNrhYsFlJUmgoAykcZ4w=="*)(*]BB*)}(*||*),(*||*){(*VB[*)(FrontEndRef["09dd1aaf-1abc-48e4-8407-2d04783abf9c"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKG1impBgmJqbpGiYmJeuaWKSa6FqYGJjrGqUYmJhbGCcmpVkmAwCN5RX9"*)(*]VB*)}}(*]GB*)

x6 time faster

We gained a bit of speed; we will see more difference once the grid is expanded and more iterations are needed.

Immediate Graphics Mode

This only means, that we will do all rendering of our primitives by ourself. The graphics API will only display the rendered frame we crafted on Wolfram Kernel. The last also meant we migh have to transfer much more data over WebSockets to frontend.

Let me show you a simple example

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

(*VB[*)(FrontEndRef["a38a8802-948f-46d8-9f9a-33a0bf6828dd"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKJxpbJFpYGBjpWppYpOmamKVY6FqmWSbqGhsnGiSlmVkYWaSkAAB9LBWL"*)(*]VB*)

Here we used so-called SDF method for drawing graphics. We iterate over all posiitons of pixels in reactangle, produce brightness values and then send them to Image, which copies them into a texture on GPU.

We can use all 4 (RGBA) color channels to draw our picture

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

(*VB[*)(FrontEndRef["37dbfc31-9a78-43de-bb3d-a7828a5639f6"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKG5unJKUlGxvqWiaaW+iaGKek6iYlGafoAnlGFommZsaWaWYAjRsV4g=="*)(*]VB*)

For the best performance use NumericArray, since it packs all numberic data to a more compact form

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

Dynamic image

There is a usual way on how to update our Image - using Offload technique as we did before

TIP

Real values takes more time to deserialize on frontend, we will use UnsignedInteger8

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

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

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

Now we can easily update our image by setting new values to buffer

t := AbsoluteTime[];

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

SetTimeout[TaskRemove[task], 10000];

If you made that far, congratulations ⭐️ Now you have learned how to do the software rendering.

Please, never do software rendering. It is slow and basically is a waste of resources of your CPU, which was not designed for rendering graphics. Use it only for educational purposes, small images or complex calulations, which are not possible to do using GPU.

Virtual Ink

To visualize velocity field, one can actually use another scalar field instead of 1000 probing balls. This scalar field is easy to imagine: ink ✍🏼 or dye or goo, which fell into a water and now is carried by the steams of fluid

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

We already know how to solve advection equation. Our function advect is used for modelling the momentum of fluid, which is a 2D vector field... Why not then to use two kinds of ink?

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

(*BB[*)(* transform to "byte" format *)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*)
nink = NumericArray[Map[255.0 {#[[1]], 0., #[[2]]} &, ink, {2}], "UnsignedInteger8", "ClipAndRound"];

Then we can directly visualize nink scalar field instead of drawing many many arrows and dots. In this way we will utilize fully our expensive software rendering.

TL;DR Final program

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;

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

  EventHandler[win, {"Closed" -> Function[Null,
    Delete[currentCell] (*BB[*)(* remove output cell if a notebook has been closed *)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*)
    (*BB[*)(* this will prevent the animation running uncontrollably on the next start *)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*)
  ]}];

  Graphics[{

    (*BB[*)(*attach listeners to a user's mouse to manipulate the grid*)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*)
    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,
                (*BB[*)(* accelerate the fluid *)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*)
                vgrid[[p[[1]],p[[2]]]] = Normalize[(xy - dest)];

                (*BB[*)(* add some ink *)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*)
                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;
        ] ]
    }], 

    (*BB[*)(*sync with browser's repaint cycle and update of fps label*)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*)
    AnimationFrameListener[vfps // Offload, "Event"->"vframe"], 

    (*BB[*)(*insert our Image*)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*)
    Inset[
    
      Image[nink // Offload, "Byte", Magnification->10]
    
    , {0,0}, {0,0}, {500,500}, 
    
      (*BB[*)(* do not apply any transformation. keep the original size *)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*)
      ViewMatrix->None
    ],
    
    Text[vfps // Offload, {0,0}]

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

(* subscribe to animation event *)

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

We do here pretty much the same as before; however, there are no longer arrows but a single Inset with a raster dynamic image inside. This places the Image on top of our Graphics canvas. This overlay is helpful since we can still listen to mouse positions and display FPS at a low cost. The fluid goes in the same direction as the mouse

I (Me - @JerryI) personally believe that we pushed the implementation to the limits of this toy model running on a CPU with an interpretive programming language. The next step definitely should be to utilize GPU compute shaders such as WebGPU, OpenCL, CUDA to calculate bigger fields and pipe the data directly to the canvas. However, this will be another story.

Thank you ❤️

Measure and optimize ⏱️​

Compile​

Immediate Graphics Mode​

Dynamic image​

Virtual Ink​

Measure and optimize ⏱️

Compile

Immediate Graphics Mode

Dynamic image

Virtual Ink