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 24, 2020, 16:39
Unhappy Stability issues with "modified pressure-correction" scheme
  #1
New Member
 
Join Date: Jun 2020
Posts: 27
Rep Power: 5
XorT is on a distinguished road
Hello,


I'm trying to develop some prototype for a realtime 2d fluid dynamics simulation. The goal is to be fast rather than physically exact as long as it is visually plausible (use-case is computer-graphics). However the simulation needs to deal with high density ratios as it needs to handle water-air-interface dynamics in a convincing way.

The method described in the paper "A fast pressure-correction method for incompressible two-fluid flows" (https://doddm.com/publications/2014-jcp-md-af.pdf) appears very interesting as it can employ a FFT based FastPoisson solver to solve the pressure equation, which can achieve realtime performance on high-resolution grids using GPUs.



However I experience severe stability issues (blow ups).
My simulation setting is fairly simple: it is only 2D and the initial conditions are a drop of water (1000kg/mē) in the center with initial velocity in the y-direction on a background of initially still air (1kg/mē). As quantities I keep only cell-centered densities and staggered velocities. I use the MacCormack scheme for advection.



Without advecting the density- and velocity field the simulation is stable and the solved divergence-free velocitiy field stabalizes to a nice (static) approximation (in which mostly the air gets accelerated to match the water' velocity, while the water is barely slowed by the initially still air, due to the massive density ratio).
However if i do advect, even for really small velocities (far less than a grid cell), some velocities will exponentially grow and blow up over multiple iterations. In contrary this does not happen for a density ratio of 1 - in this case I get a nice fluid motion, but of course it behaves as one unified fluid due to equal densities.


Since the described system uses feedback it is quite hard to analyse what causes this instability, so I wanted to ask you first if you have experienced any stability/blowup issues using this modified pressure-solve-scheme or have worked with this method at all. I checked many times my code and can not spot any error :'(

I really hope you can help me out a little to solve this stability issue as this method seems to be quite ideally suited for my use-case and solutions which permit realtime performance are rare


Best regards
XorT is offline   Reply With Quote

Old   June 24, 2020, 18:30
Default
  #2
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 XorT View Post
However the simulation needs to deal with high density ratios as it needs to handle water-air-interface dynamics in a convincing way.
Dear XorT,

the paper you have given deals with a Fourier based spectral collocation method in x- and y-direction (pseudospectral). There are some drawbacks with this kind of approximation:
  • You only can use periodic boundary condtions
  • The method is based on trigonometric polynomials and therefore not useful for strong density gradients. This is due to the Gibbs phenomena.
  • You mentioned that accuracy is not important. However with this kind of approximation you even have spatial spectral order of accuracy.

There are some possibilities to stabilize the method:
  • Filtering in the Fourier modal space (Orzag 2/3 rule which removes the aliasing error)
  • A zero padding in the Fourier modal space (Orzag 3/2 rule which removes the aliasing error)
  • Energy/Entropy stable disctretization via a skew symmetric implementation

But to be honest, if you do not need accuracy try to implement a simple first and second order pressure correction method with a fast Poisson solver. In cases of high density ratios you simply apply the first order method and in smooth regions the second order method.

Regards
FMDenaro likes this.
__________________
Check out my side project:

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

Old   June 25, 2020, 03:46
Default
  #3
New Member
 
Join Date: Jun 2020
Posts: 27
Rep Power: 5
XorT is on a distinguished road
Thank you very much for your reply.
Well, the paper particularly states it is suited for high density ratios and they are using exactly the example of a drop of water in air (and vice versa) with density ratios of 1000. Also the boundary condition does not need to be periodic (and in fact isn't in my case). You can use the DCT for neumann BCs and DST for dirichlet BCs.
Given the paper claims to solve exactly my scenario I expect some error on my side, but I just can not spot it.


In regards to you suggested possibilities to stabilize the method and to your last paragraph, can you recommend some good papers which explain this method in detail ?


Best regards.
XorT is offline   Reply With Quote

Old   June 25, 2020, 09:10
Default
  #4
Senior Member
 
Eifoehn4's Avatar
 
-
Join Date: Jul 2012
Location: Germany
Posts: 184
Rep Power: 13
Eifoehn4 is on a distinguished road
I've only skimmed the paper so far. I have read something about FFT, Fourier space and periodic BCs in x- and y direction. After another brief reading i am still quite sure that the paper only deals with trigonometric basis functions, which only can handle periodic BCs. So far I assume that your implementations are completely based on the paper.

Quote:
Originally Posted by XorT View Post
You can use the DCT for neumann BCs and DST for dirichlet BCs.
That's completely right! However, i do not see where they mentioned the fast discrete cosines/sine transform. In that case the method would rather be based on a Chebyshev expansion. The keyword was never mentioned in the paper.

Quote:
Originally Posted by XorT View Post
In regards to you suggested possibilities to stabilize the method and to your last paragraph, can you recommend some good papers which explain this method in detail ?
The stabilization methods i mentioned can also be applied to spectral collocation Chebyshev methods. The only thing you have to consider is that you are playing with Chebyshev polynomials rather than a trigonometric ones. The books of "John P. Boyd" or "Canuto, et al.", something like "Spectral methods ..." is quite standard. With the stabilization you will not prevent the Gibbs oscillations, but the code won't blow up.

From my experience (see my signature --> Youtube channel), i know that these simulations are quite fast, however i think it would be more reasonable to use a low order method in your context.

Regards
__________________
Check out my side project:

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

Old   June 25, 2020, 09:20
Default
  #5
New Member
 
Join Date: Jun 2020
Posts: 27
Rep Power: 5
XorT is on a distinguished road
Yes, the use of different BC is coming from another paper.
What I don't understand in your explanation is why the use of the particular details on how I solve the poisson equation should affect the result (or blow up in my case) as I verified that the solution suffices the poisson-equation to machine precision. Also I verified that after the correction step my velocity field really is divergence free. Therefore, I would not have suspected the problem to be in the FFT (or DCT) based solver as it obtains correct results.
XorT is offline   Reply With Quote

Old   June 25, 2020, 09:39
Default
  #6
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 XorT View Post
Yes, the use of different BC is coming from another paper.
What I don't understand in your explanation is why the use of the particular details on how I solve the poisson equation should affect the result (or blow up in my case) as I verified that the solution suffices the poisson-equation to machine precision. Also I verified that after the correction step my velocity field really is divergence free. Therefore, I would not have suspected the problem to be in the FFT (or DCT) based solver as it obtains correct results.
I am not quite sure what you mean exactly. I don't say, that the solution procedure for the Poisson equation is the key problem. With the FFT (or DCT/DST) the Poisson equation is solved in a fully coupled (spectral) sense in x- and y-direction, exactly. I don't know, what the paper does in z-direction. An exact solution of the Poisson equation up to machine precision has nothing to do with stability. You can solve the system exactly, but the code will still blow up due to some other numerical reasons. One issue may be the Gibbs oscillations. They can only be prevented with a smaller stencil.

At the moment i do not have the time to read the paper in detail. Do you know how the paper handles the huge density gradients in the context of the spectral approximation in x- and y-direction?

Regards
__________________
Check out my side project:

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

Old   June 25, 2020, 09:48
Default
  #7
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 am not quite sure what you mean exactly. I don't say, that the solution procedure for the Poisson equation is the key problem. With the FFT (or DCT/DST) the Poisson equation is solved in a fully coupled (spectral) sense in x- and y-direction, exactly. I don't know, what the paper does in z-direction. An exact solution of the Poisson equation up to machine precision has nothing to do with stability. You can solve the system exactly, but the code will still blow up due to some other numerical reasons. One issue may be the Gibbs oscillations. They can only be prevented with a smaller stencil.

At the moment i do not have the time to read the paper in detail. Do you know how the paper handles the huge density gradients in the context of the spectral approximation in x- and y-direction?

I thought you were hinting that the spectral solver and the gibbs phenomenon may be to blame. Maybe I misunderstood.
Also my simulation is 2D only, both dimensions solved using FFT. I use the smallest stencil possible - 5-point.




In the paper the important part (for me) is section 2.2.
It explains how they change the variable coefficient poisson equation (which is usually solved) to a constant coefficient poisson equation such that more efficient solvers such as FastPoissonSolvers can be employed. It basically consists of splitting the variable poisson equation into a constant part (which is actually solved) and a variable part (which is only approximated from previous timesteps).
This results in a tight coupling of the solved and the approximated part of the solution. This makes it very difficult to analyse for me where it goes wrong as it becomes a feedback system with a state (what the call p^).
XorT is offline   Reply With Quote

Old   June 25, 2020, 11:14
Default
  #8
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
Yes, the use of different BC is coming from another paper.
What I don't understand in your explanation is why the use of the particular details on how I solve the poisson equation should affect the result (or blow up in my case) as I verified that the solution suffices the poisson-equation to machine precision. Also I verified that after the correction step my velocity field really is divergence free. Therefore, I would not have suspected the problem to be in the FFT (or DCT) based solver as it obtains correct results.



How do you verify that the final velocity field is divergence-free? Provide the exact details.
FMDenaro is online now   Reply With Quote

Old   June 25, 2020, 11:15
Default
  #9
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 thought you were hinting that the spectral solver and the gibbs phenomenon may be to blame. Maybe I misunderstood.
Also my simulation is 2D only, both dimensions solved using FFT. I use the smallest stencil possible - 5-point.






This sentence has no meaning for me...
FMDenaro is online now   Reply With Quote

Old   June 25, 2020, 11:17
Default
  #10
New Member
 
Join Date: Jun 2020
Posts: 27
Rep Power: 5
XorT is on a distinguished road
Quote:
Originally Posted by FMDenaro View Post
How do you verify that the final velocity field is divergence-free? Provide the exact details.

Well, I numerically compute the divergence of the velocity field after the pressure correction.
XorT is offline   Reply With Quote

Old   June 25, 2020, 11:19
Default
  #11
New Member
 
Join Date: Jun 2020
Posts: 27
Rep Power: 5
XorT is on a distinguished road
Quote:
Originally Posted by FMDenaro View Post
This sentence has no meaning for me...
First part is about each dimension is solved using FFT (since it is separable) I use 2 1D FFTs. Second part relates to a comment made by Eifoehn4 to use a smaller stencil to avoid gibbs phenomenon.
XorT is offline   Reply With Quote

Old   June 25, 2020, 11:19
Default
  #12
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 numerically compute the divergence of the velocity field after the pressure correction.



That means nothing if you do not specify what you are really doing numerically. What kind of discretization are you using for this check? Be aware, the pressure equation is nothing but the divergence-free constraint, properly rewritten. Therefore you must check the div v constraint using the same discretization.
FMDenaro is online now   Reply With Quote

Old   June 25, 2020, 11:21
Default
  #13
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
First part is about each dimension is solved using FFT (since it is separable) I use 2 1D FFTs. Second part relates to a comment made by Eifoehn4 to use a smaller stencil to avoid gibbs phenomenon.

A spectral method uses a global stencil, you cannot think about the classic second order accurate central FD discretization of the laplacian operator.
FMDenaro is online now   Reply With Quote

Old   June 25, 2020, 11:23
Default
  #14
New Member
 
Join Date: Jun 2020
Posts: 27
Rep Power: 5
XorT is on a distinguished road
Quote:
Originally Posted by FMDenaro View Post
That means nothing if you do not specify what you are really doing numerically. What kind of discretization are you using for this check? Be aware, the pressure equation is nothing but the divergence-free constraint, properly rewritten. Therefore you must check the div v constraint using the same discretization.

Since I use a staggered grid for velocities and divergence is defined at the center I use:


d[x,y] = u[x+1/2,y] - u[x-1/2,y] + v[x,y+1/2] - v[x,y-1/2]


grid spacing is 1. pretty standard, i'd say.
XorT is offline   Reply With Quote

Old   June 25, 2020, 11:29
Default
  #15
New Member
 
Join Date: Jun 2020
Posts: 27
Rep Power: 5
XorT is on a distinguished road
Quote:
Originally Posted by FMDenaro View Post
A spectral method uses a global stencil, you cannot think about the classic second order accurate central FD discretization of the laplacian operator.

Ah sorry, then I misinterpreted.

In that case my full domain is included in the FFT. No windowing or anything like this. So the result should be exact if I'm not mistaken (and numerically this seems to be the case).
XorT is offline   Reply With Quote

Old   June 25, 2020, 11:31
Default
  #16
New Member
 
Join Date: Jun 2020
Posts: 27
Rep Power: 5
XorT is on a distinguished road
I don't think I got these fundamentals wrong, since the fluid behaves nicely and as expected if the density is constant.
It is only when non-constant density comes into play that the blow-up occurs, so I expect that feedback approximation scheme used to be the culprit here.
XorT is offline   Reply With Quote

Old   June 25, 2020, 16:11
Default
  #17
Senior Member
 
Eifoehn4's Avatar
 
-
Join Date: Jul 2012
Location: Germany
Posts: 184
Rep Power: 13
Eifoehn4 is on a distinguished road
I took a closer look at the paper and unfortunately some things are not quite clear to me. To be honest, i am not deep in the incompressible world.

The paper deals with some kind of constant/variable density splitting resulting in a pressure Poisson equation of following form:

\nabla^2 p^{n+1} = \nabla \cdot \left[  \left( 1 -  \frac{ \rho_{0} }{ \rho^{n+1}} \right) \nabla \hat{p} \right] + \frac{ \rho_0 }{ \Delta t } \nabla \cdot \boldsymbol{u}^{*}
  • \rho_0 is a constant density
  • \hat{p} is a pressure approximation for the new time step n+1
  • \boldsymbol{u}^{*} is a approximate velocity field the new time step n+1
  • \rho is a variable density

Without detailed explanations, the authors state that a fast Poisson solver via FFT is applied. The pressure field obtained via FFT is based on a spectral polynomial expansion. The spatial approximation for all other variables is based on a different discretization, equidistant first/second order stencils:
  • Does a right hand side with a density ratio of 1000/1 does not cause any problems?
  • Does a different discretization for the Poisson equation does not cause any consistency issues?
  • @XorT: Do you apply the DCT for the Poisson equation on an equidistant mesh?

Regards
__________________
Check out my side project:

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

Old   June 25, 2020, 16:18
Default
  #18
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 took a closer look at the paper and unfortunately some things are not quite clear to me. To be honest, i am not deep in the incompressible world.

The paper deals with some kind of constant/variable density splitting resulting in a pressure Poisson equation of following form:

\nabla^2 p^{n+1} = \nabla \cdot \left[  \left( 1 -  \frac{ \rho_{0} }{ \rho^{n+1}} \right) \nabla \hat{p} \right] + \frac{ \rho_0 }{ \Delta t } \nabla \cdot \boldsymbol{u}^{*}
  • \rho_0 is a constant density
  • \hat{p} is a pressure approximation for the new time step n+1
  • \boldsymbol{u}^{*} is a approximate velocity field the new time step n+1
  • \rho is a variable density

Without detailed explanations, the authors state that a fast Poisson solver via FFT is applied. The pressure field obtained via FFT is based on a spectral polynomial expansion. The spatial approximation for all other variables is based on a different discretization, equidistant first/second order stencils:
  • Does a right hand side with a density ratio of 1000/1 does not cause any problems?
  • Does a different discretization for the Poisson equation does not cause any consistency issues?
  • @XorT: Do you apply the DCT for the Poisson equation on an equidistant mesh?
Yes, I do the DCT on an equidistant mesh (not sure how i could do otherwise tbh). The grid is completely regular.
I don't really understand why this could be a problem to use a FFT based solver. It solves the constant poisson equation exactly and therefore should not be inferior to any other direct solver.

The paper claims, and the experiments show, it can simulate density ratio of even 10000. This is indeed very impressive, that's what caught my attention

I'm not sure where you see that there is a different discretization for the poisson equation. All variables (density, pressure, divergence) are defined on the same grid.
XorT is offline   Reply With Quote

Old   June 25, 2020, 16:45
Default
  #19
Senior Member
 
Eifoehn4's Avatar
 
-
Join Date: Jul 2012
Location: Germany
Posts: 184
Rep Power: 13
Eifoehn4 is on a distinguished road
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?
__________________
Check out my side project:

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

Old   June 25, 2020, 16:47
Default
  #20
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
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?
FMDenaro is online now   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 07:46.