# 3D stable fluids algorithm to simulate transition from laminar to turbulent flow

This algorithm is 3D extension of our 2D algorithm published on this page and here.
We suppose that with this code we can simulate transition from laminar to turbulent flow. In this example we compute viscous flow around cuboid at Reynolds number $$Re=6250$$. Note, that code has been tested up to $$Re=10^6$$.

dif = 1/6250; pec = .72; U0 = 1.; V0 = 0.; W0 = 0.; dn0 = 1.; n = 80;
n1 = n + 1; sm = 500; r = 20; den = ConstantArray[dn0 , {n1, n1, n1}];
u0 = ConstantArray[U0, {n1, n1, n1}];
v0 = ConstantArray[V0, {n1, n1, n1}]; w0 =
ConstantArray[W0, {n1, n1, n1}]; Do[u0[[i, j, k]] = 0;
v0[[i, j, k]] = 0;
w0[[i, j, k]] = 0;, {i, 10, 20, 1}, {j, Round[n/2 - 5],
Round[n/2 + 5], 1}, {k, Round[n/2 - 5], Round[n/2 + 5], 1}];

periodic[n_, up_, ud_, ul_, ur_, ub_] :=
Module[{bd = ub},
Do[bd[[n + 1, i, j]] = bd[[2, i, j]];
bd[[1, i, j]] = bd[[n, i, j]];, {i, 2, n}, {j, 2, n}];
Do[bd[[i, 1, j]] = ud;
bd[[i, n + 1, j]] = up; bd[[i, j, 1]] = ul;
bd[[i, j, n + 1]] = ur;, {i, 1, n + 1}, {j, 1, n + 1}];
bd];

diffuse[n_, r_, a_, c_, c0_] :=
Module[{c1 = c},
Do[Do[Do[
Do[c1[[i, j,
k]] = (c0[[i, j, k]] +
a (c1[[i - 1, j, k]] + c1[[i + 1, j, k]] +
c1[[i, j - 1, k]] + c1[[i, j + 1, k]] +
c1[[i, j, k - 1]] + c1[[i, j, k + 1]]))/(1 +
6 a);, {k, 2, n}];, {j, 2, n}];, {i, 2, n}];
Do[c1[[i, j, k]] = 0;, {i, 10, 20, 1}, {j, Round[n/2 - 5],
Round[n/2 + 5], 1}, {k, Round[n/2 - 5], Round[n/2 + 5],
1}];, {k1, 0, r}];
c1];

advect[n_, d0_, u_, v_, w_, dt_] :=
Module[{x, y, z, d1, dt0, i0, i1, j0, j1, k0, k1, s0, s1, t0, t1,
p1, p0, d000, d100, d010, d001, d110, d011, d101, d111},
d1 = ConstantArray[0, {n + 1, n + 1, n + 1}]; dt0 = dt n;
Do[Do[Do[x = i - dt0 u[[i, j, k]]; y = j - dt0 v[[i, j, k]];
z = k - dt0 w[[i, j, k]];
i0 = Which[x <= 1, 1, 1 < x < n, Floor[x], True, n];
i1 = i0 + 1;
j0 = Which[y <= 1, 1, 1 < y < n, Floor[y], True, n];
j1 = j0 + 1;
k0 = Which[z <= 1, 1, 1 < z < n, Floor[z], True, n];
k1 = k0 + 1;
d000 = (d0[[i0, j0,
k0]] + (x - i0) (d0[[i1, j0, k0]] -
d0[[i0, j0, k0]]) + (y - j0) (d0[[i0, j1, k0]] -
d0[[i0, j0, k0]]) + (z - k0) (d0[[i0, j0, k1]] -
d0[[i0, j0, k0]]));
d100 = (d0[[i1, j0,
k0]] + (x - i1) (d0[[i1, j0, k0]] -
d0[[i0, j0, k0]]) + (y - j0) (d0[[i1, j1, k0]] -
d0[[i1, j0, k0]]) + (z - k0) (d0[[i1, j0, k1]] -
d0[[i1, j0, k0]]));
d010 = (d0[[i0, j1,
k0]] + (x - i0) (d0[[i1, j1, k0]] -
d0[[i0, j1, k0]]) + (y - j1) (d0[[i0, j1, k0]] -
d0[[i0, j0, k0]]) + (z - k0) (d0[[i0, j1, k1]] -
d0[[i0, j1, k0]]));
d001 = (d0[[i0, j0,
k1]] + (x - i0) (d0[[i1, j0, k1]] -
d0[[i0, j0, k1]]) + (y - j0) (d0[[i0, j1, k1]] -
d0[[i0, j0, k1]]) + (z - k1) (d0[[i0, j0, k1]] -
d0[[i0, j0, k0]]));
d110 = (d0[[i1, j1,
k0]] + (x - i1) (d0[[i1, j1, k0]] -
d0[[i0, j1, k0]]) + (y - j1) (d0[[i1, j1, k0]] -
d0[[i1, j0, k0]]) + (z - k0) (d0[[i1, j1, k1]] -
d0[[i1, j1, k0]]));
d011 = (d0[[i0, j1,
k1]] + (x - i0) (d0[[i1, j1, k1]] -
d0[[i0, j1, k1]]) + (y - j1) (d0[[i0, j1, k1]] -
d0[[i0, j0, k1]]) + (z - k1) (d0[[i0, j1, k1]] -
d0[[i0, j1, k0]]));
d101 = (d0[[i1, j0,
k1]] + (x - i1) (d0[[i1, j0, k1]] -
d0[[i0, j0, k1]]) + (y - j0) (d0[[i1, j1, k1]] -
d0[[i1, j0, k1]]) + (z - k1) (d0[[i1, j0, k1]] -
d0[[i1, j0, k0]]));
d111 = (d0[[i1, j1,
k1]] + (x - i1) (d0[[i1, j1, k1]] -
d0[[i0, j1, k1]]) + (y - j1) (d0[[i1, j1, k1]] -
d0[[i1, j0, k1]]) + (z - k1) (d0[[i1, j1, k1]] -
d0[[i1, j1, k0]]));
d1[[i, j,
k]] = (d000 + d100 + d010 + d001 + d110 + d011 + d101 +
d111)/8;, {k, 2, n}];, {j, 2, n}];, {i, 1, n + 1}]; d1];

project[n_, r_, u0_, v0_, w0_, u_, v_, w_] :=
Module[{ux = u, vy = v, wz = w, div, p},
p = ConstantArray[0, {n + 1, n + 1, n + 1}];
div = ConstantArray[0, {n + 1, n + 1, n + 1}];
ux = ConstantArray[0, {n + 1, n + 1, n + 1}];
vy = ConstantArray[0, {n + 1, n + 1, n + 1}];
wz = ConstantArray[0, {n + 1, n + 1, n + 1}];
Do[div[[i, j,
k]] = .5/
n (u0[[i + 1, j, k]] - u0[[i - 1, j, k]] + v0[[i, 1 + j, k]] -
v0[[i, j - 1, k]] + w0[[i, j, k + 1]] -
w0[[i, j, k - 1]]);, {i, 2, n}, {j, 2, n}, {k, 2, n}];
Do[div[[i, j, k]] = 0;, {i, 10, 20, 1}, {j, Round[n/2 - 5],
Round[n/2 + 5], 1}, {k, Round[n/2 - 5], Round[n/2 + 5], 1}];
div = periodic[n, 0, 0, 0, 0, div];
Do[Do[Do[
Do[p[[i, j,
k]] = (div[[i, j,
k]] + (p[[i - 1, j, k]] + p[[i + 1, j, k]] +
p[[i, j - 1, k]] + p[[i, j + 1, k]] +
p[[i, j, k - 1]] + p[[i, j, k + 1]]))/6;, {k, 2,
n}];, {j, 2, n}], {i, 2, n}];, {k1, 0, r}];
Do[p[[i, j, k]] = 0;, {i, 10, 20, 1}, {j, Round[n/2 - 5],
Round[n/2 + 5], 1}, {k, Round[n/2 - 5], Round[n/2 + 5], 1}];
p = periodic[n, 0, 0, 0, 0, p];
Do[ux[[i, j, k]] =
u0[[i, j, k]] + .5 n (p[[i + 1, j, k]] - p[[i - 1, j, k]]);
vy[[i, j, k]] =
v0[[i, j, k]] + .5 n (p[[i, j + 1, k]] - p[[i, j - 1, k]]);
wz[[i, j, k]] =
w0[[i, j, k]] + .5 n (p[[i, j, k + 1]] - p[[i, j, k - 1]]);, {i,
2, n}, {j, 2, n}, {k, 2, n}];
Do[ux[[i, j, k]] = 0; vy[[i, j, k]] = 0;
wz[[i, j, k]] = 0;, {i, 10, 20, 1}, {j, Round[n/2 - 5],
Round[n/2 + 5], 1}, {k, Round[n/2 - 5], Round[n/2 + 5], 1}]; {ux,
vy, wz}];

onestep[n_, step_, r_, a_, uin_, vin_, win_, dt_, c_] :=
Module[{u1, v1, w1, u, v, w, u0, v0, w0},
u0 = ConstantArray[0., {n + 1, n + 1, n + 1}];
v0 = ConstantArray[0., {n + 1, n + 1, n + 1}];
w0 = ConstantArray[0., {n + 1, n + 1, n + 1}];
u = ConstantArray[0., {n + 1, n + 1, n + 1}];
v = ConstantArray[0., {n + 1, n + 1, n + 1}];
w = ConstantArray[0., {n + 1, n + 1, n + 1}];
u1 = ConstantArray[0., {n + 1, n + 1, n + 1}];
v1 = ConstantArray[0., {n + 1, n + 1, n + 1}];
w1 = ConstantArray[0., {n + 1, n + 1, n + 1}]; u0 = uin; v0 = vin;
w0 = win;
Do[u0[[i, j, k]] = 0; v0[[i, j, k]] = 0;
w0[[i, j, k]] = 0;, {i, 10, 20, 1}, {j, Round[n/2 - 5],
Round[n/2 + 5], 1}, {k, Round[n/2 - 5], Round[n/2 + 5], 1}];
u0 = advect[n, u0, u0, v0, w0, dt];
v0 = advect[n, v0, u0, v0, w0, dt];
w0 = advect[n, w0, u0, v0, w0, dt];
Do[u0[[i, j, k]] = 0;
v0[[i, j, k]] = 0;, {i, 10, 20, 1}, {j, Round[n/2 - 5],
Round[n/2 + 5], 1}, {k, Round[n/2 - 5], Round[n/2 + 5], 1}];
u0 = periodic[n, 0, 0, 0, 0, u0]; v0 = periodic[n, 0, 0, 0, 0, v0];
w0 = periodic[n, 0, 0, 0, 0, w0];
u0 = diffuse[n, r, a, c, u0]; v0 = diffuse[n, r, a, c, v0];
w0 = diffuse[n, r, a, c, w0];
Do[u0[[i, j, k]] = 0; v0[[i, j, k]] = 0;
w0[[i, j, k]] = 0;, {i, 10, 20, 1}, {j, Round[n/2 - 5],
Round[n/2 + 5], 1}, {k, Round[n/2 - 5], Round[n/2 + 5], 1}];
u0 = periodic[n, 0, 0, 0, 0, u0]; v0 = periodic[n, 0, 0, 0, 0, v0];
w0 = periodic[n, 0, 0, 0, 0, w0];
{u0, v0, w0} = project[n, r, u0, v0, w0, u, v, w];
Do[u0[[i, j, k]] = 0; v0[[i, j, k]] = 0;
w0[[i, j, k]] = 0;, {i, 10, 20, 1}, {j, Round[n/2 - 5],
Round[n/2 + 5], 1}, {k, Round[n/2 - 5], Round[n/2 + 5], 1}];
u0 = periodic[n, 0, 0, 0, 0, u0]; v0 = periodic[n, 0, 0, 0, 0, v0];
w0 = periodic[n, 0, 0, 0, 0, w0]; {u0, v0, w0}]

cf = With[{cg = CompileGetElement, hp = HoldPattern,
dv = DownValues},
Hold@Compile[{{u0argu, _Real, 3}, {v0argu, _Real,
3}, {w0argu, _Real, 3}, {denargu, _Real,
3}, {sm, _Integer}, {n, _Integer}, {r, _Integer}, dif,
pec}, Module[{u0 = u0argu, v0 = v0argu, w0 = w0argu, uu,
vv, ww, dd, den = denargu,
c = Table[0., {n + 1}, {n + 1}, {n + 1}], dt = 40./n^2,
a, dnup = den[[1, n + 1, 1]], dnd = den[[1, 1, 1]],
dnl = den[[1, 1, 1]], dnr = den[[1, 1, n + 1]]},
a = dt dif n n;

uu = vv =
ww = dd =
Table[0., {sm + 1}, {n + 1}, {n + 1}, {n + 1}];

Do[{u0, v0, w0} =
onestep[n, step, r, a, u0, v0, w0, dt, c];
uu[[step + 1]] = u0;
vv[[step + 1]] = v0; ww[[step + 1]] = w0;
den = diffuse[n, r, a/pec, c, den];
Do[den[[i, j, k]] = 0;, {i, 10, 20, 1}, {j,

Round[n/2 - 5], Round[n/2 + 5], 1}, {k,
Round[n/2 - 5], Round[n/2 + 5], 1}];
den = periodic[n, dnup, dnd, dnl, dnr, den];
den = advect[n, den, u0, v0, w0, dt];
Do[den[[i, j, k]] = 0;, {i, 10, 20, 1}, {j,
Round[n/2 - 5], Round[n/2 + 5], 1}, {k,
Round[n/2 - 5], Round[n/2 + 5], 1}];
den = periodic[n, dnup, dnd, dnl, dnr, den];
dd[[step + 1]] = den;, {step, 0, sm}]; {uu, vv, ww,
dd}], CompilationTarget -> "C",
RuntimeOptions -> "Speed"] /. dv@onestep /.
Flatten[dv /@ {advect, diffuse, periodic, project}] /.
hp@ConstantArray[c_, {i_, j_, kc_}] :>
Table[0., {i}, {j}, {kc}] /. hp@Part[a__] :> cg[a] /.
hp[cg[a__] = rhs_] :> (Part[a] = rhs) //
ReleaseHold];

rst = cf[u0, v0, w0, den, sm, n, r, dif, pec];


Visualization

Do[lst11[s] =
Table[{{(i - 1)/n, (j - 1)/n, (k - 1)/n}, rst[[1, s, i, j, k]]}, {i,
n1}, {j, n1}, {k, n1}];
lst12[s] =
Table[{{(i - 1)/n, (j - 1)/n, (k - 1)/n}, rst[[2, s, i, j, k]]}, {i,
n1}, {j, n1}, {k, n1}];
lst13[s] =
Table[{{(i - 1)/n, (j - 1)/n, (k - 1)/n}, rst[[3, s, i, j, k]]}, {i,
n1}, {j, n1}, {k, n1}];, {s, 25, sm, 25}]

Do[su1[i] =
Interpolation[Flatten[lst11[i], 2], InterpolationOrder -> 3];
sv2[i] = Interpolation[Flatten[lst12[i], 2],
InterpolationOrder -> 3];
sw3[i] = Interpolation[Flatten[lst13[i], 2],
InterpolationOrder -> 3];, {i, 25, sm, 25}]

Table[Show[
DensityPlot3D[
Norm[{su1[i][x, y, z], sv2[i][x, y, z], sw3[i][x, y, z]}], {x, 0,
1}, {y, 0.44, 1}, {z, 0, 1}, BoxRatios -> Automatic,
ColorFunction -> "BlueGreenYellow", PlotLabel -> i],
VectorPlot3D[{su1[i][x, y, z], sv2[i][x, y, z],
sw3[i][x, y, z]}, {x, 0, 1}, {y, 0, .5}, {z, 0, 1},
VectorPoints -> Fine, VectorMarkers -> "Arrow"],
Graphics3D[{{Blue,
Cuboid[{10, Round[n/2 - 5] - 1/2,
Round[n/2 - 5] - 1/2}/(n + 1), {20, Round[n/2 + 5] - 1/2,
Round[n/2 + 5] - 1/2}/(n + 1)]}}]], {i, 50, sm, 50}]


Animation

Do[lstu1[k] =
Flatten[Table[{{(i - 1)/n, (j - 1)/n},
rst[[1, k, i, Round[(n + 1)/2], j]]}, {i, n1}, {j, n1}], 1];
lstw1[k] =
Flatten[Table[{{(i - 1)/n, (j - 1)/n},
rst[[2, k, i, Round[(n + 1)/2], j]]}, {i, n1}, {j, n1}],
1];, {k, sm}];
Do[Uvel1[i] = Interpolation[lstu1[i], InterpolationOrder -> 3];, {i,
1, sm}];
Do[Wvel1[i] = Interpolation[lstw1[i], InterpolationOrder -> 3];, {i,
1, sm}];
frame = Table[
Show[DensityPlot[
Norm[{Uvel[m][x, y], Vvel[m][x, y]}], {x, 0, 1}, {y, 0, 1},
PlotRange -> All, ColorFunction -> "BlueGreenYellow",
Frame -> False, ImageSize -> Tiny, PlotLabel -> m,
PlotPoints -> 50],
Graphics[{Gray,
Rectangle[{10, Round[n/2 - 5] - 1/2}/(n + 1), {20,
Round[n/2 + 5] - 1/2}/(n + 1)]}]], {m, 5, sm, 3}];


The question is about code improvement. How can we define parameter r to solve Laplace and Poison equations in this code? Note, r is number of iterations in diffuse and project module (we use Gauss-Seidel relaxation to solve Laplace and Poison's equations).

Update 1. We can reduce r from r=20 as above to r=5 as below. Globally there is no difference in two animations.

But if we compare velocity in some point like $$x=0.65, y=0.5$$ then we have good agreement for $$r=5$$ (red points) and $$r=20$$ (gray points) in laminar flow at $$t<0.5$$, but a big difference in the turbulent flow at $$t>0.5$$.

Update 1. To test the Gauss-Seidel relaxation method itself we can compute viscous flow in a plane channel with external force (gravity) as follows

ClearAll["Global*"]

dif = 1/80; pec = 1; U0 = 0; V0 = 0; n = 11; n1 = n + 1; dt =
2./n^2; sm = 3000; r = 7; n2 = Round[n/2]; a = dt dif n n; c =
ConstantArray[0, {n1, n1}]; c0 = ConstantArray[0, {n1, n1}];; u0 =
ConstantArray[0, {n1, n1}]; v0 = ConstantArray[0, {n1, n1}]; u =
ConstantArray[0, {n1, n1}]; v = ConstantArray[0, {n1, n1}];

periodic[n_, up_, ud_, ub_] :=
Module[{bd = ub},
Do[bd[[1, i]] = .5 (bd[[n, i]] + bd[[2, i]]);
bd[[n + 1, i]] = bd[[1, i]]; bd[[i, 1]] = ud;
bd[[i, n + 1]] = up;, {i, 2, n}]; bd[[1, 1]] = ud;
bd[[n + 1, n + 1]] = up; bd[[n + 1, 1]] = ud; bd[[1, n + 1]] = up;
bd];

diffuse[n_, r_, a_, c_, c0_] :=
Module[{c1 = c}, c1 = ConstantArray[0, {n + 1, n + 1}];
Do[Do[Do[c1[[i,
j]] = (c0[[i, j]] +
a (c1[[i - 1, j]] + c1[[i + 1, j]] + c1[[i, j - 1]] +
c1[[i, j + 1]]))/(1 + 4 a);, {j, 2, n}];, {i, 2, n}];
c1 = periodic[n, 0, 0, c1];, {k, 0, r}]; c1];
advect[n_, b_, d_, d0_, u_, v_, dt_] :=
Module[{x, y}, d1 = ConstantArray[0, {n + 1, n + 1}]; dt0 = dt n;
Do[Do[x = i - dt0 u[[i, j]]; y = j - dt0 v[[i, j]];
i0 = Which[x <= 1, 1, 1 < x < n, Floor[x], x >= n, n];
i1 = i0 + 1;
j0 = Which[y <= 1, 1, 1 < y < n, Floor[y], y >= n, n];
j1 = j0 + 1; s1 = x - i0; s0 = 1 - s1; t1 = y - j0; t0 = 1 - t1;
d1[[i, j]] =
s0 (t0 d0[[i0, j0]] + t1 d0[[i0, j1]]) +
s1 (t0 d0[[i1, j0]] + t1 d0[[i1, j1]]);, {j, 1, n + 1}];, {i,
1, n + 1}]; d1];
project[n_, r_, u0_, v0_, u_, v_] :=
Module[{ux = u, vy = v}, p = ConstantArray[0, {n1, n1}];
div = ConstantArray[0, {n1, n1}]; ux = ConstantArray[0, {n1, n1}];
vy = ConstantArray[0, {n1, n1}];
Do[div[[i,
j]] = -.5 /
n (u0[[i + 1, j]] - u0[[i - 1, j]] + v0[[i, 1 + j]] -
v0[[i, j - 1]]);, {i, 2, n}, {j, 2, n}];
Do[Do[Do[
p[[i, j]] = (div[[i,
j]] + (p[[i - 1, j]] + p[[i + 1, j]] + p[[i, j - 1]] +
p[[i, j + 1]]))/4;, {j, 2, n}], {i, 2, n}];, {k, 0, r}];
Do[ux[[i, j]] = u0[[i, j]] - .5 n (p[[i + 1, j]] - p[[i - 1, j]]);
vy[[i, j]] =
v0[[i, j]] - .5 n (p[[i, j + 1]] - p[[i, j - 1]]);, {i, 2,
n}, {j, 2, n}]; {ux, vy}];
(*force*)
Fx[t_, x_, y_] := 1/10; Fy[t_, x_, y_] := 0;
f1 = ConstantArray[0, {n1, n1}]; f2 = ConstantArray[0, {n1, n1}];

Do[u0 = advect[n, 0, c, u0, u0, v0, dt];
v0 = advect[n, 0, c, v0, u0, v0, dt];
u0 = periodic[n, 0, 0, u0]; v0 = periodic[n, 0, 0, v0];
u0 = diffuse[n, r, a, c, u0]; v0 = diffuse[n, r, a, c, v0];
u0 = periodic[n, 0, 0, u0]; v0 = periodic[n, 0, 0, v0];
Do[f1[[i, j]] = Fx[dt (step + .5), (i - 1)/n, (j - 1)/n];
f2[[i, j]] = Fy[dt (step + .5), (i - 1)/n, (j - 1)/n];, {i, 2,
n1}, {j, 2, n1}]; u0 += f1 dt;
v0 += f2 dt; {u0, v0} = project[n, r, u0, v0, u, v];
u0 = periodic[n, 0, 0, u0]; v0 = periodic[n, 0, 0, v0];
uu[step] = u0; vv[step] = v0;, {step, 0, sm}];


Visualization of flow velocity and absolute error

Do[lstu[k] =
Flatten[Table[{{(i - 1)/n, (j - 1)/n}, uu[k][[i, j]]}, {i, n1}, {j,
n1}], 1];
lstv[k] =
Flatten[Table[{{(i - 1)/n, (j - 1)/n}, vv[k][[i, j]]}, {i, n1}, {j,
n1}], 1];, {k, 0, sm}];
Do[Uvel[i] = Interpolation[lstu[i], InterpolationOrder -> 3];, {i, 1,
sm}]
Do[Vvel[i] = Interpolation[lstv[i], InterpolationOrder -> 3];, {i, 1,
sm}];

{StreamDensityPlot[{Uvel[sm][x, y], Vvel[sm][x, y]}, {x, 0, 1}, {y, 0,
1}, ColorFunction -> "RoseColors", Frame -> False,
PlotLegends -> Automatic],
Plot3D[Uvel[sm][x, y], {x, 0, 1}, {y, 0, 1},
ColorFunction -> "RoseColors"]}

{ListPlot[Table[{i dt, Max[uu[i]]}, {i, sm}],
AxesLabel -> {"t", "Umax"}],
Plot[{-4 (-y + y^2), Uvel[sm][.5, y]}, {y, 0, 1},
PlotStyle -> {Thick, {Red, Dashed}}, Frame -> True,
FrameLabel -> {"y", "u"},
PlotLegends -> {"Exact solution", "Numeric solution"}],
Plot[Abs[-Uvel[sm][.5, y] - 4 (-y + y^2)], {y, 0, 1}, Frame -> True,
FrameLabel -> {"y", "Error"}]}


Actually for r=7 we have a good convergence to the exact solution on the grid of 12 points with error of $$10^{-2}$$. For r=5 the error is about 0.2 and for r=10 error is about $$10^{-2}$$, therefore r=7 is optimal value.

Update 2. In the module advect we can use trilinear interpolation as follows

advect[n_, d0_, u_, v_, w_, dt_] :=
Module[{x, y, z, d1, dt0, i0, i1, j0, j1, k0, k1, s0, s1, t0, t1,
p1, p0, d00, d10, d01, d11, cd0, cd1, xd, yd, zd},
d1 = ConstantArray[0, {n + 1, n + 1, n + 1}]; dt0 = dt n;
Do[Do[Do[x = i - dt0 u[[i, j, k]]; y = j - dt0 v[[i, j, k]];
z = k - dt0 w[[i, j, k]];
i0 = Which[x <= 1, 1, 1 < x < n, Floor[x], True, n];
i1 = i0 + 1;
j0 = Which[y <= 1, 1, 1 < y < n, Floor[y], True, n];
j1 = j0 + 1;
k0 = Which[z <= 1, 1, 1 < z < n, Floor[z], True, n];
k1 = k0 + 1;(*Trilinear interpolation*)xd = x - i0;
yd = y - j0; zd = z - k0;
d00 = d0[[i0, j0, k0]] (1 - xd) + d0[[i1, j0, k0]] xd;
d01 = d0[[i0, j0, k1]] (1 - xd) + d0[[i1, j0, k1]] xd;
d10 = d0[[i0, j1, k0]] (1 - xd) + d0[[i1, j1, k0]] xd;
d11 = d0[[i0, j1, k1]] (1 - xd) + d0[[i1, j1, k1]] xd;
cd0 = d00 (1 - yd) + d10 yd; cd1 = d01 (1 - yd) + d11 yd;
d1[[i, j, k]] = cd0 (1 - zd) + cd1 zd;
, {k, 2, n}];, {j, 2, n}];, {i, 1, n + 1}]; d1];

• What is the meaning of r? Reynolds number? I think it is not a good idea to name a parameter r because it is nearly impossible to find it in the code. How do you want to improve your code: in readability (sure you can), speed, memory use? Your code looks like Fortran in Mathematica. Why? You may gain a lot by using Mma expressive power. Commented Dec 25, 2021 at 15:06
• @PierreALBARÈDE r is number of iterations in diffuse and project module (we use Gauss-Seidel relaxation to solve Laplace and Poison's equations). Actually this code optimized for C and basically it coming from C code published in this paper dgp.toronto.edu/public_user/stam/reality/Research/pdf/ns.pdf Commented Dec 25, 2021 at 15:15
• Don't you have some convergence criterion on this r? How did you choose r=20? Commented Dec 25, 2021 at 15:28
• Er… so you mean you're looking for a systematic way to determine the proper value of r? Commented Dec 25, 2021 at 15:29
• @PierreALBARÈDE For real time simulations in NVIDIA applications they usually use $r \ge 10$ for 2D and $r \ge 4$ for 3D flow. My choose r=20 is based on several tests and comparisons with other solvers. But, I think, that r=20 is too high for this kind of problems. Commented Dec 25, 2021 at 15:45

To test projection step itself we can use Mathematica FEM and exact benchmark solution from the paper EXACT FULLY 3D NAVIER-STOKES SOLUTIONS FOR BENCHMARKING by C. ROSS ETHIER AND D. A. STEINMAN as follows. Note, that we can consider projection step as an implementation of the predictor-corrector algorithm $$\frac{u-u_n}{\tau}+(u.\nabla)u-\nu\nabla^2 u=0$$ $$\frac{u_{n+1}-u}{\tau}+\nabla p =0$$ here $$\tau$$ is time step, $$u_n, u, u_{n+1}$$ if velocity field on previous, intermediate and next step consequently, and $$p$$ is a pressure. We suppose that $$\nabla.u_{n+1}=0$$, and therefore $$\nabla^2p-\frac{\nabla.u}{\tau}=0$$ We solve these equations in the unit cuboid with Dirichlet condition using exact solution and FEM as follows

Clear["Global*"]
Needs["NDSolveFEM"]
reg = Cuboid[];
mesh = ToElementMesh[reg, MaxCellMeasure -> 0.001];
(*Exact solution*)
U[x_, y_, z_,
t_] := -a Exp[-d^2 t] (Exp[a x] Sin[a y + d z] +
Exp[a z] Cos[a x + d y]);
V[x_, y_, z_,
t_] := -a Exp[-d^2 t] (Exp[a y] Sin[a z + d x] +
Exp[a x] Cos[a y + d z]);
W[x_, y_, z_,
t_] := -a Exp[-d^2 t] (Exp[a z] Sin[a x + d y] +
Exp[a y] Cos[a z + d x]);
P[x_, y_, z_,
t_] := -a ^2/
2 Exp[-2 d^2 t] (Exp[2 a x] + Exp[2 a y] + Exp[2 a z] +
2 Sin[a x + d y] Exp[a (y + z)] Cos[a z + d x] +
2 Sin[a y + d z] Exp[a (x + z)] Cos[a x + d y] +
2 Sin[a z + d x] Exp[a (y + x)] Cos[
a y + d z]);
(*t0 is time step, nn is number of iterations *)
a = 1; d = 1; t0 = 1/400; nn = 200; \[Nu] = 1;
(*FEM implementation of predictor-corrector algorithm*)
UX[0] = U[x, y, z, 0];
VY[0] = V[x, y, z, 0]; WZ[0] = W[x, y, z, 0];
P0[0] = P[x, y, z, 0];
Do[
{UX[i], VY[i], WZ[i], P0[i]} =
NDSolveValue[{{-\[Nu]*
Laplacian[
u[x, y, z], {x, y, z}] + (u[x, y, z] - UX[i - 1])/t0 +
UX[i - 1]*D[u[x, y, z], x] + VY[i - 1]*D[u[x, y, z], y] +
WZ[i - 1]*D[u[x, y, z], z], -\[Nu]*
Laplacian[
v[x, y, z], {x, y, z}] + (v[x, y, z] - VY[i - 1])/t0 +
UX[i - 1]*D[v[x, y, z], x] + VY[i - 1]*D[v[x, y, z], y] +
WZ[i - 1]*D[v[x, y, z], z], -\[Nu]*
Laplacian[
w[x, y, z], {x, y, z}] + (w[x, y, z] - WZ[i - 1])/t0 +
UX[i - 1]*D[w[x, y, z], x] + VY[i - 1]*D[w[x, y, z], y] +
WZ[i - 1]*D[w[x, y, z], z],
Laplacian[p[x, y, z], {x, y, z}] - (
\!$$\*SuperscriptBox[\(u$$,
TagBox[
RowBox[{"(",
RowBox[{"1", ",", "0", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y, z] +
\!$$\*SuperscriptBox[\(v$$,
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y, z] +
\!$$\*SuperscriptBox[\(w$$,
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y, z])/t0} == {0, 0, 0, 0}, {
DirichletCondition[{u[x, y, z] ==
U[x, y, z, i t0 - t0/2] + t0 D[P[x, y, z, i t0], x],
v[x, y, z] ==
V[x, y, z, i t0 - t0/2] + t0 D[P[x, y, z, i t0], y],
w[x, y, z] ==
W[x, y, z, i t0 - t0/2] + t0 D[P[x, y, z, i t0], z],
p[x, y, z] == P[x, y, z, i t0 ]}, True]}}, {u[x, y, z] - t0
\!$$\*SuperscriptBox[\(p$$,
TagBox[
RowBox[{"(",
RowBox[{"1", ",", "0", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y, z], v[x, y, z] - t0
\!$$\*SuperscriptBox[\(p$$,
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y, z], w[x, y, z] - t0
\!$$\*SuperscriptBox[\(p$$,
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y, z],
p[x, y, z]}, {x, y, z} \[Element] mesh,
Method -> {"FiniteElement",
"InterpolationOrder" -> {u -> 2, v -> 2, w -> 2,
p -> 2}}];, {i, 1, nn}];


Visualization of error on every 5 steps for $$u_n=(UX[n],VY[n],WZ[n])$$, and $$p=P0[n]$$ on the line $$x=1/2, z=1/2$$

Table[{k t0 // N,
Plot[Evaluate[{U[x, y, z, k t0]/UX[k] - 1} /. {x -> 1/2,
z -> 1/2}], {y, 0, 1}, PlotRange -> All],
Plot[Evaluate[{V[x, y, z, k t0]/VY[k] - 1} /. {x -> 1/2,
z -> 1/2}], {y, 0, 1}, PlotRange -> All],
Plot[Evaluate[{W[x, y, z, k t0]/WZ[k] - 1} /. {x -> 1/2,
z -> 1/2}], {y, 0, 1}, PlotRange -> All],
Plot[Evaluate[{P[x, y, z, k t0]/P0[k] - 1} /. {x -> 1/2,
z -> 1/2}], {y, 0, 1}, PlotRange -> All]}, {k, 5, nn, 5}]
`

Please note, that in the picture shown last steps only. The relative error for pressure is about $$3\times 10^{-3}$$, and for velocity $$1.25\times10^{-3}$$ on the grid $$10\times10\times 10$$ on 200 steps in time with time step $$\tau =1/400$$.

We also can tested our FDM implementation with using exact solution.