c – 流体模拟“吹起”

以下流体模拟是 paper by Stam的翻译.发生了一件真正可怕的事情.每次程序以低DIFF = 0.01运行时,该值从小开始,然后迅速扩大或“爆炸”.我仔细检查了数学程序.由于代码从一个0.5开始,在数学上它是乘以并添加一堆零,所以最终结果应该接近零密度和其他向量.

代码很长,所以我把它分成多个块,并删除了额外的代码.减号所有的开始和SDL代码只有大约120行.我花了几个小时尝试更改无效,所以非常感谢帮助.

经过一些实验,我相信当DIFF设置得太低时可能会有一些浮点错误.当该值从0.01增加到0.02时,这些值不会爆炸.我不认为这是整个问题.

要清楚,1201ProgramAlarm和vidstige的当前答案不能解决问题.

大胆的部分是重要的部分,其余部分是完整的.

开始的东西,跳过

#include <SDL2/SDl.h>
#include <stdio.h>
#include <iostream>
#include <algorithm>


#define IX(i,j) ((i)+(N+2)*(j))
using namespace std;

// Constants
const int SCREEN_WIDTH = 600;
const int SCREEN_HEIGHT = 600;  // Should match SCREEN_WIDTH
const int n = 20;               // Grid size
const int SIM_LEN = 1000;
const int DELAY_LENGTH = 40;    // ms

const float VISC = 0.01;
const float dt = 0.1;
const float DIFF = 0.01;

const bool DISPLAY_CONSOLE = false; // Console or graphics
const bool DRAW_GRID = false; // implement later

const int nsize = (N+2)*(N+2);

数学例程漫射程序除以1 4 * a.这是否意味着密度必须<= 1？

void seT_Bnd(int n,int b,vector<float> &X)
{
    // removed
}


inline void lin_solve(int n,vector<float> &x,vector<float> &x0,float a,float C)
{
    for (int k=0; k<20; k++)
    {
        for (int i=1; i<=N; i++)
        {
            for (int j=1; j<=N; j++)
            {
                x[IX(i,j)] = (x0[IX(i,j)] + a*(x[IX(i-1,j)]+x[IX(i+1,j)]+x[IX(i,j-1)]+x[IX(i,j+1)])) / c;
            }
        }
        seT_Bnd ( N,b,x );
    }
}

// Add forces
void add_source(vector<float> &x,vector<float> &s,float dt)
{
    for (int i=0; i<nsize; i++) x[i] += dt*s[i];
}

// Diffusion with Gauss-SEIDel relaxation
void diffuse(int n,float diff,float dt)
{
    float a = dt*diff*N*N;
    lin_solve(N,x,x0,a,1+4*a);

}

// @R_696_11442@kWARDs advection
void advect(int n,vector<float> &d,vector<float> &d0,vector<float> &u,vector<float> &v,float dt)
{
    float dt0 = dt*N;
        for (int i=1; i<=N; i++)
        {
            for (int j=1; j<=N; j++)
            {
                float x = i - dt0*u[IX(i,j)];
                float y = j - dt0*v[IX(i,j)];
                if (x<0.5) x=0.5; if (x>N+0.5) x=N+0.5;
                int i0=(int)x; int i1=i0+1;
                if (y<0.5) y=0.5; if (y>N+0.5) y=N+0.5;
                int j0=(int)y; int j1=j0+1;

                float s1 = x-i0; float s0 = 1-s1; float t1 = y-j0; float t0 = 1-t1;
                d[IX(i,j)] = s0*(t0*d0[IX(i0,j0)] + t1*d0[IX(i0,j1)]) +
                             s1*(t0*d0[IX(i1,j0)] + t1*d0[IX(i1,j1)]);
            }
        }
        seT_Bnd(N,d);
    }
}

void project(int n,vector<float> &p,vector<float> &div)
{
    float h = 1.0/N;
    for (int i=1; i<=N; i++)
    {
        for (int j=1; j<=N; j++)
        {
            div[IX(i,j)] = -0.5*h*(u[IX(i+1,j)] - u[IX(i-1,j)] +
                                   v[IX(i,j+1)] - v[IX(i,j-1)]);
            p[IX(i,j)] = 0;
        }
    }
    seT_Bnd(N,div); seT_Bnd(N,p);

    lin_solve(N,p,div,1,4);

    for (int i=1; i<=N; i++)
    {
        for (int j=1; j<=N; j++)
        {
            u[IX(i,j)] -= 0.5*(p[IX(i+1,j)] - p[IX(i-1,j)])/h;
            v[IX(i,j)] -= 0.5*(p[IX(i,j+1)] - p[IX(i,j-1)])/h;
        }
    }
    seT_Bnd(N,u); seT_Bnd(N,2,v);
}

密度和速度求解器

// Density solver
void dens_step(int n,float dt)
{
    add_source(x,dt);
    swap(x0,X); diffuse(N,diff,X); advect(N,u,v,dt);
}

// VeLocity solver: addition of forces,viscous diffusion,self-advection
void vel_step(int n,vector<float> &u0,vector<float> &v0,float visc,float dt)
{
    add_source(u,u0,dt); add_source(v,v0,dt);
    swap(u0,u); diffuse(N,visc,dt);
    swap(v0,v); diffuse(N,dt);
    project(N,v0);
    swap(u0,u); swap(v0,v);
    advect(N,dt); advect(N,v0);
}

我考虑过floating-point inconsistencies,但是在用-ffloat-store编译之后,问题仍然存在.

x[i] = (x[i]+dt*s[i])/(1.0f+dt);