想了很久，还是觉得这个“问题”不太适合发在 stackexchange ，姑且在中文社区试试水吧。
大约八年前，我在这个网站（需爬墙）看到了一段有问题的卡门涡街模拟，代码如下：
(* 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}];
如注释所说，这段代码有问题。在过去八年时间里，每隔一两年，我就会心血来潮地试图修复这段代码，然后以失败告终。大家有什么主意吗？
照例说说我的失败尝试。按网页上的说法，这段代码参考的是 Jos Stam 的论文《A Simple Fluid Solver based on the FFT》及其中给出的C代码。为方便阅读，论文及代码我已在本站传了一份：
论文《A Simple Fluid Solver based on the FFT》
solver.c
但这篇文章的细节并不十分清楚，我的C语言水平也很有限（我甚至搞不清楚文章中所用的 fftw 库要怎么在 windows 下“安装”），所以看完了照样不知道该怎么修。目前，我能看出来的只有2点：
1. 这段C代码的快速傅立叶变换基于的是 fftw-2.1.5 ，从相应的文档来看，该程序包所含的快速傅立叶变换大概相当于 Fourier[…, FourierParameters -> {1, -1}] ，但把 FourierParameters 加上似乎并不足以修复上述代码。
2. ListInterpolation 应该是得强制搞成周期性的，但这依旧不足以修复上述代码。
顺便，网页上附了一个模拟结果：
据网页上的解说，此模拟用的是C++，但具体代码未给出。
