11

The following fluid simulation is a translation of a paper by Stam. Something truly terrible has happened. Each time the program is run with a low DIFF=0.01, the values start off small and then rapidly expand, or "blow up". I have checked the math routines carefully. Since the code starts off with one 0.5, mathematically it is multiplying and adding a bunch of zeros, so the end result should be close to zero density and other vectors.

The code is quite long, so I've separated it into chunks and removed extra code. Minus all the beginning and SDL code there are only about 120 lines. I have spent a few hours trying changes to no avail, so help is greatly appreciated.

After some experimentation I believe there may be some floating-point error when DIFF is set too low. When the value is increased from 0.01 to 0.02, the values don't blow up. I don't think this is the entire issue, though.

To be clear, the current answers by 1201ProgramAlarm and vidstige do not resolve the problem.

Sections in bold are important parts, the rest is for completeness.

Beginning stuff, skip

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

Math routines Diffuse routines divide by 1+4*a. Does this imply density must be <= 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 and velocity solver

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

I considered floating-point inconsistencies, but after compiling with -ffloat-store the problem still persisted.

Community
  • 1
  • 1
qwr
  • 6,786
  • 3
  • 42
  • 72
  • 1
    Try running the program under `valgrind`. It might automatically tell you what is wrong. I assume you're on Linux.... – John Zwinck Nov 27 '15 at 04:42
  • Open two separate instances in your debugger and step side-by-side and see where the two diverge and why. Also drop that macro please. An inline function will do. – Neil Kirk Nov 27 '15 at 04:44
  • @JohnZwinck Ok, I am new to debugging and will try this – qwr Nov 27 '15 at 04:46
  • 1 - Which C++ version are you using? 2 - *... which do (or should) not introduce new variables...*; do you mean *... not modifying any variables...*? If yes, then keep an eye on the reference parameters. You might be modifying them without knowing. 3 - `swap` does what it says: swapping two elements. It uses move semantics, but you can ignore that. It simply swaps, period. – Ely Nov 27 '15 at 04:58
  • 4 - In regard to *static* see http://stackoverflow.com/a/3698179/1566187 5 - If you can provide the code somewhere maybe we could play around. You would need to let us know what is the expected output for a certain input (or better many input/output pairs)... – Ely Nov 27 '15 at 05:11
  • @Elyasin The question has been greatly modified. 1. C++11 gcc 4.8.3 no optimizations 2. does not apply any more 5. all the code needed to run the program is included. Since all the math routines are adding and multiplying a bunch of zeros and one value, the expected output should tend towards zero. – qwr Nov 27 '15 at 06:02
  • 1
    What happens if you change from `float` to `double`? Calculations where floating-point errors are suspected should always be done in `double` in my experience. – SirGuy Dec 16 '15 at 19:59
  • 2
    It's related to lack of normalization in `add_source()`. When your density becomes sufficiently stationary (`x0` very similar to `x`), then `add_source()` effectively multiplies `x` by `1.0+dt`, leading to your blowup. High values of `DIFF` mask this effect by weighing `x` more heavily in `lin_solve()`, meaning that the effective multiplier is brought closer to 1 (but is still above 1). The fix is to normalize in `add_source`: `x[i] = (x[i]+dt*s[i])/(1.0f+dt);`, or to compute `c = 1+4*a+dt;`. If you compute the average of your density grid, you'll see that at some point the average [....] – Iwillnotexist Idonotexist Dec 16 '15 at 22:04
  • 1
    [...] starts increasing by a factor of `1 + (dt / (4*a))`; With your given settings (`dt=0.1`, `a=0.1*0.01*20*20=0.4`), this is `1+0.1/1.6 ~ 1.06`: Ergo, once the density becomes sufficiently stationary, it starts increasing in mass at a rate of 6% increase per iteration. – Iwillnotexist Idonotexist Dec 16 '15 at 22:08
  • @IwillnotexistIdonotexist I will be able to test it in a few hours, but wow! please make it an answer. I did not expect something so innocuous. – qwr Dec 16 '15 at 22:33

3 Answers3

6

The problem is related to a lack of normalization in add_source().

When your density becomes sufficiently stationary (x0 very similar in distribution to x, up to a scale factor), then add_source() effectively multiplies x by about 1+dt, leading to your exponential blowup. High values of DIFF mask this effect by weighing x more heavily over x0 in lin_solve(), meaning that the effective multiplier is brought closer to 1, but is still above 1.

The effect, then is that with every iteration, more mass is added. If it cannot "spread out" fast enough at the edges, it will start piling up. Once the density becomes perfectly stationary, it will increase in mass at an exponential rate determined by 1+dt/(4a).

With your given settings (dt=0.1, a=0.1*0.01*20*20=0.4), this is 1+0.1/1.6 ~ 1.06.

The fix is to normalize in add_source:

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

, or to compute the c argument to lin_solve() as 1+4*a+dt. Either will force the mass to drop.

qwr
  • 6,786
  • 3
  • 42
  • 72
Iwillnotexist Idonotexist
  • 12,747
  • 3
  • 39
  • 62
3

One source of trouble is in lin_solve. Your i and j loops start at zero, but you reference IX(i-1,j), which will access the out of bounds array element x[-1].

1201ProgramAlarm
  • 30,320
  • 7
  • 40
  • 49
2

Seeing this I immediately felt I had to answer. I read this article way back when it was published. I've implemented his stuff on Android and just love it. I even met the fellow when speaking at Umeå in the early 2000s, and he's a very friendly fellow. And tall. :)

So to the problem. You are not doing a velocity propagation step, I think this is critical for not "blowing up" if I remember correctly.

vidstige
  • 11,161
  • 7
  • 58
  • 95
  • That was purposely left out to shorten this SO question, as values still "blow up" with the propagation step. I can put those steps back in when I get around to it. But your help would be much appreciated. – qwr Dec 16 '15 at 18:41