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

<p>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 <span class="math-container">$\mathbf{x}$</span>, and do
<span class="math-container">$$
\begin{align}
\mathbf{x}' &amp;= (\mathbf{I}-\mathbf{\hat{u}}\mathbf{\hat{u}}^T)\mathbf{x}\\
&amp;= \mathbf{x} - (\mathbf{\hat{u}} \cdot \mathbf{x} )\mathbf{\hat{u}}
\end{align}
$$</span>
where <span class="math-container">$\mathbf{\hat{u}}$</span> is the normalized null space of gradient operator:
<span class="math-container">$\mathbf{u} = (1,1,\dots,1)$</span>, <span class="math-container">$\mathbf{\hat{u}}=\frac{\mathbf{u}}{\lVert \mathbf{u} \rVert}$</span>.</p>

<p>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 <code>hasLiquid</code> is not enough), and apply correct discretization for them (see below).</p>

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

<hr>

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

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

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

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

CC BY-SA 4.0
2019-02-21T03:08:51.373Z
2019-02-21T03:08:51.373Z
TheBusyTypist
<p>I am trying to implement a 2D version of Foster and Fedkiw's paper, "Practical Animation of Liquids" here: <a href="http://physbam.stanford.edu/~fedkiw/papers/stanford2001-02.pdf" rel="nofollow noreferrer">http://physbam.stanford.edu/~fedkiw/papers/stanford2001-02.pdf</a></p>

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

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

<p>Here are my steps for generating the matrix A:</p>

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

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

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

<p>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?</p>

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

<pre><code>-1   0   1

 0  -1   1

 1   1  -2
</code></pre>

<hr>

<p><strong>Edit</strong></p>

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

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

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

<p>In <a href="http://www.math.toronto.edu/mpugh/Teaching/Mat1062/notes2.pdf" rel="nofollow noreferrer">these notes</a>, 3 different ways are given to apply boundary conditions, to the best of my understanding:</p>

<ol>
<li><p>Dirichlet - specifies absolute values at the boundaries.</p></li>
<li><p>Neummann - specifies the derivative at the boundaries.</p></li>
<li><p>Robin - specifies some kind of linear combination of the absolute value and the derivative at the boundaries.</p></li>
</ol>

<p>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."</p>

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

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

<p><a href="http://www.math.toronto.edu/mpugh/Teaching/Mat1062/notes2.pdf" rel="nofollow noreferrer">Notes on boundary conditions for Poisson Matrices</a></p>

<p><a href="https://scicomp.stackexchange.com/questions/5355/writing-the-poisson-equation-finite-difference-matrix-with-neumann-boundary-cond">Computational Science StackExchange post on Neumann boundary conditions</a></p>

<p><a href="https://scicomp.stackexchange.com/questions/19485/2d-poisson-solver-for-taylor-green-vortex-problem">Computational Science StackExchange post on Poisson Solver</a></p>

<p><a href="http://physbam.stanford.edu/~mlentine/project.html#water" rel="nofollow noreferrer">Water Physbam Implementation</a></p>

<hr>

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

<pre><code>for (int i = 0; i &lt; cells.length; i++) {
  for (int j = 0; j &lt; 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
        }
      }
    }
</code></pre>

CC BY-SA 3.0
2016-03-03T17:13:24.183Z
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?
trichoplax is on Codidact now
CC BY-SA 3.0
2016-03-03T17:14:50.840Z
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](http://stackoverflow.com/help/mcve)
trichoplax is on Codidact now
CC BY-SA 3.0
2016-03-03T17:57:59.703Z
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).
Jared Counts
CC BY-SA 3.0
2016-03-03T17:58:37.703Z
I will edit an MCVE in tonight when I am near my computer with the source.
Jared Counts
CC BY-SA 3.0
2016-03-03T18:00:50.033Z
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?
Jared Counts
CC BY-SA 3.0
2016-03-03T19:49:29.360Z
I edited the question to include my code for generating the matrix. It is in Processing/Java.
Jared Counts
CC BY-SA 3.0
2016-03-04T18:02:45.977Z
I added more information to the post.
Jared Counts
<simulation><fluid-sim>
Conserving Mass in Liquid Simulation
Jared Counts
# Conserving Mass in Liquid Simulation

(Asked by: username_0 on 2016-03-02)

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](http://www.math.toronto.edu/mpugh/Teaching/Mat1062/notes2.pdf), 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](http://www.math.toronto.edu/mpugh/Teaching/Mat1062/notes2.pdf)

[Computational Science StackExchange post on Neumann boundary conditions](https://scicomp.stackexchange.com/questions/5355/writing-the-poisson-equation-finite-difference-matrix-with-neumann-boundary-cond)

[Computational Science StackExchange post on Poisson Solver](https://scicomp.stackexchange.com/questions/19485/2d-poisson-solver-for-taylor-green-vortex-problem)

[Water Physbam Implementation](http://physbam.stanford.edu/~mlentine/project.html#water)

---

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

```

**Tags:** simulation, fluid-sim

> **Comments:**
> **username_1**: 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?  
> **username_1**: 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](http://stackoverflow.com/help/mcve)  
> **username_0**: Hey username_1. 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).  
> **username_0**: I will edit an MCVE in tonight when I am near my computer with the source.  
> **username_0**: 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?  
> **username_0**: I edited the question to include my code for generating the matrix. It is in Processing/Java.  
> **username_0**: I added more information to the post.  

---

## 1 answers

### Answer by username_2 on 2019-02-21. Score: 2

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

```


---

