基于栅格的流体模拟错误

所以我正在写一个基于http://jamie-wong.com/2016/08/05/webgl-fluid-simulation/.

这是用C++编写的。

class Grid只是一个具有以下内容的类：

Grid.data[][] : float[][] - floating point data
Grid.copy(): Grid -  creates a new grid with identical elements
Grid.set(int x,int y, float f): void - sets the data at (x,y)
Grid.operator() (int x,int y): float - gets the data at (x,y)

这是一个单独的类的原因是Grid类可以处理坐标的环绕(我对边界条件太懒了(。

平流：

void advect(Grid &from, Grid &to, Grid &vx, Grid &vy);

我确信平流函数是完美的。

这些是压力和发散函数：

void calcDivergence(Grid &density, Grid &divergence, Grid &vx, Grid &vy){
for(int x=0;x<GRIDSIZE;x++){
for(int y=0;y<GRIDSIZE;y++){
divergence.set(x,y,  (-2/dt)*(vx(x+1,y)-vx(x-1,y)+vy(x,y+1)-vy(x,y-1))   );
}
}
}
void solvePressure(Grid &divergence, Grid &pressure){
pressure.init();//set to 0
Grid pressure0=pressure.copy();
for(int i=0;i<PRESSURE_SOLVE_ITERS;i++){
for(int x=0;x<GRIDSIZE;x++){
for(int y=0;y<GRIDSIZE;y++){
pressure.set(x,y,  (divergence(x,y)+pressure0(x-1,y)+pressure0(x+1,y)+pressure0(x,y-1)+pressure0(x,y+1))/4.0f   );
}
}
pressure0=pressure.copy();
}
}
void fixDivergence(Grid &pressure, Grid&vx, Grid &vy){
for(int x=0;x<GRIDSIZE;x++){
for(int y=0;y<GRIDSIZE;y++){
float gradX=dt*0.5*(pressure(x+1,y)-pressure(x-1,y));
float gradY=dt*0.5*(pressure(x,y+1)-pressure(x,y-1));
//            printf("%f,%fn",gradX,gradY);//why is this zero?  solvePressure must be broken, check swap code
vx.set(x,y,  vx(x,y)-gradX);
vy.set(x,y,  vy(x,y)-gradY);
}
}
}

用于确定压力求解器的精度：

float pressureAccuracy(Grid &density, Grid &pressure){
float f=0;
for(int x=0;x<GRIDSIZE;x++){
for(int y=0;y<GRIDSIZE;y++){
float realValue=4*pressure(x,y)-pressure(x-1,y)-pressure(x+1,y)-pressure(x,y-1)-pressure(x,y+1);
f+=abs(realValue-density(x,y));
}
}
return f/GRIDSIZE/GRIDSIZE;
}

最后，更新功能：

void update(Grid &vx, Grid &vy, Grid &density, Grid &pressure, Grid &divergence){
Grid vx0=vx.copy();
Grid vy0=vy.copy();
advect(vx0,vx,vx0,vy0);
advect(vy0,vy,vx0,vy0);
calcDivergence(density, divergence, vx, vy);
solvePressure(divergence, pressure);
fixDivergence(pressure, vx, vy);
Grid density0=density.copy();
advect(density0,density,vx,vy);
}

网格初始化如下：

vx[x][y]=cos(4pi*y/GRIDSIZE)
vy[x][y]=sin(4pi*x/GRIDSIZE)
pressure[x][y]=0
density[x][y]=1 if (x/100)%2+(y/100)%2==1, otherwise 0
divergence[x][y]=0

网格大小为200。

这是屏幕截图：https://ibb.co/dmi66K

出了什么问题？

编辑1:澄清：当你看到形状奇怪的密度等值线时，这些等值线随机振荡，但永远保持在同一位置。

编辑2:过了一段时间，网格平均密度均匀，没有振荡。