Stability issues with "modified pressure-correction" scheme

XorT · June 24, 2020, 17:39

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

Eifoehn4 · June 24, 2020, 19:30

Quote:

Originally Posted by XorT

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

XorT · June 25, 2020, 04:46

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.

Eifoehn4 · June 25, 2020, 10:10

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

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

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

XorT · June 25, 2020, 10:20

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.

Eifoehn4 · June 25, 2020, 10:39

Quote:

Originally Posted by XorT

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

XorT · June 25, 2020, 10:48

Quote:

Originally Posted by Eifoehn4

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^).

FMDenaro · June 25, 2020, 12:14

Quote:

Originally Posted by XorT

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 · June 25, 2020, 12:15

Quote:

Originally Posted by XorT

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...

XorT · June 25, 2020, 12:17

Quote:

Originally Posted by FMDenaro

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 · June 25, 2020, 12:19

Quote:

Originally Posted by FMDenaro

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.

FMDenaro · June 25, 2020, 12:19

Quote:

Originally Posted by XorT

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 · June 25, 2020, 12:21

Quote:

Originally Posted by XorT

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.

XorT · June 25, 2020, 12:23

Quote:

Originally Posted by FMDenaro

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 · June 25, 2020, 12:29

Quote:

Originally Posted by FMDenaro

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 · June 25, 2020, 12:31

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.

Eifoehn4 · June 25, 2020, 17:11

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

XorT · June 25, 2020, 17:18

Quote:

Originally Posted by Eifoehn4

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.

Eifoehn4 · June 25, 2020, 17:45

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?

FMDenaro · June 25, 2020, 17:47

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?

June 24, 2020, 17:39	Stability issues with "modified pressure-correction" scheme	#1
XorT New Member Join Date: Jun 2020 Posts: 27 Rep Power: 6	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

June 25, 2020, 17:11		#17
Eifoehn4 Senior Member - Join Date: Jul 2012 Location: Germany Posts: 184 Rep Power: 14	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.

June 25, 2020, 17:45		#19
Eifoehn4 Senior Member - Join Date: Jul 2012 Location: Germany Posts: 184 Rep Power: 14	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.

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 02:33
Wind tunnel Boundary Conditions in Fluent	metmet	FLUENT	6	October 30, 2019 13:23
outlet pressure Boundary settings -velocity streamline under ambient temp.conditions	Vishnu_bharathi	CFX	12	November 21, 2017 07:56
Setting up the pressure variation due to tornado in a duct(UDF)+animation	guillaume1990	FLUENT	0	March 3, 2014 12:59
Pressure BC for combustion chamber	Giuki	FLUENT	1	July 19, 2011 12:35

June 25, 2020, 04:46		#3
XorT New Member Join Date: Jun 2020 Posts: 27 Rep Power: 6	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.

June 25, 2020, 10:20		#5
XorT New Member Join Date: Jun 2020 Posts: 27 Rep Power: 6	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.

June 25, 2020, 12:31		#16
XorT New Member Join Date: Jun 2020 Posts: 27 Rep Power: 6	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.

June 25, 2020, 17:47		#20
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,882 Rep Power: 73	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?