8
$\begingroup$

I am trying to implement a 2D version of Foster and Fedkiw's paper, "Practical Animation of Liquids" here: http://physbam.stanford.edu/~fedkiw/papers/stanford2001-02.pdf

Mostly everything works, except for section 8: "Conservation of Mass." There, we set up a matrix of equations to compute the pressures needed to make the liquid divergent free.

I believe my code matches the paper, however I am getting an unsolvable matrix during the conservation of mass step.

Here are my steps for generating the matrix A:

  1. Set the diagonal entries $A_{i,i}$ to the negative of number of adjacent liquid cells to cell i.
  2. Set the entries $A_{i,j}$ and $A_{j,i}$ to 1 if both cells i and j have liquid.

Note that, in my implementation, cell $i$,$j$ in the liquid grid corresponds to row $i + $gridWidth$ * j$ in the matrix.

The paper mentions, "Static object and empty cells don’t disrupt this structure. In that case pressure and velocity terms can disappear from both sides", so I delete the columns and rows for cells that have no liquid.

So my question is: Why is my matrix singular? Am I missing some kind of boundary condition in some other place in the paper? Is it the fact that my implementation is 2D?

Here is an example matrix from my implementation for a 2x2 grid where the cell at 0,0 has no liquid:

-1   0   1

 0  -1   1

 1   1  -2

Edit

My research has led me to believe that I'm not properly handling the boundary conditions.

First of all, at this point I can say that my matrix represents the discrete pressure Poisson equation. It is the discrete analog of applying the Laplacian operator coupling local pressure changes to cell divergence.

As far as I can understand, since we're dealing with pressure differences, boundary conditions are needed to "anchor" the pressures to an absolute reference value. Otherwise there may be an infinite number of solutions to the set of equations.

In these notes, 3 different ways are given to apply boundary conditions, to the best of my understanding:

  1. Dirichlet - specifies absolute values at the boundaries.

  2. Neummann - specifies the derivative at the boundaries.

  3. Robin - specifies some kind of linear combination of the absolute value and the derivative at the boundaries.

Foster and Fedki's paper does not mention any of these, but I believe that they enforce Dirichlet boundary conditions, notable because of this statement at the end of 7.1.2, "The pressure in a surface cell is set to atmospheric pressure."

I've read the notes I linked a few times and still don't quite understand the math going on. How exactly do we enforce these boundary conditions? Looking at other implementations, there seems to be some kind of notion of a "Ghost" cells that lie at the boundary.

Here I've linked to a few sources that may be helpful to others reading this.

Notes on boundary conditions for Poisson Matrices

Computational Science StackExchange post on Neumann boundary conditions

Computational Science StackExchange post on Poisson Solver

Water Physbam Implementation


Here is the code I use to generate the matrix. Note that, instead of explicitly deleting columns and rows, I generate and use a map from liquid cell indices to the final matrix columns/rows.

for (int i = 0; i < cells.length; i++) {
  for (int j = 0; j < cells[i].length; j++) {
    FluidGridCell cell = cells[i][j];

    if (!cell.hasLiquid)
      continue;

    // get indices for the grid and matrix
    int gridIndex = i + cells.length * j;
    int matrixIndex = gridIndexToMatrixIndex.get((Integer)gridIndex);

    // count the number of adjacent liquid cells
    int adjacentLiquidCellCount = 0;
    if (i != 0) {
      if (cells[i-1][j].hasLiquid)
        adjacentLiquidCellCount++;
    }
    if (i != cells.length-1) {
      if (cells[i+1][j].hasLiquid)
        adjacentLiquidCellCount++;
    }
    if (j != 0) {
      if (cells[i][j-1].hasLiquid)
      adjacentLiquidCellCount++;
    }
    if (j != cells[0].length-1) {
      if (cells[i][j+1].hasLiquid)
        adjacentLiquidCellCount++;
    }

    // the diagonal entries are the negative count of liquid cells
    liquidMatrix.setEntry(matrixIndex, // column
                          matrixIndex, // row
                          -adjacentLiquidCellCount); // value

    // set off-diagonal values of the pressure matrix
    if (cell.hasLiquid) {
      if (i != 0) {
        if (cells[i-1][j].hasLiquid) {
          int adjacentGridIndex = (i-1) + j * cells.length;
          int adjacentMatrixIndex = gridIndexToMatrixIndex.get((Integer)adjacentGridIndex);
          liquidMatrix.setEntry(matrixIndex, // column
                                adjacentMatrixIndex, // row
                                1.0); // value
          liquidMatrix.setEntry(adjacentMatrixIndex, // column
                                matrixIndex, // row
                                1.0); // value
        }
      }
      if (i != cells.length-1) {
        if (cells[i+1][j].hasLiquid) {
          int adjacentGridIndex = (i+1) + j * cells.length;
          int adjacentMatrixIndex = gridIndexToMatrixIndex.get((Integer)adjacentGridIndex);
          liquidMatrix.setEntry(matrixIndex, // column
                                adjacentMatrixIndex, // row
                                1.0); // value
          liquidMatrix.setEntry(adjacentMatrixIndex, // column
                                matrixIndex, // row
                                1.0); // value
        }
      }
      if (j != 0) {
        if (cells[i][j-1].hasLiquid) {
          int adjacentGridIndex = i + (j-1) * cells.length;
          int adjacentMatrixIndex = gridIndexToMatrixIndex.get((Integer)adjacentGridIndex);
          liquidMatrix.setEntry(matrixIndex, // column
                                adjacentMatrixIndex, // row
                                1.0); // value
          liquidMatrix.setEntry(adjacentMatrixIndex, // column
                                matrixIndex, // row
                                1.0); // value
        }
      }
      if (j != cells[0].length-1) {
        if (cells[i][j+1].hasLiquid) {
          int adjacentGridIndex = i + (j+1) * cells.length;
          int adjacentMatrixIndex = gridIndexToMatrixIndex.get((Integer)adjacentGridIndex);
          liquidMatrix.setEntry(matrixIndex, // column
                                adjacentMatrixIndex, // row
                                1.0); // value
          liquidMatrix.setEntry(adjacentMatrixIndex, // column
                                matrixIndex, // row
                                1.0); // value
        }
      }
    }
$\endgroup$
  • $\begingroup$ It isn't clear why static objects and empty cells should allow the deletion of rows and columns. Are you setting these rows and columns to zero or removing them altogether to give a smaller matrix? $\endgroup$ – trichoplax Mar 3 '16 at 17:13
  • $\begingroup$ In case the problem is somewhere other than where you guess, it would help to see the code, if this is something you are happy to share. Ideally an MCVE $\endgroup$ – trichoplax Mar 3 '16 at 17:14
  • $\begingroup$ Hey trichoplax. A matrix with an all zero row or column would be singular, as far as I know, so I instead remove them from the matrix to make a smaller matrix (as well as their corresponding entries in the b vector). $\endgroup$ – Jared Counts Mar 3 '16 at 17:57
  • $\begingroup$ I will edit an MCVE in tonight when I am near my computer with the source. $\endgroup$ – Jared Counts Mar 3 '16 at 17:58
  • $\begingroup$ I also suspected that I was maybe making a wrong assumption somewhere else in the code, however this only pertains to the matrix structure (and whether or not it's singular). The only thing I can think of is what qualifies as a "surface cell" vs an air cell or a liquid cell. If this is a liquid cell adjacent to an air cell, is there something different that I should be doing with its corresponding columns/rows? $\endgroup$ – Jared Counts Mar 3 '16 at 18:00

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Browse other questions tagged or ask your own question.