Fluid simulation boundaries and advect

Question

This fluid simulation is based off of a paper by Stam. On page 7, he describes the basic idea behind advection:

Start with two grids: one that contains the density values from the previous time step and one that will contain the new values. For each grid cell of the latter we trace the cell’s center position backwards through the velocity field. We then linearly interpolate from the grid of previous density values and assign this value to the current grid cell.

Advect code. The two density grids are d and d0, u and v are velocity components, dt is the time step, N (global) is grid size, b can be ignored:

void advect(int b, vfloat &d, const vfloat &d0, const vfloat &u, const vfloat &v, float dt, std::vector<bool> &bound)
{
    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(b, d, bound);
}

This method is concise and works well enough, but implementing object boundaries is tricky for me to figure out because values are traced backwards and interpolated. My current solution is to simply push density out of boundaries if there is an empty space (or spaces) next to it, but that is inaccurate and causes density to build up, especially on corners and areas with diagonal velocity. only visually accurate. I'm looking for "correctness" now.

Relevant parts of my boundary code:

void set_bnd(const int b, vfloat &x, std::vector<bool> &bound)
{
    //...
    for (int i=1; i<=N; i++)
    {
        for (int j=1; j<=N; j++)
        {
            if (bound[IX(i,j)])
            {
                //...
                else if (b==0)
                {
                    // Distribute density from bound to surrounding cells
                    int nearby_count = !bound[IX(i+1,j)] + !bound[IX(i-1,j)] + !bound[IX(i,j+1)] + !bound[IX(i,j-1)];
                    if (!nearby_count) x[IX(i,j)] = 0;
                    else
                        x[IX(i,j)] = ((bound[IX(i+1,j)] ? 0 : x[IX(i+1,j)]) +
                                      (bound[IX(i-1,j)] ? 0 : x[IX(i-1,j)]) +
                                      (bound[IX(i,j+1)] ? 0 : x[IX(i,j+1)]) +
                                      (bound[IX(i,j-1)] ? 0 : x[IX(i,j-1)])) / surround;
                }
            }
        }
    }
}

bound is a vector of bools with rows and columns 0 to N+1. Boundary objects are set up before the main loop by setting cell coordinates in bound to 1.

The paper vaguely states "Then we simply have to add some code to the set_bnd() routine to fill in values for the occupied cells from the values of their direct neighbors", which is sort of what I'm doing. I am looking for a way to implement boundaries more accurately, that is having non-fluid solid boundaries and maybe eventually supporting boundaries for multiple fluids. Visual quality is much more important than physics correctness.

Answer 1

Your answer comes from physics rather than simulation. Since you're dealing with boundaries, your velocity field needs to meet the Prandtl no-slip boundary condition, which says simply that the velocity at the boundary must be zero. See https://en.wikipedia.org/wiki/Boundary_layer for (a lot) more information. If your velocity field does not meet this criterion, you'll have the difficulties you describe, including advecting mass back across a boundary, which is a pretty basic violation of the model.

You should also be aware that this advection code does not conserve density (by design) and that the conservation law is fixed up at the end. You'll need to pay attention to that step, since the Hodge decomposition of the vector field also has applicable boundary conditions.

Answer 2

You may be interested in "The Art of Fluid Animation" by Jos Stam (Sept. 2015). Around page 69 he discusses boundary conditions in some detail..

Perhaps also of interest: https://software.intel.com/en-us/articles/fluid-simulation-for-video-games-part-1/.

"The Perfect Storm" was a while ago so now your fluid sim has to be either very big, very fast, or very detailed. Preferably all three. Some might use a GPU if their use case allows.

Maybe it helps..