基于网格的流体模拟误差

所以我正在写一个基于 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)

这是一个单独的类,这样网格类就可以处理坐标的包装(对于边界条件,我太懒了)。

平流:

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,%f\n",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: 过了一会儿,网格平均到均匀密度,没有振荡。