Numerical Hodge decomposition with boundaries

Question

For a fluid simulation, I'm using the algorithm proposed by Jos Stam (http://www.intpowertechcorp.com/GDC03.pdf).

One step of the algorithm is a routine that projects the velocity vector field so that it becomes an incompressible field (div F = 0). To achieve that, it uses the Hodge decomposition to decompose the field into the sum an incompressible field (div F = 0) and a gradient field (curl F = 0).

Now we'd like to introduce boundaries that act like walls in the fluid. The fluid on one side of the wall must not influence the fluid on the other side.

This is achieved easily in the advect(), setBnd() and diffuse() steps of the solver. However, I'm not sure how to introduce this boundary separation in the project() step. My first approach was to just cut all the influences of those 'cell neighbors' that are on the other side of the boundary (see image above: the orange arrows represent neighbour cells across the boundary).

But this leads to an unstable solution. The fluid begins to move across the boundary without any force present.

Is there an easy solution to this problem?

Here's the code of the project() routine from Jos Stam's paper:

void project ( int N, float * u, float * v, float * p, float * div )
{
    int i, j, k;
    float h;
    h = 1.0/N;

    // calculate div
    for ( i=1 ; i<=N ; i++ ) {
        for ( 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 boundaries
    set_bnd ( N, 0, div ); set_bnd ( N, 0, p );

    // solve linear equation system using Gauss Seidel
    for ( k=0 ; k<20 ; k++ ) {
        for ( i=1 ; i<=N ; i++ ) {
            for ( j=1 ; j<=N ; j++ ) {
                p[IX(i,j)] = (div[IX(i,j)]+
                    p[IX(i-1,j)]+p[IX(i+1,j)]+
                    p[IX(i,j-1)]+p[IX(i,j+1)])/4;
            }
        }
        set_bnd ( N, 0, p );
    }

    // subtract the gradient field from our vector field ==> result has zero divergence
    for ( i=1 ; i<=N ; i++ ) {
        for ( 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 );
}

Answer 1

Solved it using the following approach, which works great. It also seems to make some sense:

• First, mark the 'arrows' which cross a boundary. In the image above, they are the orange arrows.
• In the calculate div step of the code, check which arrows are orange. For all 'orange neighbours', replace their values (for example u[IX(i+1,j)] for the right neighbour) by the respective value of the cell itself (for example u[IX(i,j)]).
• In the solve linear equation step, remove the values of the 'orange neighbours' from the inner sum, and set the divisor to 4 minus the number of orange neighbours.
• In the subtract gradient field step, do the same as in the calculate div step: Replace the neighbour values by the value of the cell itself, if the neighbour is across the boundary.