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

Stability issues with "modified pressure-correction" scheme

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

Like Tree4Likes

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   June 25, 2020, 16:53
Default
  #21
New Member
 
Join Date: Jun 2020
Posts: 27
Rep Power: 5
XorT is on a distinguished road
Quote:
Originally Posted by Eifoehn4 View Post
I ask because:
  • the collocations points of the trigonometric Fourier basis lives on a equidistant grid
  • the Chebyshev basis and the related DCT lives on a Chebyshev grid

Can you briefly explain how you handle this?

Hmm, I'm affraid I don't know enough about it to answer this. However, what I did verify is that after applying the pressure correction, the divergence of the velocity field is in fact zero everywhere. Regardless of this, though, the velocities grow indefinetely every timestep.
XorT is offline   Reply With Quote

Old   June 25, 2020, 16:55
Default
  #22
Senior Member
 
Eifoehn4's Avatar
 
-
Join Date: Jul 2012
Location: Germany
Posts: 184
Rep Power: 13
Eifoehn4 is on a distinguished road
Quote:
Originally Posted by FMDenaro View Post
To be honest, the paper is not very clear...
That's true.
__________________
Check out my side project:

A multiphysics discontinuous Galerkin framework: Youtube, Gitlab.
Eifoehn4 is offline   Reply With Quote

Old   June 25, 2020, 16:57
Default
  #23
New Member
 
Join Date: Jun 2020
Posts: 27
Rep Power: 5
XorT is on a distinguished road
Quote:
Originally Posted by FMDenaro View Post
To be honest, the paper is not very clear and, in my opinion, some issues appear...for example, consider Eq.(39) and (40) and compute the div v^n+1 on the staggered grid. If you use the second order central discretization on both sides, the pressure contribution appears in a discretized second order laplacian for p^n+1. But the pressure solution is obtained by a spectral representation, therefore an incongruence appears. That should converte the procedure in an approximate projection rather than in an exact projection. As a consequence a non-vanishing (I mean at machine precision) divergence will appear in the velocity. That would be of some small magnitude order but can cumulate in time.

However, several other reasons for a numerical instability could be possible, a wrong time step, a bug in the code...


PS: I am not sure, you are using h=1 as computational step size?

I think you are mistaken. The constant poisson equation is solved directly and accurately. The solution exactly solves the equation. Also the velocity field has zero divergence after the correction.


PS.: yes, i assume h=1 for simplicity
XorT is offline   Reply With Quote

Old   June 25, 2020, 17:01
Default
  #24
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,793
Rep Power: 71
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
Quote:
Originally Posted by XorT View Post
I think you are mistaken. The constant poisson equation is solved directly and accurately. The solution exactly solves the equation. Also the velocity field has zero divergence after the correction.


PS.: yes, i assume h=1 for simplicity



I just considered the equations in the paper. Try yourself by hand, build the discrete div v on the staggered grid and check the result. You will see that the discrete 5-point FD Poisson equation appears. But if you use the spectral solution that is not verified exactly.
FMDenaro is offline   Reply With Quote

Old   June 25, 2020, 17:19
Default
  #25
New Member
 
Join Date: Jun 2020
Posts: 27
Rep Power: 5
XorT is on a distinguished road
Quote:
Originally Posted by FMDenaro View Post
I just considered the equations in the paper. Try yourself by hand, build the discrete div v on the staggered grid and check the result. You will see that the discrete 5-point FD Poisson equation appears. But if you use the spectral solution that is not verified exactly.

I don't know what you mean exactly. Are you saying that the use of the poisson equation is wrong ? Because the kind of solver used to solve it should be completely irrelevant. Either it solves the equation or it doesn't (in which case it is not really a solver, but merely an approximation of some kind).
I'm not getting why we're stuck on the FFT based solver topic.
XorT is offline   Reply With Quote

Old   June 25, 2020, 17:25
Default
  #26
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,793
Rep Power: 71
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
Quote:
Originally Posted by XorT View Post
I don't know what you mean exactly. Are you saying that the use of the poisson equation is wrong ? Because the kind of solver used to solve it should be completely irrelevant. Either it solves the equation or it doesn't (in which case it is not really a solver, but merely an approximation of some kind).
I'm not getting why we're stuck on the FFT based solver topic.



Maybe you are not aware of the basic principle about approximate and exact projection methods. One of the fundamental is the use of the same discrete operators in each step. Using the FFT means you are adopting a spectral interpolation for the pressure. Then the pressure gradient should be computed from the spectral polynomial not from the discret gradient shown in Eqs. (39) and (40). That is not clear to me, how can be ensured the divergence-free field up to machine precision?

Again, use (39) and (40) to build the Poisson equation. You get a discrete FD equation while the paper addressed the spectral interpolation for the continuous Poisson problem.
FMDenaro is offline   Reply With Quote

Old   June 25, 2020, 17:33
Default
  #27
New Member
 
Join Date: Jun 2020
Posts: 27
Rep Power: 5
XorT is on a distinguished road
Quote:
Originally Posted by FMDenaro View Post
Maybe you are not aware of the basic principle about approximate and exact projection methods. One of the fundamental is the use of the same discrete operators in each step. Using the FFT means you are adopting a spectral interpolation for the pressure. Then the pressure gradient should be computed from the spectral polynomial not from the discret gradient shown in Eqs. (39) and (40). That is not clear to me, how can be ensured the divergence-free field up to machine precision?

Again, use (39) and (40) to build the Poisson equation. You get a discrete FD equation while the paper addressed the spectral interpolation for the continuous Poisson problem.

Isn't the whole paper using the discretized version ? Isn't the discretized poisson equation solved exactly (by the solver of your choice). I don't see this dependent on a spectral solver at all. Where do you see this statement ? You can use whatever solver you want as long as it solves the fully discretized poisson equation. They even employ MG methods as an alternative.
I really don't get it...
XorT is offline   Reply With Quote

Old   June 25, 2020, 17:40
Default
  #28
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,793
Rep Power: 71
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
Quote:
Originally Posted by XorT View Post
Isn't the whole paper using the discretized version ? Isn't the discretized poisson equation solved exactly (by the solver of your choice). I don't see this dependent on a spectral solver at all. Where do you see this statement ? You can use whatever solver you want as long as it solves the fully discretized poisson equation. They even employ MG methods as an alternative.
I really don't get it...



The FFT is not just a solver! It implies a functional representation of the function p by means of trigonometric functions. If you start from Lap p that is a continuous representation of the operator, equivalent to Div Grad. In projection methods you have to discretize Grad and Div in a congruent way but that is not necessarily equivalent to the discrete Lap operator.

That is the basic theory of the projection methods you can find in the original Chorin's paper as well as in the textbooks of CFD.
The paper addressed (39) and (40) as if the Div Grad operator was discretized by second order central formula. But, actually, it seem that Lap p is directly solved by the FFT.
FMDenaro is offline   Reply With Quote

Old   June 26, 2020, 03:55
Default
  #29
New Member
 
Join Date: Jun 2020
Posts: 27
Rep Power: 5
XorT is on a distinguished road
Quote:
Originally Posted by FMDenaro View Post
The FFT is not just a solver! It implies a functional representation of the function p by means of trigonometric functions. If you start from Lap p that is a continuous representation of the operator, equivalent to Div Grad. In projection methods you have to discretize Grad and Div in a congruent way but that is not necessarily equivalent to the discrete Lap operator.

That is the basic theory of the projection methods you can find in the original Chorin's paper as well as in the textbooks of CFD.
The paper addressed (39) and (40) as if the Div Grad operator was discretized by second order central formula. But, actually, it seem that Lap p is directly solved by the FFT.

Well, I'm not expert on this, but as I said, I verified the results. Here is some proof for two test-cases (once I hard transition, once a random field for the RHS). Computed the pressure, computed the gradient, computed the 5-point laplacian and computed the sum of squared differences to the RHS as error (see above laplacian image). The error is in the real of machine precision I'd say, at least it is so tiny that this is not the problem here, I'm quite sure. (These are the DST versions for Dirichlet BC, btw)



https://ibb.co/LtQHKcY
https://ibb.co/HVwdCBF
XorT is offline   Reply With Quote

Old   June 26, 2020, 04:33
Default
  #30
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,793
Rep Power: 71
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
Quote:
Originally Posted by XorT View Post
Well, I'm not expert on this, but as I said, I verified the results. Here is some proof for two test-cases (once I hard transition, once a random field for the RHS). Computed the pressure, computed the gradient, computed the 5-point laplacian and computed the sum of squared differences to the RHS as error (see above laplacian image). The error is in the real of machine precision I'd say, at least it is so tiny that this is not the problem here, I'm quite sure. (These are the DST versions for Dirichlet BC, btw)



https://ibb.co/LtQHKcY
https://ibb.co/HVwdCBF

From your figures I don't see anything I am talking about. The problem is not the solution of the Laplacian but the computation of div v^n+1 in each node.

Could you post the field of the result ?
Again, this is a mathematical conclusion from Eq.(39) and (40). Take the discrete divergence of the LHS and set to zero. Then take the same divergence of the RHS and develop the discrete pressure equation. This is the FD equation that should be solved, not the spectral representation of the laplacian of the pressure. The article is not clear and I do not understand what is done.
FMDenaro is offline   Reply With Quote

Old   June 26, 2020, 05:02
Default
  #31
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,793
Rep Power: 71
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
Maybe the key is in the source term, it is modified to include the FD laplacian contribution?
FMDenaro is offline   Reply With Quote

Old   June 26, 2020, 05:28
Default
  #32
New Member
 
Join Date: Jun 2020
Posts: 27
Rep Power: 5
XorT is on a distinguished road
Well, we seem to run in circles.
Maybe, for reference, check the standard approach to pressure correction used in fluid dynamics for computer graphics here: http://users.encs.concordia.ca/~grog...cs/fluid-5.pdf (section 5.3.3). They do exactly the same thing. Compute the divergence (cell centered) from the staggered velocities (at cell boundaries), use the divergence (at cell center) to build the RHS of the poisson equation which is then solved to yield pressure (at cell center). Then the gradient of the pressure (which is at cell boundary) is applied to the staggered velocities (also at cell boundary).

Are you saying this a problem by itself already ?
XorT is offline   Reply With Quote

Old   June 26, 2020, 05:46
Default
  #33
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,793
Rep Power: 71
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
I am stating that you cannot mix different discrete operatore in these steps!
FMDenaro is offline   Reply With Quote

Old   June 26, 2020, 06:03
Default
  #34
Senior Member
 
Eifoehn4's Avatar
 
-
Join Date: Jul 2012
Location: Germany
Posts: 184
Rep Power: 13
Eifoehn4 is on a distinguished road
Quote:
Originally Posted by FMDenaro View Post
I am stating that you cannot mix different discrete operatore in these steps!
Maybe we should look at this more from the point of view of a computer-graphic scientist.

Quote:
Originally Posted by XorT View Post
Well, I'm not expert on this, but as I said, I verified the results. Here is some proof for two test-cases (once I hard transition, once a random field for the RHS). Computed the pressure, computed the gradient, computed the 5-point laplacian and computed the sum of squared differences to the RHS as error (see above laplacian image). The error is in the real of machine precision I'd say, at least it is so tiny that this is not the problem here, I'm quite sure. (These are the DST versions for Dirichlet BC, btw)

https://ibb.co/LtQHKcY
https://ibb.co/HVwdCBF
Did you also perform a convergence test, e.g. with a manufactured solutions? This would be a first step to ensure that you at least converge to the right solution no matter what order. Moreover what did you do so far to tackle the instability issue?
__________________
Check out my side project:

A multiphysics discontinuous Galerkin framework: Youtube, Gitlab.
Eifoehn4 is offline   Reply With Quote

Old   June 26, 2020, 06:14
Default
  #35
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,793
Rep Power: 71
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
Quote:
Originally Posted by Eifoehn4 View Post
Maybe we should look at this more from the point of view of a computer-graphic scientist.



Did you also perform a convergence test, e.g. with a manufactured solutions? This would be a first step to ensure that you at least converge to the right solution no matter what order.
But that jcp paper is about a CFD topic... I would assess if the finally velocity is discretely divergence-free at machine accuracy or at local truncation error accuracy.
From a point of a computer graphics I can accept also a first order stable solution but as here the problem is the instability ...
FMDenaro is offline   Reply With Quote

Old   June 26, 2020, 06:23
Default
  #36
Senior Member
 
Eifoehn4's Avatar
 
-
Join Date: Jul 2012
Location: Germany
Posts: 184
Rep Power: 13
Eifoehn4 is on a distinguished road
Quote:
Originally Posted by FMDenaro View Post
But that jcp paper is about a CFD topic... I would assess if the finally velocity is discretely divergence-free at machine accuracy or at local truncation error accuracy.
From a point of a computer graphics I can accept also a first order stable solution but as here the problem is the instability ...
Yes that's true. Such inconsistencies may also influence the stability.
__________________
Check out my side project:

A multiphysics discontinuous Galerkin framework: Youtube, Gitlab.
Eifoehn4 is offline   Reply With Quote

Old   June 26, 2020, 06:48
Default
  #37
New Member
 
Join Date: Jun 2020
Posts: 27
Rep Power: 5
XorT is on a distinguished road
Quote:
Originally Posted by FMDenaro View Post
I am stating that you cannot mix different discrete operatore in these steps!

so http://users.encs.concordia.ca/~grog...cs/fluid-5.pdf (section 5.3.3) does it wrong as well ? because this is well established standard way of doing it when it comes to CFD in CG.



Quote:
Originally Posted by FMDenaro View Post
But that jcp paper is about a CFD topic... I would assess if the finally velocity is discretely divergence-free at machine accuracy or at local truncation error accuracy.

Both papers are about a CFD topic. And as I said before, I did assess that the final velocity is discretely divergence-free at machine accuracy.
XorT is offline   Reply With Quote

Old   June 26, 2020, 06:53
Default
  #38
New Member
 
Join Date: Jun 2020
Posts: 27
Rep Power: 5
XorT is on a distinguished road
Quote:
Originally Posted by Eifoehn4 View Post
Did you also perform a convergence test, e.g. with a manufactured solutions? This would be a first step to ensure that you at least converge to the right solution no matter what order. Moreover what did you do so far to tackle the instability issue?

I'm trying to simplify the simulation and scenario to home in on the source of the error numerically. I apply no forces, use really small timesteps, but this does not prevent the blow-up UNLESS i use constant density (in which case everything works as expected). The latter also being the reason to doubt, that this discussion is going in the right direction. It is clear to me that the problem is in the "hack" of approximating the pressure for the variable density from past pressures and density gradients. Not what we're discussing about here for days now
XorT is offline   Reply With Quote

Old   June 26, 2020, 07:01
Default
  #39
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,793
Rep Power: 71
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
Quote:
Originally Posted by XorT View Post
so http://users.encs.concordia.ca/~grog...cs/fluid-5.pdf (section 5.3.3) does it wrong as well ? because this is well established standard way of doing it when it comes to CFD in CG.






Both papers are about a CFD topic. And as I said before, I did assess that the final velocity is discretely divergence-free at machine accuracy.

Look, exactly what I am stating, the pressure equation is not solved by a spectra interpolation
Attached Images
File Type: jpg 1D439959-F73F-48D5-A217-4240848644C8.jpg (155.7 KB, 2 views)
FMDenaro is offline   Reply With Quote

Old   June 26, 2020, 07:05
Default
  #40
New Member
 
Join Date: Jun 2020
Posts: 27
Rep Power: 5
XorT is on a distinguished road
Quote:
Originally Posted by FMDenaro View Post
Look, exactly what I am stating, the pressure equation is not solved by a spectra interpolation

yes, because it is variable coefficient, which the spectral solver is not able to cope with. for constant coefficient it would not matter. any solver could be used, including FFT based ones. They are direct, precise solvers. They solve the given discrete poisson equation. And as such be as good as any other solver which solves the given discrete poisson equation. The algorithm is agnostic to the solver used (as long as it solves the equation).
XorT is offline   Reply With Quote

Reply

Tags
modified solver, pressure based solver

Thread Tools Search this Thread
Search this Thread:

Advanced Search
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 Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Pressure fields in FOAM, p field, total pressure, etc. Tobi OpenFOAM Post-Processing 9 March 25, 2022 01:33
Wind tunnel Boundary Conditions in Fluent metmet FLUENT 6 October 30, 2019 12:23
outlet pressure Boundary settings -velocity streamline under ambient temp.conditions Vishnu_bharathi CFX 12 November 21, 2017 06:56
Setting up the pressure variation due to tornado in a duct(UDF)+animation guillaume1990 FLUENT 0 March 3, 2014 11:59
Pressure BC for combustion chamber Giuki FLUENT 1 July 19, 2011 11:35


All times are GMT -4. The time now is 00:21.