CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Main CFD Forum

Poisson equation,Neumann BCs

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree3Likes

Reply
 
LinkBack Thread Tools Display Modes
Old   May 1, 2011, 23:42
Default Poisson equation,Neumann BCs
  #1
New Member
 
vivien
Join Date: Jan 2011
Posts: 8
Rep Power: 6
zhangweisnoopy is on a distinguished road
Hi~all
I need to solve a pressure Poisson equation with only Neumann boundaries with F.D. method. Unfortunately, it leads to the sparse linear system equation Ax = b, where A is singular (because the Neumann BC in all boundaries).And I use matlab ....Any advise with possibility to solve singular linear systems? Thank you~
zhangweisnoopy is offline   Reply With Quote

Old   May 2, 2011, 02:35
Default
  #2
Senior Member
 
Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 369
Rep Power: 10
arjun is on a distinguished road
Quote:
Originally Posted by zhangweisnoopy View Post
Hi~all
I need to solve a pressure Poisson equation with only Neumann boundaries with F.D. method. Unfortunately, it leads to the sparse linear system equation Ax = b, where A is singular (because the Neumann BC in all boundaries).And I use matlab ....Any advise with possibility to solve singular linear systems? Thank you~

pick any boundary cell treat it as fixed value of say 0. So all the boundary faces have neumann condition except for 1. if you do so, your equation will no longer be singular. Then use whatever matrix solver you want to solve the system.
arjun is offline   Reply With Quote

Old   May 2, 2011, 20:09
Default
  #3
Senior Member
 
Julien de Charentenay
Join Date: Jun 2009
Location: Australia
Posts: 229
Rep Power: 9
julien.decharentenay is on a distinguished road
Send a message via Skype™ to julien.decharentenay
I think that you can fix any point at the boundary (as mentioned by arjun) or anywhere within the computational domain.
__________________
---
Julien de Charentenay
julien.decharentenay is offline   Reply With Quote

Old   May 3, 2011, 10:29
Default
  #4
New Member
 
vivien
Join Date: Jan 2011
Posts: 8
Rep Power: 6
zhangweisnoopy is on a distinguished road
Hi
I've tried it ,but it seems that it's still hard to solve the equations.
In fact my problem is ∂p/∂x+∂p/∂y= ∂u(x,y)/∂x+∂v(x,y)/∂y
u(x,y) and v(x,y) are known in all points with 1 Dirichlet BC on west and 3 N BCs on the others.
(uu(x,y),vv(x,y))=( u(x,y), v(x,y) )-(∂p/∂x,∂p/∂y)
uu(x,y) and vv(x,y) have the same boudary with u(x,y) and v(x,y).
And p(x,y) is N BCs according to the paper .
But I don't know how to figure out the boundary conditions of p ...any advise?Thank you all~
zhangweisnoopy is offline   Reply With Quote

Old   May 3, 2011, 20:30
Default
  #5
Senior Member
 
Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 369
Rep Power: 10
arjun is on a distinguished road
Quote:
Originally Posted by zhangweisnoopy View Post
Hi
I've tried it ,but it seems that it's still hard to solve the equations.
In fact my problem is ∂p/∂x+∂p/∂y= ∂u(x,y)/∂x+∂v(x,y)/∂y
u(x,y) and v(x,y) are known in all points with 1 Dirichlet BC on west and 3 N BCs on the others.
(uu(x,y),vv(x,y))=( u(x,y), v(x,y) )-(∂p/∂x,∂p/∂y)
uu(x,y) and vv(x,y) have the same boudary with u(x,y) and v(x,y).
And p(x,y) is N BCs according to the paper .
But I don't know how to figure out the boundary conditions of p ...any advise?Thank you all~

If you tried and it does not work then there are many possiblities:

1. You made mistake in implementing what is advised.

2. If matlab is using direct solver to solve that matrix system then point (1) is your only possiblity.

3. If matlab is using iterative method then it shall be noted that not all iterative solvers can solve all neumann problem (it also depends on size of your problem). As the size increases the difficulty in solving increases.

4. Usually though (3) only reduces rate of convergence that means you would observe some convergence but it would be converging very slowly. When the problem size increase after some point there will be virtually no convergence.



PS: In the end if you are using finite difference and you mesh is cartesian mesh then have a look at fishpack and its routine called blktri .
(But i think it can also not handle all neumann problem, but have a look to make sure).
arjun is offline   Reply With Quote

Old   January 27, 2014, 12:00
Default Implementing pure Neumann Boundary conditions for Poisson equations with FD method
  #6
New Member
 
Patricio Cumsille
Join Date: Jan 2014
Posts: 4
Rep Power: 3
pcumsill is on a distinguished road
Quote:
Originally Posted by arjun View Post
pick any boundary cell treat it as fixed value of say 0. So all the boundary faces have neumann condition except for 1. if you do so, your equation will no longer be singular. Then use whatever matrix solver you want to solve the system.
I have a question about the implementation of what you have said: I have to fix the value of the numerical solution at only one grid point (e.g. specifying the solution at one corner) or I have to fix the values of the solution at the four grid points of the boundary cell?

Thank you very much in advance

Best regards,
Patricio Cumsille
pcumsill is offline   Reply With Quote

Old   January 27, 2014, 14:56
Default
  #7
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 1,612
Rep Power: 23
FMDenaro will become famous soon enough
Quote:
Originally Posted by pcumsill View Post
I have a question about the implementation of what you have said: I have to fix the value of the numerical solution at only one grid point (e.g. specifying the solution at one corner) or I have to fix the values of the solution at the four grid points of the boundary cell?

Thank you very much in advance

Best regards,
Patricio Cumsille

You can fix the value of the solution in anyone of the grid points.
However, this is not necessary, the system admits infinite solutions and during a linear solver, your solution automatically sets a constant reference value.
FMDenaro is offline   Reply With Quote

Old   January 27, 2014, 19:23
Default
  #8
Member
 
adrin
Join Date: Mar 2009
Posts: 42
Rep Power: 8
adrin is on a distinguished road
Not setting a Dirichlet BC for at least one point will not lead to a converged solution unless the implementation is noisy enough to remove the matrix singularity just enough (for convergence to take place).

On the other hand, setting a node value to a constant Dirichlet BC, as proposed in this thread, is a poor choice/strategy as well. Yes, it allows for the matrix to converge, but it yields a spike at, and in the immediate neighborhood of, the point the Dirichlet BC is applied. A few years ago, I'd developed (and published) a highly accurate and inexpensive solution, without a spike, in a boundary element methods setting. I've just begun searching for similarly accurate methods in a finite-difference/volume approach, but I haven't yet found an easy-to-implement and/or inexpensive strategy!

Adrin
sbaffini likes this.
adrin is offline   Reply With Quote

Old   January 27, 2014, 21:45
Default
  #9
New Member
 
Patricio Cumsille
Join Date: Jan 2014
Posts: 4
Rep Power: 3
pcumsill is on a distinguished road
Thank you very much for the responses.

Adrin: Please, could you send to me your publication? I am interested! My email is pcumsill@gmail.com

Thank you again!

Best regards,
PC
pcumsill is offline   Reply With Quote

Old   January 28, 2014, 03:54
Default
  #10
Senior Member
 
Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 369
Rep Power: 10
arjun is on a distinguished road
Quote:
Originally Posted by adrin View Post
Not setting a Dirichlet BC for at least one point will not lead to a converged solution unless the implementation is noisy enough to remove the matrix singularity just enough (for convergence to take place).

On the other hand, setting a node value to a constant Dirichlet BC, as proposed in this thread, is a poor choice/strategy as well. Yes, it allows for the matrix to converge, but it yields a spike at, and in the immediate neighborhood of, the point the Dirichlet BC is applied. A few years ago, I'd developed (and published) a highly accurate and inexpensive solution, without a spike, in a boundary element methods setting. I've just begun searching for similarly accurate methods in a finite-difference/volume approach, but I haven't yet found an easy-to-implement and/or inexpensive strategy!

Adrin
you are correct about the spikes part. But could you point to the paper you published, it would be good read. May be it could be applied to FV/FD scenario make things better.
arjun is offline   Reply With Quote

Old   January 28, 2014, 04:36
Default
  #11
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 1,612
Rep Power: 23
FMDenaro will become famous soon enough
Quote:
Originally Posted by adrin View Post
Not setting a Dirichlet BC for at least one point will not lead to a converged solution unless the implementation is noisy enough to remove the matrix singularity just enough (for convergence to take place).

On the other hand, setting a node value to a constant Dirichlet BC, as proposed in this thread, is a poor choice/strategy as well. Yes, it allows for the matrix to converge, but it yields a spike at, and in the immediate neighborhood of, the point the Dirichlet BC is applied. A few years ago, I'd developed (and published) a highly accurate and inexpensive solution, without a spike, in a boundary element methods setting. I've just begun searching for similarly accurate methods in a finite-difference/volume approach, but I haven't yet found an easy-to-implement and/or inexpensive strategy!

Adrin

well, this is not an issue of CFD but a mathematical property... The Poisson equation Div Grad (phi) = q with Neumann BC.s admits a solution (apart a constant) provided that the compatibility relation is satisfied. This does not require complicate implementation of the BC.s but only the fact that:

Int [S] d phi/dn dS = Int [S] q dS

The CFD task consists in fulfill such relation in discrete sense.

It is a common experience that fixing an arbitrary value leads to a slower convergence
FMDenaro is offline   Reply With Quote

Old   January 28, 2014, 06:03
Default
  #12
Senior Member
 
sbaffini's Avatar
 
Paolo Lampitella
Join Date: Mar 2009
Location: Italy
Posts: 531
Blog Entries: 14
Rep Power: 17
sbaffini will become famous soon enough
My experience agrees with that of adrin, for scale resolving simulations (e.g., LES/DNS) pressure should not be fixed, otherwise you cannot get the right statistics from it. Still, there should be no problem for the momentum equations
sbaffini is offline   Reply With Quote

Old   January 28, 2014, 17:55
Default
  #13
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 1,612
Rep Power: 23
FMDenaro will become famous soon enough
Let me say that the correct way to fix a value must accord with the correct BCs otherwise you get convergence towards some uncorrect solution
Try to fix the value not in the interior but by entering it by means of the Neumann BCs, it should work correctly without spike. However, as I said, the convergence rate is slower.
FMDenaro is offline   Reply With Quote

Old   January 28, 2014, 18:12
Default
  #14
Member
 
adrin
Join Date: Mar 2009
Posts: 42
Rep Power: 8
adrin is on a distinguished road
>>>it should work correctly without spike. However, as I said, the convergence rate is slower.

Well, I happen to be working on this problem these days, and I can say that a pure neumann BC without any other modifications will _not_ work! I'm using a multigrid preconditioned krylov solver (PCG), which converges in ~10 iterations (to error of order 1.E-8) for the 3D poisson problem that I've tried. With a pure neumann BC the solution does not converge to even order 1.E-3 for up to 3000 iterations, which is the upper limit I set. In contrast, setting one of the boundary nodes to zero (dirichlet BC) leads to convergence in 7 iterations! I agree that convergence doesn't necessarily mean convergence to a correct solution (but in this case it seems it does)

Adrin
adrin is offline   Reply With Quote

Old   January 28, 2014, 18:18
Default
  #15
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 1,612
Rep Power: 23
FMDenaro will become famous soon enough
Quote:
Originally Posted by adrin View Post
>>>it should work correctly without spike. However, as I said, the convergence rate is slower.

Well, I happen to be working on this problem these days, and I can say that a pure neumann BC without any other modifications will _not_ work! I'm using a multigrid preconditioned krylov solver (PCG), which converges in ~10 iterations (to error of order 1.E-8) for the 3D poisson problem that I've tried. With a pure neumann BC the solution does not converge to even order 1.E-3 for up to 3000 iterations, which is the upper limit I set. In contrast, setting one of the boundary nodes to zero (dirichlet BC) leads to convergence in 7 iterations! I agree that convergence doesn't necessarily mean convergence to a correct solution (but in this case it seems it does)

Adrin

Adrin, something in the BC implementation could be wrong... I assume you are working on the pressure problem derived from the continuity equation in which you substitute the Hodge decomposition. The same decomposition, projected along the normal direction to the boundaries provides the correct Neumann BC.s that fulfill the compatibility relation ensuring the convergence.
I work usually with Neumann BC.s and it works....
Please, check if your BC.s satisfy numerically the relation

Int [S] d phi/dn dS - Int [S] q.n dS=0
FMDenaro is offline   Reply With Quote

Old   January 28, 2014, 22:16
Default
  #16
Member
 
adrin
Join Date: Mar 2009
Posts: 42
Rep Power: 8
adrin is on a distinguished road
>> could you point to the paper you published

A. Gharakhani and A. F. Ghoniem,"BEM Solution of the 3D Internal Neumann Problem and A Regularized Formulation for the Potential Velocity Gradients," International Journal of Numerical Methods in Fluids, Vol. 24, No. 1, pp. 81-100, 1997.
adrin is offline   Reply With Quote

Old   January 29, 2014, 09:24
Default
  #17
Senior Member
 
sbaffini's Avatar
 
Paolo Lampitella
Join Date: Mar 2009
Location: Italy
Posts: 531
Blog Entries: 14
Rep Power: 17
sbaffini will become famous soon enough
I remember having a possibly identical problem with a Vortex panel method applied to closed surfaces with lift (in practice the outside potential is well defined by its behavior at infinity, but the inside one is not). However, there the problem can be easily solved by assigning the circulation value for a panel as, just like the pressure, the solution is in terms of differences of circulations among adjacent panels and not their absolute values. Still, i never had the chance to apply the method to multiple closed surfaces to see if it really doesn't matter in general...

Actually, Adrin himself helped me in solving my issue here on CFD-ONLINE (3D PANEL METHOD ).
I should read your paper to check if the matter is exactly the same.

Coming to the general issue, i guess that you both (Adrin and Filippo) are saying exactly the same thing, one of the equations should be switched to a global constraint instead of fixing it arbitrarily.

However, there might possibly be differences in the single approaches. I guess that if pressure has infinitely many solutions, the solution should simply converge to one of them just like for the pressure checkerboard case the solution converge to a specific checkerboard pattern. My experience is that point Gauss-Seidel iterations can converge to machine accuracy (not so fast actually) if the global b.c. constraint is preserved.

However, this global constraint implies a full coupling among the equations which might be difficult to solve in the context of an otherwise banded matrix.
sbaffini is offline   Reply With Quote

Old   February 4, 2014, 18:02
Default
  #18
New Member
 
Patricio Cumsille
Join Date: Jan 2014
Posts: 4
Rep Power: 3
pcumsill is on a distinguished road
Quote:
Originally Posted by FMDenaro View Post
Let me say that the correct way to fix a value must accord with the correct BCs otherwise you get convergence towards some uncorrect solution
Try to fix the value not in the interior but by entering it by means of the Neumann BCs, it should work correctly without spike. However, as I said, the convergence rate is slower.

The problem I see is that the convergence rate is really slow. I have carried out some numerical tests and I got a convergence rate even less than 1!!!

Or even worst, for certain problems (for certain data) I did not get convergence!

Thus, the technique of fixing a value at a grid point is not good in general!!!
pcumsill is offline   Reply With Quote

Old   February 6, 2014, 03:27
Default
  #19
Senior Member
 
Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 369
Rep Power: 10
arjun is on a distinguished road
Quote:
Originally Posted by pcumsill View Post
The problem I see is that the convergence rate is really slow. I have carried out some numerical tests and I got a convergence rate even less than 1!!!

Or even worst, for certain problems (for certain data) I did not get convergence!

Thus, the technique of fixing a value at a grid point is not good in general!!!

Yes, this was the main problem I faced when I tried to fix a single point for neumann. Even the additive corrective AMG was not that effective specially when case sizes were large.
If you could use, try smoothed aggregation AMG. I ended up taking out additive corrective and replacing it with smoothed aggregation that worked well.



Arjun
arjun is offline   Reply With Quote

Old   March 25, 2014, 17:01
Default Revisiting Neumann BC
  #20
Member
 
adrin
Join Date: Mar 2009
Posts: 42
Rep Power: 8
adrin is on a distinguished road
It turns out that FMDenaro is correct. So long as the Neumann BC compatibility is satisfied _explicitly_ one need not set a Dirichlet BC at one node, and one can still get a solution. In my previous experiments in the finite volume formulation I was using a simple manufactured problem for benchmarking. It turns out that although compatibility was theoretically ensured it was _not_ satisfied numerically (in discrete form). So, a simple correction led to "a" solution. That solution shifts by a constant as a function of grid resolution, but that's not an issue. I still maintain the only reason any solution may be obtained with such an approach is sufficient numerical noise to remove the linear dependence of two arbitrary equations - this approach would fail in an "infinite accuracy" solution.

So, in summary, if the _discrete_ surface integral of the fluxes is not equal to the volume integral of the Poisson source term, perform the following pre-processing before assembling the matrix and RHS. Find the difference between the aforementioned two terms, and either (a) divide that difference by the total volume and then add the volume averaged error to the source term for every node, or (b) divide the difference by the total surface area and then add the surface averaged error to the fluxes on all surface nodes. Make sure the signs are such that compatibility is now satisfied. This yields a quick, converged solution!

Adrin
FMDenaro likes this.
adrin is offline   Reply With Quote

Reply

Tags
neumann bcs, poisson equation

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Dealing with BC's in OF 1.6 vkrastev OpenFOAM Running, Solving & CFD 5 September 4, 2012 11:58
BCs for Pressure Correction Equation (SIMPLE) Bharath Somayaji Main CFD Forum 1 March 1, 2006 07:12
Recommendation of a good poisson solver Quarkz Main CFD Forum 2 December 2, 2005 10:12
Poisson Equation in CFD Maciej Matyka Main CFD Forum 9 November 10, 2004 12:30
Poisson BC's Peter Main CFD Forum 2 May 18, 2001 04:40


All times are GMT -4. The time now is 06:13.