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.

Creatinine

Can you post the setup code for bound? My guess is that boundary cells are row 0, row N+1, column 0, and column N+1 [from diagram in paper]? Also, set_bnd has else if (b==0) and the if for that else [above] is missing. Can you add a bit more? I do see some problems [problem similar to video frame boundaries and motion vectors], but before I attempt an answer I'd like to have a bit more info. Also, when is b zero/non-zero? – Casease

@CraigEstey I've updated my question. Some code is removed to make the question more concise; the commented out parts of set_bnd deal with edge conditions that aren't needed for the question. A full version of my code is here – Creatinine

From what I see in the paper, you want to just use continuity on the density, not 'reflect' it from the boundary. This should mean (for simple cases as described by the paper) setting the boundary cell to be the same as the non-boundary cell it shares a face with. – Selinski

I've build it and have been running it [adding single step, debug, etc.]. For the bounds rect, the inner areas are obviously solid, but how do you want the edges? (e.g. Is [i15,j20] solid or fluid?). Or, from the left going right, what is the last fluid x value 14 or 15? I presume 14? This affects things. Should the edge of the solid have any non-zero [interpolated] density [or velocity] values? I presume not? How do you quantify errors? Negative density values? Anything else? How to spot "badness" visually or numerically from the console dump? – Casease

@CraigEstey Currently, edges are treated as fluid, but ideally they would be solid. I think I have already taken care of velocity, so density is the issue right now. "Badness" is largely density disappearing (not being pushed away) and areas producing density where they shouldn't. From the start this simulation has only gone for visual accuracy, so as long as it looks realistic that's enough. – Creatinine

Noticed [solid] edges as fluid--hence my Q. The way would be to up dims on rect by +1 each dir [i.e. now it's off-by-one]. To me, velocity should reflect at angle of incidence against wall? I've been considering replacing bounds vector with something more complex (e.g. linked-list of shapes) and suppressing calculations on inner points in solid [wasteful/inaccurate?]). What I noticed was gray accretion on the SW rect corner--is that the issue? Q: Why is your handling in set_bnd of b==1/b==2 different from Stam's [in loop above]? – Casease

@CraigEstey Currently I apply Stam's method for boundary velocity for my objects: by making horizontal and vertical velocity negative of what they would be, so the net effect is zero. The solid edges are fluid, just because that's how they're described in the paper. I don't think replacing the bounds vector will get much performance increase and the vector is the simplest way. The gray accretion was the issue, but that has been fixed by what lrm29 mentioned. Now the question isn't as important as I thought it would be - it's about correctness now. – Creatinine

What did you change to implement Irm29's fix [seems like I'm "downrev"]? I had been suspicious of the use of nearby_count--Is that what you changed? To me, top edge of rect should work like bottom edge of outer box, etc. The alternate way for bounds may just be list of points. Then, forall bpt in bound: bpt.i,bpt.j,bpt.xxx (e.g. precalc more stuff, so fewer if's in loop). To see better, I did the grid at bottom of render loop: SDL_SetRenderDrawColor(renderer,100,100,0,0); SDL_RenderDrawRect(renderer,&r); but some lines are double wide due to rounding – Casease

@CraigEstey I've updated the question. I hope you haven't spent too much time messing with my code (or if you have, at least it was worthwhile for you); I'm not sure what you're doing with bounds now but I'm just asking for a change in how advect is done – Creatinine

No worries. For me, got an intro to SDL and may add it to a gfx wrapper lib I have that does GTK3/GDK3 at present [install of latter is like 20 pkgs]. Stam's advect seems very similar to what is done in [H.264] video compression [motion vectors and macroblocks and convolution kernels, etc.]. Alternate inspiration for you? At this point, I'm going to assume you're on your own, but, just curious, what's wrong with the present advect? – Casease

@CraigEstey I guess I'm just looking for how to implement advect to work with collisions instead of what I'm doing now. The current code looks visually acceptable with stationary boundaries, but an extension like moving boundaries such as two fluids (fire and air, water and air) would need a "real" advect solution. – Creatinine

@Creatinine What do you mean by the 'correct way of advecting'? All you want is to ensure that there is no advection through the boundary, and you do this by making the gradient across it equal to zero. If you're really interested I'd recommend trying out OpenFOAM, a proper CFD code. There is a tutorial for transport of scalar quantities that you may find useful. – Selinski

@Creatinine You may get faster and better help over at cfd-online. – Selinski

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.

Mia

This answer is helpful but ultimately does not provide a solution, if one exists. – Creatinine

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..

Ramey