流体模拟“爆炸”答案

【问题标题】：Fluid Simulation "Blows Up"流体模拟“爆炸”
【发布时间】：2016-03-01 05:46:38
【问题描述】：

以下流体模拟是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。这是否意味着密度必须

void set_bnd(int N, int b, vector<float> &x)
{
    // removed
}


inline void lin_solve(int N, int b, 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, int b, vector<float> &x, vector<float> &x0, float diff, float dt)
{
    float a = dt*diff*N*N;
    lin_solve(N, b, x, x0, a, 1+4*a);

}

// Backwards advection
void advect(int N, int b, 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, b, d);
    }
}

void project(int N, vector<float> &u, vector<float> &v, 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, 0, div); set_bnd(N, 0, p);

    lin_solve(N, 0, 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, 1, u); set_bnd(N, 2, v);
}

密度和速度求解器

// Density solver
void dens_step(int N, vector<float> &x, vector<float> &x0, vector<float> &u, vector<float> &v, float diff, float dt)
{
    add_source(x, x0, dt);
    swap(x0, x); diffuse(N, 0, x, x0, diff, dt);
    swap(x0, x); advect(N, 0, x, x0, u, v, dt);
}

// Velocity solver: addition of forces, viscous diffusion, self-advection
void vel_step(int N, vector<float> &u, vector<float> &v, 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, 1, u, u0, visc, dt);
    swap(v0, v); diffuse(N, 2, v, v0, visc, dt);
    project(N, u, v, u0, v0);
    swap(u0, u); swap(v0, v);
    advect(N, 1, u, u0, u0, v0, dt); advect(N, 2, v, v0, u0, v0, dt);
    project(N, u, v, u0, v0);
}

我考虑过floating-point inconsistencies，但在使用-ffloat-store 编译后问题仍然存在。

【问题讨论】：

尝试在valgrind下运行程序。它可能会自动告诉您出了什么问题。我假设您使用的是 Linux....
在您的调试器中打开两个单独的实例，并肩并肩地查看两者的不同之处以及原因。也请放下那个宏。内联函数就可以了。
如果从float 更改为double 会发生什么？根据我的经验，怀疑浮点错误的计算应始终在double 中完成。
这与add_source()中缺乏规范化有关。当你的密度变得足够稳定时（x0 与x 非常相似），然后add_source() 有效地将x 乘以1.0+dt，导致你的爆炸。 DIFF 的高值通过在 lin_solve() 中更重地加权 x 来掩盖这种影响，这意味着有效乘数更接近 1（但仍高于 1）。修复方法是在add_source:x[i] = (x[i]+dt*s[i])/(1.0f+dt); 中标准化，或计算c = 1+4*a+dt;。如果您计算密度网格的平均值，您会发现在某个点上，平均值 [....]
[...] 开始增加1 + (dt / (4*a)) 的系数；使用给定的设置（dt=0.1、a=0.1*0.01*20*20=0.4），这是1+0.1/1.6 ~ 1.06：因此，一旦密度变得足够稳定，它的质量就会开始以每次迭代 6% 的速度增加。

标签： c++ sdl game-physics fluid-dynamics

【解决方案1】：

麻烦的根源之一是lin_solve。您的i 和j 循环从零开始，但您引用IX(i-1,j)，它将访问越界数组元素x[-1]。

【讨论】：

@qwr 但您仍然在 lin_solve 中引用元素 x[-1]。
另外，后面 i 和 j 值的 x 计算取决于早期 i 和 j 值的 x 计算。这是原始论文的一个特点，还是一个错误？
@JamesPicone 我也怀疑过，但它在原始论文中
@JamesPicone 是的，我相信这是一个跨越式的集成。所以故意的。

【解决方案2】：

看到这里，我立刻觉得我必须回答。我在这篇文章发表时读过这篇文章。我已经在 Android 上实现了他的东西并且很喜欢它。我什至在 2000 年代初在于默奥演讲时遇到了这个人，他是一个非常友好的人。而且个子高。 :)

所以问题。您没有进行速度传播步骤，如果我没记错的话，我认为这对于不“炸毁”至关重要。

【讨论】：

这是为了缩短这个 SO 问题而故意省略的，因为传播步骤中的值仍然“爆炸”。当我开始处理它时，我可以把这些步骤放回去。但您的帮助将不胜感激。

【解决方案3】：

问题与add_source() 中缺乏规范化有关。

当您的密度变得足够稳定时（x0 的分布与x 非常相似，达到比例因子），然后add_source() 有效地将x 乘以大约1+dt，导致您的指数爆炸。 DIFF 的高值通过在 lin_solve() 中对 x0 的权重更重于 x0 来掩盖这种影响，这意味着有效乘数更接近 1，但仍高于 1。

效果就是每次迭代都会增加更多的质量。如果它不能在边缘足够快地“散开”，它将开始堆积。一旦密度变得完全静止，它的质量将以1+dt/(4a) 确定的指数速率增加。

根据您给定的设置 (dt=0.1, a=0.1*0.01*20*20=0.4)，这是1+0.1/1.6 ~ 1.06。

解决方法是在 add_source 中规范化：

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

，或将 lin_solve() 的 c 参数计算为 1+4*a+dt。两者都会迫使质量下降。

【讨论】：