# What's wrong with this FFT-based Von Kármán vortex street simulation?

About 9 years ago, I came across this interesting website, and found the following paragraph with a broken Mathematica code sample:

When fluid passes an object, it can leave a trail of vortices called a Von Kármán Vortex Street. This animation shows the vorticity where blue is clockwise and red is counter-clockwise. This simulation assumes unsteady, incompressible, viscid, laminar flow. For solenoidal flows, mass conservation can be achieved by taking the Fast Fourier Transform (FFT) of the velocity and then removing the radial component of the wave number vectors. The code was adapted from Jos Stam’s C Source Code and paper which is patented by Alias.

(* Note: Something is wrong with this Mathematica code.
Please tell me if you find out what it is. runtime: 14 seconds *)

n = 64; dt = 0.3; mu = 0.001; v = Table[{0, 0}, {n}, {n}];

Do[
Do[If[i < 5, v[[i, j]] = {0.1, 0}];
If[(i - n/4)^2 + (j -n/2)^2 < 4^2, v[[i, j]] = {0, 0}],
{i, 1, n}, {j, 1, n}];

ui = ListInterpolation[v[[All, All, 1]]];
vi = ListInterpolation[v[[All, All, 2]]];
v = Table[{i2, j2} = {i, j} - n dt v[[i, j]];
{ui[i2, j2], vi[i2, j2]}, {i, 1, n}, {j, 1, n}];

v = Transpose[Map[Fourier[v[[All, All, #]]] &, {1, 2}], {3, 1, 2}];
v = Table[x = Mod[i + n/2, n] - n/2;
y = Mod[j + n/2, n] - n/2;
k = x^2 + y^2;
Exp[-k dt mu] If[k > 0, (v[[i, j]].{-y, x}/k){-y, x}, v[[i, j]]],
{i, 1, n}, {j, 1, n}];
v = Transpose[Map[Re[InverseFourier[v[[All, All, #]]]] &,{1,2}], {3, 1, 2}];

ListDensityPlot[Table[((v[[i + 1, j, 2]] - v[[i - 1, j, 2]]) -
(v[[i, j + 1, 1]] - v[[i, j - 1, 1]]))/2,
{j, 2, n - 1}, {i, 2, n - 1}],
Mesh -> False,Frame->False,ColorFunction -> (Hue[2#/3] &)],
{t, 1, 25}];


Since then, every a year or two, I would be seized by a whim and try to fix the broken code, and fail then. And yes, I just failed one more time :D . So I think it may be the time to cry out loud: Why doesn't this sample work? How can it be fixed? By "doesn't work", I mean the obtained simulation result isn't even close to the GIF shown above. As an example, the following is the ListDensityPlot for t == 25:

I know there exist multiple ways to simulate vortex street, but in this question I'm particularly interested in fixing this code sample.

Several possible issues I can spot:

1. The ListDensityPlot is inside a Do loop, this is probably because the code is written in a early version of Mathematica, and *Plot is automatically Printed at that time, as discussed in this post. This isn't a big problem, we just need to e.g. use Table instead of the outermost Do.

2. ui and vi should be periodic, this can be easily fixed with e.g.

{ui, vi} = Module[{ui = #}, ui[[-1]] = ui[[1]]; ui[[All, -1]] = ui[[All, 1]];
ListInterpolation[ui, InterpolationOrder -> 1,
PeriodicInterpolation -> True]] & /@ Transpose[v, {2, 3, 1}]


but fixing this won't fix the simulation.

3. The original C code normalizes the velocity by dividing it by n^2 as the last step, but adding this step to Mathematica code doesn't fix the code. this is not needed for Mathematica code, because according to the document of fftw, the convention of fftw is amount to FourierParameters -> {-1, 1}. The default setting of Fourier is FourierParameters -> {0, 1} i.e. the normalization is already done.

4. The grid used for simulation seems to be too coarse, but adjusting various parameters in the code doesn't seem to help. (I admit my adjustion is somewhat blind, though. )

For your convenience, the following is the code sample aiming at fixing the mentioned issues above. Of course, it still doesn't work:

n = 128; dt = 0.3; mu = 0.001; nt = 50;
v = Table[{0., 0.}, {n}, {n}];
separate = Transpose[#, {2, 3, 1}] &;
combine = Transpose[#, {3, 1, 2}] &;

vlst = Table[Do[If[i < 1 + n/16, v[[i, j]] = {0.1, 0.}];
If[(i - n/4)^2 + (j - n/2)^2 < (n/16)^2,
v[[i, j]] = {0., 0.}], {i, n}, {j, n}];
{ui, vi} =
Module[{umat = #}, umat[[-1]] = umat[[1]];
umat[[All, -1]] = umat[[All, 1]];
ListInterpolation[umat, InterpolationOrder -> 1,
PeriodicInterpolation -> True]] & /@ separate@v;
v = Table[{i2, j2} = {i, j} - n dt v[[i, j]];
{ui[i2, j2], vi[i2, j2]}, {i, n}, {j, n}];
v = Fourier /@ separate@v // combine;
v = Table[x = Mod[i + n/2, n] - n/2;
y = Mod[j + n/2, n] - n/2;
k = x^2 + y^2;
Exp[-k dt mu] If[k > 0, (v[[i, j]] . {-y, x}/k) {-y, x},
v[[i, j]]],
{i, n}, {j, n}];
v = InverseFourier /@ separate@v // combine // Re,
{t, nt}]; // AbsoluteTiming
(* {24.4562, Null} *)
vortex = Compile[{{v, _Real, 3}},
With[{n = Length@v},
Table[((v[[i + 1, j, 2]] - v[[i - 1, j, 2]]) -
(v[[i, j + 1, 1]] - v[[i, j - 1, 1]]))/2,
{j, 2, n - 1}, {i, 2, n - 1}]]];
arrayplot =
ArrayPlot[#, DataReversed -> True, ColorFunction -> "Rainbow"] &;
vortex@vlst[[-1]] // arrayplot


## Update

This is a bit embrassing: there exists another simple mistake in the implementation, if it's fixed, the vortex street will be obtained with minimal parameter adjustion. In order not to waste the bounty, I'd like not to make the answer public for the moment. The bounty will be awarded to the first answerer figuring out the mistake and obtaining the vortex street (not necessarily exactly the same as the GIF above).

## Update 2: Hint

Since the bounty will be expired in 23 hours but nobody has figured out the answer so far, let me give a hint: moving the cursor with mouse, the code can be fixed within 5 keystrokes.

• Maybe it would help if you make your post self-contained. For instance, what equations you are trying to implement? What exactly is the problem (apart from it wouldn't work)? Notice that attached paper is not very informative, not even a single equation there. It is hard to fix something that is not even written out. Oct 16, 2022 at 8:08
• How you tried to run the c++ code to see if it is producing the right picture? Oct 16, 2022 at 8:12
• @yarchik According to my understanding, this algorithm is related to NS equation for incompressible flow (together with the mass continuity equation), but it doesn't solve the PDE in a rigorous manner, for example, it simulates the viscosity of flow by multiplying a Exp[-k dt mu] term to the velocity in frequency domain. The problem is, in principle we should get a simulation result similar to the GIF in the question, but currently the density plot of the calculated vorticity is vastly different from it. Oct 16, 2022 at 8:21
• I doubt that the code matches the visualization. The paper can handle only periodic BCs. Are you familiar with Fourier Pseudospectral methods? Oct 16, 2022 at 8:48
• @yarchik I'd argue the information is enough. The algorithm is published and in public. The question is why the implementation doesn't give a reasonable result? I admit the question is far beyond ideal. (That's why I didn't post it in past 9 years! ) But it's already an answerable question, it's just hard. Oct 16, 2022 at 9:27

I found the index should start from 0 instead of 1. This affects the (squared) wave number k.

x = Mod[i + n/2, n] - n/2;
y = Mod[j + n/2, n] - n/2;


Minus one and modify parameters, I got this

n = 64; dt = 0.3; mu = 0.0001; nt = 300;

vortex = Compile[{{v, _Real, 3}},
With[{n = Length@v},
Table[((v[[i + 1, j, 2]] - v[[i - 1, j, 2]]) -
(v[[i, j + 1, 1]] - v[[i, j - 1, 1]]))/2
, {j, 2, n - 1}, {i, 2, n - 1}]]];
arrayplot = ArrayPlot[#, DataReversed -> True, ColorFunction -> "Rainbow"] &;

vlst = Table[
Do[
Which[
i <= 1 + n/16,
v[[i, j]] = {0.1, 0.},
(i - n/4)^2 + (j - n/2)^2 <= (n/16)^2,
v[[i, j]] = {0., 0.}
]
, {i, n}, {j, n}];
{ui, vi} = Module[{umat = #},
umat[[-1]] = umat[[1]];
umat[[All, -1]] = umat[[All, 1]];
ListInterpolation[umat, InterpolationOrder -> 1, PeriodicInterpolation -> True]
] & /@ separate@v;
v = Table[
{i2, j2} = {i, j} - n dt v[[i, j]];
{ui[i2, j2], vi[i2, j2]}
, {i, n}, {j, n}];
v = Fourier /@ separate@v // combine;
v = Table[
x = Mod[i + n/2-1, n] - n/2;
y = Mod[j + n/2-1, n] - n/2;
k = x^2 + y^2;
Exp[-k dt mu] If[k > 0,
(v[[i, j]] . {-y, x}/k) {-y, x},
v[[i, j]]
]
, {i, n}, {j, n}];
v = InverseFourier /@ separate@v // combine //Re
, {t, nt}]; // AbsoluteTiming


Is that all?

• Yes, this is the correct answer. Will award the bounty tomorrow morning. (You can set nt to 350 so we're within 5 keystrokes :) . ) Oct 24, 2022 at 15:09
• It is a nice solution (+1). Oct 25, 2022 at 5:27

Finally, rnotlnglgq finds the correct answer in time, so let me share my answer, too :) .

First I'd like to explain a bit more about the theoretical basis of the mass conservation correction in the algorithm, given the original paper is a bit brief. In 2D Cartesian coordinates, the continuity equation of incompressible flow is $$\newcommand{\f}{\frac}$$ $$\newcommand{\p}{\partial}$$ $$\f{\p u(x,y)}{\p x}+\f{\p v(x,y)}{\p y}=0$$

Now, if the $$(u(x,y),v(x,y))$$ field is well-behaved, and is periodic or defined on an infinite domain, then we can apply Fourier transform or finite Fourier transform to the equation and obtain

$$\omega _x \mathcal{F}_{x,y}[u(x,y)]\left(\omega _x,\omega _y\right) +\omega _y \mathcal{F}_{x,y}[v(x,y)]\left(\omega _x,\omega _y\right)=0$$

where $$\mathcal{F}_{x,y}[u(x,y)]\left(\omega _x,\omega _y\right)$$ stands for 2D Fourier transform/finite Fourier transform of $$u(x,y)$$. Clearly it's equivalent to

$$(\omega _x,\omega _y)\cdot (\mathcal{F}_{x,y}[u(x,y)]\left(\omega _x,\omega _y\right),\mathcal{F}_{x,y}[v(x,y)]\left(\omega _x,\omega _y\right))=0$$

i.e. the projection of 2D Fourier transform/finite Fourier transform of $$(u,v)$$ in $$(\omega _x,\omega _y)$$ direction should be zero, otherwise the continuity equation won't hold.

In the algorithm Fourier has been used for numerically calculate Fourier transform/finite Fourier transform, but the output of Fourier is arranged in a somewhat special way which is discussed in detail in e.g. this post, so we need to shift the frequency. This is incorrectly done in original code as

…
x = Mod[i + n/2, n] - n/2;
y = Mod[j + n/2, n] - n/2;
…


The correct one, as shown by rnotlnglgq, should be

…
x = Mod[i + n/2-1, n] - n/2;
y = Mod[j + n/2-1, n] - n/2;
…


or equivalently:

…
x = Mod[i - 1, nx, -nx/2];
y = Mod[j - 1, ny, -ny/2];
…


One may find this suspicious at first glance. "Yeah, the original code is incorrect, but is it that bad? It's almost correct, and I'd expect it producing an almost correct result." This is the surprising part, the algorithm is very sensitive to the correctness of frequency shift. Such a small mistake ruins everything. (BTW, avoiding the singularity at k == 0 with e.g. k = x^2 + y^2 + 10^-16. ruins the algorithm, too. )

Now that the algorithm is fixed, let me show an optimized implementation:

SetAttributes[compile, HoldAll];
compile[argu__] :=
With[{cg = CompileGetElement},
Hold@Compile[argu, RuntimeOptions -> "Speed", CompilationTarget -> "C"] /.
Part -> cg /. HoldPattern[cg[a__] = rhs_] :> (Part[a] = rhs)] // ReleaseHold //
Last

vortex = compile[{{v, _Real, 3}}, Module[{nx, ny}, {nx, ny} = Dimensions@v[[1]];
Table[((v[[2, i + 1, j]] - v[[2, i - 1, j]]) - (v[[1, i, j + 1]] -
v[[1, i, j - 1]]))/2, {j, 2, ny - 1}, {i, 2, nx - 1}]]];

interrule = inter[valueL_, valueR_, scale_] :> (1 - scale) valueL + scale valueR;

Hold@compile[{{arg, _Real, 3}, {dt, _Real}},
Module[{u, v, unew, vnew, nx, ny, inew, jnew, iL, jL, iR, jR},
{u, v} = arg;
{nx, ny} = Dimensions@v;
unew = vnew = Table[0., {nx}, {ny}];
(*Alternative force:*)
(*Do[v[[1,j]]={.1,0.},{j,ny}];*)
Do[If[i < 1 + nx/16, u[[i, j]] = 0.1; v[[i, j]] = 0.];
If[(i - nx/4)^2 + (j - ny/2)^2 < (nx/16)^2, u[[i, j]] = 0.;
v[[i, j]] = 0.], {i, nx}, {j, ny}];
Do[
inew = Mod[i - nx dt u[[i, j]], nx - 1, 1];
jnew = Mod[j - nx dt v[[i, j]], ny - 1, 1];
iL = Floor@inew; jL = Floor@jnew;
iR = iL + 1; jR = jL + 1;
unew[[i, j]] = interfunc[u];
vnew[[i, j]] = interfunc[v];
, {i, nx}, {j, ny}]; {unew, vnew}]] /.
interfunc[v_] :>
inter[inter[v[[iL, jL]], v[[iL, jR]], jnew - jL],
inter[v[[iR, jL]], v[[iR, jR]], jnew - jL], inew - iL] //. interrule //
ReleaseHold;
viscosity…and…conservation =
compile[{{arg, _Complex, 3}, dt, mu}, Module[{x, y, nx, ny, k, u, v},
{u, v} = arg; {nx, ny} = Dimensions@v;
Do[x = Mod[i - 1, nx, -nx/2];
y = Mod[j - 1, ny, -ny/2];
k = x^2 + y^2;
With[{coef = Exp[-k dt mu], mid = -y u[[i, j]] + x v[[i, j]]},
If[k > 0,
u[[i, j]] = -y mid/k coef; v[[i, j]] = x mid/k coef, 0.]], {i, nx}, {j, ny}];
{u, v}]];

nx = 128; ny = nx; dt = 0.3; mu = 0.001; nt = 800;
v = Table[0., {2}, {nx}, {ny}];
separate = Transpose[#, {2, 3, 1}] &;
combine = Transpose[#, {3, 1, 2}] &;
vlst = Table[v = force…and…advection[v, dt];
v = Fourier /@ v;
v = viscosity…and…conservation[v, dt, mu];
v = InverseFourier /@ v // Re, {t, nt}]; // AbsoluteTiming
(* {3.28788, Null} *)

cut = nx/16;
disk = Graphics[{Opacity[1/2], Disk[{nx/4 - cut, ny/2}, nx/16]}];
arrayplot =
Show[ArrayPlot[#, DataReversed -> True, ColorFunction -> "Rainbow",
ImageSize -> nx, PlotRange -> {0, 1}], disk] &;

piclst = arrayplot /@
Rescale[vortex@separate@combine[#][[Round@cut ;;]] & /@
vlst[[;; ;; 10]]]; // AbsoluteTiming
(* {2.35815, Null} *)

piclst // ListAnimate


## Remark

1. If you don't have a C compiler installed, please take away the CompilationTarget -> "C" option and the // Last in definition of compile. I do recommend you to install one, though.

2. The color of argu in definition of compile is red, but don't worry.

3. You may have noticed the obtained flow is different from those obtained with the fixed original code. This isn't surprising because the implementations for periodic linear interpolation aren't exactly the same. Since the algorithm isn't a rigorous PDE solver, but more of a creator for something looks like flow, I'd like to leave the discrepency alone.

4. If nx is small, ListDensityPlot isn't a bad choice for visualization. The following is another possible definition of arrayplot using ListDensityPlot:

arrayplot =
ListDensityPlot[#, ColorFunction -> "Rainbow", ImageSize -> 4 nx,
PlotRange -> {0, 1}, AspectRatio -> Automatic, Epilog -> disk[[1]],
Mesh -> False, Frame -> False] &;

5. You may have noticed I've avoided arithmetic operation on lists in the code. For example, x = Mod[i - 1, nx, -nx/2]; y = Mod[j - 1, ny, -ny/2] can be shortened to {x, y} = Mod[{i, j} - 1, {nx, ny}, -{nx, ny}/2]. I don't do this because it significantly slows down the code when compiling to C.

6. I've used nx = 128, but even nx = 32 is enough to generate observable vortex street.

7. Viscosity turns out to be unnecessary to produce vortex street i.e. the Exp[-k dt mu] in code can actually be taken away. The resulting pattern looks more beautiful (at least in my view):

8. A large enough ny is necessary, if we cut ny to ny = nx/2:

As we can see, the obtained pattern isn't that great.

## A Version 3 Compatible Implementation

Just for fun, I further modified the code so it works from v3.0 to v13.1:

Off[General::spell]
Off[General::spell1]

SetAttributes[compile, HoldAll];

If[$VersionNumber < 8, compile[argu__] := Compile[argu], compile[argu__] := With[{cg = CompileGetElement}, Hold@Compile[argu, RuntimeOptions -> "Speed", CompilationTarget -> "C"] /. Part -> cg /. HoldPattern[cg[a__] = rhs_] :> (Part[a] = rhs)] // ReleaseHold(*// Last*)]; vortex = compile[{{v, _Real, 3}}, Module[{nx, ny}, {nx, ny} = Dimensions@v[[1]]; Table[((v[[2, i + 1, j]] - v[[2, i - 1, j]]) - (v[[1, i, j + 1]] - v[[1, i, j - 1]]))/2, {j, 2, ny - 1}, {i, 2, nx - 1}]]]; interrule = inter[valueL_, valueR_, scale_] :> (1 - scale) valueL + scale valueR; force…and…advection = Hold@compile[{{arg, _Real, 3}, {dt, _Real}}, Module[{u, v, unew, vnew, nx, ny, inew, jnew, iL, jL, iR, jR}, {u, v} = arg; {nx, ny} = Dimensions@v; unew = vnew = Table[0., {nx}, {ny}]; (*Alternative force:*)(*Do[v[[1,j]]={.1,0.},{j,ny}];*) Do[If[i < 1 + nx/16, u[[i, j]] = 0.1; v[[i, j]] = 0., 0.]; If[(i - nx/4)^2 + (j - ny/2)^2 < (nx/16)^2, u[[i, j]] = 0.;v[[i, j]] = 0., 0.], {i, nx}, {j, ny}]; Do[inew = Mod[i - nx dt u[[i, j]] - 1, nx - 1] + 1; jnew = Mod[j - nx dt v[[i, j]] - 1, ny - 1] + 1; iL = Floor@inew; jL = Floor@jnew; iR = iL + 1; jR = jL + 1; unew[[i, j]] = interfunc[u]; vnew[[i, j]] = interfunc[v];, {i, nx}, {j, ny}]; {unew, vnew}]] /. interfunc[v_] :> inter[inter[v[[iL, jL]], v[[iL, jR]], jnew - jL], inter[v[[iR, jL]], v[[iR, jR]], jnew - jL], inew - iL] //. interrule // ReleaseHold; viscosity…and…conservation = compile[{{arg, _Complex, 3}, dt, mu}, Module[{x, y, nx, ny, k, u, v}, {u, v} = arg; {nx, ny} = Dimensions@v; Do[ x = Mod[i - 1 + nx/2, nx] - nx/2; y = Mod[j - 1 + ny/2, ny] - ny/2; k = x^2 + y^2; With[{coef = Exp[-k dt mu], mid = -y u[[i, j]] + x v[[i, j]]}, If[k > 0, u[[i, j]] = -y mid/k coef; v[[i, j]] = x mid/k coef, 0.]], {i, nx}, {j, ny}]; {u, v}]]; (*Obtained with FindFormula:*) rainbowcore = compile[s, {0.4812174314044116 - 3.7442543302272293 s + 21.92055029805457 s^2. - 63.09750453225312 s^3. + 108.54708507988678 s^4.002496025211754 - 94.17958853869598 s^5. + 30.931639079029107 s^6., 0.12036767137989492 - 1.2465284793071336 s + 26.535159362530283 s^2 - 94.55369722465234 s^3 + 155.26914677786502 s^4 - 124.36433714169121 s^5 + 38.37235660032842 s^6, 0.5273478922125036 + 2.1548191214387797 s + 2.632115374325022 s^2 - 50.21875613715874 s^3 + 112.94351047860702 s^4 - 98.6320053106276 s^5 + 30.726738288661593 s^6}]; rainbowfunc = RGBColor @@ rainbowcore@# &; If[$VersionNumber < 5.1, Rescale = (# - Min[#])/(Max[#] - Min[#]) &];
nx = 64; ny = nx; dt = 0.3; mu = 0.001; nt = 800;
v = Table[0., {2}, {nx}, {ny}];
separate = Transpose[#, {2, 3, 1}] &;
combine = Transpose[#, {3, 1, 2}] &;
vlst = Table[v = force…and…advection[v, dt];
v = Fourier /@ v;
v = viscosity…and…conservation[v, dt, mu];
v = InverseFourier /@ v // Re, {t, nt}]; // Timing

cut = nx/16;
If[$VersionNumber < 6, disk = Disk[{nx/4 - cut, ny/2}, nx/16]; arrayplot = ListDensityPlot[#, ColorFunction -> rainbowfunc, ImageSize -> 4 nx, PlotRange -> {0, 1}, AspectRatio -> Automatic, Epilog -> disk, Mesh -> False] &, disk = {Opacity[1/2], Disk[{nx/4 - cut, ny/2}, nx/16]}; arrayplot = ArrayPlot[#, DataReversed -> True, ColorFunction -> "Rainbow", Epilog -> disk, ImageSize -> 4 nx, PlotRange -> {0, 1}] &;]; piclst = arrayplot /@ Rescale[vortex@ separate@Drop[combine[#1], Round[cut]] & /@ (Transpose[ Partition[vlst, 10]][[1]])]; // Timing If[$VersionNumber < 6, "You can use Ctrl+Y to create animation. ",
piclst // ListAnimate]


Warning: since packed array isn't yet introduced in v3.0, the code is rather slow therein. For nx = 64; ny = nx/2, the timing is about 461 seconds on 2GHz machine.

The timing will reduce to about 17 seconds in v4, and about 15 seconds in v5.

• It is very good solution (+1). Oct 25, 2022 at 5:26
• The name for module viscosity\[Ellipsis]and\[Ellipsis]conservation is not right. We don't use conservation, not solve equation $\nabla.\vec{v}$ in this algorithm. The last step of predictor-corrector algorithm is projection. For this step we can use FFT. Jos Stam combined diffusion step and projection step in one step. So, the right name could be viscosity\[Ellipsis]and\[Ellipsis]projection. Oct 26, 2022 at 3:02
• @Alex Jam Stam uses the term "Conservation of Mass" in his paper…: i.stack.imgur.com/idgyj.png Oct 26, 2022 at 3:45
• Yes, but he used name Project for projection module in the paper dgp.toronto.edu/public_user/stam/reality/Research/pdf/ns.pdf . The advance of the projection step is that we compute pressure gradient and divergent free velocity field in one step. See my answer on mathematica.stackexchange.com/questions/261185/… Oct 26, 2022 at 4:17
• @AlexTrounev I think it's not surprising. The paper you linked is a rigorous PDE solver, and the projection is imposed on the whole RHS $-(\mathbf{u}\cdot\nabla)\mathbf{u}+\nu \nabla^2 \mathbf{u}+\mathbf{f}$, but the algorithm discussed in this question is, as mentioned above, more of a creator of something looks like a flow, and only the speed field is considered, so it's OK to call/explain this step as "force our field to be mass conserving". Oct 26, 2022 at 6:08

As well known the stably fluids algorithm is some kind of predictor corrector algorithm - see my answer here. This algorithm includes 3 steps - advection, diffusion and projection. In Fourie space the diffusion and projection can be combine in one step as follows

$$(u-u_n)/dt+(u.\nabla) u=0, (u_{n+1}-u)/dt+\nabla p-\mu \nabla^2 u=0$$

Apply $$\nabla .$$ to the last equation and use $$\nabla .u_n=0$$, then we have

$$\nabla^2 p-\nabla. u/dt-\mu \nabla^2( \nabla .u)=0$$

Using FFT we can transform last 2 linear equations to the system of algebraic equations and express $$u_{n+1}$$ Fourie image as

$$\vec{u}_{n+1}=\vec{u}-(\vec{k}.\vec{u})\vec{k}/k^2+\mu dt k^2((\vec{k}.\vec{u})\vec{k}/k^2-\vec{u})$$

where $$\vec{k}=(k_x,k_y),k^2=k_x^2+k_y^2$$.

Using Mathematica FFT, we should know how $$k_x, k_y$$ look like. These vectors are given by FFT matrix in a form

mat = Table[E^(2 \[Pi] I (r - 1) (s - 1)/n), {r, 1, n}, {s, 1, n}];r = (Log[Flatten[mat]]/I) // DeleteDuplicates;


Consequently, $$k_x=r, k_y=r$$ since FFT for matrix A in 2D is mat.A.mat (the normalization depends on FourierParameters option). The final code is given by

n = 64; dt = 0.3; mu = 0.0001; nt = 300; mat =
Table[E^(
2 \[Pi] I (r - 1) (s - 1)/n), {r, 1, n}, {s, 1,
n}]; r = (Log[Flatten[mat]]/I) // DeleteDuplicates;
v = Table[{0., 0.}, {n}, {n}];
separate = Transpose[#, {2, 3, 1}] &;
combine = Transpose[#, {3, 1, 2}] &;

g = Graphics[{Blue, Disk[{n/4, n/2}, n/16]}];

Do[Do[If[i < 1 + n/16, v[[i, j]] = {0.1, 0.}];
If[(i - n/4)^2 + (j - n/2)^2 < (n/16)^2,
v[[i, j]] = {0., 0.}], {i, n}, {j, n}];
{ui, vi} = Module[{umat = #}, umat[[-1]] = umat[[1]];
umat[[All, -1]] = umat[[All, 1]];
ListInterpolation[umat, InterpolationOrder -> 1,
PeriodicInterpolation -> True]] & /@ separate@v;
v = Table[{i2, j2} = {i, j} - n dt v[[i, j]];
{ui[i2, j2], vi[i2, j2]}, {i, n}, {j, n}];
v = Fourier /@ separate@v // combine;
v = Table[x = r[[i]];
y = r[[j]];
k = x^2 + y^2;
If[k > 0, (v[[i, j]] - (v[[i, j]] . {x, y}) {x, y}/k) +
mu dt ((v[[i, j]] . {x, y}) {x, y} - k v[[i, j]]),
v[[i, j]]], {i, n}, {j, n}];
v = (InverseFourier /@ separate@v // combine // Re);
vs[t] = v;, {t, nt}]; // AbsoluteTiming


It takes about one minute on my laptop. Visualization

lst = Table[
Show[ListDensityPlot[
Table[((vs[t][[i + 1, j, 2]] -
vs[t][[i - 1, j, 2]]) - (vs[t][[i, j + 1, 1]] -
vs[t][[i, j - 1, 1]]))/2, {j, 2, n - 1}, {i, 2, n - 1}],
Frame -> False, ColorFunction -> Hue, PlotRange -> All,
ImageSize -> Tiny], g], {t, 2, nt, 2}];
ListAnimate[lst]


• I am not quite sure about your implementation. However, if you want to implement an simple spectral code using FFT with non-periodic boundary conditions the easiest way would be to consider a Chebyshev expansion in x-direction and a Fourier expansion in y-direction. Or both directions using Chebyshev. Oct 21, 2022 at 13:12
• @ConvexHull You are right. We can't use FFT to solve this problem since domain has a hole. Oct 21, 2022 at 15:55
• @xzczd Yes, you are right. If we use same normalization at k=0, then Max[Abs[(vAlex - v) ]] is about 10^-15 for your test as well. I see, that my normalization is not good for some fields. Oct 25, 2022 at 5:09
• Use p = pk Exp[I k1 x + I k2 y]; u = {uk, vk} Exp[I k1 x + I k2 y]; u1 = {u1k, v1k} Exp[I k1 x + I k2 y]; and solve equations eq1 = Laplacian[p, {x, y}] - Div[u, {x, y}]/dt - mu Laplacian[Div[u, {x, y}], {x, y}]==0; eq2 = (u1 - u)/dt + Grad[p, {x, y}] - mu Laplacian[u, {x, y}] to express pk, u1k, v1k. Answer for velocity {{u1k -> ( k2^2 uk - dt k1^2 k2^2 mu uk - dt k2^4 mu uk - k1 k2 vk + dt k1^3 k2 mu vk + dt k1 k2^3 mu vk)/(k1^2 + k2^2), v1k -> (-k1 k2 uk + dt k1^3 k2 mu uk + dt k1 k2^3 mu uk + k1^2 vk - dt k1^4 mu vk - dt k1^2 k2^2 mu vk)/(k1^2 + k2^2)}} Oct 27, 2022 at 14:55
• Therefore, equation for Fourie velocity image does not depends on pressure while pressure is given by pk -> (I (-k1 uk + dt k1^3 mu uk + dt k1 k2^2 mu uk - k2 vk + dt k1^2 k2 mu vk + dt k2^3 mu vk))/(dt (k1^2 + k2^2)) Oct 27, 2022 at 15:02

You must tell ListDensityPlot that it should print the image, otherwise it is absorbed in the Do loop.

From the following code I get 25 images. Due to space, I only show the first and last. I am not sure if it is what you want:

n = 64; dt = 0.3; mu = 0.001; v = Table[{0, 0}, {n}, {n}];

Do[

Do[If[i < 5, v[[i, j]] = {0.1, 0}];
If[(i - n/4)^2 + (j - n/2)^2 < 4^2, v[[i, j]] = {0, 0}], {i, 1,
n}, {j, 1, n}];

ui = ListInterpolation[v[[All, All, 1]]];
vi = ListInterpolation[v[[All, All, 2]]];

{ui, vi} =
Module[{ui = #}, ui[[-1]] = ui[[1]]; ui[[All, -1]] = ui[[All, 1]];
ListInterpolation[ui, InterpolationOrder -> 1,
PeriodicInterpolation -> True]] & /@ Transpose[v, {2, 3, 1}];

v = Table[{i2, j2} = {i, j} - n dt v[[i, j]];
{ui[i2, j2], vi[i2, j2]}, {i, 1, n}, {j, 1, n}];

v = Transpose[Map[Fourier[v[[All, All, #]]] &, {1, 2}], {3, 1, 2}];
v = Table[x = Mod[i + n/2, n] - n/2;
y = Mod[j + n/2, n] - n/2;
k = x^2 + y^2;
Exp[-k dt mu] If[k > 0, (v[[i, j]] . {-y, x}/k) {-y, x},
v[[i, j]]], {i, 1, n}, {j, 1, n}];
v = Transpose[
Map[Re[InverseFourier[v[[All, All, #]]]] &, {1, 2}], {3, 1, 2}];
ListDensityPlot[
Table[((v[[i + 1, j, 2]] - v[[i - 1, j, 2]]) - (v[[i, j + 1, 1]] -
v[[i, j - 1, 1]]))/2, {j, 2, n - 1}, {i, 2, n - 1}],
Mesh -> False, Frame -> False, ColorFunction -> (Hue[2 #/3] &)] //
Print, {t, 1, 25}];


• As already mentioned in 1st possible issue of question, this isn't a big problem, the real problem is, the obtained picture is vastly different from the GIF. Oct 16, 2022 at 8:26