Real-time Fluid Simulation Part 2

· 9 min read

Maintainer

Using Wolfram Language and WLJS libraries

In this notebook, we will continue to explore a simple technique for simulating 2D incompressible fluids for visual effects. This work is mostly based on Jos Stam. Stable Fluids SIGGRAPH 1999 as well as a tutorial by Karl Sims

Download original notebook

Advection

From mathematica point of view it describes a general process of moving something on a scalar or vector field in the direction given by another vector field.

For simplicity's sake, we will first start with scalar fields.

grid = Table[{0,0}, {5}, {5}]; grid[[3,3]] = {0,1};
u = Table[1., {5}, {5}]; (*BB[*)(*our scalar field we will move*)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*)

One can think about u as if it were a mass of something, which will be carried by the stream grid.

Naive approach

Looking at the borders

One can consider interpolating the velocity at borders and then move the mass between them. We can write it as follows:

\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}

Warning: this is wrong

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

Now if we test it on our u field

arrow = Graphics[Arrow[{{0,0}, {1,0}}], ImagePadding->None, ImageSize->30];

Table[MatrixPlot[Nest[advect[grid, #]&, u, n], ImageSize->150], {n, {0, 1, 100}}];
Riffle[%, arrow] // Row 

(*GB[*){{(*VB[*)(FrontEndRef["e1451c50-456f-4076-9ed4-1aeca5bf0e88"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKpxqamBommxrompiapemaGJib6VqmppjoGiamJieaJqUZpFpYAAB7ihWk"*)(*]VB*)(*|*),(*|*)(*VB[*)(Graphics[Arrow[{{0, 0}, {1, 0}}], ImagePadding -> None, ImageSize -> 30])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KWnMIB4HkHAvSizIyEwuhsizAgnHoqL88jQmmHKfzOISVF4mkGYAE2jijJjiQaU5qcU8QIZnbmJ6akBiSkpmXjpYxi8/LxVNHSdMXXBmVWqmHJAHALiLJE4="*)(*]VB*)(*|*),(*|*)(*VB[*)(FrontEndRef["b238a126-d81c-40fe-8c34-f1d04780de66"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKJxkZWyQaGpnpplgYJuuaGKSl6lokG5vophmmGJiYWxikpJqZAQCAdxVw"*)(*]VB*)(*|*),(*|*)(*VB[*)(Graphics[Arrow[{{0, 0}, {1, 0}}], ImagePadding -> None, ImageSize -> 30])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KWnMIB4HkHAvSizIyEwuhsizAgnHoqL88jQmmHKfzOISVF4mkGYAE2jijJjiQaU5qcU8QIZnbmJ6akBiSkpmXjpYxi8/LxVNHSdMXXBmVWqmHJAHALiLJE4="*)(*]VB*)(*|*),(*|*)(*VB[*)(FrontEndRef["9b540d0e-6b20-4b41-91d4-cb4acf1fb194"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKWyaZmhikGKTqmiUZGeiaJJkY6loappjoJieZJCanGaYlGVqaAACB3RXB"*)(*]VB*)}}(*]GB*)

After 100 iterations, the space in the middle becomes negative. In our code

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

If org (which stands for the current cell) equals 0, it still does not mean that it cannot be reduced even more if we have a velocity vector to the nearest cell.

Donor-cell advection

The way to fix it is to use the Donor-cell method described well in Heidelberg University's lecture on numerical fluid dynamics. There we again consider velocities at borders, but based on the relative sign (outtake or intake) we decide which cell is going to be a donor.

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

Using these notations one can write the donor-cell method as follows:

\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}

Warning: this is correct

Adding the projections from other sides, our final function will be

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

Now if we run the same test, the result should be correct

arrow = Graphics[Arrow[{{0,0}, {1,0}}], ImagePadding->None, ImageSize->30];

Table[MatrixPlot[Nest[advect[grid, #]&, u, n], ImageSize->150], {n, {0, 1, 100}}];
Riffle[%, arrow] // Row 

(*GB[*){{(*VB[*)(FrontEndRef["e8fdf197-1b9a-4d20-bb58-7d80a23f0253"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKp1qkpaQZWprrGiZZJuqapBgZ6CYlmVromqdYGCQaGacZGJkaAwCNkhWi"*)(*]VB*)(*|*),(*|*)(*VB[*)(Graphics[Arrow[{{0, 0}, {1, 0}}], ImagePadding -> None, ImageSize -> 30])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KWnMIB4HkHAvSizIyEwuhsizAgnHoqL88jQmmHKfzOISVF4mkGYAE2jijJjiQaU5qcU8QIZnbmJ6akBiSkpmXjpYxi8/LxVNHSdMXXBmVWqmHJAHALiLJE4="*)(*]VB*)(*|*),(*|*)(*VB[*)(FrontEndRef["f41c04c8-2eef-469c-bf5d-4f8fe8dd9191"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKp5kYJhuYJFvoGqWmpumamFkm6yalmabomqRZpKVapKRYGloaAgCQbhZK"*)(*]VB*)(*|*),(*|*)(*VB[*)(Graphics[Arrow[{{0, 0}, {1, 0}}], ImagePadding -> None, ImageSize -> 30])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KWnMIB4HkHAvSizIyEwuhsizAgnHoqL88jQmmHKfzOISVF4mkGYAE2jijJjiQaU5qcU8QIZnbmJ6akBiSkpmXjpYxi8/LxVNHSdMXXBmVWqmHJAHALiLJE4="*)(*]VB*)(*|*),(*|*)(*VB[*)(FrontEndRef["086251cd-1038-4ff8-989d-cbc828c24c25"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKG1iYGZkaJqfoGhoYW+iapKVZ6FpaWKboJiclWxhZJBuZJBuZAgB5KRVP"*)(*]VB*)}}(*]GB*)

We can check it numerically as well

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

((*GB[*){{1.`(*|*),(*|*)1.`(*|*),(*|*)1.`(*|*),(*|*)1.`(*|*),(*|*)1.`}(*||*),(*||*){1.`(*|*),(*|*)1.`(*|*),(*|*)1.`(*|*),(*|*)1.`(*|*),(*|*)1.`}(*||*),(*||*){1.`(*|*),(*|*)0.005920529220334025`(*|*),(*|*)0.0370812093273552`(*|*),(*|*)2.9569982614523105`(*|*),(*|*)1.`}(*||*),(*||*){1.`(*|*),(*|*)1.`(*|*),(*|*)1.`(*|*),(*|*)1.`(*|*),(*|*)1.`}(*||*),(*||*){1.`(*|*),(*|*)1.`(*|*),(*|*)1.`(*|*),(*|*)1.`(*|*),(*|*)1.`}}(*]GB*))

or in a more creative way

Map[Function[cell, RGBColor[Tanh[cell], 0,-Tanh[cell]]], Nest[advect[grid, #]&, u, 100], {2}];
% // MatrixForm

((*GB[*){{(*VB[*)(RGBColor[0.7615941559557649, 0, -0.7615941559557649])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeGJAIcndyzs/JLwouTyxJzghJzS3ISSxJTWMGyXMgyRdN+Syy9VfcC/siBiiAiewHAHfzHEw="*)(*]VB*)(*|*),(*|*)(*VB[*)(RGBColor[0.7615941559557649, 0, -0.7615941559557649])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeGJAIcndyzs/JLwouTyxJzghJzS3ISSxJTWMGyXMgyRdN+Syy9VfcC/siBiiAiewHAHfzHEw="*)(*]VB*)(*|*),(*|*)(*VB[*)(RGBColor[0.7615941559557649, 0, -0.7615941559557649])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeGJAIcndyzs/JLwouTyxJzghJzS3ISSxJTWMGyXMgyRdN+Syy9VfcC/siBiiAiewHAHfzHEw="*)(*]VB*)(*|*),(*|*)(*VB[*)(RGBColor[0.7615941559557649, 0, -0.7615941559557649])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeGJAIcndyzs/JLwouTyxJzghJzS3ISSxJTWMGyXMgyRdN+Syy9VfcC/siBiiAiewHAHfzHEw="*)(*]VB*)(*|*),(*|*)(*VB[*)(RGBColor[0.7615941559557649, 0, -0.7615941559557649])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeGJAIcndyzs/JLwouTyxJzghJzS3ISSxJTWMGyXMgyRdN+Syy9VfcC/siBiiAiewHAHfzHEw="*)(*]VB*)}(*||*),(*||*){(*VB[*)(RGBColor[0.7615941559557649, 0, -0.7615941559557649])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeGJAIcndyzs/JLwouTyxJzghJzS3ISSxJTWMGyXMgyRdN+Syy9VfcC/siBiiAiewHAHfzHEw="*)(*]VB*)(*|*),(*|*)(*VB[*)(RGBColor[0.7615941559557649, 0, -0.7615941559557649])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeGJAIcndyzs/JLwouTyxJzghJzS3ISSxJTWMGyXMgyRdN+Syy9VfcC/siBiiAiewHAHfzHEw="*)(*]VB*)(*|*),(*|*)(*VB[*)(RGBColor[0.7615941559557649, 0, -0.7615941559557649])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeGJAIcndyzs/JLwouTyxJzghJzS3ISSxJTWMGyXMgyRdN+Syy9VfcC/siBiiAiewHAHfzHEw="*)(*]VB*)(*|*),(*|*)(*VB[*)(RGBColor[0.7615941559557649, 0, -0.7615941559557649])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeGJAIcndyzs/JLwouTyxJzghJzS3ISSxJTWMGyXMgyRdN+Syy9VfcC/siBiiAiewHAHfzHEw="*)(*]VB*)(*|*),(*|*)(*VB[*)(RGBColor[0.7615941559557649, 0, -0.7615941559557649])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeGJAIcndyzs/JLwouTyxJzghJzS3ISSxJTWMGyXMgyRdN+Syy9VfcC/siBiiAiewHAHfzHEw="*)(*]VB*)}(*||*),(*||*){(*VB[*)(RGBColor[0.7615941559557649, 0, -0.7615941559557649])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeGJAIcndyzs/JLwouTyxJzghJzS3ISSxJTWMGyXMgyRdN+Syy9VfcC/siBiiAiewHAHfzHEw="*)(*]VB*)(*|*),(*|*)(*VB[*)(RGBColor[0.005920460044525685, 0, -0.005920460044525685])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeGJAIcndyzs/JLwouTyxJzghJzS3ISSxJTWMGyXMgyRedvOydwutQYV/EAAUwkf0AVBEZTA=="*)(*]VB*)(*|*),(*|*)(*VB[*)(RGBColor[0.037064222916834776, 0, -0.037064222916834776])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeGJAIcndyzs/JLwouTyxJzghJzS3ISSxJTWMGyXMgyRc5i/Z/Evm1yL6IAQpgIvsBV8caPg=="*)(*]VB*)(*|*),(*|*)(*VB[*)(RGBColor[0.9946118170687553, 0, -0.9946118170687553])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeGJAIcndyzs/JLwouTyxJzghJzS3ISSxJTWMGyXMgyRelXgjQvHP5vX0RAxTARPYDAG2TG8Q="*)(*]VB*)(*|*),(*|*)(*VB[*)(RGBColor[0.7615941559557649, 0, -0.7615941559557649])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeGJAIcndyzs/JLwouTyxJzghJzS3ISSxJTWMGyXMgyRdN+Syy9VfcC/siBiiAiewHAHfzHEw="*)(*]VB*)}(*||*),(*||*){(*VB[*)(RGBColor[0.7615941559557649, 0, -0.7615941559557649])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeGJAIcndyzs/JLwouTyxJzghJzS3ISSxJTWMGyXMgyRdN+Syy9VfcC/siBiiAiewHAHfzHEw="*)(*]VB*)(*|*),(*|*)(*VB[*)(RGBColor[0.7615941559557649, 0, -0.7615941559557649])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeGJAIcndyzs/JLwouTyxJzghJzS3ISSxJTWMGyXMgyRdN+Syy9VfcC/siBiiAiewHAHfzHEw="*)(*]VB*)(*|*),(*|*)(*VB[*)(RGBColor[0.7615941559557649, 0, -0.7615941559557649])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeGJAIcndyzs/JLwouTyxJzghJzS3ISSxJTWMGyXMgyRdN+Syy9VfcC/siBiiAiewHAHfzHEw="*)(*]VB*)(*|*),(*|*)(*VB[*)(RGBColor[0.7615941559557649, 0, -0.7615941559557649])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeGJAIcndyzs/JLwouTyxJzghJzS3ISSxJTWMGyXMgyRdN+Syy9VfcC/siBiiAiewHAHfzHEw="*)(*]VB*)(*|*),(*|*)(*VB[*)(RGBColor[0.7615941559557649, 0, -0.7615941559557649])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeGJAIcndyzs/JLwouTyxJzghJzS3ISSxJTWMGyXMgyRdN+Syy9VfcC/siBiiAiewHAHfzHEw="*)(*]VB*)}(*||*),(*||*){(*VB[*)(RGBColor[0.7615941559557649, 0, -0.7615941559557649])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeGJAIcndyzs/JLwouTyxJzghJzS3ISSxJTWMGyXMgyRdN+Syy9VfcC/siBiiAiewHAHfzHEw="*)(*]VB*)(*|*),(*|*)(*VB[*)(RGBColor[0.7615941559557649, 0, -0.7615941559557649])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeGJAIcndyzs/JLwouTyxJzghJzS3ISSxJTWMGyXMgyRdN+Syy9VfcC/siBiiAiewHAHfzHEw="*)(*]VB*)(*|*),(*|*)(*VB[*)(RGBColor[0.7615941559557649, 0, -0.7615941559557649])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeGJAIcndyzs/JLwouTyxJzghJzS3ISSxJTWMGyXMgyRdN+Syy9VfcC/siBiiAiewHAHfzHEw="*)(*]VB*)(*|*),(*|*)(*VB[*)(RGBColor[0.7615941559557649, 0, -0.7615941559557649])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeGJAIcndyzs/JLwouTyxJzghJzS3ISSxJTWMGyXMgyRdN+Syy9VfcC/siBiiAiewHAHfzHEw="*)(*]VB*)(*|*),(*|*)(*VB[*)(RGBColor[0.7615941559557649, 0, -0.7615941559557649])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeGJAIcndyzs/JLwouTyxJzghJzS3ISSxJTWMGyXMgyRdN+Syy9VfcC/siBiiAiewHAHfzHEw="*)(*]VB*)}}(*]GB*))

Now we are ready to apply this knowledge to our fluid simulation. But how?

Advection differential equation

I would like to draw your attention to the differential equation of advection

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

If one tries to descritize it on 1D grid

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

where $f_{q \pm 1/2}$ we define flux at the interface of $q$ th cell. The exact form if a flux is given by the algorythm used. For above-mentioned donor-cell we find

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}

After substituting it back, one finds, that it matches our previous equation we found for a problem of moving "mass" on a grid. ==We solved this exact differnetial advection equation while implementing advect function.==

Fluid Momentum

Newton's second law in the Navier-Stokes equation is defined by the following term.

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

Does it remind you anything?.. The conservation of momentum of our fluid is actual self advection.

To simulate fluid momentum, the flow field can simply flow itself. Flow vectors are updated by reading interpolated values from upstream grid cells in the same way the tracer image is advected above. New vector values are pulled from the negative direction of flow.

To test it, one can construct the same sandbox used in Part 1, but change the update loop slightly.

vgrid = Table[{0.,0.}, {i,15}, {j,15}];
vcolors = Table[1.0, {Length[vgrid]}, {Length[vgrid]}];
vfps = 0;

add mouse handler

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]

Subscribe for events

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

Assemble the widget

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

In the end it should be similar to what is shown on the animation

Test particles

In order to better visualize the fluid, we need something to be moved by our stream. The simplest solution is to spread hundreds of tracer particles. The velocity field can push them in the direction interpolated bilinearly from the grid.

Bilinear interpolation

This is the most basic thing to implement. This probably is not the fastest version, but good enough for educational purposes.

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

One can check if it works correctly by scaling up some raster images

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 

(*GB[*){{(*VB[*)(FrontEndRef["a2a331ed-cc9c-4b80-9377-d58dd2533045"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKJxolGhsbpqboJidbJuuaJFkY6Foam5vrpphapKQYmRobG5iYAgCKcRVw"*)(*]VB*)(*|*),(*|*)(*VB[*)(FrontEndRef["c2d79718-ca37-454d-8357-b2c9a1224b84"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKJxulmFuaG1roJicam+uamJqk6FoYm5rrJhklWyYaGhmZJFmYAAB9KBUa"*)(*]VB*)}}(*]GB*)

On the right, there is our interpolated array, while on the left, the original one

Uniform distribution

There is another short problem to solve: how to distribute particles evenly over a defined region. However, there is a solution out of the box in the Wolfram Language standard library.

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

(*VB[*)(FrontEndRef["1b5a0b27-3686-4aa2-a7a8-ce20c1a34c56"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKGyaZJhokGZnrGptZmOmaJCYa6SaaJ1roJqcaGSQbJhqbJJuaAQCBHhWL"*)(*]VB*)

The last thing is to implement a function, which updates all probe particles

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

Complete Real-Time fluid simulation

Wrapping everything up together with the divergence solver from the previous part, we can finally assemble the model

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

TL;DR Here is the final code

Evaluate the cell below in the notebook. Use your mouse to draw velocity vectors.

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


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

Use your mouse to draw velocity vectors

There is still a lot that can be improved further. Especially advection and the code for removing divergence can be optimized and compiled as well.

Optimizations

Since all our main functions used in the animation loop are dealing only with numerics, one can apply Compile to them. With some minor changes made (which unfortunately make the code a bit ugly, because not all functions are compilable) we have

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

and

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

After that the main bottleneck is still the way how we draw and update SVG graphics (as well as WebSocket connection latency). However, these changes improves average FPS from 17-24 up to 30-35 on my Mac Air M1.

Advection​

Naive approach​

Donor-cell advection​

Advection differential equation​

Fluid Momentum​

Test particles​

Bilinear interpolation​

Uniform distribution​

Complete Real-Time fluid simulation​

Optimizations​