Следующее моделирование жидкости является переводом статья Стама. Произошло нечто действительно ужасное. Каждый раз, когда программа запускается с низким 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, 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);
}
я считал несоответствия с плавающей точкой, но после компиляции с -ffloat-store
проблема все еще сохраняется.
Проблема связана с отсутствием нормализации в add_source()
,
Когда ваша плотность становится достаточно постоянной (x0
очень похожа по распределению на x
до коэффициента масштабирования), затем add_source()
эффективно умножает x
примерно 1+dt
, приводя к вашему экспоненциальному взрыву. Высокие значения DIFF
маскировать этот эффект путем взвешивания x
более тяжело над x0
в lin_solve()
Это означает, что эффективный множитель приближается к 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);
, или вычислить c
аргумент lin_solve()
как 1+4*a+dt
, Либо заставит массу падать.
Один из источников проблем в lin_solve
, Ваш i
а также j
Циклы начинаются с нуля, но вы ссылаетесь IX(i-1,j)
, который будет обращаться к элементу массива вне границ x[-1]
,
Увидев это, я сразу почувствовал, что должен ответить. Я прочитал эту статью еще тогда, когда она была опубликована. Я реализовал его вещи на Android и просто в восторге. Я даже встретил этого парня, когда говорил в Умео в начале 2000-х, и он очень дружелюбный парень. И высокий. 🙂
Так что к проблеме. Вы не делаете шаг распространения скорости, я думаю, что это важно для того, чтобы не «взорваться», если я правильно помню.