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:
- Set the diagonal entries $A_{i,i}$ to the negative of number of adjacent liquid cells to cell i.
- 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:
Dirichlet - specifies absolute values at the boundaries.
Neummann - specifies the derivative at the boundaries.
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
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
}
}
}