Sıvı Simülasyonu "Üfleme"

Question

Aşağıdaki sıvı simülasyonu, paper by Stam'un bir çevirisi. Gerçekten korkunç bir şey oldu. Program her zaman düşük bir DIFF=0.01 ile çalıştırıldığında, değerler küçükten başlar ve sonra hızla genişler veya "patlatır". Matematik rutinlerini dikkatli bir şekilde kontrol ettim. Kod bir 0.5 ile başladığından, matematiksel olarak çarparak ve bir grup sıfır ekleyerek, sonuç sıfır sıfıra yakın ve diğer vektörlere yakın olmalıdır.

Kod oldukça uzun, ben de bunları parçalara ayırdım ve fazladan kodları kaldırdım. Tüm başlangıç ​​ve SDL kodları eksi olarak sadece 120 satır vardır. Birkaç saat boyunca hiçbir değişiklik yapmadan değişiklikler yapmaya çalıştım, bu yüzden yardım çok takdir ediliyor.

Bazı deneylerden sonra, DIFF çok düşük ayarlandığında bazı kayan nokta hatası olabileceğine inanıyorum. Değer 0.01'dan 0.02'a artırıldığında, değerler patlamaz. Bununla birlikte, bunun tüm sorun olduğunu düşünmüyorum.

Net olmak gerekirse, 1201ProgramAlarm ve vidstige tarafından geçerli olan yanıtlar sorunu çözmez.

kalın kısımlarındaki bölümleri önemli parçalardır, gerisi tamlık içindir.

---

başlayarak şeyler,

#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); 

Matematik Diffüz rutinleri 1+4*a bölün rutinleri atlayın. Bu ima yoğunluğu < = 1 olmalıdır?

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); 
} 

Yoğunluk ve hız çözücü

// 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); 
} 

Ben floating-point inconsistencies kabul ama -ffloat-store ile derleme sonra sorun hala devam etti.

Answer 1

Sorun, add_source() numaralı telefondaki normalleştirme eksikliğiyle ilgilidir.

Eğer yoğunluk yeterince sabit (bir ölçek faktörüne kadar, x dağılımda x0 çok benzer) olur, daha sonra add_source() etkili üstel blowup giden yaklaşık 1+dt ile x çarpar.DIFF yüksek değerler etkili çarpan yakın 1 getirilir, yani lin_solve() daha ağır x0 üzerinde x tartılarak bu etkiyi maske, ancak 1.

etkisi üzerinde hala, daha sonra her yineleme ile, daha kitlesel eklendi. Kenarlarda yeterince hızlı "yayılmazsa", birikmeye başlayacaktır. Yoğunluk mükemmel bir şekilde sabitlendiğinde, 1+dt/(4a) tarafından belirlenen bir üstel oranda kütle içinde artacaktır. size verilen ayarlarla (dt=0.1, a=0.1*0.01*20*20=0.4) ile

, bu 1+0.1/1.6 ~ 1.06 olduğunu.

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

veya 1+4*a+dt olarak lin_solve() için c argüman hesaplamak için:

düzeltme add_source olarak normalize etmektir. Ya kütleyi düşürmeye zorlar.

Answer 2

Sorunların bir tanesi lin_solve'dadır. i ve j döngüleriniz sıfırdan başlar, ancak x[-1] sınır dizisi elemanına erişen IX(i-1,j) referansını kullanırsınız.

Answer 3

Bunu görmek, hemen cevap vermem gerektiğini hissettim. Yayınlandığı zaman bu makaleyi okudum. Eşyalarını Android'e uyguladım ve sadece sevdim. 2000'li yılların başında Umeå'da konuştuğumla tanıştım ve çok arkadaş canlısı biri. Ve uzun. Bu nedenle problem için. Bir hız yayılımı adımı atmıyorsunuz, bence doğru hatırlıyorsam bu "havaya uçurmak" için kritik bir durum.