代码很长,所以我把它分成多个块,并删除了额外的代码.减号所有的开始和SDL代码只有大约120行.我花了几个小时尝试更改无效,所以非常感谢帮助.@H_404_3@
经过一些实验,我相信当DIFF设置得太低时可能会有一些浮点错误.当该值从0.01增加到0.02时,这些值不会爆炸.我不认为这是整个问题.@H_404_3@
要清楚,1201ProgramAlarm和vidstige的当前答案不能解决问题.@H_404_3@
大胆的部分是重要的部分,其余部分是完整的.@H_404_3@
开始的东西,跳过@H_404_3@
#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? @H_404_3@
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); } // Backwards 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); }
密度和速度求解器@H_404_3@
// 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编译之后,问题仍然存在.@H_404_3@
解决方法
当您的密度变得足够稳定(x0在x分布非常相似时,直到比例因子),那么add_source()有效地将x乘以约1 dt,导致您的指数爆发. DIFF的高值通过在lin_solve()中称重x超过x0来掩盖此效果,这意味着有效乘数接近1,但仍高于1.@H_404_3@
那么效果就在于每一次迭代都会增加更多的质量.如果它不能在边缘快速“展开”,它将开始堆积.一旦密度变得完全静止,它将以由1 dt /(4a)确定的指数速率增加质量.@H_404_3@
给定的设置(dt = 0.1,a = 0.1 * 0.01 * 20 * 20 = 0.4),这是1 0.1 / 1.6〜1.06.@H_404_3@
修正是在add_source中进行规范化:@H_404_3@
x[i] = (x[i]+dt*s[i])/(1.0f+dt);
,或者将lin_solve()的c参数计算为1 4 * a dt.要么会迫使群众掉下来.@H_404_3@