# Stability issues with "modified pressure-correction" scheme

 User Name Remember Me Password
 Register Blogs Members List Search Today's Posts Mark Forums Read

June 25, 2020, 16:53
#21
New Member

Join Date: Jun 2020
Posts: 27
Rep Power: 5
Quote:
 Originally Posted by Eifoehn4 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.

June 25, 2020, 16:55
#22
Senior Member

-
Join Date: Jul 2012
Location: Germany
Posts: 184
Rep Power: 13
Quote:
 Originally Posted by FMDenaro To be honest, the paper is not very clear...
That's true.
__________________
Check out my side project:

A multiphysics discontinuous Galerkin framework: Youtube, Gitlab.

June 25, 2020, 16:57
#23
New Member

Join Date: Jun 2020
Posts: 27
Rep Power: 5
Quote:
 Originally Posted by FMDenaro 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

June 25, 2020, 17:01
#24
Senior Member

Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,793
Rep Power: 71
Quote:
 Originally Posted by XorT 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.

June 25, 2020, 17:19
#25
New Member

Join Date: Jun 2020
Posts: 27
Rep Power: 5
Quote:
 Originally Posted by FMDenaro 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.

June 25, 2020, 17:25
#26
Senior Member

Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,793
Rep Power: 71
Quote:
 Originally Posted by XorT 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.

June 25, 2020, 17:33
#27
New Member

Join Date: Jun 2020
Posts: 27
Rep Power: 5
Quote:
 Originally Posted by FMDenaro 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...

June 25, 2020, 17:40
#28
Senior Member

Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,793
Rep Power: 71
Quote:
 Originally Posted by XorT 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.

June 26, 2020, 03:55
#29
New Member

Join Date: Jun 2020
Posts: 27
Rep Power: 5
Quote:
 Originally Posted by FMDenaro 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

June 26, 2020, 04:33
#30
Senior Member

Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,793
Rep Power: 71
Quote:
 Originally Posted by XorT 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.

 June 26, 2020, 05:02 #31 Senior Member   Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,793 Rep Power: 71 Maybe the key is in the source term, it is modified to include the FD laplacian contribution?

 June 26, 2020, 05:28 #32 New Member   Join Date: Jun 2020 Posts: 27 Rep Power: 5 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 ?

 June 26, 2020, 05:46 #33 Senior Member   Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,793 Rep Power: 71 I am stating that you cannot mix different discrete operatore in these steps!

June 26, 2020, 06:03
#34
Senior Member

-
Join Date: Jul 2012
Location: Germany
Posts: 184
Rep Power: 13
Quote:
 Originally Posted by FMDenaro 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 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.

June 26, 2020, 06:14
#35
Senior Member

Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,793
Rep Power: 71
Quote:
 Originally Posted by Eifoehn4 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 ...

June 26, 2020, 06:23
#36
Senior Member

-
Join Date: Jul 2012
Location: Germany
Posts: 184
Rep Power: 13
Quote:
 Originally Posted by FMDenaro 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.

June 26, 2020, 06:48
#37
New Member

Join Date: Jun 2020
Posts: 27
Rep Power: 5
Quote:
 Originally Posted by FMDenaro 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 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.

June 26, 2020, 06:53
#38
New Member

Join Date: Jun 2020
Posts: 27
Rep Power: 5
Quote:
 Originally Posted by Eifoehn4 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

June 26, 2020, 07:01
#39
Senior Member

Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,793
Rep Power: 71
Quote:
 Originally Posted by XorT 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
 1D439959-F73F-48D5-A217-4240848644C8.jpg (155.7 KB, 2 views)

June 26, 2020, 07:05
#40
New Member

Join Date: Jun 2020
Posts: 27
Rep Power: 5
Quote:
 Originally Posted by FMDenaro 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).

 Tags modified solver, pressure based solver

 Thread Tools Search this Thread Search this Thread: Advanced Search Display Modes Linear Mode

 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 OffTrackbacks are Off Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post Tobi OpenFOAM Post-Processing 9 March 25, 2022 01:33 metmet FLUENT 6 October 30, 2019 12:23 Vishnu_bharathi CFX 12 November 21, 2017 06:56 guillaume1990 FLUENT 0 March 3, 2014 11:59 Giuki FLUENT 1 July 19, 2011 11:35

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

 Contact Us - CFD Online - Privacy Statement - Top