Real-time Fluid Simulation Part 1

· 7 min read

Maintainer

Using Wolfram Language and WLJS libraries

In this notebook, we will 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

Working with a Grid

Here we will use the Euler method to represent a fluid, i.e., by storing velocity as a vector field discretized to a uniform grid.

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

((*GB[*){{((*GB[*){{1}(*||*),(*||*){-1}}(*]GB*))(*|*),(*|*)((*GB[*){{2}(*||*),(*||*){-1}}(*]GB*))(*|*),(*|*)((*GB[*){{3}(*||*),(*||*){-1}}(*]GB*))(*|*),(*|*)((*GB[*){{4}(*||*),(*||*){-1}}(*]GB*))(*|*),(*|*)((*GB[*){{5}(*||*),(*||*){-1}}(*]GB*))}(*||*),(*||*){((*GB[*){{1}(*||*),(*||*){-2}}(*]GB*))(*|*),(*|*)((*GB[*){{2}(*||*),(*||*){-2}}(*]GB*))(*|*),(*|*)((*GB[*){{3}(*||*),(*||*){-2}}(*]GB*))(*|*),(*|*)((*GB[*){{4}(*||*),(*||*){-2}}(*]GB*))(*|*),(*|*)((*GB[*){{5}(*||*),(*||*){-2}}(*]GB*))}(*||*),(*||*){((*GB[*){{1}(*||*),(*||*){-3}}(*]GB*))(*|*),(*|*)((*GB[*){{2}(*||*),(*||*){-3}}(*]GB*))(*|*),(*|*)((*GB[*){{3}(*||*),(*||*){-3}}(*]GB*))(*|*),(*|*)((*GB[*){{4}(*||*),(*||*){-3}}(*]GB*))(*|*),(*|*)((*GB[*){{5}(*||*),(*||*){-3}}(*]GB*))}(*||*),(*||*){((*GB[*){{1}(*||*),(*||*){-4}}(*]GB*))(*|*),(*|*)((*GB[*){{2}(*||*),(*||*){-4}}(*]GB*))(*|*),(*|*)((*GB[*){{3}(*||*),(*||*){-4}}(*]GB*))(*|*),(*|*)((*GB[*){{4}(*||*),(*||*){-4}}(*]GB*))(*|*),(*|*)((*GB[*){{5}(*||*),(*||*){-4}}(*]GB*))}(*||*),(*||*){((*GB[*){{1}(*||*),(*||*){-5}}(*]GB*))(*|*),(*|*)((*GB[*){{2}(*||*),(*||*){-5}}(*]GB*))(*|*),(*|*)((*GB[*){{3}(*||*),(*||*){-5}}(*]GB*))(*|*),(*|*)((*GB[*){{4}(*||*),(*||*){-5}}(*]GB*))(*|*),(*|*)((*GB[*){{5}(*||*),(*||*){-5}}(*]GB*))}}(*]GB*))

Then one can easily visualize it as a vector field:

grid // ListVectorPlot

(*VB[*)(FrontEndRef["7fa3d105-973d-43df-aa2e-88d286b02c82"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKm6clGqcYGpjqWpobp+iaGKek6SYmGqXqWlikGFmYJRkYJVsYAQCGjBWk"*)(*]VB*)

For the best performance, we will use Map or MapIndexed with a pure function inside. Then Wolfram Kernel is able to use JIT compilation. For example:

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

{{{1,1},{1,2},{1,3},{1,4},{1,5}},{{2,1},{2,2},{2,3},{2,4},{2,5}},{{3,1},{3,2},{3,3},{3,4},{3,5}},{{4,1},{4,2},{4,3},{4,4},{4,5}},{{5,1},{5,2},{5,3},{5,4},{5,5}}}

(*VB[*)(FrontEndRef["9d29d953-9bbb-46aa-805e-391ff5054311"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKW6YYWaZYmhrrWiYlJemamCUm6loYmKbqGlsapqWZGpiaGBsaAgCEmxVJ"*)(*]VB*)

Custom Syntax Sugar

Another possibility to visualize using pure syntax sugar would be to define our own output form for a single velocity vector.

The most efficient way is to define a decorator function in Javascript. Let it be a small arrow pointing to a desired direction.

.js
core.ArrowGuide = async (args, env) => {
  const dir = await interpretate(args[0], env);
  const mag = dir.map((el)=>el*el).reduce((c,a) => c+a);

  const angle = Math.round(180 * Math.atan2(dir[0], dir[1]) / Math.PI);
  
  console.log(angle);
  env.element.style.transform = `rotate(${angle}deg)`;
  
  if (mag < 0.01) { //do not display if it is too small
    env.element.style.opacity = 0.5;
    env.element.innerHTML = ".";
  } else {
    env.element.innerHTML = "&uarr;";
  }
}

core.ArrowGuide = async (args, env) => {
  const dir = await interpretate(args[0], env);
  const mag = dir.map((el)=>el*el).reduce((c,a) => c+a);

  const angle = Math.round(180 * Math.atan2(dir[0], dir[1]) / Math.PI);
  
  console.log(angle);
  env.element.style.transform = `rotate(${angle}deg)`;
  
  if (mag < 0.01) { //do not display if it is too small
    env.element.style.opacity = 0.5;
    env.element.innerHTML = ".";
  } else {
    env.element.innerHTML = "&uarr;";
  }
}

Now we define a standard form ArrowGuide, which will use the defined Javascript function to display itself inside a cell.

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

The first argument will be an actual expression used as input; First will break our wrapper and substitute an original vector.

ArrowGuide[{1,1}]

(*VB[*)({1, 1})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBKkAEwCJmgvW"*)(*]VB*)

Now one can do some cool tricks like this one:

(*VB[*)({1, 1})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBKkAEwCJmgvW"*)(*]VB*) . (*VB[*)({1, 0})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBKnIBGIGAImWC9U="*)(*]VB*)

Now we can basically apply this decoration wrapper to our matrix using multidimensional Map function and then again output it as a matrix

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

((*GB[*){{(*VB[*)({1, -1})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBKnI/A8EAJOMD9E="*)(*]VB*)(*|*),(*|*)(*VB[*)({2, -1})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBPEy/wMBAJOVD9I="*)(*]VB*)(*|*),(*|*)(*VB[*)({3, -1})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JZAYyMv8DAQCTng/T"*)(*]VB*)(*|*),(*|*)(*VB[*)({4, -1})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBDEy/wMBAJOnD9Q="*)(*]VB*)(*|*),(*|*)(*VB[*)({5, -1})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JZAUyMv8DAQCTsA/V"*)(*]VB*)}(*||*),(*||*){(*VB[*)({1, -2})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBKnI/Pf//38Ak4gP0A=="*)(*]VB*)(*|*),(*|*)(*VB[*)({2, -2})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBPEy//3//x8Ak5EP0Q=="*)(*]VB*)(*|*),(*|*)(*VB[*)({3, -2})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JZAYyMv/9//8fAJOaD9I="*)(*]VB*)(*|*),(*|*)(*VB[*)({4, -2})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBDEy//3//x8Ak6MP0w=="*)(*]VB*)(*|*),(*|*)(*VB[*)({5, -2})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JZAUyMv/9//8fAJOsD9Q="*)(*]VB*)}(*||*),(*||*){(*VB[*)({1, -3})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBKnI/Pv//38Ak4QPzw=="*)(*]VB*)(*|*),(*|*)(*VB[*)({2, -3})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBPEy//7//x8Ak40P0A=="*)(*]VB*)(*|*),(*|*)(*VB[*)({3, -3})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JZAYyMv/+//8fAJOWD9E="*)(*]VB*)(*|*),(*|*)(*VB[*)({4, -3})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBDEy//7//x8Ak58P0g=="*)(*]VB*)(*|*),(*|*)(*VB[*)({5, -3})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JZAUyMv/+//8fAJOoD9M="*)(*]VB*)}(*||*),(*||*){(*VB[*)({1, -4})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBKnI/PP//38Ak4APzg=="*)(*]VB*)(*|*),(*|*)(*VB[*)({2, -4})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBPEy//z//x8Ak4kPzw=="*)(*]VB*)(*|*),(*|*)(*VB[*)({3, -4})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JZAYyMv/8//8fAJOSD9A="*)(*]VB*)(*|*),(*|*)(*VB[*)({4, -4})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBDEy//z//x8Ak5sP0Q=="*)(*]VB*)(*|*),(*|*)(*VB[*)({5, -4})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JZAUyMv/8//8fAJOkD9I="*)(*]VB*)}(*||*),(*||*){(*VB[*)({1, -5})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBKnI/P3//38Ak3wPzQ=="*)(*]VB*)(*|*),(*|*)(*VB[*)({2, -5})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBPEyf////x8Ak4UPzg=="*)(*]VB*)(*|*),(*|*)(*VB[*)({3, -5})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JZAYyMn////8fAJOOD88="*)(*]VB*)(*|*),(*|*)(*VB[*)({4, -5})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBDEyf////x8Ak5cP0A=="*)(*]VB*)(*|*),(*|*)(*VB[*)({5, -5})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JZAUyMn////8fAJOgD9E="*)(*]VB*)}}(*]GB*))

Let's define a helper function

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

Divergence

Conservation of mass for incompressible fluid dictates

div ~\mathbf{v} = 0

where $\mathbf{v}$ is velocity vector field we played with earlier - grid.

Example

Let us start with a simples example

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

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

And put a single vector in the middle

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

To visualize the divergence, one can use DensityPlot

{intX, intY} = {
    ListInterpolation[grid[[All,All,1]]],
    ListInterpolation[grid[[All,All,2]]]
}; (*BB[*)(*interpolate data for the convenience*)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*)

CoolColor[ z_ ] := RGBColor[z, 0.5, 1 - z];

(*BB[*)(*find divergence*)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*)
Div[{intX[x,y], intY[x,y]}, {x,y}]; 
ContourPlot[%, {x,1,5}, {y,1,5}, ColorFunction->CoolColor]

(*VB[*)(FrontEndRef["15f4029b-3d6d-46e2-8893-95deb4f9c83e"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKG5qmmRgYWSbpGqeYpeiamKUa6VpYWBrrWpqmpCaZpFkmWxinAgB8URWN"*)(*]VB*)

From a physical point of view, it means that there is a source node at the bottom and a drain at the top, which should not be a feature of a closed system with incompressible fluid.

Removing the Divergence

One can try to solve this problem iteratively by taking a field with non-zero divergency and apply artificial tranformation on it to make it satisfy the equation.

There is a clever algorithm for removing the divergence of an arbitrary 2D vector field, which acts like a cellular automaton

\begin{align*} v^\prime_{(x=0,y=0)} &= v_{(0,0)} + \frac{1}{8} \Big\{ \\ &\Big((v_{(-1,-1)} + v_{(1,1)})\cdot \begin{bmatrix} 1 \\ 1 \end{bmatrix} \Big) \begin{bmatrix} 1 \\ 1 \end{bmatrix} \Big\} + \\ &\Big((v_{(-1,1)} + v_{(1,-1)})\cdot \begin{bmatrix} 1 \\ -1 \end{bmatrix} \Big) \begin{bmatrix} 1 \\ -1 \end{bmatrix} \Big\} + \\ \begin{bmatrix} 2 & 0\\ 0 & -2 \end{bmatrix} &\cdot (v_{(-1,0)} + v_{(1,0)} - v_{(0,-1)} - v_{(0,1)}) + \\ & (-4) v_{(0,0)} \Big\}, \end{align*}

where $v^{\prime}_{(0,0)}$ is a new velocity vector at coordinates $(x,y)$ and others $v_{(i,j)}$ are ones from the previous state at $(x + i, y + j)$ .

This is done by adding flow away from high pressure converging areas, and toward low pressure diverging areas.

Here is our implementation

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

Let us test it.

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

(*GB[*){{((*GB[*){{(*VB[*)({0, 0})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMYAIAiY0L1A=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0, 0})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMYAIAiY0L1A=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0, 0})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMYAIAiY0L1A=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0, 0})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMYAIAiY0L1A=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0, 0})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMYAIAiY0L1A=="*)(*]VB*)}(*||*),(*||*){(*VB[*)({0, 0})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMYAIAiY0L1A=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0, 0})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMYAIAiY0L1A=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0, 0})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMYAIAiY0L1A=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0, 0})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMYAIAiY0L1A=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0, 0})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMYAIAiY0L1A=="*)(*]VB*)}(*||*),(*||*){(*VB[*)({0, 0})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMYAIAiY0L1A=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0, 0})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMYAIAiY0L1A=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0, 1})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMmSBlAImRC9U="*)(*]VB*)(*|*),(*|*)(*VB[*)({0, 0})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMYAIAiY0L1A=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0, 0})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMYAIAiY0L1A=="*)(*]VB*)}(*||*),(*||*){(*VB[*)({0, 0})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMYAIAiY0L1A=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0, 0})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMYAIAiY0L1A=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0, 0})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMYAIAiY0L1A=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0, 0})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMYAIAiY0L1A=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0, 0})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMYAIAiY0L1A=="*)(*]VB*)}(*||*),(*||*){(*VB[*)({0, 0})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMYAIAiY0L1A=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0, 0})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMYAIAiY0L1A=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0, 0})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMYAIAiY0L1A=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0, 0})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMYAIAiY0L1A=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0, 0})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC7JBNIMYAIAiY0L1A=="*)(*]VB*)}}(*]GB*))(*|*),(*|*)((*GB[*){{(*VB[*)({0., 0.})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYoACOAMA53wL5g=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0., 0.})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYoACOAMA53wL5g=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0., 0.})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYoACOAMA53wL5g=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0., 0.})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYoACOAMA53wL5g=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0., 0.})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYoACOAMA53wL5g=="*)(*]VB*)}(*||*),(*||*){(*VB[*)({0., 0.})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYoACOAMA53wL5g=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0.125, 0.125})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYgCDA/ZwBgDz8Q3k"*)(*]VB*)(*|*),(*|*)(*VB[*)({0., -0.25})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYoACKOPCfgDp2w11"*)(*]VB*)(*|*),(*|*)(*VB[*)({-0.125, 0.125})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYgCDA/thDHsA+PEOZA=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0., 0.})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYoACOAMA53wL5g=="*)(*]VB*)}(*||*),(*||*){(*VB[*)({0., 0.})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYoACOAMA53wL5g=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0., 0.25})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYoACKOOCPQDpWwz1"*)(*]VB*)(*|*),(*|*)(*VB[*)({0., 0.5})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYoACKOOBPQDpew0F"*)(*]VB*)(*|*),(*|*)(*VB[*)({0., 0.25})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYoACKOOCPQDpWwz1"*)(*]VB*)(*|*),(*|*)(*VB[*)({0., 0.})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYoACOAMA53wL5g=="*)(*]VB*)}(*||*),(*||*){(*VB[*)({0., 0.})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYoACOAMA53wL5g=="*)(*]VB*)(*|*),(*|*)(*VB[*)({-0.125, 0.125})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYgCDA/thDHsA+PEOZA=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0., -0.25})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYoACKOPCfgDp2w11"*)(*]VB*)(*|*),(*|*)(*VB[*)({0.125, 0.125})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYgCDA/ZwBgDz8Q3k"*)(*]VB*)(*|*),(*|*)(*VB[*)({0., 0.})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYoACOAMA53wL5g=="*)(*]VB*)}(*||*),(*||*){(*VB[*)({0., 0.})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYoACOAMA53wL5g=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0., 0.})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYoACOAMA53wL5g=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0., 0.})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYoACOAMA53wL5g=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0., 0.})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYoACOAMA53wL5g=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0., 0.})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYoACOAMA53wL5g=="*)(*]VB*)}}(*]GB*))(*|*),(*|*)((*GB[*){{(*VB[*)({0.03125, 0.03125})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYgCDBfZwBgDyUQ2k"*)(*]VB*)(*|*),(*|*)(*VB[*)({0., -0.0625})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYoACKGPDfgDpmw1V"*)(*]VB*)(*|*),(*|*)(*VB[*)({0., 0.0625})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYoACKGODPQDpGwzV"*)(*]VB*)(*|*),(*|*)(*VB[*)({0., -0.0625})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYoACKGPDfgDpmw1V"*)(*]VB*)(*|*),(*|*)(*VB[*)({-0.03125, 0.03125})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYgCDBfthDHsA91EOJA=="*)(*]VB*)}(*||*),(*||*){(*VB[*)({0., 0.0625})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYoACKGODPQDpGwzV"*)(*]VB*)(*|*),(*|*)(*VB[*)({0.125, 0.})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYgCDA/ZQBgMDAPIyDOU="*)(*]VB*)(*|*),(*|*)(*VB[*)({0., -0.125})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYoACKOPAfgDpuw1l"*)(*]VB*)(*|*),(*|*)(*VB[*)({-0.125, 0.})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYgCDA/uhDAYGAPcyDWU="*)(*]VB*)(*|*),(*|*)(*VB[*)({0., 0.0625})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYoACKGODPQDpGwzV"*)(*]VB*)}(*||*),(*||*){(*VB[*)({0., 0.0625})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYoACKGODPQDpGwzV"*)(*]VB*)(*|*),(*|*)(*VB[*)({0., 0.125})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYoACKOOAPQDpOwzl"*)(*]VB*)(*|*),(*|*)(*VB[*)({0., 0.625})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYoACKOOJPQDpgw0J"*)(*]VB*)(*|*),(*|*)(*VB[*)({0., 0.125})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYoACKOOAPQDpOwzl"*)(*]VB*)(*|*),(*|*)(*VB[*)({0., 0.0625})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYoACKGODPQDpGwzV"*)(*]VB*)}(*||*),(*||*){(*VB[*)({0., 0.0625})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYoACKGODPQDpGwzV"*)(*]VB*)(*|*),(*|*)(*VB[*)({-0.125, 0.})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYgCDA/uhDAYGAPcyDWU="*)(*]VB*)(*|*),(*|*)(*VB[*)({0., -0.125})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYoACKOPAfgDpuw1l"*)(*]VB*)(*|*),(*|*)(*VB[*)({0.125, 0.})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYgCDA/ZQBgMDAPIyDOU="*)(*]VB*)(*|*),(*|*)(*VB[*)({0., 0.0625})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYoACKGODPQDpGwzV"*)(*]VB*)}(*||*),(*||*){(*VB[*)({-0.03125, 0.03125})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYgCDBfthDHsA91EOJA=="*)(*]VB*)(*|*),(*|*)(*VB[*)({0., -0.0625})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYoACKGPDfgDpmw1V"*)(*]VB*)(*|*),(*|*)(*VB[*)({0., 0.0625})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYoACKGODPQDpGwzV"*)(*]VB*)(*|*),(*|*)(*VB[*)({0., -0.0625})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYoACKGPDfgDpmw1V"*)(*]VB*)(*|*),(*|*)(*VB[*)({0.03125, 0.03125})(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeF5BwLCrKL3cvzUxJTWOCqfDJLC4pYgCDBfZwBgDyUQ2k"*)(*]VB*)}}(*]GB*))}(*||*),(*||*)}(*]GB*)

As you can see, if we have a stream of fluid from A to B, there must be some other streams made from B to A. The applied procedure produces the missing one spread over the whole grid if removeDivergence is applied multiple times.

We can check it directly by calculating the divergence.

CoolColor[ z_ ] := RGBColor[z, 0.5, 1 - z];

plotDiv[grid_] := Module[{intX, intY, x, y},

  {intX, intY} = {
      ListInterpolation[grid[[All,All,1]]],
      ListInterpolation[grid[[All,All,2]]]
  }; (*BB[*)(*interpolate data for the convenience*)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*)

  With[{div = Div[{intX[x,y], intY[x,y]}, {x,y}]},
    ContourPlot[
      div, {x,1,5}, {y,1,5}, 
      ColorFunctionScaling -> False, (*BB[*)(*make sure to have the same scale*)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*)
      ColorFunction -> ColorData[{"ThermometerColors", {-1, 1}}]
    ]
  ]
]

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

(*GB[*){{(*BB[*)("Original")(*,*)(*"1:eJxTTMoPSmNiYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCRYg4ZGfkwLhcQKJ4JKizLx0p/yKTD4gD6IdpCqoNCc1mA2uACwWUlSaCgCn0xcV"*)(*]BB*)(*|*),(*|*)(*BB[*)("8 iterations")(*,*)(*"1:eJxTTMoPSmNiYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCRYg4ZGfkwLhcQKJ4JKizLx0p/yKTD4gD6IdpCqoNCc1mA2uACwWUlSaCgCn0xcV"*)(*]BB*)}(*||*),(*||*){(*VB[*)(FrontEndRef["62b4aa00-b857-4b7d-8530-a563afd164ec"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKmxklmSQmGhjoJlmYmuuaJJmn6FqYGhvoJpqaGSempRiamaQmAwB/dRWX"*)(*]VB*)(*|*),(*|*)(*VB[*)(FrontEndRef["b99a3f97-955e-4d38-b5cb-4ec0903c77e3"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKJ1laJhqnWZrrWpqapuqapBhb6CaZJifpmqQmG1gaGCebm6caAwCG1hW5"*)(*]VB*)}}(*]GB*)

Apart from obvious artifacts caused by interpolation on a small grid, this algorithm definitely works. However, it is still far from modeling a fluid.

Interactive Example

We can try to solve this in real-time and see how it reacts live. However, the current approach to visualizing will be quite slow when it comes to dynamics.

If we just try to reevaluate ListVectorPlot directly, it will be horribly slow. In general, there is no need to reevaluate the whole expression for graphics. The grid is the same, but vectors differ. We can go either fully raster using Image or use SVG graphics. The latter is the easiest way. Let us start from the ground up.

grid = grid // removeDivergence;

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

(*VB[*)(FrontEndRef["46f8ea9d-0f01-4681-82d9-35d40f2fba91"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKm5ilWaQmWqboGqQZGOqamFkY6loYpVjqGpummBikGaUlJVoaAgCBGxV6"*)(*]VB*)

What about a color you ask. One can calculate Hue based on the polar angle of our vector

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

(*VB[*)(FrontEndRef["5e1eca05-6839-466c-9b1d-9ac89c5ca036"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKm6YapiYnGpjqmlkYW+qamJkl61omGaboWiYmW1gmmwKljM0AhKkVrA=="*)(*]VB*)

Now we basically recreated VectorPlot function in our primitive way. The next step will be to implement dynamic updates.

Dynamic Evaluation

The key feature to fast dynamic evaluation is to minimize the amount of data transferred between the Kernel and frontend as well as the number of things updated.

We need to couple the Arrow primitive to a specific element of our array grid. To do this, one can use the Offload wrapper, which acts similarly to Hold

grid = Table[{0.,0.}, {i,5}, {j,5}];
grid[[3,3]] = {0,1.0}; (*BB[*)(* initial divergence *)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*)

Graphics[{
  Table[With[{i=i, j=j}, (*BB[*)(* to substitute numbers, not symbols *)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*)
    Arrow[{
      {i,j}, 
      {i,j} + 1.5 grid[[i]][[j]]
    } // Offload ]], {i, 5}, {j, 5}]
}, PlotRange->{{0.5,5.5}, {0.5,5.5}}, ImagePadding->None]

Now try to evaluate the next cell multiple times.

grid = grid // removeDivergence;

Let us increase the resolution and automate the process using AnimationFrameListener. In general, one can also add dynamic color to our arrows as well.

TL;DR

bgrid = Table[{0.,0.}, {i,15}, {j,15}];
bcolors = Table[1.0, {Length[bgrid]}, {Length[bgrid]}];
bfps = 0;

mouseHandler[bgrid_, win_, canvas_] = Module[{
  arrow = False,
  dest = {0,0},
  start = {0,0}
}, {
      "mouseup" -> Function[xy,
        With[{v = -Normalize[start - xy]},
          Do[ 
            With[{p = Clip[#, {1,15}] &/@ Round /@ ((xy - start) a + start)},
              bgrid[[p[[1]],p[[2]]]] = {v[[1]], v[[2]]};
            ];
          , {a, 0, 1,0.1}]; (* 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]

(*BB[*)(*define handler for the animation*)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*)
btime = AbsoluteTime[];

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


With[{
  win = CurrentWindow[], (*BB[*)(* save the current window to append graphics *)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*)
  currentCell = ResultCell[],
  canvasId = CreateUUID[]
},

  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[{
    Table[With[{i=i, j=j}, 
      (*BB[*)(* now we have dynamic Hue and dynamic Arrow *)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*)
      Offload[{ 
        Hue[bcolors[[i]][[j]]],
        Arrow[{{i,j}, {i,j} +  bgrid[[i]][[j]]}]
      }] 
    
    ], {i,15}, {j,15}],

    EventHandler[Graphics`Canvas[], mouseHandler[bgrid, win, MetaMarker[canvasId]]], 

    (*BB[*)(*sync with browser's repaint cycle and update of fps label*)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*)
    AnimationFrameListener[bfps // Offload, "Event"->"bframe"], 
    (*BB[*)(*mark this instance of Graphics with uid to append new elements*)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*)
    MetaMarker[canvasId], 
    Text[bfps // Offload, {0,0}]
  }, 
    Controls->False, 
    ImageSize->500, 
    PlotRange->{{-0.5,15.5}, {-0.5,15.5}}, 
    ImagePadding->None, 
    TransitionType->None
  ]
]

The code can be divided into two parts: visualization using Graphics and Arrow and calculation loop.

As we discussed before, we draw each point from a grid individually with a pair of Hue and Arrow. At the same place, we attach listeners to the Graphics canvas to capture mouse clicks and movements to modify the velocity field.

Animation is done using AnimationFrameListener coupled to fpsLabel, which works as a trigger to a new frame. It serves two goals:

Sync animation with the browser's repaint cycle.
Do not start a new frame until the data is ready, which is indicated by fpsLabel. It can also be grid, but since fpsLabel is the last symbol to be updated before a new frame, we have chosen it to be a trigger.

Then, before a new frame, it fires an event object with "bframe" uid, which is captured later by EventHandler. In the animation loop, it removes divergence, recalculates colors, and then also updates the FPS counter just for the statistics.

An expected result is shown on a GIF below

It is quite far from an actual fluid simulation. We should also consider the momentum of it.

See you in Part 2 🦄

Working with a Grid​

Custom Syntax Sugar​

Divergence​

Example​

Removing the Divergence​

Interactive Example​

Dynamic Evaluation​