流体シミュレーション "ブローアップ"

Question

以下の流体シミュレーションは、paper by Stamの翻訳です。本当に恐ろしいことが起きました。プログラムが低いDIFF=0.01で実行されるたびに、値は小さくなり始め、次に急速に拡大したり、 "爆発"します。私は慎重に数学ルーチンをチェックしました。コードは1つの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, 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でコンパイルした後に問題がまだ持続しました。

Answer 1

add_source()に正規化の不足が問題です。

あなたの密度が十分に固定（スケールファクタまで、xに分布x0非常に類似している）になり、その後add_source()効果的に指数関数的爆発につながる、約1+dtによってxを乗算します。DIFFの高い値は、効果的な乗数が近い1にもたらしたが、その効果1.

の上にまだある、そしてそれは、反復ごとに複数の大容量であることを意味し、lin_solve()でx0上より重くxの重量を測定することによって、この効果をマスクしますが加えられる。エッジで十分速く「広がらない」ことができない場合、それは積み重なるようになります。密度が完全に静止すると、1+dt/(4a)によって決定される指数関数的な速度で質量が増加します。

指定された設定（dt=0.1, a=0.1*0.01*20*20=0.4）の場合、これは1+0.1/1.6 ~ 1.06です。

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

、または1+4*a+dtようlin_solve()からc引数を計算するために：

修正はadd_sourceで正規化することです。どちらの場合でも強制的に質量が落ちます。

Answer 2

問題の原因の1つはlin_solveです。 iとjループはゼロから開始しますが、範囲外の配列要素x[-1]にアクセスするIX(i-1,j)を参照しています。

Answer 3

これを見てすぐに私は答えなければならないと感じました。私はこの記事を、それが出版されたときに戻って読んだ。私はAndroidで彼のものを実装して、それを愛しています。私は2000年代初めにウメオで話すときに仲間に会いました。彼はとてもフレンドリーです。そして背の高い。 :)

だから問題に。あなたは速度の伝播のステップを行っていない、私は正しく覚えている場合、これは "爆破"しないために重要だと思う。