
[Sponsors] 
Stability issues with "modified pressurecorrection" scheme 

LinkBack  Thread Tools  Search this Thread  Display Modes 
June 24, 2020, 16:39 
Stability issues with "modified pressurecorrection" scheme

#1 
New Member
Join Date: Jun 2020
Posts: 27
Rep Power: 5 
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 (usecase is computergraphics). However the simulation needs to deal with high density ratios as it needs to handle waterairinterface dynamics in a convincing way. The method described in the paper "A fast pressurecorrection method for incompressible twofluid flows" (https://doddm.com/publications/2014jcpmdaf.pdf) appears very interesting as it can employ a FFT based FastPoisson solver to solve the pressure equation, which can achieve realtime performance on highresolution 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 ydirection on a background of initially still air (1kg/mē). As quantities I keep only cellcentered 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 divergencefree 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 pressuresolvescheme 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 usecase and solutions which permit realtime performance are rare Best regards 

June 24, 2020, 18:30 

#2  
Senior Member

Join Date: Jul 2012
Location: Germany
Posts: 184
Rep Power: 13 
Quote:
the paper you have given deals with a Fourier based spectral collocation method in x and ydirection (pseudospectral). There are some drawbacks with this kind of approximation:
There are some possibilities to stabilize the method:
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 

June 25, 2020, 03:46 

#3 
New Member
Join Date: Jun 2020
Posts: 27
Rep Power: 5 
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, 09:10 

#4  
Senior Member

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

June 25, 2020, 09:20 

#5 
New Member
Join Date: Jun 2020
Posts: 27
Rep Power: 5 
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 poissonequation 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, 09:39 

#6  
Senior Member

Join Date: Jul 2012
Location: Germany
Posts: 184
Rep Power: 13 
Quote:
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 ydirection? Regards 

June 25, 2020, 09:48 

#7  
New Member
Join Date: Jun 2020
Posts: 27
Rep Power: 5 
Quote:
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  5point. 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^). 

June 25, 2020, 11:14 

#8  
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,793
Rep Power: 71 
Quote:
How do you verify that the final velocity field is divergencefree? Provide the exact details. 

June 25, 2020, 11:15 

#9  
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,793
Rep Power: 71 
Quote:
This sentence has no meaning for me... 

June 25, 2020, 11:17 

#10 
New Member
Join Date: Jun 2020
Posts: 27
Rep Power: 5 

June 25, 2020, 11:19 

#11 
New Member
Join Date: Jun 2020
Posts: 27
Rep Power: 5 

June 25, 2020, 11:19 

#12  
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,793
Rep Power: 71 
Quote:
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 divergencefree constraint, properly rewritten. Therefore you must check the div v constraint using the same discretization. 

June 25, 2020, 11:21 

#13  
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,793
Rep Power: 71 
Quote:
A spectral method uses a global stencil, you cannot think about the classic second order accurate central FD discretization of the laplacian operator. 

June 25, 2020, 11:23 

#14  
New Member
Join Date: Jun 2020
Posts: 27
Rep Power: 5 
Quote:
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[x1/2,y] + v[x,y+1/2]  v[x,y1/2] grid spacing is 1. pretty standard, i'd say. 

June 25, 2020, 11:29 

#15  
New Member
Join Date: Jun 2020
Posts: 27
Rep Power: 5 
Quote:
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). 

June 25, 2020, 11:31 

#16 
New Member
Join Date: Jun 2020
Posts: 27
Rep Power: 5 
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 nonconstant density comes into play that the blowup occurs, so I expect that feedback approximation scheme used to be the culprit here. 

June 25, 2020, 16:11 

#17 
Senior Member

Join Date: Jul 2012
Location: Germany
Posts: 184
Rep Power: 13 
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:
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:
Regards 

June 25, 2020, 16:18 

#18 
New Member
Join Date: Jun 2020
Posts: 27
Rep Power: 5 
Quote:
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. 

June 25, 2020, 16:45 

#19 
Senior Member

Join Date: Jul 2012
Location: Germany
Posts: 184
Rep Power: 13 
I ask because:
Can you briefly explain how you handle this? 

June 25, 2020, 16:47 

#20 
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,793
Rep Power: 71 
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 nonvanishing (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? 

Tags 
modified solver, pressure based solver 
Thread Tools  Search this Thread 
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Pressure fields in FOAM, p field, total pressure, etc.  Tobi  OpenFOAM PostProcessing  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 