Conserving Mass in Liquid Simulation

Question

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

Answer 1

From your code snippet and your result for 2x2 example, I can see that you are actually simulating a domain with only Neumann boundary conditions (slip wall). In this case, the system contains a null space and your matrix is singular.

If this is the simulation configuration you want (i.e. no Dirichlet(pressure) BC), you will need to project out the null space from your solution. This is straightforward if you are using the conjugate gradient (CG) as suggested in that paper. In each iteration of your CG iteration, just take the current solution vector $\mathbf{x}$, and do $$ \begin{align} \mathbf{x}' &= (\mathbf{I}-\mathbf{\hat{u}}\mathbf{\hat{u}}^T)\mathbf{x}\\ &= \mathbf{x} - (\mathbf{\hat{u}} \cdot \mathbf{x} )\mathbf{\hat{u}} \end{align} $$ where $\mathbf{\hat{u}}$ is the normalized null space of gradient operator: $\mathbf{u} = (1,1,\dots,1)$, $\mathbf{\hat{u}}=\frac{\mathbf{u}}{\lVert \mathbf{u} \rVert}$.

Otherwise if you do want to simulate air (free boundary or Dirichlet BC), you will need to distinguish a wall and an air cell (i.e. having a boolean hasLiquid is not enough), and apply correct discretization for them (see below).

As a final note, your diagonal entries are negative. You may want to flip the signs so that the CG method works.

---

Below I would like to show more details. Consider the pressure projection process. Denote the velocity before the pressure projection as $\mathbf{v}^*$. It could be divergent, so we compute pressure to correct it and get the divergence free velocity $\mathbf{v}^{n+1}$. That is, $$ \begin{align} \mathbf{v}^{n+1} &= \mathbf{v}^* - \frac{\Delta t}{\rho} \nabla P \end{align} $$ Take divergence of it, and since $\mathbf{v}^{n+1}$ is divergence free, $$ \begin{align} \nabla \cdot \nabla P &= \nabla \cdot \mathbf{v}^* \end{align} $$ Suppose that there is no pressure Dirichlet BC presesnt and we have one solution $P_0$, then $P_0+c$ for any constant $c$ is also a solution because $\nabla \cdot \nabla (P_0+c) = \nabla \cdot \nabla P_0 = \nabla \cdot \mathbf{v}^*$. $c$ is the null space we want to project out.

To handle the Dirichlet BC, lets consider the 1D case as an example. Assume you are using staggered grid, where pressures $p_i$ are located at grid centers and velocities $v_{i+1/2}$ are located at the faces between nodes $i$ and $i+1$. Then the general discretization for one cell is $$ \frac{p_{i+1}-p_i-(p_i-p_{i-1})}{\Delta x^2} = \text{rhs} $$ Suppose that $p_{i+1}$ is an air cell, that is, its pressure has been specified, then the $p_{i+1}$ term is moved to the right hand side and disappears from the matrix. Note that the count of diagonal term $p_i$ is still two. That is why I said your 2x2 example did not contain a Dirichlet BC.

With either Dirichlet or Neumann BC the matrix is always symmetric positive definite. That is why the authors said

Static object and empty cells don’t disrupt this structure.
In that case pressure and velocity terms can disappear from both sides