<p>First of all, I would like to state that handling the boundary conditions correctly is challenging, both from mathematical perspective and from an implementation perspective. Additionally, I think most of the papers don't explain the handling of boundary conditions sufficiently. </p>

<p>Furthermore, it is worth noticing that you have to apply boundary conditions to both the velocity field and the pressure field (and they are linked to each other). Additionally, you have to distinguish between is the pressure solve (or projection step) and the diffusion step, which must be treated separately. </p>

<p>In the very beginning you need to decide what kind of physical behavior you want to have at the boundaries. The standard example is a square surrounded by solid walls where no fluid can enter or leave. In this case you constrain the normal component of the velocity field to be zero ( a dirichlet boundary condition for the diffusion step). As we don't want changes for the normal component of the velocity in the projection step,  we also need to constrain the pressure. If you check the pressure update routine (which makes the velocity field divergency free) you can see that the derivative of the pressure field in direction of the normal must be zero, I.e. a Neumann boundary condition for the pressure field. </p>

<p>Let's assume you are talking about the boundary conditions in the pressure solve/ projection step. Stam states that for pressure there is a Neumann boundary condition. As written above this condition implies that the velocity won't be modified at the last stage of the projection step (the pressure update). In order to achieve the dirichlet boundaries for the velocity as bridson states it, you require to have homogenous Neumann conditions. So for pressure, both authors have equivalent boundary conditions. </p>

<p>With these considerations you need to set up your system of linear equations. For this it is worth to note the equations for some cells (I.e. one row of the system of linear equations) and check how the matrix and the right hand side must be modified to account for the boundary conditions. In a 2d scenario we have a 5 point Laplace stencil that results in a 4 on the diagonal and -1 for each of the adjacent cells (times some constants). In cells that are adjacent to the boundary there is at least one value that is outside of the simulation domain and care has to be taken what we do with this value. For homogenous dirichlet conditions this value is simply zero and no changes are required for the matrix and the right hand side. For Neumann boundaries the constraint is that there is no change across the boundary,  I.e. the value should be identical to the one of the cell. This reduces the value on the diagonal of the matrix by -1, I.e. 3 for a cell with one Neumann boundary. These considerations result in the linear system that you have to solve. </p>

<p>The set bnd method (that doesn't appear in the paper by the way) is a trick to make the computations without the necessity to modify the matrix. This trick only works with gauss seidel or Jacobi iterations (As commented on your other question). The set bnd method is part of the code that stam published. Personally I find it difficult to read and understand because of the variable names and the "b" that changes the behavior significantly and that seems to confuse you as well (see your other question).</p>

<p>Edit:
At first glance it looks like stam and bridson apply different boundary conditions, but they are actually the same. If we assume a solid body that is not moving and the velocity is zero before the pressure correction then the difference between them is zero. The term in parentheses in (4.23) will then be evaluated to -p_{i,j} which is exactly a  homogeneous Neumann boundary condition for the pressure. I've got the second edition of bridsons book,  where  this is explained in equation 5.3 and at the bottom of the page in more detail. </p>

<p>Edit2: up to know we discussed the Neumann boundary conditions for the pressure. No slip/slip boundary conditions are applied to the velocity and for both (with non-moving walls) you need homogenized Neumann boundary conditions for the pressure. With stam's approach of discretization, I.e. storing all velocity components at the cell center it is easier to realize no slip boundary conditions as the boundary for all velocity components are at cell faces. So to realize no slip boundary conditions in the diffusion step (mathematically you can't realize no slip conditions with the Euler equations ( no diffusion)) you need to use homogeneous dirichlet boundary conditions for all velocity components.</p>

<p>Edit3: It is not straight forward to use the set_bnd method on a staggered grid for no slip boundary conditions. The storage for the tangential (to the boundary) components is not ideal. </p>

<p>Edit 4: A Neumann boundary condition for the pressure implies that there are no changes for the (normal) velocity at this boundary. If you subtract the pressure gradient at such a boundary, the pressure gradient is simply zero due to the boundary condition. </p>

<p>Edit 5: Slip / no slip boundary conditions (for the velocity) are not(!) realized in the pressure correction step. For non-moving solid walls there is always a homogeneous Neumann boundary condition for the pressure no matter if you impose a slip or no-slip boundary condition. </p>

<p>Edit 6: I don't "have" facebook ;)</p>

<p>There are more subtleties like singular systems for Neumann- only boundary conditions, central differences in staggered or standard grids, distance to the boundary (at cell border or at cell center outside three domain) etc., that I could address if you detail your question. </p>

CC BY-SA 3.0
2017-06-05T01:30:43.107Z
Why can't I use the same technique in the PCG Method? I mean the set_bnd
Anas Alaa
CC BY-SA 3.0
2017-06-05T01:48:09.417Z
I edited my answer w.r.t. the equivalence between stam and bridson
dweber
CC BY-SA 3.0
2017-06-05T01:49:28.723Z
... If you have a more detailed question, please add page numbers / equation numbers, so I can find the right spots.
dweber
CC BY-SA 3.0
2017-06-05T01:52:59.430Z
Did you try the set bnd method for pcg? I think it is unclear to which variables it needs to be applied to. You can give it a try and report about your results.
dweber
CC BY-SA 3.0
2017-06-05T01:58:31.000Z
Firstly, thank you for your detailed answer. Stam applies a no slip boundary condition for the velocity and Bridson applies slip boundary condition. I don't have a problem with the velocity right now. I am sorry I didn't put the code because i thought it's in the paper. The problem is at the pressure boundary condition which is indeed a neuman (the normal derivative of p is just zero) and it's applied in the projection step but stam applies it in a different way from Bridson. If you see the code of stam (I will modify the question to add it),  TBC in the next comment
Anas Alaa
CC BY-SA 3.0
2017-06-05T01:59:46.137Z
If you see Stam code, you will see that he applied the boundaries as mentioned here in GPU GEMS Equation 17 and 18.    http://developer.download.nvidia.com/books/HTML/gpugems/gpugems_ch38.html.    But Bridson derives a formula to calculate the pressure of the boundary cells dependant on the velocity which is in the SIGGRAPH paper i mentioned above in chapter 4 equation number 4.10
Anas Alaa
CC BY-SA 3.0
2017-06-05T02:05:46.487Z
Actually stam doesn't imply the boundaries in the matrix of the pcg like the math. He just finds the values of the pressure and the velocity that applies the boundaries and put them in the matrix. Actually I see not the code of the no slip of Stam which I will try to understand (the set_bnd and how it works according to the paper of cpu gems) but I didn't find a way for the free slip or no stick condition in which the fluid moves freely in the tangential direction. Thanks for your answers again and again. I would like to chat with you at any time.
Anas Alaa
CC BY-SA 3.0
2017-06-05T02:19:41.760Z
Sorry, Bridson changes the equations so that they fit the boundary conditions for the pressure and the velocity. Now i understand this. But how can we implement the free slip boundaries?
Anas Alaa
CC BY-SA 3.0
2017-06-05T10:25:14.187Z
The gpu gems code is implementing the boundary conditions the same way as stam but again "only" for Jacobi iterations,  not for a cg solver.
dweber
CC BY-SA 3.0
2017-06-05T14:00:47.980Z
Now I understand how stam realizes the velocity boundary condition on an allocation grid. Can we use his implementation in a staggered grid? I still don't understand the relation between forcing the normal derivative to be zero in the way stam is doing and the equation that Bridson derives of the pressure at boundaries. Bridson says that when we subtract the pressure gradient with the velocity of the slip boundary conditions it will just enforce the pressure boundary but he didn't explain why. One time I asked someone and he told me that it's prescribed in the momentum and I didn't understand
Anas Alaa
CC BY-SA 3.0
2017-06-05T14:01:50.823Z
Also there is no code implementing the slip boundaries of velocity like the code of stam is implementing the no slip. I added this in an answer by mistake lol
Anas Alaa
CC BY-SA 3.0
2017-06-05T15:05:28.133Z
Another thing please, Stam says that the normal derivative of pressure should be zero in stick condition but Bridson says a difference thing in equ in the secondary edition book you told me about. Is the normal derivative of pressure different in the two cases? I would like to add you on Facebook or something
Anas Alaa
CC BY-SA 3.0
2017-06-06T02:48:18.047Z
I still don't understand everything. Why can't I use set_bnd in staggered grids? What's the tangential component and how can I code some function like set_bnd to enforce the slip boundary? (yes or no) Does the velocity boundary condition whether it's slip or no slip affects the normal derivative of pressure? I mean that if slip, the normal derivative will be like Bridson stated and in no slip, it will be just zero?
Anas Alaa
CC BY-SA 3.0
2017-06-06T02:51:40.587Z
Where can I read more about this? I don't find any paper say something clear about what is a boundary condition!!!
Anas Alaa
CC BY-SA 3.0
2017-06-06T02:54:52.727Z
Why stam uses the normal derivative equals zero and Bridson uses the normal derivative of pressure to be different in page 83 in his book
Anas Alaa
CC BY-SA 3.0
2017-06-06T02:57:49.383Z
I really want to contact you in private even Skype chat or anything you want please. I want to hit a wall with my confused head.
Anas Alaa
CC BY-SA 3.0
2017-06-04T21:12:34.727Z
Really thanks for this detailed answer. Can you read Bridson paper who has a different way of boundaries from Stam's one and I can't relate them. Another thing, can we just evaluate the values of the velocity and pressure at the boundaries then put the values directly in the matrix?
Anas Alaa
CC BY-SA 3.0
2017-06-04T20:32:50.323Z
2017-06-05T20:29:23.180Z
2017-06-05T20:29:23.180Z
dweber
<p>The answers of my questions is that Stam uses set_bnd to implement both the pressure boundary condition and the velocity. Bridson uses slip boundaries and he derives an equation for the pressure at boundaries. When you subtract the pressure gradient in the chorin's projection method he uses, the boundary condition of the pressure is forced. </p>

CC BY-SA 3.0
2017-06-05T14:00:28.983Z
2017-06-05T14:00:28.983Z
Anas Alaa
<p>I read <a href="https://www.autodeskresearch.com/pdf/ns.pdf" rel="nofollow noreferrer">Stable Fluids</a> and <a href="https://www.autodeskresearch.com/sites/default/files/TheArtOfFluidAnimationCode.zip" rel="nofollow noreferrer">the code</a> by Jos Stam and <a href="http://www.cs.uu.nl/docs/vakken/mgp/2013-2014/literature/Fluid%20Simulation.pdf" rel="nofollow noreferrer">FLUID SIMULATION SIGGRAPH 2007 Course Notes</a> by Robert Bridson. The boundary conditions in the two papers are different in my sight. I don't understand them. How are these two ways related to solving a linear system of PDEs with boundary conditions? </p>

<p>In Stam's paper, he uses a function called <code>set_bnd</code> to set boundaries, so is this function just setting the boundaries that I need in my matrix of linear systems of Pressure Poisson and Diffuse Equation? In Numerical Analysis, we don't modify the matrix with the boundary values, but we put the boundary condition either it's a Neumann or Dirichlet. I don't really understand what kind of boundaries is this <code>set_bnd</code> function. My thinking now is that it sets the boundary values which will be used to solve linear systems but what kind of boundaries?</p>

<p>In Bridson's paper, I don't find any implementation about boundaries except for speaking about the normal components of velocity. Actually, I can't find any relation between this and Stam's way. In the section on making fluid incompressible, he says something about boundaries but I don't understand. He just satisfies boundary conditions on the boundary cells then he solves the system or he put the boundaries in the matrix?
Actually I don't know how to satisfy <code>u*n = u_solid * n</code> which I think is the no stick boundaries for interacting between a fluid and solid boundaries.</p>

<fluid-sim>
Rigorous and simple explanation for satisfying boundary conditions in Stam and Bridson papers
Anas Alaa
# Rigorous and simple explanation for satisfying boundary conditions in Stam and Bridson papers

(Asked by: username_0 on 2017-05-31)

I read [Stable Fluids](https://www.autodeskresearch.com/pdf/ns.pdf) and [the code](https://www.autodeskresearch.com/sites/default/files/TheArtOfFluidAnimationCode.zip) by Jos Stam and [FLUID SIMULATION SIGGRAPH 2007 Course Notes](http://www.cs.uu.nl/docs/vakken/mgp/2013-2014/literature/Fluid%20Simulation.pdf) by Robert Bridson. The boundary conditions in the two papers are different in my sight. I don't understand them. How are these two ways related to solving a linear system of PDEs with boundary conditions?

In Stam's paper, he uses a function called `set_bnd` to set boundaries, so is this function just setting the boundaries that I need in my matrix of linear systems of Pressure Poisson and Diffuse Equation? In Numerical Analysis, we don't modify the matrix with the boundary values, but we put the boundary condition either it's a Neumann or Dirichlet. I don't really understand what kind of boundaries is this `set_bnd` function. My thinking now is that it sets the boundary values which will be used to solve linear systems but what kind of boundaries?

In Bridson's paper, I don't find any implementation about boundaries except for speaking about the normal components of velocity. Actually, I can't find any relation between this and Stam's way. In the section on making fluid incompressible, he says something about boundaries but I don't understand. He just satisfies boundary conditions on the boundary cells then he solves the system or he put the boundaries in the matrix?
Actually I don't know how to satisfy `u*n = u_solid * n` which I think is the no stick boundaries for interacting between a fluid and solid boundaries.

**Tags:** fluid-sim


---

## 2 answers

### Answer by username_1 on 2017-06-04. Score: 3 | Accepted answer

First of all, I would like to state that handling the boundary conditions correctly is challenging, both from mathematical perspective and from an implementation perspective. Additionally, I think most of the papers don't explain the handling of boundary conditions sufficiently.

Furthermore, it is worth noticing that you have to apply boundary conditions to both the velocity field and the pressure field (and they are linked to each other). Additionally, you have to distinguish between is the pressure solve (or projection step) and the diffusion step, which must be treated separately.

In the very beginning you need to decide what kind of physical behavior you want to have at the boundaries. The standard example is a square surrounded by solid walls where no fluid can enter or leave. In this case you constrain the normal component of the velocity field to be zero ( a dirichlet boundary condition for the diffusion step). As we don't want changes for the normal component of the velocity in the projection step, we also need to constrain the pressure. If you check the pressure update routine (which makes the velocity field divergency free) you can see that the derivative of the pressure field in direction of the normal must be zero, I.e. a Neumann boundary condition for the pressure field.

Let's assume you are talking about the boundary conditions in the pressure solve/ projection step. Stam states that for pressure there is a Neumann boundary condition. As written above this condition implies that the velocity won't be modified at the last stage of the projection step (the pressure update). In order to achieve the dirichlet boundaries for the velocity as bridson states it, you require to have homogenous Neumann conditions. So for pressure, both authors have equivalent boundary conditions.

With these considerations you need to set up your system of linear equations. For this it is worth to note the equations for some cells (I.e. one row of the system of linear equations) and check how the matrix and the right hand side must be modified to account for the boundary conditions. In a 2d scenario we have a 5 point Laplace stencil that results in a 4 on the diagonal and -1 for each of the adjacent cells (times some constants). In cells that are adjacent to the boundary there is at least one value that is outside of the simulation domain and care has to be taken what we do with this value. For homogenous dirichlet conditions this value is simply zero and no changes are required for the matrix and the right hand side. For Neumann boundaries the constraint is that there is no change across the boundary, I.e. the value should be identical to the one of the cell. This reduces the value on the diagonal of the matrix by -1, I.e. 3 for a cell with one Neumann boundary. These considerations result in the linear system that you have to solve.

The set bnd method (that doesn't appear in the paper by the way) is a trick to make the computations without the necessity to modify the matrix. This trick only works with gauss seidel or Jacobi iterations (As commented on your other question). The set bnd method is part of the code that stam published. Personally I find it difficult to read and understand because of the variable names and the "b" that changes the behavior significantly and that seems to confuse you as well (see your other question).

Edit:
At first glance it looks like stam and bridson apply different boundary conditions, but they are actually the same. If we assume a solid body that is not moving and the velocity is zero before the pressure correction then the difference between them is zero. The term in parentheses in (4.23) will then be evaluated to -p\_{i,j} which is exactly a homogeneous Neumann boundary condition for the pressure. I've got the second edition of bridsons book, where this is explained in equation 5.3 and at the bottom of the page in more detail.

Edit2: up to know we discussed the Neumann boundary conditions for the pressure. No slip/slip boundary conditions are applied to the velocity and for both (with non-moving walls) you need homogenized Neumann boundary conditions for the pressure. With stam's approach of discretization, I.e. storing all velocity components at the cell center it is easier to realize no slip boundary conditions as the boundary for all velocity components are at cell faces. So to realize no slip boundary conditions in the diffusion step (mathematically you can't realize no slip conditions with the Euler equations ( no diffusion)) you need to use homogeneous dirichlet boundary conditions for all velocity components.

Edit3: It is not straight forward to use the set\_bnd method on a staggered grid for no slip boundary conditions. The storage for the tangential (to the boundary) components is not ideal.

Edit 4: A Neumann boundary condition for the pressure implies that there are no changes for the (normal) velocity at this boundary. If you subtract the pressure gradient at such a boundary, the pressure gradient is simply zero due to the boundary condition.

Edit 5: Slip / no slip boundary conditions (for the velocity) are not(!) realized in the pressure correction step. For non-moving solid walls there is always a homogeneous Neumann boundary condition for the pressure no matter if you impose a slip or no-slip boundary condition.

Edit 6: I don't "have" facebook ;)

There are more subtleties like singular systems for Neumann- only boundary conditions, central differences in staggered or standard grids, distance to the boundary (at cell border or at cell center outside three domain) etc., that I could address if you detail your question.

> **Comments:**
> **username_0**: Really thanks for this detailed answer. Can you read Bridson paper who has a different way of boundaries from Stam's one and I can't relate them. Another thing, can we just evaluate the values of the velocity and pressure at the boundaries then put the values directly in the matrix?  
> **username_0**: Why can't I use the same technique in the PCG Method? I mean the set_bnd  
> **username_1**: I edited my answer w.r.t. the equivalence between stam and bridson  
> **username_1**: ... If you have a more detailed question, please add page numbers / equation numbers, so I can find the right spots.  
> **username_1**: Did you try the set bnd method for pcg? I think it is unclear to which variables it needs to be applied to. You can give it a try and report about your results.  
> **username_0**: Firstly, thank you for your detailed answer. Stam applies a no slip boundary condition for the velocity and Bridson applies slip boundary condition. I don't have a problem with the velocity right now. I am sorry I didn't put the code because i thought it's in the paper. The problem is at the pressure boundary condition which is indeed a neuman (the normal derivative of p is just zero) and it's applied in the projection step but stam applies it in a different way from Bridson. If you see the code of stam (I will modify the question to add it),  TBC in the next comment  
> **username_0**: If you see Stam code, you will see that he applied the boundaries as mentioned here in GPU GEMS Equation 17 and 18.    http://developer.download.nvidia.com/books/HTML/gpugems/gpugems_ch38.html.    But Bridson derives a formula to calculate the pressure of the boundary cells dependant on the velocity which is in the SIGGRAPH paper i mentioned above in chapter 4 equation number 4.10  
> **username_0**: Actually stam doesn't imply the boundaries in the matrix of the pcg like the math. He just finds the values of the pressure and the velocity that applies the boundaries and put them in the matrix. Actually I see not the code of the no slip of Stam which I will try to understand (the set_bnd and how it works according to the paper of cpu gems) but I didn't find a way for the free slip or no stick condition in which the fluid moves freely in the tangential direction. Thanks for your answers again and again. I would like to chat with you at any time.  
> **username_0**: Sorry, Bridson changes the equations so that they fit the boundary conditions for the pressure and the velocity. Now i understand this. But how can we implement the free slip boundaries?  
> **username_1**: The gpu gems code is implementing the boundary conditions the same way as stam but again "only" for Jacobi iterations,  not for a cg solver.  
> **username_0**: Now I understand how stam realizes the velocity boundary condition on an allocation grid. Can we use his implementation in a staggered grid? I still don't understand the relation between forcing the normal derivative to be zero in the way stam is doing and the equation that Bridson derives of the pressure at boundaries. Bridson says that when we subtract the pressure gradient with the velocity of the slip boundary conditions it will just enforce the pressure boundary but he didn't explain why. One time I asked someone and he told me that it's prescribed in the momentum and I didn't understand  
> **username_0**: Also there is no code implementing the slip boundaries of velocity like the code of stam is implementing the no slip. I added this in an answer by mistake lol  
> **username_0**: Another thing please, Stam says that the normal derivative of pressure should be zero in stick condition but Bridson says a difference thing in equ in the secondary edition book you told me about. Is the normal derivative of pressure different in the two cases? I would like to add you on Facebook or something  
> **username_0**: I still don't understand everything. Why can't I use set_bnd in staggered grids? What's the tangential component and how can I code some function like set_bnd to enforce the slip boundary? (yes or no) Does the velocity boundary condition whether it's slip or no slip affects the normal derivative of pressure? I mean that if slip, the normal derivative will be like Bridson stated and in no slip, it will be just zero?  
> **username_0**: Where can I read more about this? I don't find any paper say something clear about what is a boundary condition!!!  
> **username_0**: Why stam uses the normal derivative equals zero and Bridson uses the normal derivative of pressure to be different in page 83 in his book  
> **username_0**: I really want to contact you in private even Skype chat or anything you want please. I want to hit a wall with my confused head.  

---

### Answer by username_0 on 2017-06-05. Score: 0

The answers of my questions is that Stam uses set\_bnd to implement both the pressure boundary condition and the velocity. Bridson uses slip boundaries and he derives an equation for the pressure at boundaries. When you subtract the pressure gradient in the chorin's projection method he uses, the boundary condition of the pressure is forced.


---

